www.mooseframework.org
ComputeMultipleInelasticStress.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 
12 #include "StressUpdateBase.h"
13 #include "MooseException.h"
14 
16 
17 template <>
18 InputParameters
20 {
21  InputParameters params = validParams<ComputeFiniteStrainElasticStress>();
22  params.addClassDescription("Compute state (stress and internal parameters such as plastic "
23  "strains and internal parameters) using an iterative process. "
24  "Combinations of creep models and plastic models may be used.");
25  params.addParam<unsigned int>("max_iterations",
26  30,
27  "Maximum number of the stress update "
28  "iterations over the stress change after all "
29  "update materials are called");
30  params.addParam<Real>("relative_tolerance",
31  1e-5,
32  "Relative convergence tolerance for the stress "
33  "update iterations over the stress change "
34  "after all update materials are called");
35  params.addParam<Real>("absolute_tolerance",
36  1e-5,
37  "Absolute convergence tolerance for the stress "
38  "update iterations over the stress change "
39  "after all update materials are called");
40  params.addParam<bool>(
41  "internal_solve_full_iteration_history",
42  false,
43  "Set to true to output stress update iteration information over the stress change");
44  params.addParam<bool>("perform_finite_strain_rotations",
45  true,
46  "Tensors are correctly rotated in "
47  "finite-strain simulations. For "
48  "optimal performance you can set "
49  "this to 'false' if you are only "
50  "ever using small strains");
51  MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
52  params.addParam<MooseEnum>(
53  "tangent_operator",
54  tangent_operator,
55  "Type of tangent operator to return. 'elastic': return the "
56  "elasticity tensor. 'nonlinear': return the full, general consistent tangent "
57  "operator.");
58  params.addRequiredParam<std::vector<MaterialName>>(
59  "inelastic_models",
60  "The material objects to use to calculate stress and inelastic strains. "
61  "Note: specify creep models first and plasticity models second.");
62  params.addParam<std::vector<Real>>("combined_inelastic_strain_weights",
63  "The combined_inelastic_strain Material Property is a "
64  "weighted sum of the model inelastic strains. This parameter "
65  "is a vector of weights, of the same length as "
66  "inelastic_models. Default = '1 1 ... 1'. This "
67  "parameter is set to 1 if the number of models = 1");
68  params.addParam<bool>(
69  "cycle_models", false, "At timestep N use only inelastic model N % num_models.");
70  return params;
71 }
72 
75  _max_iterations(parameters.get<unsigned int>("max_iterations")),
76  _relative_tolerance(parameters.get<Real>("relative_tolerance")),
77  _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
78  _internal_solve_full_iteration_history(getParam<bool>("internal_solve_full_iteration_history")),
79  _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
80  _elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "elastic_strain")),
81  _strain_increment(getMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
82  _inelastic_strain(declareProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
83  _inelastic_strain_old(
84  getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
85  _tangent_operator_type(getParam<MooseEnum>("tangent_operator").getEnum<TangentOperatorEnum>()),
86  _num_models(getParam<std::vector<MaterialName>>("inelastic_models").size()),
87  _tangent_computation_flag(_num_models, false),
88  _tangent_calculation_method(TangentCalculationMethod::ELASTIC),
89  _inelastic_weights(isParamValid("combined_inelastic_strain_weights")
90  ? getParam<std::vector<Real>>("combined_inelastic_strain_weights")
91  : std::vector<Real>(_num_models, true)),
92  _consistent_tangent_operator(_num_models),
93  _cycle_models(getParam<bool>("cycle_models")),
94  _matl_timestep_limit(declareProperty<Real>("matl_timestep_limit")),
95  _identity_symmetric_four(RankFourTensor::initIdentitySymmetricFour)
96 {
97  if (_inelastic_weights.size() != _num_models)
98  mooseError(
99  "ComputeMultipleInelasticStress: combined_inelastic_strain_weights must contain the same "
100  "number of entries as inelastic_models ",
101  _inelastic_weights.size(),
102  " ",
103  _num_models);
104 }
105 
106 void
108 {
110  _inelastic_strain[_qp].zero();
111 }
112 
113 void
115 {
118 
119  std::vector<MaterialName> models = getParam<std::vector<MaterialName>>("inelastic_models");
120 
121  for (unsigned int i = 0; i < _num_models; ++i)
122  {
123  StressUpdateBase * rrr = dynamic_cast<StressUpdateBase *>(&getMaterialByName(models[i]));
124 
125  if (rrr)
126  {
127  _models.push_back(rrr);
129  mooseError("Model " + models[i] +
130  " requires an isotropic elasticity tensor, but the one supplied is not "
131  "guaranteed isotropic");
132  }
133  else
134  mooseError("Model " + models[i] + " is not compatible with ComputeMultipleInelasticStress");
135  }
136 
137  // Check if tangent calculation methods are consistent. If all models have
138  // TangentOperatorEnum::ELASTIC or tangent_operator is set by the user as elasic, then the tangent
139  // is never calculated: J_tot = C. If PARTIAL and NONE models are present, utilize PARTIAL
140  // formulation: J_tot = (I + J_1 + ... J_N)^-1 C. If FULL and NONE models are present, utilize
141  // FULL formulation: J_tot = J_1 * C^-1 * J_2 * C^-1 * ... J_N * C. If PARTIAL and FULL models are
142  // present, error out.
143 
145  {
146  bool full_present = false;
147  bool partial_present = false;
148  for (unsigned int i = 0; i < _num_models; ++i)
149  {
150  if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::FULL)
151  {
152  full_present = true;
153  _tangent_computation_flag[i] = true;
155  }
156  else if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::PARTIAL)
157  {
158  partial_present = true;
159  _tangent_computation_flag[i] = true;
161  }
162  }
163  if (full_present && partial_present)
164  mooseError("In ",
165  _name,
166  ": Models that calculate the full tangent operator and the partial tangent "
167  "operator are being combined. Either set tangent_operator to elastic, implement "
168  "the corrent tangent formulations, or use different models.");
169  }
170 }
171 
172 void
174 {
178 }
179 
180 void
182 {
183  RankTwoTensor elastic_strain_increment;
184  RankTwoTensor combined_inelastic_strain_increment;
185 
186  if (_num_models == 0)
187  {
189 
190  // If the elasticity tensor values have changed and the tensor is isotropic,
191  // use the old strain to calculate the old stress
194  else
195  _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * _strain_increment[_qp];
196 
197  if (_fe_problem.currentlyComputingJacobian())
199 
200  _matl_timestep_limit[_qp] = std::numeric_limits<Real>::max();
201  }
202  else
203  {
204  if (_num_models == 1 || _cycle_models)
205  updateQpStateSingleModel((_t_step - 1) % _num_models,
206  elastic_strain_increment,
207  combined_inelastic_strain_increment);
208  else
209  updateQpState(elastic_strain_increment, combined_inelastic_strain_increment);
210 
211  _elastic_strain[_qp] = _elastic_strain_old[_qp] + elastic_strain_increment;
212  _inelastic_strain[_qp] = _inelastic_strain_old[_qp] + combined_inelastic_strain_increment;
213  }
214 }
215 
216 void
217 ComputeMultipleInelasticStress::finiteStrainRotation(const bool force_elasticity_rotation)
218 {
219  _elastic_strain[_qp] =
220  _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
221  _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
222  _inelastic_strain[_qp] =
223  _rotation_increment[_qp] * _inelastic_strain[_qp] * _rotation_increment[_qp].transpose();
224  if (force_elasticity_rotation ||
227  _Jacobian_mult[_qp].rotate(_rotation_increment[_qp]);
228 }
229 
230 void
232  RankTwoTensor & combined_inelastic_strain_increment)
233 {
235  {
236  _console << std::endl
237  << "iteration output for ComputeMultipleInelasticStress solve:"
238  << " time=" << _t << " int_pt=" << _qp << std::endl;
239  }
240  Real l2norm_delta_stress;
241  Real first_l2norm_delta_stress = 1.0;
242  unsigned int counter = 0;
243 
244  std::vector<RankTwoTensor> inelastic_strain_increment;
245  inelastic_strain_increment.resize(_num_models);
246 
247  for (unsigned i_rmm = 0; i_rmm < _models.size(); ++i_rmm)
248  inelastic_strain_increment[i_rmm].zero();
249 
250  RankTwoTensor stress_max, stress_min;
251 
252  do
253  {
254  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
255  {
256  _models[i_rmm]->setQp(_qp);
257 
258  // initially assume the strain is completely elastic
259  elastic_strain_increment = _strain_increment[_qp];
260  // and subtract off all inelastic strain increments calculated so far
261  // except the one that we're about to calculate
262  for (unsigned j_rmm = 0; j_rmm < _num_models; ++j_rmm)
263  if (i_rmm != j_rmm)
264  elastic_strain_increment -= inelastic_strain_increment[j_rmm];
265 
266  // form the trial stress, with the check for changed elasticity constants
268  _stress[_qp] =
269  _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
270  else
271  _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
272 
273  // given a trial stress (_stress[_qp]) and a strain increment (elastic_strain_increment)
274  // let the i^th model produce an admissible stress (as _stress[_qp]), and decompose
275  // the strain increment into an elastic part (elastic_strain_increment) and an
276  // inelastic part (inelastic_strain_increment[i_rmm])
278  elastic_strain_increment,
279  inelastic_strain_increment[i_rmm],
281 
282  if (i_rmm == 0)
283  {
284  stress_max = _stress[_qp];
285  stress_min = _stress[_qp];
286  }
287  else
288  {
289  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
290  {
291  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
292  {
293  if (_stress[_qp](i, j) > stress_max(i, j))
294  stress_max(i, j) = _stress[_qp](i, j);
295  else if (stress_min(i, j) > _stress[_qp](i, j))
296  stress_min(i, j) = _stress[_qp](i, j);
297  }
298  }
299  }
300  }
301 
302  // now check convergence in the stress:
303  // once the change in stress is within tolerance after each recompute material
304  // consider the stress to be converged
305  l2norm_delta_stress = (stress_max - stress_min).L2norm();
306  if (counter == 0 && l2norm_delta_stress > 0.0)
307  first_l2norm_delta_stress = l2norm_delta_stress;
308 
310  {
311  _console << "stress iteration number = " << counter << "\n"
312  << " relative l2 norm delta stress = "
313  << (0 == first_l2norm_delta_stress ? 0
314  : l2norm_delta_stress / first_l2norm_delta_stress)
315  << "\n"
316  << " stress convergence relative tolerance = " << _relative_tolerance << "\n"
317  << " absolute l2 norm delta stress = " << l2norm_delta_stress << "\n"
318  << " stress convergence absolute tolerance = " << _absolute_tolerance << std::endl;
319  }
320  ++counter;
321  } while (counter < _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
322  (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance &&
323  _num_models != 1);
324 
325  if (counter == _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
326  (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance)
327  throw MooseException("Max stress iteration hit during ComputeMultipleInelasticStress solve!");
328 
329  combined_inelastic_strain_increment.zero();
330  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
331  combined_inelastic_strain_increment +=
332  _inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
333 
334  if (_fe_problem.currentlyComputingJacobian())
336 
337  _matl_timestep_limit[_qp] = 0.0;
338  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
339  _matl_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
340 
341  if (MooseUtils::absoluteFuzzyEqual(_matl_timestep_limit[_qp], 0.0))
342  _matl_timestep_limit[_qp] = std::numeric_limits<Real>::max();
343  else
344  _matl_timestep_limit[_qp] = 1.0 / _matl_timestep_limit[_qp];
345 }
346 
347 void
349 {
353  {
355  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
356  A += _consistent_tangent_operator[i_rmm];
357  mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
358  _Jacobian_mult[_qp] = A.invSymm() * _elasticity_tensor[_qp];
359  }
360  else
361  {
362  const RankFourTensor E_inv = _elasticity_tensor[_qp].invSymm();
364  for (unsigned i_rmm = 1; i_rmm < _num_models; ++i_rmm)
365  _Jacobian_mult[_qp] = _consistent_tangent_operator[i_rmm] * E_inv * _Jacobian_mult[_qp];
366  }
367 }
368 
369 void
371  unsigned model_number,
372  RankTwoTensor & elastic_strain_increment,
373  RankTwoTensor & combined_inelastic_strain_increment)
374 {
375  for (auto model : _models)
376  model->setQp(_qp);
377 
378  elastic_strain_increment = _strain_increment[_qp];
379 
380  // If the elasticity tensor values have changed and the tensor is isotropic,
381  // use the old strain to calculate the old stress
383  _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
384  else
385  _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
386 
387  computeAdmissibleState(model_number,
388  elastic_strain_increment,
389  combined_inelastic_strain_increment,
391 
392  if (_fe_problem.currentlyComputingJacobian())
393  {
397  {
399  mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
400  _Jacobian_mult[_qp] = A.invSymm() * _elasticity_tensor[_qp];
401  }
402  else
404  }
405 
406  _matl_timestep_limit[_qp] = _models[0]->computeTimeStepLimit();
407 
408  /* propagate internal variables, etc, to this timestep for those inelastic models where
409  * "updateState" is not called */
410  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
411  if (i_rmm != model_number)
412  _models[i_rmm]->propagateQpStatefulProperties();
413 }
414 
415 void
417  RankTwoTensor & elastic_strain_increment,
418  RankTwoTensor & inelastic_strain_increment,
419  RankFourTensor & consistent_tangent_operator)
420 {
421  const bool jac = _fe_problem.currentlyComputingJacobian();
422  _models[model_number]->updateState(elastic_strain_increment,
423  inelastic_strain_increment,
424  _rotation_increment[_qp],
425  _stress[_qp],
426  _stress_old[_qp],
427  _elasticity_tensor[_qp],
428  _elastic_strain_old[_qp],
429  (jac && _tangent_computation_flag[model_number]),
430  consistent_tangent_operator);
431 
432  if (jac && !_tangent_computation_flag[model_number])
433  {
435  consistent_tangent_operator.zero();
436  else
437  consistent_tangent_operator = _elasticity_tensor[_qp];
438  }
439 }
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)
std::vector< RankFourTensor > _consistent_tangent_operator
the consistent tangent operators computed by each plastic model
ComputeMultipleInelasticStress computes the stress, the consistent tangent operator (or an approximat...
const MaterialProperty< RankTwoTensor > & _rotation_increment
Rotation increment material property.
const MaterialProperty< RankTwoTensor > & _elastic_strain_old
Strain tensors.
virtual void computeQpStress() override
Compute the stress and store it in the _stress material property for the current quadrature point...
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 MaterialProperty< RankFourTensor > & _elasticity_tensor
Elasticity tensor material property.
const std::vector< Real > _inelastic_weights
_inelastic_strain = sum_i (_inelastic_weights_i * inelastic_strain_from_model_i)
MaterialProperty< RankTwoTensor > & _stress
Stress material property.
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 ...
registerMooseObject("TensorMechanicsApp", ComputeMultipleInelasticStress)
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...
const std::string _elasticity_tensor_name
Name of the elasticity tensor material property.
std::vector< bool > _tangent_computation_flag
Flags to compute tangent during updateState call.
TangentOperatorEnum
what sort of Tangent operator to calculate
InputParameters validParams< ComputeFiniteStrainElasticStress >()
const RankFourTensor _identity_symmetric_four
Rank four symmetric identity tensor.
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
TangentCalculationMethod _tangent_calculation_method
Calculation method for the tangent modulus.
const unsigned _num_models
number of plastic models
ComputeFiniteStrainElasticStress computes the stress following elasticity theory for finite strains...
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
Old state of the stress tensor material property.
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.
Real L2norm(const RankTwoTensor &r2tensor)
const bool _perform_finite_strain_rotations
after updateQpState, rotate the stress, elastic_strain, inelastic_strain and Jacobian_mult using _rot...
MaterialProperty< RankTwoTensor > & _elastic_strain
Elastic strain material property.
bool hasGuaranteedMaterialProperty(const MaterialPropertyName &prop, Guarantee guarantee)
TangentCalculationMethod
TangentCalculationMethod is an enum that determines the calculation method for the tangent operator...
InputParameters validParams< ComputeMultipleInelasticStress >()
virtual void computeQpJacobianMult()
Using _elasticity_tensor[_qp] and the consistent tangent operators, _consistent_tangent_operator[...] computed by the inelastic models, compute _Jacobian_mult[_qp].
static unsigned int counter
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...