www.mooseframework.org
ComputeMultipleInelasticStress.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 
9 #include "StressUpdateBase.h"
10 #include "MooseException.h"
11 
12 template <>
13 InputParameters
15 {
16  InputParameters params = validParams<ComputeFiniteStrainElasticStress>();
17  params.addClassDescription("Compute state (stress and internal parameters such as plastic "
18  "strains and internal parameters) using an iterative process. "
19  "Combinations of creep models and plastic models may be used");
20  params.addParam<unsigned int>("max_iterations",
21  30,
22  "Maximum number of the stress update "
23  "iterations over the stress change after all "
24  "update materials are called");
25  params.addParam<Real>("relative_tolerance",
26  1e-5,
27  "Relative convergence tolerance for the stress "
28  "update iterations over the stress change "
29  "after all update materials are called");
30  params.addParam<Real>("absolute_tolerance",
31  1e-5,
32  "Absolute convergence tolerance for the stress "
33  "update iterations over the stress change "
34  "after all update materials are called");
35  params.addParam<bool>(
36  "output_iteration_info",
37  false,
38  "Set to true to output stress update iteration information over the stress change");
39  params.addParam<bool>("perform_finite_strain_rotations",
40  true,
41  "Tensors are correctly rotated in "
42  "finite-strain simulations. For "
43  "optimal performance you can set "
44  "this to 'false' if you are only "
45  "ever using small strains");
46  MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
47  params.addParam<MooseEnum>(
48  "tangent_operator",
49  tangent_operator,
50  "Type of tangent operator to return. 'elastic': return the "
51  "elasticity tensor. 'nonlinear': return the full, general consistent tangent "
52  "operator.");
53  params.addRequiredParam<std::vector<MaterialName>>(
54  "inelastic_models",
55  "The material objects to use to calculate stress and inelastic strains. "
56  "Note: specify creep models first and plasticity models second.");
57  params.addParam<std::vector<Real>>("combined_inelastic_strain_weights",
58  "The combined_inelastic_strain Material Property is a "
59  "weighted sum of the model inelastic strains. This parameter "
60  "is a vector of weights, of the same length as "
61  "inelastic_models. Default = '1 1 ... 1'. This "
62  "parameter is set to 1 if the number of models = 1");
63  params.addParam<bool>(
64  "cycle_models", false, "At timestep N use only inelastic model N % num_models.");
65  return params;
66 }
67 
70  _max_iterations(parameters.get<unsigned int>("max_iterations")),
71  _relative_tolerance(parameters.get<Real>("relative_tolerance")),
72  _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
73  _output_iteration_info(getParam<bool>("output_iteration_info")),
74  _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
75  _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_base_name + "elasticity_tensor")),
76  _elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "elastic_strain")),
77  _strain_increment(getMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
78  _inelastic_strain(declareProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
79  _inelastic_strain_old(
80  getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
81  _tangent_operator_type(getParam<MooseEnum>("tangent_operator").getEnum<TangentOperatorEnum>()),
82  _num_models(getParam<std::vector<MaterialName>>("inelastic_models").size()),
83  _inelastic_weights(isParamValid("combined_inelastic_strain_weights")
84  ? getParam<std::vector<Real>>("combined_inelastic_strain_weights")
85  : std::vector<Real>(_num_models, true)),
86  _consistent_tangent_operator(_num_models),
87  _cycle_models(getParam<bool>("cycle_models")),
88  _matl_timestep_limit(declareProperty<Real>("matl_timestep_limit"))
89 {
90  if (_inelastic_weights.size() != _num_models)
91  mooseError(
92  "ComputeMultipleInelasticStress: combined_inelastic_strain_weights must contain the same "
93  "number of entries as inelastic_models ",
94  _inelastic_weights.size(),
95  " ",
96  _num_models);
97 }
98 
99 void
101 {
103  _inelastic_strain[_qp].zero();
104 }
105 
106 void
108 {
111 
112  std::vector<MaterialName> models = getParam<std::vector<MaterialName>>("inelastic_models");
113  for (unsigned int i = 0; i < _num_models; ++i)
114  {
115  StressUpdateBase * rrr = dynamic_cast<StressUpdateBase *>(&getMaterialByName(models[i]));
116  if (rrr)
117  {
118  _models.push_back(rrr);
120  mooseError("Model " + models[i] +
121  " requires an isotropic elasticity tensor, but the one supplied is not "
122  "guaranteed isotropic");
123  }
124  else
125  mooseError("Model " + models[i] + " is not compatible with ComputeMultipleInelasticStress");
126  }
127 }
128 
129 void
131 {
135 }
136 
137 void
139 {
140  RankTwoTensor elastic_strain_increment;
141  RankTwoTensor combined_inelastic_strain_increment;
142 
143  if (_num_models == 0)
144  {
146 
147  // If the elasticity tensor values have changed and the tensor is isotropic,
148  // use the old strain to calculate the old stress
150  {
151  _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + _strain_increment[_qp]);
152  // InitialStress Deprecation: remove these lines
156  }
157  else
158  _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * _strain_increment[_qp];
159 
160  if (_fe_problem.currentlyComputingJacobian())
162  }
163  else
164  {
165  if (_num_models == 1 || _cycle_models)
166  updateQpStateSingleModel((_t_step - 1) % _num_models,
167  elastic_strain_increment,
168  combined_inelastic_strain_increment);
169  else
170  updateQpState(elastic_strain_increment, combined_inelastic_strain_increment);
171 
172  _elastic_strain[_qp] = _elastic_strain_old[_qp] + elastic_strain_increment;
173  _inelastic_strain[_qp] = _inelastic_strain_old[_qp] + combined_inelastic_strain_increment;
174  }
175 }
176 
177 void
178 ComputeMultipleInelasticStress::finiteStrainRotation(const bool force_elasticity_rotation)
179 {
180  _elastic_strain[_qp] =
181  _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
182  _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
183  _inelastic_strain[_qp] =
184  _rotation_increment[_qp] * _inelastic_strain[_qp] * _rotation_increment[_qp].transpose();
185  if (force_elasticity_rotation ||
188  _Jacobian_mult[_qp].rotate(_rotation_increment[_qp]);
189 }
190 
191 void
192 ComputeMultipleInelasticStress::updateQpState(RankTwoTensor & elastic_strain_increment,
193  RankTwoTensor & combined_inelastic_strain_increment)
194 {
195  if (_output_iteration_info == true)
196  {
197  _console << std::endl
198  << "iteration output for ComputeMultipleInelasticStress solve:"
199  << " time=" << _t << " int_pt=" << _qp << std::endl;
200  }
201  Real l2norm_delta_stress, first_l2norm_delta_stress;
202  unsigned int counter = 0;
203 
204  std::vector<RankTwoTensor> inelastic_strain_increment;
205  inelastic_strain_increment.resize(_num_models);
206 
207  for (unsigned i_rmm = 0; i_rmm < _models.size(); ++i_rmm)
208  inelastic_strain_increment[i_rmm].zero();
209 
210  RankTwoTensor stress_max, stress_min;
211 
212  do
213  {
214  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
215  {
216  _models[i_rmm]->setQp(_qp);
217 
218  // initially assume the strain is completely elastic
219  elastic_strain_increment = _strain_increment[_qp];
220  // and subtract off all inelastic strain increments calculated so far
221  // except the one that we're about to calculate
222  for (unsigned j_rmm = 0; j_rmm < _num_models; ++j_rmm)
223  if (i_rmm != j_rmm)
224  elastic_strain_increment -= inelastic_strain_increment[j_rmm];
225 
226  // form the trial stress, with the check for changed elasticity constants
228  {
229  _stress[_qp] =
230  _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
231  // InitialStress Deprecation: remove these lines
235  }
236  else
237  _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
238 
239  // given a trial stress (_stress[_qp]) and a strain increment (elastic_strain_increment)
240  // let the i^th model produce an admissible stress (as _stress[_qp]), and decompose
241  // the strain increment into an elastic part (elastic_strain_increment) and an
242  // inelastic part (inelastic_strain_increment[i_rmm])
244  elastic_strain_increment,
245  inelastic_strain_increment[i_rmm],
247 
248  if (i_rmm == 0)
249  {
250  stress_max = _stress[_qp];
251  stress_min = _stress[_qp];
252  }
253  else
254  {
255  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
256  {
257  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
258  {
259  if (_stress[_qp](i, j) > stress_max(i, j))
260  stress_max(i, j) = _stress[_qp](i, j);
261  else if (stress_min(i, j) > _stress[_qp](i, j))
262  stress_min(i, j) = _stress[_qp](i, j);
263  }
264  }
265  }
266  }
267 
268  // now check convergence in the stress:
269  // once the change in stress is within tolerance after each recompute material
270  // consider the stress to be converged
271  RankTwoTensor delta_stress(stress_max - stress_min);
272  l2norm_delta_stress = delta_stress.L2norm();
273  if (counter == 0)
274  first_l2norm_delta_stress = l2norm_delta_stress;
275 
276  if (_output_iteration_info == true)
277  {
278  _console << "stress iteration number = " << counter << "\n"
279  << " relative l2 norm delta stress = "
280  << (0 == first_l2norm_delta_stress ? 0
281  : l2norm_delta_stress / first_l2norm_delta_stress)
282  << "\n"
283  << " stress convergence relative tolerance = " << _relative_tolerance << "\n"
284  << " absolute l2 norm delta stress = " << l2norm_delta_stress << "\n"
285  << " stress convergence absolute tolerance = " << _absolute_tolerance << std::endl;
286  }
287  ++counter;
288  } while (counter < _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
289  (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance &&
290  _num_models != 1);
291 
292  if (counter == _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
293  (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance)
294  throw MooseException("Max stress iteration hit during ComputeMultipleInelasticStress solve!");
295 
296  combined_inelastic_strain_increment.zero();
297  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
298  combined_inelastic_strain_increment +=
299  _inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
300 
301  if (_fe_problem.currentlyComputingJacobian())
303 
304  _matl_timestep_limit[_qp] = 0.0;
305  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
306  _matl_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
307 
308  _matl_timestep_limit[_qp] = 1.0 / _matl_timestep_limit[_qp];
309 }
310 
311 void
313 {
316  else
317  {
318  const RankFourTensor E_inv = _elasticity_tensor[_qp].invSymm();
320  for (unsigned i_rmm = 1; i_rmm < _num_models; ++i_rmm)
321  _Jacobian_mult[_qp] = _consistent_tangent_operator[i_rmm] * E_inv * _Jacobian_mult[_qp];
322  }
323 }
324 
325 void
327  unsigned model_number,
328  RankTwoTensor & elastic_strain_increment,
329  RankTwoTensor & combined_inelastic_strain_increment)
330 {
331  for (auto model : _models)
332  model->setQp(_qp);
333 
334  elastic_strain_increment = _strain_increment[_qp];
335 
336  // If the elasticity tensor values have changed and the tensor is isotropic,
337  // use the old strain to calculate the old stress
339  {
340  _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
341  // InitialStress Deprecation: remove these lines
345  }
346  else
347  _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
348 
349  computeAdmissibleState(model_number,
350  elastic_strain_increment,
351  combined_inelastic_strain_increment,
352  _Jacobian_mult[_qp]);
353 
354  _matl_timestep_limit[_qp] = _models[0]->computeTimeStepLimit();
355 
356  /* propagate internal variables, etc, to this timestep for those inelastic models where
357  * "updateState" is not called */
358  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
359  if (i_rmm != model_number)
360  _models[i_rmm]->propagateQpStatefulProperties();
361 }
362 
363 void
365  RankTwoTensor & elastic_strain_increment,
366  RankTwoTensor & inelastic_strain_increment,
367  RankFourTensor & consistent_tangent_operator)
368 {
369  _models[model_number]->updateState(elastic_strain_increment,
370  inelastic_strain_increment,
371  _rotation_increment[_qp],
372  _stress[_qp],
373  _stress_old[_qp],
374  _elasticity_tensor[_qp],
375  _elastic_strain_old[_qp],
377  consistent_tangent_operator);
378 }
bool _is_elasticity_tensor_guaranteed_isotropic
is the elasticity tensor guaranteed to be isotropic?
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
const std::string _elasticity_tensor_name
std::vector< RankFourTensor > _consistent_tangent_operator
the consistent tangent operators computed by each plastic model
const MaterialProperty< RankTwoTensor > & _rotation_increment
const MaterialProperty< RankTwoTensor > & _elastic_strain_old
MaterialProperty< Real > & _matl_timestep_limit
const MaterialProperty< RankTwoTensor > & _strain_increment
virtual void updateQpState(RankTwoTensor &elastic_strain_increment, RankTwoTensor &combined_inelastic_strain_increment)
Given the _strain_increment[_qp], iterate over all of the user-specified recompute materials in order...
StressUpdateBase is a material that is not called by MOOSE because of the compute=false flag set in t...
const std::vector< Real > _inelastic_weights
_inelastic_strain = sum_i (_inelastic_weights_i * inelastic_strain_from_model_i)
MaterialProperty< RankTwoTensor > & _stress
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 ...
virtual void updateQpStateSingleModel(unsigned model_number, RankTwoTensor &elastic_strain_increment, RankTwoTensor &combined_inelastic_strain_increment)
An optimised version of updateQpState that gets used when the number of plastic models is unity...
TangentOperatorEnum
what sort of Tangent operator to calculate
virtual bool requiresIsotropicTensor()=0
Does the model require the elasticity tensor to be isotropic?
const bool _cycle_models
whether to cycle through the models, using only one model per timestep
const unsigned _num_models
number of plastic models
const MaterialProperty< RankFourTensor > & _elasticity_tensor
Rank-4 and Rank-2 elasticity and elastic strain tensors.
ComputeFiniteStrainElasticStress computes the stress following elasticity theory for finite strains...
void addQpInitialStress()
InitialStress Deprecation: remove this method.
virtual void computeAdmissibleState(unsigned model_number, RankTwoTensor &elastic_strain_increment, RankTwoTensor &inelastic_strain_increment, RankFourTensor &consistent_tangent_operator)
Given a trial stress (_stress[_qp]) and a strain increment (elastic_strain_increment) let the model_n...
enum ComputeMultipleInelasticStress::TangentOperatorEnum _tangent_operator_type
const MaterialProperty< RankTwoTensor > & _stress_old
virtual void finiteStrainRotation(const bool force_elasticity_rotation=false)
Rotate _elastic_strain, _stress, _inelastic_strain, and _Jacobian_mult to the new configuration...
virtual void initQpStatefulProperties() override
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...
MaterialProperty< RankTwoTensor > & _elastic_strain
bool hasGuaranteedMaterialProperty(const MaterialPropertyName &prop, Guarantee guarantee)
InputParameters validParams< ComputeMultipleInelasticStress >()
virtual void rotateQpInitialStress()
InitialStress Deprecation: remove this method.
virtual void computeQpJacobianMult()
Using _elasticity_tensor[_qp] and the consistent tangent operators, _comsistent_tangent_operator[...] computed by the inelastic models, compute _Jacobian_mult[_qp].
static unsigned int counter
InputParameters validParams< ComputeFiniteStrainElasticStress >()
ComputeMultipleInelasticStress(const InputParameters &parameters)
MaterialProperty< RankTwoTensor > & _inelastic_strain
The sum of the inelastic strains that come from the plastic models.
const unsigned int _max_iterations
Input parameters associated with the recompute iteration to return the stress state to the yield surf...