www.mooseframework.org
ADComputeMultipleInelasticStress.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 "ADStressUpdateBase.h"
12 #include "MooseException.h"
13 
15 
19  params.addClassDescription("Compute state (stress and internal parameters such as plastic "
20  "strains and internal parameters) using an iterative process. "
21  "Combinations of creep models and plastic models may be used.");
22  params.addParam<unsigned int>("max_iterations",
23  30,
24  "Maximum number of the stress update "
25  "iterations over the stress change after all "
26  "update materials are called");
27  params.addParam<Real>("relative_tolerance",
28  1e-5,
29  "Relative convergence tolerance for the stress "
30  "update iterations over the stress change "
31  "after all update materials are called");
32  params.addParam<Real>("absolute_tolerance",
33  1e-5,
34  "Absolute convergence tolerance for the stress "
35  "update iterations over the stress change "
36  "after all update materials are called");
37  params.addParam<bool>(
38  "internal_solve_full_iteration_history",
39  false,
40  "Set to true to output stress update iteration information over the stress change");
41  params.addParam<bool>("perform_finite_strain_rotations",
42  true,
43  "Tensors are correctly rotated in "
44  "finite-strain simulations. For "
45  "optimal performance you can set "
46  "this to 'false' if you are only "
47  "ever using small strains");
48  params.addRequiredParam<std::vector<MaterialName>>(
49  "inelastic_models",
50  "The material objects to use to calculate stress and inelastic strains. "
51  "Note: specify creep models first and plasticity models second.");
52  params.addParam<std::vector<Real>>(
53  "combined_inelastic_strain_weights",
54  "The combined_inelastic_strain Material Property is a "
55  "weighted sum of the model inelastic strains. This parameter "
56  "is a vector of weights, of the same length as "
57  "inelastic_models. Default = '1 1 ... 1'. This "
58  "parameter is set to 1 if the number of models = 1");
59  params.addParam<bool>("cycle_models",
60  false,
61  "At time step N use only inelastic model N % num_models."););
62 
63 template <ComputeStage compute_stage>
65  const InputParameters & parameters)
66  : ADComputeFiniteStrainElasticStress<compute_stage>(parameters),
67  _max_iterations(parameters.get<unsigned int>("max_iterations")),
68  _relative_tolerance(parameters.get<Real>("relative_tolerance")),
69  _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
70  _internal_solve_full_iteration_history(
71  adGetParam<bool>("internal_solve_full_iteration_history")),
72  _perform_finite_strain_rotations(adGetParam<bool>("perform_finite_strain_rotations")),
73  _inelastic_strain(adDeclareADProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
74  _inelastic_strain_old(
75  adGetMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
76  _num_models(adGetParam<std::vector<MaterialName>>("inelastic_models").size()),
77  _inelastic_weights(isParamValid("combined_inelastic_strain_weights")
78  ? adGetParam<std::vector<Real>>("combined_inelastic_strain_weights")
79  : std::vector<Real>(_num_models, true)),
80  _cycle_models(adGetParam<bool>("cycle_models")),
81  _matl_timestep_limit(adDeclareProperty<Real>("matl_timestep_limit"))
82 {
83  if (_inelastic_weights.size() != _num_models)
84  paramError("combined_inelastic_strain_weights",
85  "must contain the same number of entries as inelastic_models ",
86  _inelastic_weights.size(),
87  " vs. ",
88  _num_models);
89 }
90 
91 template <ComputeStage compute_stage>
92 void
94 {
96  _inelastic_strain[_qp].zero();
97 }
98 
99 template <ComputeStage compute_stage>
100 void
102 {
103  _is_elasticity_tensor_guaranteed_isotropic =
104  this->hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC);
105 
106  std::vector<MaterialName> models = adGetParam<std::vector<MaterialName>>("inelastic_models");
107 
108  for (unsigned int i = 0; i < _num_models; ++i)
109  {
111  &this->template getMaterialByName<compute_stage>(models[i]));
112 
113  if (rrr)
114  {
115  _models.push_back(rrr);
116  if (rrr->requiresIsotropicTensor() && !_is_elasticity_tensor_guaranteed_isotropic)
117  mooseError("Model " + models[i] +
118  " requires an isotropic elasticity tensor, but the one supplied is not "
119  "guaranteed isotropic");
120  }
121  else
122  mooseError("Model " + models[i] + " is not compatible with ADComputeMultipleInelasticStress");
123  }
124 }
125 
126 template <ComputeStage compute_stage>
127 void
129 {
130  computeQpStressIntermediateConfiguration();
131  if (_perform_finite_strain_rotations)
132  finiteStrainRotation();
133 }
134 
135 template <ComputeStage compute_stage>
136 void
138 {
139  ADRankTwoTensor elastic_strain_increment;
140  ADRankTwoTensor combined_inelastic_strain_increment;
141 
142  if (_num_models == 0)
143  {
144  _elastic_strain[_qp] = _elastic_strain_old[_qp] + _strain_increment[_qp];
145 
146  // If the elasticity tensor values have changed and the tensor is isotropic,
147  // use the old strain to calculate the old stress
148  if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
149  _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + _strain_increment[_qp]);
150  else
151  _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * _strain_increment[_qp];
152  }
153  else
154  {
155  if (_num_models == 1 || _cycle_models)
156  updateQpStateSingleModel((_t_step - 1) % _num_models,
157  elastic_strain_increment,
158  combined_inelastic_strain_increment);
159  else
160  updateQpState(elastic_strain_increment, combined_inelastic_strain_increment);
161 
162  _elastic_strain[_qp] = _elastic_strain_old[_qp] + elastic_strain_increment;
163  _inelastic_strain[_qp] = _inelastic_strain_old[_qp] + combined_inelastic_strain_increment;
164  }
165 }
166 
167 template <ComputeStage compute_stage>
168 void
170 {
171  _elastic_strain[_qp] =
172  _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
173  _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
174  _inelastic_strain[_qp] =
175  _rotation_increment[_qp] * _inelastic_strain[_qp] * _rotation_increment[_qp].transpose();
176 }
177 
178 template <ComputeStage compute_stage>
179 void
181  ADRankTwoTensor & elastic_strain_increment,
182  ADRankTwoTensor & combined_inelastic_strain_increment)
183 {
184  if (_internal_solve_full_iteration_history == true)
185  {
186  _console << std::endl
187  << "iteration output for ADComputeMultipleInelasticStress solve:"
188  << " time=" << _t << " int_pt=" << _qp << std::endl;
189  }
190  Real l2norm_delta_stress;
191  Real first_l2norm_delta_stress = 1.0;
192  unsigned int counter = 0;
193 
194  std::vector<ADRankTwoTensor> inelastic_strain_increment;
195  inelastic_strain_increment.resize(_num_models);
196 
197  for (unsigned i_rmm = 0; i_rmm < _models.size(); ++i_rmm)
198  inelastic_strain_increment[i_rmm].zero();
199 
200  ADRankTwoTensor stress_max, stress_min;
201 
202  do
203  {
204  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
205  {
206  _models[i_rmm]->setQp(_qp);
207 
208  // initially assume the strain is completely elastic
209  elastic_strain_increment = _strain_increment[_qp];
210  // and subtract off all inelastic strain increments calculated so far
211  // except the one that we're about to calculate
212  for (unsigned j_rmm = 0; j_rmm < _num_models; ++j_rmm)
213  if (i_rmm != j_rmm)
214  elastic_strain_increment -= inelastic_strain_increment[j_rmm];
215 
216  // form the trial stress, with the check for changed elasticity constants
217  if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
218  _stress[_qp] =
219  _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
220  else
221  _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
222 
223  // given a trial stress (_stress[_qp]) and a strain increment (elastic_strain_increment)
224  // let the i^th model produce an admissible stress (as _stress[_qp]), and decompose
225  // the strain increment into an elastic part (elastic_strain_increment) and an
226  // inelastic part (inelastic_strain_increment[i_rmm])
227  computeAdmissibleState(i_rmm, elastic_strain_increment, inelastic_strain_increment[i_rmm]);
228 
229  if (i_rmm == 0)
230  {
231  stress_max = _stress[_qp];
232  stress_min = _stress[_qp];
233  }
234  else
235  {
236  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
237  {
238  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
239  {
240  if (_stress[_qp](i, j) > stress_max(i, j))
241  stress_max(i, j) = _stress[_qp](i, j);
242  else if (stress_min(i, j) > _stress[_qp](i, j))
243  stress_min(i, j) = _stress[_qp](i, j);
244  }
245  }
246  }
247  }
248 
249  // now check convergence in the stress:
250  // once the change in stress is within tolerance after each recompute material
251  // consider the stress to be converged
252  l2norm_delta_stress = MetaPhysicL::raw_value((stress_max - stress_min).L2norm());
253  if (counter == 0 && l2norm_delta_stress > 0.0)
254  first_l2norm_delta_stress = l2norm_delta_stress;
255 
256  if (_internal_solve_full_iteration_history == true)
257  {
258  _console << "stress iteration number = " << counter << "\n"
259  << " relative l2 norm delta stress = "
260  << (0 == first_l2norm_delta_stress ? 0
261  : l2norm_delta_stress / first_l2norm_delta_stress)
262  << "\n"
263  << " stress convergence relative tolerance = " << _relative_tolerance << "\n"
264  << " absolute l2 norm delta stress = " << l2norm_delta_stress << "\n"
265  << " stress convergence absolute tolerance = " << _absolute_tolerance << std::endl;
266  }
267  ++counter;
268  } while (counter < _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
269  (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance &&
270  _num_models != 1);
271 
272  if (counter == _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
273  (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance)
274  throw MooseException("Max stress iteration hit during ADComputeMultipleInelasticStress solve!");
275 
276  combined_inelastic_strain_increment.zero();
277  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
278  combined_inelastic_strain_increment +=
279  _inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
280 
281  _matl_timestep_limit[_qp] = 0.0;
282  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
283  _matl_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
284 
285  if (MooseUtils::absoluteFuzzyEqual(_matl_timestep_limit[_qp], 0.0))
286  _matl_timestep_limit[_qp] = std::numeric_limits<Real>::max();
287  else
288  _matl_timestep_limit[_qp] = 1.0 / _matl_timestep_limit[_qp];
289 }
290 
291 template <ComputeStage compute_stage>
292 void
294  unsigned model_number,
295  ADRankTwoTensor & elastic_strain_increment,
296  ADRankTwoTensor & combined_inelastic_strain_increment)
297 {
298  for (auto model : _models)
299  model->setQp(_qp);
300 
301  elastic_strain_increment = _strain_increment[_qp];
302 
303  // If the elasticity tensor values have changed and the tensor is isotropic,
304  // use the old strain to calculate the old stress
305  if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
306  _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
307  else
308  _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
309 
310  computeAdmissibleState(
311  model_number, elastic_strain_increment, combined_inelastic_strain_increment);
312 
313  _matl_timestep_limit[_qp] = _models[0]->computeTimeStepLimit();
314 
315  /* propagate internal variables, etc, to this timestep for those inelastic models where
316  * "updateState" is not called */
317  for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
318  if (i_rmm != model_number)
319  _models[i_rmm]->propagateQpStatefulProperties();
320 }
321 
322 template <ComputeStage compute_stage>
323 void
325  unsigned model_number,
326  ADRankTwoTensor & elastic_strain_increment,
327  ADRankTwoTensor & inelastic_strain_increment)
328 {
329  _models[model_number]->updateState(elastic_strain_increment,
330  inelastic_strain_increment,
331  _rotation_increment[_qp],
332  _stress[_qp],
333  _stress_old[_qp],
334  _elasticity_tensor[_qp],
335  _elastic_strain_old[_qp]);
336 }
virtual void computeAdmissibleState(unsigned model_number, ADRankTwoTensor &elastic_strain_increment, ADRankTwoTensor &inelastic_strain_increment)
Given a trial stress (_stress[_qp]) and a strain increment (elastic_strain_increment) let the model_n...
const unsigned _num_models
number of plastic models
registerADMooseObject("TensorMechanicsApp", ADComputeMultipleInelasticStress)
defineADValidParams(ADComputeMultipleInelasticStress, ADComputeFiniteStrainElasticStress, params.addClassDescription("Compute state (stress and internal parameters such as plastic " "strains and internal parameters) using an iterative process. " "Combinations of creep models and plastic models may be used.");params.addParam< unsigned int >("max_iterations", 30, "Maximum number of the stress update " "iterations over the stress change after all " "update materials are called");params.addParam< Real >("relative_tolerance", 1e-5, "Relative convergence tolerance for the stress " "update iterations over the stress change " "after all update materials are called");params.addParam< Real >("absolute_tolerance", 1e-5, "Absolute convergence tolerance for the stress " "update iterations over the stress change " "after all update materials are called");params.addParam< bool >("internal_solve_full_iteration_history", false, "Set to true to output stress update iteration information over the stress change");params.addParam< bool >("perform_finite_strain_rotations", true, "Tensors are correctly rotated in " "finite-strain simulations. For " "optimal performance you can set " "this to 'false' if you are only " "ever using small strains");params.addRequiredParam< std::vector< MaterialName >>("inelastic_models", "The material objects to use to calculate stress and inelastic strains. " "Note: specify creep models first and plasticity models second.");params.addParam< std::vector< Real >>("combined_inelastic_strain_weights", "The combined_inelastic_strain Material Property is a " "weighted sum of the model inelastic strains. This parameter " "is a vector of weights, of the same length as " "inelastic_models. Default = '1 1 ... 1'. This " "parameter is set to 1 if the number of models = 1");params.addParam< bool >("cycle_models", false, "At time step N use only inelastic model N % num_models.");)
ADComputeMultipleInelasticStress computes the stress and a decomposition of the strain into elastic a...
ADStressUpdateBase is a material that is not called by MOOSE because of the compute=false flag set in...
virtual void updateQpState(ADRankTwoTensor &elastic_strain_increment, ADRankTwoTensor &combined_inelastic_strain_increment)
Given the _strain_increment[_qp], iterate over all of the user-specified recompute materials in order...
virtual void initQpStatefulProperties() override
virtual void finiteStrainRotation()
Rotate _elastic_strain, _stress, and _inelastic_strain to the new configuration.
virtual void computeQpStressIntermediateConfiguration()
Compute the stress for the current QP, but do not rotate tensors from the intermediate configuration ...
ADComputeMultipleInelasticStress(const InputParameters &parameters)
const std::vector< Real > _inelastic_weights
_inelastic_strain = sum_i (_inelastic_weights_i * inelastic_strain_from_model_i)
Real L2norm(const RankTwoTensor &r2tensor)
static unsigned int counter
virtual void updateQpStateSingleModel(unsigned model_number, ADRankTwoTensor &elastic_strain_increment, ADRankTwoTensor &combined_inelastic_strain_increment)
An optimised version of updateQpState that gets used when the number of plastic models is unity...
virtual bool requiresIsotropicTensor()=0
Does the model require the elasticity tensor to be isotropic?
ADComputeFiniteStrainElasticStress computes the stress following elasticity theory for finite strains...