Line data Source code
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 :
10 : #include "ADComputeMultipleInelasticStress.h"
11 : #include "MooseException.h"
12 :
13 : registerMooseObject("TensorMechanicsApp", ADComputeMultipleInelasticStress);
14 :
15 : InputParameters
16 957 : ADComputeMultipleInelasticStress::validParams()
17 : {
18 957 : InputParameters params = ADComputeFiniteStrainElasticStress::validParams();
19 957 : 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 1914 : params.addParam<unsigned int>("max_iterations",
23 1914 : 30,
24 : "Maximum number of the stress update "
25 : "iterations over the stress change after all "
26 : "update materials are called");
27 1914 : params.addParam<Real>("relative_tolerance",
28 1914 : 1e-5,
29 : "Relative convergence tolerance for the stress "
30 : "update iterations over the stress change "
31 : "after all update materials are called");
32 1914 : params.addParam<Real>("absolute_tolerance",
33 1914 : 1e-5,
34 : "Absolute convergence tolerance for the stress "
35 : "update iterations over the stress change "
36 : "after all update materials are called");
37 1914 : params.addParam<bool>(
38 : "internal_solve_full_iteration_history",
39 1914 : false,
40 : "Set to true to output stress update iteration information over the stress change");
41 1914 : params.addParam<bool>("perform_finite_strain_rotations",
42 1914 : 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 1914 : 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 1914 : params.addParam<std::vector<Real>>("combined_inelastic_strain_weights",
53 : "The combined_inelastic_strain Material Property is a "
54 : "weighted sum of the model inelastic strains. This parameter "
55 : "is a vector of weights, of the same length as "
56 : "inelastic_models. Default = '1 1 ... 1'. This "
57 : "parameter is set to 1 if the number of models = 1");
58 1914 : params.addParam<bool>(
59 1914 : "cycle_models", false, "At time step N use only inelastic model N % num_models.");
60 1914 : params.addParam<MaterialName>("damage_model", "Name of the damage model");
61 :
62 957 : return params;
63 0 : }
64 :
65 718 : ADComputeMultipleInelasticStress::ADComputeMultipleInelasticStress(
66 718 : const InputParameters & parameters)
67 : : ADComputeFiniteStrainElasticStress(parameters),
68 718 : _max_iterations(parameters.get<unsigned int>("max_iterations")),
69 718 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
70 718 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
71 1436 : _internal_solve_full_iteration_history(getParam<bool>("internal_solve_full_iteration_history")),
72 1436 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
73 718 : _inelastic_strain(declareADProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
74 718 : _inelastic_strain_old(
75 718 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
76 1436 : _num_models(getParam<std::vector<MaterialName>>("inelastic_models").size()),
77 1436 : _inelastic_weights(isParamValid("combined_inelastic_strain_weights")
78 718 : ? getParam<std::vector<Real>>("combined_inelastic_strain_weights")
79 718 : : std::vector<Real>(_num_models, true)),
80 1436 : _cycle_models(getParam<bool>("cycle_models")),
81 718 : _material_timestep_limit(declareProperty<Real>(_base_name + "material_timestep_limit")),
82 1436 : _is_elasticity_tensor_guaranteed_isotropic(false)
83 : {
84 718 : if (_inelastic_weights.size() != _num_models)
85 0 : paramError("combined_inelastic_strain_weights",
86 : "must contain the same number of entries as inelastic_models ",
87 : _inelastic_weights.size(),
88 : " vs. ",
89 : _num_models);
90 718 : }
91 :
92 : void
93 17816 : ADComputeMultipleInelasticStress::initQpStatefulProperties()
94 : {
95 17816 : ADComputeFiniteStrainElasticStress::initQpStatefulProperties();
96 17816 : _inelastic_strain[_qp].zero();
97 17816 : }
98 :
99 : void
100 710 : ADComputeMultipleInelasticStress::initialSetup()
101 : {
102 710 : _damage_model = isParamValid("damage_model")
103 720 : ? dynamic_cast<DamageBaseTempl<true> *>(&getMaterial("damage_model"))
104 : : nullptr;
105 :
106 2130 : if (isParamValid("damage_model") && !_damage_model)
107 1 : paramError("damage_model",
108 1 : "Damage Model " + getMaterial("damage_model").name() +
109 : " is not compatible with ADComputeMultipleInelasticStress");
110 :
111 709 : _is_elasticity_tensor_guaranteed_isotropic =
112 709 : this->hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC);
113 :
114 1418 : std::vector<MaterialName> models = getParam<std::vector<MaterialName>>("inelastic_models");
115 :
116 1657 : for (unsigned int i = 0; i < _num_models; ++i)
117 : {
118 : ADStressUpdateBase * rrr =
119 948 : dynamic_cast<ADStressUpdateBase *>(&this->getMaterialByName(models[i]));
120 :
121 948 : if (rrr)
122 : {
123 948 : _models.push_back(rrr);
124 948 : if (rrr->requiresIsotropicTensor() && !_is_elasticity_tensor_guaranteed_isotropic)
125 0 : mooseError("Model " + models[i] +
126 : " requires an isotropic elasticity tensor, but the one supplied is not "
127 : "guaranteed isotropic");
128 : }
129 : else
130 0 : mooseError("Model " + models[i] + " is not compatible with ADComputeMultipleInelasticStress");
131 : }
132 709 : }
133 :
134 : void
135 4285214 : ADComputeMultipleInelasticStress::computeQpStress()
136 : {
137 4285214 : if (_damage_model)
138 : {
139 1984 : _undamaged_stress_old = _stress_old[_qp];
140 1984 : _damage_model->setQp(_qp);
141 1984 : _damage_model->computeUndamagedOldStress(_undamaged_stress_old);
142 : }
143 :
144 4285214 : computeQpStressIntermediateConfiguration();
145 :
146 4285179 : if (_damage_model)
147 : {
148 1984 : _damage_model->setQp(_qp);
149 1984 : _damage_model->updateDamage();
150 1984 : _damage_model->updateStressForDamage(_stress[_qp]);
151 1984 : _damage_model->finiteStrainRotation(_rotation_increment[_qp]);
152 :
153 1984 : const Real damage_timestep_limit = _damage_model->computeTimeStepLimit();
154 1984 : if (_material_timestep_limit[_qp] > damage_timestep_limit)
155 1792 : _material_timestep_limit[_qp] = damage_timestep_limit;
156 : }
157 :
158 4285179 : if (_perform_finite_strain_rotations)
159 4285179 : finiteStrainRotation();
160 4285179 : }
161 :
162 : void
163 4311902 : ADComputeMultipleInelasticStress::computeQpStressIntermediateConfiguration()
164 : {
165 4311902 : ADRankTwoTensor elastic_strain_increment;
166 4311902 : ADRankTwoTensor combined_inelastic_strain_increment;
167 :
168 4311902 : if (_num_models == 0)
169 : {
170 26688 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + _strain_increment[_qp];
171 :
172 : // If the elasticity tensor values have changed and the tensor is isotropic,
173 : // use the old strain to calculate the old stress
174 26688 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
175 26688 : _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + _strain_increment[_qp]);
176 : else
177 : {
178 0 : if (_damage_model)
179 0 : paramError(
180 : "damage_model",
181 : "Damage models cannot be used with inelastic models and elastic anisotropic behavior");
182 :
183 0 : ADRankFourTensor elasticity_tensor_rotated = _elasticity_tensor[_qp];
184 0 : elasticity_tensor_rotated.rotate(_rotation_total_old[_qp]);
185 :
186 0 : _stress[_qp] =
187 0 : elasticity_tensor_rotated * (_elastic_strain_old[_qp] + _strain_increment[_qp]);
188 :
189 : // Update current total rotation matrix to be used in next step
190 0 : _rotation_total[_qp] = _rotation_increment[_qp] * _rotation_total_old[_qp];
191 : }
192 : }
193 : else
194 : {
195 4285214 : if (_num_models == 1 || _cycle_models)
196 3582194 : updateQpStateSingleModel((_t_step - 1) % _num_models,
197 : elastic_strain_increment,
198 : combined_inelastic_strain_increment);
199 : else
200 703020 : updateQpState(elastic_strain_increment, combined_inelastic_strain_increment);
201 :
202 4285179 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + elastic_strain_increment;
203 4285179 : _inelastic_strain[_qp] = _inelastic_strain_old[_qp] + combined_inelastic_strain_increment;
204 : }
205 4311867 : }
206 :
207 : void
208 4700651 : ADComputeMultipleInelasticStress::finiteStrainRotation()
209 : {
210 4700651 : _elastic_strain[_qp] =
211 4700651 : _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
212 4700651 : _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
213 4700651 : _inelastic_strain[_qp] =
214 4700651 : _rotation_increment[_qp] * _inelastic_strain[_qp] * _rotation_increment[_qp].transpose();
215 4700651 : }
216 :
217 : void
218 703020 : ADComputeMultipleInelasticStress::updateQpState(
219 : ADRankTwoTensor & elastic_strain_increment,
220 : ADRankTwoTensor & combined_inelastic_strain_increment)
221 : {
222 703020 : if (_internal_solve_full_iteration_history == true)
223 : {
224 0 : _console << std::endl
225 0 : << "iteration output for ADComputeMultipleInelasticStress solve:"
226 0 : << " time=" << _t << " int_pt=" << _qp << std::endl;
227 : }
228 : Real l2norm_delta_stress;
229 : Real first_l2norm_delta_stress = 1.0;
230 703020 : unsigned int counter = 0;
231 :
232 : std::vector<ADRankTwoTensor> inelastic_strain_increment;
233 703020 : inelastic_strain_increment.resize(_num_models);
234 :
235 2251780 : for (unsigned i_rmm = 0; i_rmm < _models.size(); ++i_rmm)
236 1548760 : inelastic_strain_increment[i_rmm].zero();
237 :
238 703020 : ADRankTwoTensor stress_max, stress_min;
239 :
240 : do
241 : {
242 15966562 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
243 : {
244 10691948 : _models[i_rmm]->setQp(_qp);
245 :
246 : // initially assume the strain is completely elastic
247 10691948 : elastic_strain_increment = _strain_increment[_qp];
248 : // and subtract off all inelastic strain increments calculated so far
249 : // except the one that we're about to calculate
250 33074884 : for (unsigned j_rmm = 0; j_rmm < _num_models; ++j_rmm)
251 22382936 : if (i_rmm != j_rmm)
252 11690988 : elastic_strain_increment -= inelastic_strain_increment[j_rmm];
253 :
254 : // form the trial stress, with the check for changed elasticity constants
255 10691948 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
256 10651196 : _stress[_qp] =
257 10651196 : _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
258 : else
259 : {
260 40752 : if (_damage_model)
261 0 : paramError("damage_model",
262 : "Damage models cannot be used with inelastic models and elastic anisotropic "
263 : "behavior");
264 :
265 40752 : ADRankFourTensor elasticity_tensor_rotated = _elasticity_tensor[_qp];
266 40752 : elasticity_tensor_rotated.rotate(_rotation_total_old[_qp]);
267 :
268 81504 : _stress[_qp] =
269 40752 : elasticity_tensor_rotated * (_elastic_strain_old[_qp] + elastic_strain_increment);
270 :
271 : // Update current total rotation matrix to be used in next step
272 40752 : _rotation_total[_qp] = _rotation_increment[_qp] * _rotation_total_old[_qp];
273 : }
274 : // given a trial stress (_stress[_qp]) and a strain increment (elastic_strain_increment)
275 : // let the i^th model produce an admissible stress (as _stress[_qp]), and decompose
276 : // the strain increment into an elastic part (elastic_strain_increment) and an
277 : // inelastic part (inelastic_strain_increment[i_rmm])
278 10691948 : computeAdmissibleState(i_rmm, elastic_strain_increment, inelastic_strain_increment[i_rmm]);
279 :
280 10691948 : if (i_rmm == 0)
281 : {
282 5274614 : stress_max = _stress[_qp];
283 5274614 : stress_min = _stress[_qp];
284 : }
285 : else
286 : {
287 21669336 : for (const auto i : make_range(Moose::dim))
288 : {
289 65008008 : for (const auto j : make_range(Moose::dim))
290 : {
291 48756006 : if (_stress[_qp](i, j) > stress_max(i, j))
292 14968741 : stress_max(i, j) = _stress[_qp](i, j);
293 33787265 : else if (stress_min(i, j) > _stress[_qp](i, j))
294 9915826 : stress_min(i, j) = _stress[_qp](i, j);
295 : }
296 : }
297 : }
298 : }
299 :
300 : // now check convergence in the stress:
301 : // once the change in stress is within tolerance after each recompute material
302 : // consider the stress to be converged
303 5274614 : l2norm_delta_stress = MetaPhysicL::raw_value((stress_max - stress_min).L2norm());
304 5274614 : if (counter == 0 && l2norm_delta_stress > 0.0)
305 : first_l2norm_delta_stress = l2norm_delta_stress;
306 :
307 5274614 : if (_internal_solve_full_iteration_history == true)
308 : {
309 0 : _console << "stress iteration number = " << counter << "\n"
310 0 : << " relative l2 norm delta stress = "
311 0 : << (0 == first_l2norm_delta_stress ? 0
312 0 : : l2norm_delta_stress / first_l2norm_delta_stress)
313 0 : << "\n"
314 0 : << " stress convergence relative tolerance = " << _relative_tolerance << "\n"
315 0 : << " absolute l2 norm delta stress = " << l2norm_delta_stress << "\n"
316 0 : << " stress convergence absolute tolerance = " << _absolute_tolerance << std::endl;
317 : }
318 5274614 : ++counter;
319 :
320 5273594 : } while (counter < _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
321 10218409 : (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance &&
322 4571594 : _num_models != 1);
323 :
324 703020 : if (counter == _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
325 0 : (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance)
326 0 : mooseException(
327 : "In ", _name, ": Max stress iteration hit during ADComputeMultipleInelasticStress solve!");
328 :
329 703020 : combined_inelastic_strain_increment.zero();
330 2251780 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
331 : combined_inelastic_strain_increment +=
332 3097520 : _inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
333 :
334 703020 : _material_timestep_limit[_qp] = 0.0;
335 2251780 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
336 1548760 : _material_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
337 :
338 703020 : if (MooseUtils::absoluteFuzzyEqual(_material_timestep_limit[_qp], 0.0))
339 46050 : _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
340 : else
341 656970 : _material_timestep_limit[_qp] = 1.0 / _material_timestep_limit[_qp];
342 703020 : }
343 :
344 : void
345 3582194 : ADComputeMultipleInelasticStress::updateQpStateSingleModel(
346 : unsigned model_number,
347 : ADRankTwoTensor & elastic_strain_increment,
348 : ADRankTwoTensor & combined_inelastic_strain_increment)
349 : {
350 7417972 : for (auto model : _models)
351 3835778 : model->setQp(_qp);
352 :
353 3582194 : elastic_strain_increment = _strain_increment[_qp];
354 :
355 : // If the elasticity tensor values have changed and the tensor is isotropic,
356 : // use the old strain to calculate the old stress
357 3582194 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
358 3582194 : _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
359 : else
360 : {
361 0 : if (_damage_model)
362 0 : paramError(
363 : "damage_model",
364 : "Damage models cannot be used with inelastic models and elastic anisotropic behavior");
365 :
366 0 : ADRankFourTensor elasticity_tensor_rotated = _elasticity_tensor[_qp];
367 0 : elasticity_tensor_rotated.rotate(_rotation_total_old[_qp]);
368 :
369 0 : _stress[_qp] =
370 0 : elasticity_tensor_rotated * (_elastic_strain_old[_qp] + elastic_strain_increment);
371 :
372 : // Update current total rotation matrix to be used in next step
373 0 : _rotation_total[_qp] = _rotation_increment[_qp] * _rotation_total_old[_qp];
374 : }
375 :
376 3582194 : computeAdmissibleState(
377 : model_number, elastic_strain_increment, combined_inelastic_strain_increment);
378 :
379 3582159 : _material_timestep_limit[_qp] = _models[0]->computeTimeStepLimit();
380 :
381 : /* propagate internal variables, etc, to this timestep for those inelastic models where
382 : * "updateState" is not called */
383 7417902 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
384 3835743 : if (i_rmm != model_number)
385 253584 : _models[i_rmm]->propagateQpStatefulProperties();
386 3582159 : }
387 :
388 : void
389 14274142 : ADComputeMultipleInelasticStress::computeAdmissibleState(
390 : unsigned model_number,
391 : ADRankTwoTensor & elastic_strain_increment,
392 : ADRankTwoTensor & inelastic_strain_increment)
393 : {
394 : // Properly update material properties (necessary if substepping is employed).
395 14274142 : _models[model_number]->resetIncrementalMaterialProperties();
396 :
397 14274142 : if (_damage_model)
398 1984 : _models[model_number]->updateState(elastic_strain_increment,
399 : inelastic_strain_increment,
400 1984 : _rotation_increment[_qp],
401 1984 : _stress[_qp],
402 1984 : _undamaged_stress_old,
403 1984 : _elasticity_tensor[_qp],
404 1984 : _elastic_strain_old[_qp]);
405 14272158 : else if (_models[model_number]->substeppingCapabilityEnabled() &&
406 13085 : (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations))
407 : {
408 13085 : _models[model_number]->updateStateSubstep(elastic_strain_increment,
409 : inelastic_strain_increment,
410 13085 : _rotation_increment[_qp],
411 13085 : _stress[_qp],
412 13085 : _stress_old[_qp],
413 13085 : _elasticity_tensor[_qp],
414 13085 : _elastic_strain_old[_qp]);
415 : }
416 : else
417 14259073 : _models[model_number]->updateState(elastic_strain_increment,
418 : inelastic_strain_increment,
419 14259073 : _rotation_increment[_qp],
420 14259073 : _stress[_qp],
421 14259073 : _stress_old[_qp],
422 14259073 : _elasticity_tensor[_qp],
423 14259073 : _elastic_strain_old[_qp]);
424 14274107 : }
|
