Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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("SolidMechanicsApp", ADComputeMultipleInelasticStress);
14 :
15 : InputParameters
16 1217 : ADComputeMultipleInelasticStress::validParams()
17 : {
18 1217 : InputParameters params = ADComputeFiniteStrainElasticStress::validParams();
19 1217 : 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 2434 : params.addParam<unsigned int>("max_iterations",
23 2434 : 30,
24 : "Maximum number of the stress update "
25 : "iterations over the stress change after all "
26 : "update materials are called");
27 2434 : params.addParam<Real>("relative_tolerance",
28 2434 : 1e-5,
29 : "Relative convergence tolerance for the stress "
30 : "update iterations over the stress change "
31 : "after all update materials are called");
32 2434 : params.addParam<Real>("absolute_tolerance",
33 2434 : 1e-5,
34 : "Absolute convergence tolerance for the stress "
35 : "update iterations over the stress change "
36 : "after all update materials are called");
37 2434 : params.addParam<bool>(
38 : "internal_solve_full_iteration_history",
39 2434 : false,
40 : "Set to true to output stress update iteration information over the stress change");
41 2434 : params.addParam<bool>("perform_finite_strain_rotations",
42 2434 : 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 2434 : 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 2434 : 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 2434 : params.addParam<bool>(
59 2434 : "cycle_models", false, "At time step N use only inelastic model N % num_models.");
60 2434 : params.addParam<MaterialName>("damage_model", "Name of the damage model");
61 :
62 1217 : return params;
63 0 : }
64 :
65 913 : ADComputeMultipleInelasticStress::ADComputeMultipleInelasticStress(
66 913 : const InputParameters & parameters)
67 : : ADComputeFiniteStrainElasticStress(parameters),
68 913 : _max_iterations(parameters.get<unsigned int>("max_iterations")),
69 913 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
70 913 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
71 1826 : _internal_solve_full_iteration_history(getParam<bool>("internal_solve_full_iteration_history")),
72 1826 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
73 913 : _inelastic_strain(declareADProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
74 913 : _inelastic_strain_old(
75 913 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
76 1826 : _num_models(getParam<std::vector<MaterialName>>("inelastic_models").size()),
77 1826 : _inelastic_weights(isParamValid("combined_inelastic_strain_weights")
78 913 : ? getParam<std::vector<Real>>("combined_inelastic_strain_weights")
79 913 : : std::vector<Real>(_num_models, true)),
80 1826 : _cycle_models(getParam<bool>("cycle_models")),
81 913 : _material_timestep_limit(declareProperty<Real>(_base_name + "material_timestep_limit")),
82 1826 : _is_elasticity_tensor_guaranteed_isotropic(false)
83 : {
84 913 : 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 913 : }
91 :
92 : void
93 21536 : ADComputeMultipleInelasticStress::initQpStatefulProperties()
94 : {
95 21536 : ADComputeFiniteStrainElasticStress::initQpStatefulProperties();
96 21536 : _inelastic_strain[_qp].zero();
97 21536 : }
98 :
99 : void
100 905 : ADComputeMultipleInelasticStress::initialSetup()
101 : {
102 905 : _damage_model = isParamValid("damage_model")
103 918 : ? dynamic_cast<DamageBaseTempl<true> *>(&getMaterial("damage_model"))
104 : : nullptr;
105 :
106 2715 : 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 904 : _is_elasticity_tensor_guaranteed_isotropic =
112 904 : this->hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC);
113 :
114 1808 : std::vector<MaterialName> models = getParam<std::vector<MaterialName>>("inelastic_models");
115 :
116 2077 : for (unsigned int i = 0; i < _num_models; ++i)
117 : {
118 : ADStressUpdateBase * rrr =
119 1173 : dynamic_cast<ADStressUpdateBase *>(&this->getMaterialByName(models[i]));
120 :
121 1173 : if (rrr)
122 : {
123 1173 : _models.push_back(rrr);
124 1173 : 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 904 : }
133 :
134 : void
135 4542436 : ADComputeMultipleInelasticStress::computeQpStress()
136 : {
137 4542436 : if (_damage_model)
138 : {
139 2800 : _undamaged_stress_old = _stress_old[_qp];
140 : _damage_model->setQp(_qp);
141 2800 : _damage_model->computeUndamagedOldStress(_undamaged_stress_old);
142 : }
143 :
144 4542436 : computeQpStressIntermediateConfiguration();
145 :
146 4542428 : if (_damage_model)
147 : {
148 2800 : _damage_model->setQp(_qp);
149 2800 : _damage_model->updateDamage();
150 2800 : _damage_model->updateStressForDamage(_stress[_qp]);
151 2800 : _damage_model->finiteStrainRotation(_rotation_increment[_qp]);
152 :
153 2800 : const Real damage_timestep_limit = _damage_model->computeTimeStepLimit();
154 2800 : if (_material_timestep_limit[_qp] > damage_timestep_limit)
155 2536 : _material_timestep_limit[_qp] = damage_timestep_limit;
156 : }
157 :
158 4542428 : if (_perform_finite_strain_rotations)
159 4542428 : finiteStrainRotation();
160 4542428 : }
161 :
162 : void
163 4575244 : ADComputeMultipleInelasticStress::computeQpStressIntermediateConfiguration()
164 : {
165 4575244 : ADRankTwoTensor elastic_strain_increment;
166 4575244 : ADRankTwoTensor combined_inelastic_strain_increment;
167 :
168 4575244 : if (_num_models == 0)
169 : {
170 32808 : _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 32808 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
175 32808 : _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 4542436 : if (_num_models == 1 || _cycle_models)
196 3720884 : updateQpStateSingleModel((_t_step - 1) % _num_models,
197 : elastic_strain_increment,
198 : combined_inelastic_strain_increment);
199 : else
200 821552 : updateQpState(elastic_strain_increment, combined_inelastic_strain_increment);
201 :
202 4542428 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + elastic_strain_increment;
203 4542428 : _inelastic_strain[_qp] = _inelastic_strain_old[_qp] + combined_inelastic_strain_increment;
204 : }
205 4575236 : }
206 :
207 : void
208 5006248 : ADComputeMultipleInelasticStress::finiteStrainRotation()
209 : {
210 5006248 : _elastic_strain[_qp] =
211 5006248 : _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
212 5006248 : _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
213 5006248 : _inelastic_strain[_qp] =
214 5006248 : _rotation_increment[_qp] * _inelastic_strain[_qp] * _rotation_increment[_qp].transpose();
215 5006248 : }
216 :
217 : void
218 821552 : ADComputeMultipleInelasticStress::updateQpState(
219 : ADRankTwoTensor & elastic_strain_increment,
220 : ADRankTwoTensor & combined_inelastic_strain_increment)
221 : {
222 821552 : 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 821552 : unsigned int counter = 0;
231 :
232 : std::vector<ADRankTwoTensor> inelastic_strain_increment;
233 821552 : inelastic_strain_increment.resize(_num_models);
234 :
235 2527056 : for (unsigned i_rmm = 0; i_rmm < _models.size(); ++i_rmm)
236 1705504 : inelastic_strain_increment[i_rmm].zero();
237 :
238 821552 : ADRankTwoTensor stress_max, stress_min;
239 :
240 : do
241 : {
242 19677885 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
243 : {
244 13139390 : _models[i_rmm]->setQp(_qp);
245 :
246 : // initially assume the strain is completely elastic
247 13139390 : 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 39854970 : for (unsigned j_rmm = 0; j_rmm < _num_models; ++j_rmm)
251 26715580 : if (i_rmm != j_rmm)
252 13576190 : elastic_strain_increment -= inelastic_strain_increment[j_rmm];
253 :
254 : // form the trial stress, with the check for changed elasticity constants
255 13139390 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
256 13086686 : _stress[_qp] =
257 26173372 : _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
258 : else
259 : {
260 52704 : if (_damage_model)
261 0 : paramError("damage_model",
262 : "Damage models cannot be used with inelastic models and elastic anisotropic "
263 : "behavior");
264 :
265 52704 : ADRankFourTensor elasticity_tensor_rotated = _elasticity_tensor[_qp];
266 52704 : elasticity_tensor_rotated.rotate(_rotation_total_old[_qp]);
267 :
268 52704 : _stress[_qp] =
269 52704 : elasticity_tensor_rotated * (_elastic_strain_old[_qp] + elastic_strain_increment);
270 :
271 : // Update current total rotation matrix to be used in next step
272 52704 : _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 13139390 : computeAdmissibleState(i_rmm, elastic_strain_increment, inelastic_strain_increment[i_rmm]);
279 :
280 13139390 : if (i_rmm == 0)
281 : {
282 6538495 : stress_max = _stress[_qp];
283 6538495 : stress_min = _stress[_qp];
284 : }
285 : else
286 : {
287 26403580 : for (const auto i : make_range(Moose::dim))
288 : {
289 79210740 : for (const auto j : make_range(Moose::dim))
290 : {
291 59408055 : if (_stress[_qp](i, j) > stress_max(i, j))
292 18762788 : stress_max(i, j) = _stress[_qp](i, j);
293 40645267 : else if (stress_min(i, j) > _stress[_qp](i, j))
294 12360091 : 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 6538495 : l2norm_delta_stress = MetaPhysicL::raw_value((stress_max - stress_min).L2norm());
304 6538495 : if (counter == 0 && l2norm_delta_stress > 0.0)
305 : first_l2norm_delta_stress = l2norm_delta_stress;
306 :
307 6538495 : 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 6538495 : ++counter;
319 :
320 6537091 : } while (counter < _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
321 12716821 : (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance &&
322 5716943 : _num_models != 1);
323 :
324 821552 : 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 821552 : combined_inelastic_strain_increment.zero();
330 2527056 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
331 : combined_inelastic_strain_increment +=
332 3411008 : _inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
333 :
334 821552 : _material_timestep_limit[_qp] = 0.0;
335 2527056 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
336 1705504 : _material_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
337 :
338 821552 : if (MooseUtils::absoluteFuzzyEqual(_material_timestep_limit[_qp], 0.0))
339 31204 : _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
340 : else
341 790348 : _material_timestep_limit[_qp] = 1.0 / _material_timestep_limit[_qp];
342 821552 : }
343 :
344 : void
345 3720884 : ADComputeMultipleInelasticStress::updateQpStateSingleModel(
346 : unsigned model_number,
347 : ADRankTwoTensor & elastic_strain_increment,
348 : ADRankTwoTensor & combined_inelastic_strain_increment)
349 : {
350 7721344 : for (auto model : _models)
351 4000460 : model->setQp(_qp);
352 :
353 3720884 : 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 3720884 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
358 3720884 : _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 3720884 : computeAdmissibleState(
377 : model_number, elastic_strain_increment, combined_inelastic_strain_increment);
378 :
379 3720876 : _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 7721328 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
384 4000452 : if (i_rmm != model_number)
385 279576 : _models[i_rmm]->propagateQpStatefulProperties();
386 3720876 : }
387 :
388 : void
389 16860274 : 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 16860274 : _models[model_number]->resetIncrementalMaterialProperties();
396 :
397 16860274 : if (_damage_model)
398 2800 : _models[model_number]->updateState(elastic_strain_increment,
399 : inelastic_strain_increment,
400 2800 : _rotation_increment[_qp],
401 2800 : _stress[_qp],
402 2800 : _undamaged_stress_old,
403 2800 : _elasticity_tensor[_qp],
404 2800 : _elastic_strain_old[_qp]);
405 16857474 : else if (_models[model_number]->substeppingCapabilityEnabled() &&
406 3138 : (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations))
407 : {
408 3138 : _models[model_number]->updateStateSubstep(elastic_strain_increment,
409 : inelastic_strain_increment,
410 3138 : _rotation_increment[_qp],
411 3138 : _stress[_qp],
412 3138 : _stress_old[_qp],
413 3138 : _elasticity_tensor[_qp],
414 3138 : _elastic_strain_old[_qp]);
415 : }
416 : else
417 16854336 : _models[model_number]->updateState(elastic_strain_increment,
418 : inelastic_strain_increment,
419 16854336 : _rotation_increment[_qp],
420 16854336 : _stress[_qp],
421 16854336 : _stress_old[_qp],
422 16854336 : _elasticity_tensor[_qp],
423 16854336 : _elastic_strain_old[_qp]);
424 16860266 : }
|
