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 "ComputeMultipleInelasticStressBase.h"
11 :
12 : #include "StressUpdateBase.h"
13 : #include "MooseException.h"
14 : #include "DamageBase.h"
15 : #include "libmesh/int_range.h"
16 :
17 : InputParameters
18 2797 : ComputeMultipleInelasticStressBase::validParams()
19 : {
20 2797 : InputParameters params = ComputeFiniteStrainElasticStress::validParams();
21 2797 : params.addClassDescription("Compute state (stress and internal parameters such as plastic "
22 : "strains and internal parameters) using an iterative process. "
23 : "Combinations of creep models and plastic models may be used.");
24 5594 : params.addParam<unsigned int>("max_iterations",
25 5594 : 30,
26 : "Maximum number of the stress update "
27 : "iterations over the stress change after all "
28 : "update materials are called");
29 5594 : params.addParam<Real>("relative_tolerance",
30 5594 : 1e-5,
31 : "Relative convergence tolerance for the stress "
32 : "update iterations over the stress change "
33 : "after all update materials are called");
34 5594 : params.addParam<Real>("absolute_tolerance",
35 5594 : 1e-5,
36 : "Absolute convergence tolerance for the stress "
37 : "update iterations over the stress change "
38 : "after all update materials are called");
39 5594 : params.addParam<bool>(
40 : "internal_solve_full_iteration_history",
41 5594 : false,
42 : "Set to true to output stress update iteration information over the stress change");
43 5594 : params.addParam<bool>("perform_finite_strain_rotations",
44 5594 : true,
45 : "Tensors are correctly rotated in "
46 : "finite-strain simulations. For "
47 : "optimal performance you can set "
48 : "this to 'false' if you are only "
49 : "ever using small strains");
50 5594 : MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
51 5594 : params.addParam<MooseEnum>(
52 : "tangent_operator",
53 : tangent_operator,
54 : "Type of tangent operator to return. 'elastic': return the "
55 : "elasticity tensor. 'nonlinear': return the full, general consistent tangent "
56 : "operator.");
57 5594 : params.addParam<std::vector<Real>>("combined_inelastic_strain_weights",
58 : "The combined_inelastic_strain Material Property is a "
59 : "weighted sum of the model inelastic strains. This parameter "
60 : "is a vector of weights, of the same length as "
61 : "inelastic_models. Default = '1 1 ... 1'. This "
62 : "parameter is set to 1 if the number of models = 1");
63 5594 : params.addParam<bool>(
64 5594 : "cycle_models", false, "At timestep N use only inelastic model N % num_models.");
65 5594 : params.addParam<MaterialName>("damage_model", "Name of the damage model");
66 :
67 2797 : return params;
68 2797 : }
69 :
70 2095 : ComputeMultipleInelasticStressBase::ComputeMultipleInelasticStressBase(
71 2095 : const InputParameters & parameters)
72 : : ComputeFiniteStrainElasticStress(parameters),
73 2095 : _max_iterations(parameters.get<unsigned int>("max_iterations")),
74 2095 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
75 2095 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
76 4190 : _internal_solve_full_iteration_history(getParam<bool>("internal_solve_full_iteration_history")),
77 4190 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
78 4190 : _elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "elastic_strain")),
79 4190 : _strain_increment(getMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
80 2095 : _inelastic_strain(declareProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
81 2095 : _inelastic_strain_old(
82 2095 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
83 4190 : _tangent_operator_type(getParam<MooseEnum>("tangent_operator").getEnum<TangentOperatorEnum>()),
84 2095 : _tangent_calculation_method(TangentCalculationMethod::ELASTIC),
85 4190 : _cycle_models(getParam<bool>("cycle_models")),
86 2095 : _material_timestep_limit(declareProperty<Real>(_base_name + "material_timestep_limit")),
87 2095 : _identity_symmetric_four(RankFourTensor::initIdentitySymmetricFour),
88 4190 : _all_models_isotropic(true)
89 : {
90 2095 : }
91 :
92 : void
93 81736 : ComputeMultipleInelasticStressBase::initQpStatefulProperties()
94 : {
95 81736 : ComputeStressBase::initQpStatefulProperties();
96 81736 : _inelastic_strain[_qp].zero();
97 81736 : }
98 :
99 : void
100 2007 : ComputeMultipleInelasticStressBase::initialSetup()
101 : {
102 2007 : _damage_model = isParamValid("damage_model")
103 2025 : ? dynamic_cast<DamageBaseTempl<false> *>(&getMaterial("damage_model"))
104 : : nullptr;
105 :
106 2007 : _is_elasticity_tensor_guaranteed_isotropic =
107 2007 : hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC);
108 :
109 2007 : std::vector<MaterialName> models = getInelasticModelNames();
110 :
111 2007 : _num_models = models.size();
112 2007 : _tangent_computation_flag.resize(_num_models, false);
113 2007 : _consistent_tangent_operator.resize(_num_models);
114 :
115 2007 : _inelastic_weights = isParamValid("combined_inelastic_strain_weights")
116 4032 : ? getParam<std::vector<Real>>("combined_inelastic_strain_weights")
117 4014 : : std::vector<Real>(_num_models, 1.0);
118 :
119 2007 : if (_inelastic_weights.size() != _num_models)
120 0 : mooseError("ComputeMultipleInelasticStressBase: combined_inelastic_strain_weights must contain "
121 : "the same "
122 : "number of entries as inelastic_models ",
123 0 : _inelastic_weights.size(),
124 : " ",
125 0 : _num_models);
126 :
127 4006 : for (const auto i : make_range(_num_models))
128 : {
129 2000 : StressUpdateBase * rrr = dynamic_cast<StressUpdateBase *>(&getMaterialByName(models[i]));
130 :
131 2000 : if (rrr)
132 : {
133 2000 : _models.push_back(rrr);
134 2000 : _all_models_isotropic = _all_models_isotropic && rrr->isIsotropic();
135 2000 : if (rrr->requiresIsotropicTensor() && !_is_elasticity_tensor_guaranteed_isotropic)
136 2 : mooseError("Model " + models[i] +
137 : " requires an isotropic elasticity tensor, but the one supplied is not "
138 : "guaranteed isotropic");
139 : }
140 : else
141 0 : mooseError("Model " + models[i] +
142 : " is not compatible with ComputeMultipleInelasticStressBase");
143 : }
144 :
145 : // Check if tangent calculation methods are consistent. If all models have
146 : // TangentOperatorEnum::ELASTIC or tangent_operator is set by the user as elasic, then the tangent
147 : // is never calculated: J_tot = C. If PARTIAL and NONE models are present, utilize PARTIAL
148 : // formulation: J_tot = (I + J_1 + ... J_N)^-1 C. If FULL and NONE models are present, utilize
149 : // FULL formulation: J_tot = J_1 * C^-1 * J_2 * C^-1 * ... J_N * C. If PARTIAL and FULL models are
150 : // present, error out.
151 :
152 2006 : if (_tangent_operator_type != TangentOperatorEnum::elastic)
153 : {
154 : bool full_present = false;
155 : bool partial_present = false;
156 3607 : for (const auto i : make_range(_num_models))
157 : {
158 1773 : if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::FULL)
159 : {
160 : full_present = true;
161 : _tangent_computation_flag[i] = true;
162 1284 : _tangent_calculation_method = TangentCalculationMethod::FULL;
163 : }
164 489 : else if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::PARTIAL)
165 : {
166 : partial_present = true;
167 : _tangent_computation_flag[i] = true;
168 375 : _tangent_calculation_method = TangentCalculationMethod::PARTIAL;
169 : }
170 : }
171 1834 : if (full_present && partial_present)
172 0 : mooseError("In ",
173 0 : _name,
174 : ": Models that calculate the full tangent operator and the partial tangent "
175 : "operator are being combined. Either set tangent_operator to elastic, implement "
176 : "the corrent tangent formulations, or use different models.");
177 : }
178 :
179 6018 : if (isParamValid("damage_model") && !_damage_model)
180 0 : paramError("damage_model",
181 0 : "Damage Model " + _damage_model->name() +
182 : " is not compatible with ComputeMultipleInelasticStressBase");
183 :
184 : // This check prevents the hierarchy from silently skipping substepping without informing the user
185 4004 : for (const auto model_number : index_range(_models))
186 : {
187 1999 : const bool use_substep = _models[model_number]->substeppingCapabilityRequested();
188 1999 : if (use_substep && !_models[model_number]->substeppingCapabilityEnabled())
189 : {
190 1 : std::stringstream error_message;
191 : error_message << "Usage of substepping has been requested, but the inelastic model "
192 2 : << _models[model_number]->name() << " does not implement substepping yet.";
193 1 : mooseError(error_message.str());
194 0 : }
195 : }
196 2005 : }
197 :
198 : void
199 9664287 : ComputeMultipleInelasticStressBase::computeQpStress()
200 : {
201 9664287 : if (_damage_model)
202 : {
203 3504 : _undamaged_stress_old = _stress_old[_qp];
204 3504 : _damage_model->setQp(_qp);
205 3504 : _damage_model->computeUndamagedOldStress(_undamaged_stress_old);
206 : }
207 9664287 : computeQpStressIntermediateConfiguration();
208 :
209 9664229 : if (_damage_model)
210 : {
211 3504 : _damage_model->setQp(_qp);
212 3504 : _damage_model->updateDamage();
213 3504 : _damage_model->updateStressForDamage(_stress[_qp]);
214 3504 : _damage_model->finiteStrainRotation(_rotation_increment[_qp]);
215 3504 : _damage_model->updateJacobianMultForDamage(_Jacobian_mult[_qp]);
216 :
217 3504 : const Real damage_timestep_limit = _damage_model->computeTimeStepLimit();
218 3504 : if (_material_timestep_limit[_qp] > damage_timestep_limit)
219 3168 : _material_timestep_limit[_qp] = damage_timestep_limit;
220 : }
221 :
222 9664229 : if (_perform_finite_strain_rotations)
223 8795981 : finiteStrainRotation();
224 9664229 : }
225 :
226 : void
227 9692527 : ComputeMultipleInelasticStressBase::computeQpStressIntermediateConfiguration()
228 : {
229 9692527 : RankTwoTensor elastic_strain_increment;
230 9692527 : RankTwoTensor combined_inelastic_strain_increment;
231 :
232 9692527 : if (_num_models == 0)
233 : {
234 104272 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + _strain_increment[_qp];
235 :
236 : // If the elasticity tensor values have changed and the tensor is isotropic,
237 : // use the old strain to calculate the old stress
238 104272 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
239 104272 : _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + _strain_increment[_qp]);
240 : else
241 : {
242 0 : if (_damage_model)
243 0 : _stress[_qp] = _undamaged_stress_old + _elasticity_tensor[_qp] * _strain_increment[_qp];
244 : else
245 0 : _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * _strain_increment[_qp];
246 : }
247 104272 : if (_fe_problem.currentlyComputingJacobian())
248 22848 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
249 :
250 104272 : _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
251 : }
252 : else
253 : {
254 9588255 : if (_num_models == 1 || _cycle_models)
255 9167031 : updateQpStateSingleModel((_t_step - 1) % _num_models,
256 : elastic_strain_increment,
257 : combined_inelastic_strain_increment);
258 : else
259 421224 : updateQpState(elastic_strain_increment, combined_inelastic_strain_increment);
260 :
261 9588197 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + elastic_strain_increment;
262 9588197 : _inelastic_strain[_qp] = _inelastic_strain_old[_qp] + combined_inelastic_strain_increment;
263 : }
264 9692469 : }
265 :
266 : void
267 9343053 : ComputeMultipleInelasticStressBase::finiteStrainRotation(const bool force_elasticity_rotation)
268 : {
269 : _elastic_strain[_qp] =
270 9343053 : _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
271 9343053 : _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
272 : _inelastic_strain[_qp] =
273 9343053 : _rotation_increment[_qp] * _inelastic_strain[_qp] * _rotation_increment[_qp].transpose();
274 :
275 9343053 : if (force_elasticity_rotation ||
276 8824221 : !(_is_elasticity_tensor_guaranteed_isotropic &&
277 8819805 : (_all_models_isotropic ||
278 619496 : _tangent_calculation_method == TangentCalculationMethod::ELASTIC || _num_models == 0)))
279 538520 : _Jacobian_mult[_qp].rotate(_rotation_increment[_qp]);
280 9343053 : }
281 :
282 : void
283 66232 : ComputeMultipleInelasticStressBase::computeQpJacobianMult()
284 : {
285 66232 : if (_tangent_calculation_method == TangentCalculationMethod::ELASTIC)
286 63336 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
287 2896 : else if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
288 : {
289 2880 : RankFourTensor A = _identity_symmetric_four;
290 8640 : for (const auto i_rmm : make_range(_num_models))
291 5760 : A += _consistent_tangent_operator[i_rmm];
292 : mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
293 2880 : _Jacobian_mult[_qp] = A.invSymm() * _elasticity_tensor[_qp];
294 : }
295 : else
296 : {
297 16 : const RankFourTensor E_inv = _elasticity_tensor[_qp].invSymm();
298 16 : _Jacobian_mult[_qp] = _consistent_tangent_operator[0];
299 32 : for (const auto i_rmm : make_range(1u, _num_models))
300 16 : _Jacobian_mult[_qp] = _consistent_tangent_operator[i_rmm] * E_inv * _Jacobian_mult[_qp];
301 : }
302 66232 : }
303 :
304 : void
305 9167031 : ComputeMultipleInelasticStressBase::updateQpStateSingleModel(
306 : unsigned model_number,
307 : RankTwoTensor & elastic_strain_increment,
308 : RankTwoTensor & combined_inelastic_strain_increment)
309 : {
310 18341422 : for (auto model : _models)
311 9174391 : model->setQp(_qp);
312 :
313 9167031 : elastic_strain_increment = _strain_increment[_qp];
314 :
315 : // If the elasticity tensor values have changed and the tensor is isotropic,
316 : // use the old strain to calculate the old stress
317 9167031 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
318 9164087 : _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
319 : else
320 : {
321 2944 : if (_damage_model)
322 0 : _stress[_qp] = _undamaged_stress_old + _elasticity_tensor[_qp] * elastic_strain_increment;
323 : else
324 2944 : _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
325 : }
326 :
327 9167031 : computeAdmissibleState(model_number,
328 : elastic_strain_increment,
329 : combined_inelastic_strain_increment,
330 : _consistent_tangent_operator[0]);
331 :
332 9166973 : if (_fe_problem.currentlyComputingJacobian())
333 : {
334 256432 : if (_tangent_calculation_method == TangentCalculationMethod::ELASTIC)
335 56384 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
336 200048 : else if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
337 : {
338 52216 : RankFourTensor A = _identity_symmetric_four + _consistent_tangent_operator[0];
339 : mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
340 52216 : _Jacobian_mult[_qp] = A.invSymm() * _elasticity_tensor[_qp];
341 : }
342 : else
343 147832 : _Jacobian_mult[_qp] = _consistent_tangent_operator[0];
344 : }
345 :
346 9166973 : _material_timestep_limit[_qp] = _models[0]->computeTimeStepLimit();
347 :
348 : /* propagate internal variables, etc, to this timestep for those inelastic models where
349 : * "updateState" is not called */
350 18341306 : for (const auto i_rmm : make_range(_num_models))
351 9174333 : if (i_rmm != model_number)
352 7360 : _models[i_rmm]->propagateQpStatefulProperties();
353 9166973 : }
354 :
355 : void
356 14745705 : ComputeMultipleInelasticStressBase::computeAdmissibleState(
357 : unsigned model_number,
358 : RankTwoTensor & elastic_strain_increment,
359 : RankTwoTensor & inelastic_strain_increment,
360 : RankFourTensor & consistent_tangent_operator)
361 : {
362 : // Reset properties to the beginning of the time step (necessary if substepping is employed).
363 14745705 : _models[model_number]->resetIncrementalMaterialProperties();
364 :
365 14745705 : const bool jac = _fe_problem.currentlyComputingJacobian();
366 14745705 : if (_damage_model)
367 7008 : _models[model_number]->updateState(elastic_strain_increment,
368 : inelastic_strain_increment,
369 3504 : _rotation_increment[_qp],
370 3504 : _stress[_qp],
371 3504 : _undamaged_stress_old,
372 3504 : _elasticity_tensor[_qp],
373 3504 : _elastic_strain_old[_qp],
374 3504 : (jac && _tangent_computation_flag[model_number]),
375 : consistent_tangent_operator);
376 14742201 : else if (_models[model_number]->substeppingCapabilityEnabled() &&
377 16533 : (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations))
378 33066 : _models[model_number]->updateStateSubstep(elastic_strain_increment,
379 : inelastic_strain_increment,
380 16533 : _rotation_increment[_qp],
381 16533 : _stress[_qp],
382 16533 : _stress_old[_qp],
383 16533 : _elasticity_tensor[_qp],
384 16533 : _elastic_strain_old[_qp],
385 16533 : (jac && _tangent_computation_flag[model_number]),
386 : consistent_tangent_operator);
387 : else
388 29451336 : _models[model_number]->updateState(elastic_strain_increment,
389 : inelastic_strain_increment,
390 14725668 : _rotation_increment[_qp],
391 14725668 : _stress[_qp],
392 14725668 : _stress_old[_qp],
393 14725668 : _elasticity_tensor[_qp],
394 14725668 : _elastic_strain_old[_qp],
395 14725668 : (jac && _tangent_computation_flag[model_number]),
396 : consistent_tangent_operator);
397 :
398 14745647 : if (jac && !_tangent_computation_flag[model_number])
399 : {
400 779744 : if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
401 0 : consistent_tangent_operator.zero();
402 : else
403 779744 : consistent_tangent_operator = _elasticity_tensor[_qp];
404 : }
405 14745647 : }
|
