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 1962 : ADComputeMultipleInelasticStress::validParams()
17 : {
18 1962 : InputParameters params = ADComputeFiniteStrainElasticStress::validParams();
19 1962 : 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 3924 : params.addParam<unsigned int>("max_iterations",
23 3924 : 30,
24 : "Maximum number of the stress update "
25 : "iterations over the stress change after all "
26 : "update materials are called");
27 3924 : params.addParam<Real>("relative_tolerance",
28 3924 : 1e-5,
29 : "Relative convergence tolerance for the stress "
30 : "update iterations over the stress change "
31 : "after all update materials are called");
32 3924 : params.addParam<Real>("absolute_tolerance",
33 3924 : 1e-5,
34 : "Absolute convergence tolerance for the stress "
35 : "update iterations over the stress change "
36 : "after all update materials are called");
37 3924 : params.addParam<bool>(
38 : "internal_solve_full_iteration_history",
39 3924 : false,
40 : "Set to true to output stress update iteration information over the stress change");
41 3924 : params.addParam<bool>("perform_finite_strain_rotations",
42 3924 : 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 3924 : 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 3924 : 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 3924 : params.addParam<bool>(
59 3924 : "cycle_models", false, "At time step N use only inelastic model N % num_models.");
60 3924 : params.addParam<MaterialName>("damage_model", "Name of the damage model");
61 :
62 1962 : return params;
63 0 : }
64 :
65 1472 : ADComputeMultipleInelasticStress::ADComputeMultipleInelasticStress(
66 1472 : const InputParameters & parameters)
67 : : ADComputeFiniteStrainElasticStress(parameters),
68 1472 : _max_iterations(parameters.get<unsigned int>("max_iterations")),
69 1472 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
70 1472 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
71 2944 : _internal_solve_full_iteration_history(getParam<bool>("internal_solve_full_iteration_history")),
72 2944 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
73 1472 : _inelastic_strain(declareADProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
74 1472 : _inelastic_strain_old(
75 1472 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
76 2944 : _num_models(getParam<std::vector<MaterialName>>("inelastic_models").size()),
77 2944 : _inelastic_weights(isParamValid("combined_inelastic_strain_weights")
78 1472 : ? getParam<std::vector<Real>>("combined_inelastic_strain_weights")
79 1472 : : std::vector<Real>(_num_models, true)),
80 2944 : _cycle_models(getParam<bool>("cycle_models")),
81 1472 : _material_timestep_limit(declareProperty<Real>(_base_name + "material_timestep_limit")),
82 2944 : _is_elasticity_tensor_guaranteed_isotropic(false)
83 : {
84 1472 : 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 1472 : }
91 :
92 : void
93 35824 : ADComputeMultipleInelasticStress::initQpStatefulProperties()
94 : {
95 35824 : ADComputeFiniteStrainElasticStress::initQpStatefulProperties();
96 35824 : _inelastic_strain[_qp].zero();
97 35824 : }
98 :
99 : void
100 1456 : ADComputeMultipleInelasticStress::initialSetup()
101 : {
102 1456 : _damage_model = isParamValid("damage_model")
103 1476 : ? dynamic_cast<DamageBaseTempl<true> *>(&getMaterial("damage_model"))
104 : : nullptr;
105 :
106 4368 : if (isParamValid("damage_model") && !_damage_model)
107 2 : paramError("damage_model",
108 2 : "Damage Model " + getMaterial("damage_model").name() +
109 : " is not compatible with ADComputeMultipleInelasticStress");
110 :
111 1454 : _is_elasticity_tensor_guaranteed_isotropic =
112 1454 : this->hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC);
113 :
114 2908 : std::vector<MaterialName> models = getParam<std::vector<MaterialName>>("inelastic_models");
115 :
116 3386 : for (unsigned int i = 0; i < _num_models; ++i)
117 : {
118 : ADStressUpdateBase * rrr =
119 1932 : dynamic_cast<ADStressUpdateBase *>(&this->getMaterialByName(models[i]));
120 :
121 1932 : if (rrr)
122 : {
123 1932 : _models.push_back(rrr);
124 1932 : 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 1454 : }
133 :
134 : void
135 8390502 : ADComputeMultipleInelasticStress::computeQpStress()
136 : {
137 8390502 : if (_damage_model)
138 : {
139 3616 : _undamaged_stress_old = _stress_old[_qp];
140 : _damage_model->setQp(_qp);
141 3616 : _damage_model->computeUndamagedOldStress(_undamaged_stress_old);
142 : }
143 :
144 8390502 : computeQpStressIntermediateConfiguration();
145 :
146 8390432 : if (_damage_model)
147 : {
148 3616 : _damage_model->setQp(_qp);
149 3616 : _damage_model->updateDamage();
150 3616 : _damage_model->updateStressForDamage(_stress[_qp]);
151 3616 : _damage_model->finiteStrainRotation(_rotation_increment[_qp]);
152 :
153 3616 : const Real damage_timestep_limit = _damage_model->computeTimeStepLimit();
154 3616 : if (_material_timestep_limit[_qp] > damage_timestep_limit)
155 3264 : _material_timestep_limit[_qp] = damage_timestep_limit;
156 : }
157 :
158 8390432 : if (_perform_finite_strain_rotations)
159 8390432 : finiteStrainRotation();
160 8390432 : }
161 :
162 : void
163 8433190 : ADComputeMultipleInelasticStress::computeQpStressIntermediateConfiguration()
164 : {
165 8433190 : ADRankTwoTensor elastic_strain_increment;
166 8433190 : ADRankTwoTensor combined_inelastic_strain_increment;
167 :
168 8433190 : if (_num_models == 0)
169 : {
170 42688 : _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 42688 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
175 42688 : _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 8390502 : if (_num_models == 1 || _cycle_models)
196 6999094 : updateQpStateSingleModel((_t_step - 1) % _num_models,
197 : elastic_strain_increment,
198 : combined_inelastic_strain_increment);
199 : else
200 1391408 : updateQpState(elastic_strain_increment, combined_inelastic_strain_increment);
201 :
202 8390432 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + elastic_strain_increment;
203 8390432 : _inelastic_strain[_qp] = _inelastic_strain_old[_qp] + combined_inelastic_strain_increment;
204 : }
205 8433120 : }
206 :
207 : void
208 9124064 : ADComputeMultipleInelasticStress::finiteStrainRotation()
209 : {
210 9124064 : _elastic_strain[_qp] =
211 9124064 : _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
212 9124064 : _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
213 9124064 : _inelastic_strain[_qp] =
214 9124064 : _rotation_increment[_qp] * _inelastic_strain[_qp] * _rotation_increment[_qp].transpose();
215 9124064 : }
216 :
217 : void
218 1391408 : ADComputeMultipleInelasticStress::updateQpState(
219 : ADRankTwoTensor & elastic_strain_increment,
220 : ADRankTwoTensor & combined_inelastic_strain_increment)
221 : {
222 1391408 : 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 1391408 : unsigned int counter = 0;
231 :
232 : std::vector<ADRankTwoTensor> inelastic_strain_increment;
233 1391408 : inelastic_strain_increment.resize(_num_models);
234 :
235 4443664 : for (unsigned i_rmm = 0; i_rmm < _models.size(); ++i_rmm)
236 3052256 : inelastic_strain_increment[i_rmm].zero();
237 :
238 1391408 : ADRankTwoTensor stress_max, stress_min;
239 :
240 : do
241 : {
242 31793572 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
243 : {
244 21285528 : _models[i_rmm]->setQp(_qp);
245 :
246 : // initially assume the strain is completely elastic
247 21285528 : 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 65742664 : for (unsigned j_rmm = 0; j_rmm < _num_models; ++j_rmm)
251 44457136 : if (i_rmm != j_rmm)
252 23171608 : elastic_strain_increment -= inelastic_strain_increment[j_rmm];
253 :
254 : // form the trial stress, with the check for changed elasticity constants
255 21285528 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
256 21212760 : _stress[_qp] =
257 42425520 : _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
258 : else
259 : {
260 72768 : if (_damage_model)
261 0 : paramError("damage_model",
262 : "Damage models cannot be used with inelastic models and elastic anisotropic "
263 : "behavior");
264 :
265 72768 : ADRankFourTensor elasticity_tensor_rotated = _elasticity_tensor[_qp];
266 72768 : elasticity_tensor_rotated.rotate(_rotation_total_old[_qp]);
267 :
268 72768 : _stress[_qp] =
269 72768 : elasticity_tensor_rotated * (_elastic_strain_old[_qp] + elastic_strain_increment);
270 :
271 : // Update current total rotation matrix to be used in next step
272 72768 : _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 21285528 : computeAdmissibleState(i_rmm, elastic_strain_increment, inelastic_strain_increment[i_rmm]);
279 :
280 21285528 : if (i_rmm == 0)
281 : {
282 10508044 : stress_max = _stress[_qp];
283 10508044 : stress_min = _stress[_qp];
284 : }
285 : else
286 : {
287 43109936 : for (const auto i : make_range(Moose::dim))
288 : {
289 129329808 : for (const auto j : make_range(Moose::dim))
290 : {
291 96997356 : if (_stress[_qp](i, j) > stress_max(i, j))
292 29836014 : stress_max(i, j) = _stress[_qp](i, j);
293 67161342 : else if (stress_min(i, j) > _stress[_qp](i, j))
294 19756312 : 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 10508044 : l2norm_delta_stress = MetaPhysicL::raw_value((stress_max - stress_min).L2norm());
304 10508044 : if (counter == 0 && l2norm_delta_stress > 0.0)
305 : first_l2norm_delta_stress = l2norm_delta_stress;
306 :
307 10508044 : 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 10508044 : ++counter;
319 :
320 10506004 : } while (counter < _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
321 20364506 : (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance &&
322 9116636 : _num_models != 1);
323 :
324 1391408 : 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 1391408 : combined_inelastic_strain_increment.zero();
330 4443664 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
331 : combined_inelastic_strain_increment +=
332 6104512 : _inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
333 :
334 1391408 : _material_timestep_limit[_qp] = 0.0;
335 4443664 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
336 3052256 : _material_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
337 :
338 1391408 : if (MooseUtils::absoluteFuzzyEqual(_material_timestep_limit[_qp], 0.0))
339 85884 : _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
340 : else
341 1305524 : _material_timestep_limit[_qp] = 1.0 / _material_timestep_limit[_qp];
342 1391408 : }
343 :
344 : void
345 6999094 : ADComputeMultipleInelasticStress::updateQpStateSingleModel(
346 : unsigned model_number,
347 : ADRankTwoTensor & elastic_strain_increment,
348 : ADRankTwoTensor & combined_inelastic_strain_increment)
349 : {
350 14502476 : for (auto model : _models)
351 7503382 : model->setQp(_qp);
352 :
353 6999094 : 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 6999094 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
358 6999094 : _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 6999094 : computeAdmissibleState(
377 : model_number, elastic_strain_increment, combined_inelastic_strain_increment);
378 :
379 6999024 : _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 14502336 : for (unsigned i_rmm = 0; i_rmm < _num_models; ++i_rmm)
384 7503312 : if (i_rmm != model_number)
385 504288 : _models[i_rmm]->propagateQpStatefulProperties();
386 6999024 : }
387 :
388 : void
389 28284622 : 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 28284622 : _models[model_number]->resetIncrementalMaterialProperties();
396 :
397 28284622 : if (_damage_model)
398 3616 : _models[model_number]->updateState(elastic_strain_increment,
399 : inelastic_strain_increment,
400 3616 : _rotation_increment[_qp],
401 3616 : _stress[_qp],
402 3616 : _undamaged_stress_old,
403 3616 : _elasticity_tensor[_qp],
404 3616 : _elastic_strain_old[_qp]);
405 28281006 : else if (_models[model_number]->substeppingCapabilityEnabled() &&
406 21306 : (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations))
407 : {
408 21306 : _models[model_number]->updateStateSubstep(elastic_strain_increment,
409 : inelastic_strain_increment,
410 21306 : _rotation_increment[_qp],
411 21306 : _stress[_qp],
412 21306 : _stress_old[_qp],
413 21306 : _elasticity_tensor[_qp],
414 21306 : _elastic_strain_old[_qp]);
415 : }
416 : else
417 28259700 : _models[model_number]->updateState(elastic_strain_increment,
418 : inelastic_strain_increment,
419 28259700 : _rotation_increment[_qp],
420 28259700 : _stress[_qp],
421 28259700 : _stress_old[_qp],
422 28259700 : _elasticity_tensor[_qp],
423 28259700 : _elastic_strain_old[_qp]);
424 28284552 : }
|
