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 "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 6630 : ComputeMultipleInelasticStressBase::validParams()
19 : {
20 6630 : InputParameters params = ComputeFiniteStrainElasticStress::validParams();
21 6630 : 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 13260 : params.addParam<unsigned int>("max_iterations",
25 13260 : 30,
26 : "Maximum number of the stress update "
27 : "iterations over the stress change after all "
28 : "update materials are called");
29 13260 : params.addParam<Real>("relative_tolerance",
30 13260 : 1e-5,
31 : "Relative convergence tolerance for the stress "
32 : "update iterations over the stress change "
33 : "after all update materials are called");
34 13260 : params.addParam<Real>("absolute_tolerance",
35 13260 : 1e-5,
36 : "Absolute convergence tolerance for the stress "
37 : "update iterations over the stress change "
38 : "after all update materials are called");
39 13260 : params.addParam<bool>(
40 : "internal_solve_full_iteration_history",
41 13260 : false,
42 : "Set to true to output stress update iteration information over the stress change");
43 13260 : params.addParam<bool>("perform_finite_strain_rotations",
44 13260 : 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 13260 : MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
51 13260 : 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 13260 : 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 13260 : params.addParam<bool>(
64 13260 : "cycle_models", false, "At timestep N use only inelastic model N % num_models.");
65 13260 : params.addParam<MaterialName>("damage_model", "Name of the damage model");
66 :
67 6630 : return params;
68 6630 : }
69 :
70 4967 : ComputeMultipleInelasticStressBase::ComputeMultipleInelasticStressBase(
71 4967 : const InputParameters & parameters)
72 : : ComputeFiniteStrainElasticStress(parameters),
73 4967 : _max_iterations(parameters.get<unsigned int>("max_iterations")),
74 4967 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
75 4967 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
76 9934 : _internal_solve_full_iteration_history(getParam<bool>("internal_solve_full_iteration_history")),
77 9934 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
78 9934 : _elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "elastic_strain")),
79 9934 : _strain_increment(getMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
80 4967 : _inelastic_strain(declareProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
81 4967 : _inelastic_strain_old(
82 4967 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
83 9934 : _tangent_operator_type(getParam<MooseEnum>("tangent_operator").getEnum<TangentOperatorEnum>()),
84 4967 : _tangent_calculation_method(TangentCalculationMethod::ELASTIC),
85 9934 : _cycle_models(getParam<bool>("cycle_models")),
86 4967 : _material_timestep_limit(declareProperty<Real>(_base_name + "material_timestep_limit")),
87 4967 : _identity_symmetric_four(RankFourTensor::initIdentitySymmetricFour),
88 9934 : _all_models_isotropic(true)
89 : {
90 4967 : }
91 :
92 : void
93 113016 : ComputeMultipleInelasticStressBase::initQpStatefulProperties()
94 : {
95 113016 : ComputeStressBase::initQpStatefulProperties();
96 113016 : _inelastic_strain[_qp].zero();
97 113016 : }
98 :
99 : void
100 4752 : ComputeMultipleInelasticStressBase::initialSetup()
101 : {
102 4752 : _damage_model = isParamValid("damage_model")
103 4794 : ? dynamic_cast<DamageBaseTempl<false> *>(&getMaterial("damage_model"))
104 : : nullptr;
105 :
106 4752 : _is_elasticity_tensor_guaranteed_isotropic =
107 4752 : hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC);
108 :
109 4752 : std::vector<MaterialName> models = getInelasticModelNames();
110 :
111 4752 : _num_models = models.size();
112 4752 : _tangent_computation_flag.resize(_num_models, false);
113 4752 : _consistent_tangent_operator.resize(_num_models);
114 :
115 4752 : _inelastic_weights = isParamValid("combined_inelastic_strain_weights")
116 14298 : ? getParam<std::vector<Real>>("combined_inelastic_strain_weights")
117 9504 : : std::vector<Real>(_num_models, 1.0);
118 :
119 4752 : 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 9536 : for (const auto i : make_range(_num_models))
128 : {
129 4786 : StressUpdateBase * rrr = dynamic_cast<StressUpdateBase *>(&getMaterialByName(models[i]));
130 :
131 4786 : if (rrr)
132 : {
133 4786 : _models.push_back(rrr);
134 4786 : _all_models_isotropic = _all_models_isotropic && rrr->isIsotropic();
135 4786 : if (rrr->requiresIsotropicTensor() && !_is_elasticity_tensor_guaranteed_isotropic)
136 4 : 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 4750 : if (_tangent_operator_type != TangentOperatorEnum::elastic)
153 : {
154 : bool full_present = false;
155 : bool partial_present = false;
156 8018 : for (const auto i : make_range(_num_models))
157 : {
158 3942 : if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::FULL)
159 : {
160 : full_present = true;
161 : _tangent_computation_flag[i] = true;
162 2802 : _tangent_calculation_method = TangentCalculationMethod::FULL;
163 : }
164 1140 : else if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::PARTIAL)
165 : {
166 : partial_present = true;
167 : _tangent_computation_flag[i] = true;
168 870 : _tangent_calculation_method = TangentCalculationMethod::PARTIAL;
169 : }
170 : }
171 4076 : 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 14250 : 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 9532 : for (const auto model_number : index_range(_models))
186 : {
187 4784 : const bool use_substep = _models[model_number]->substeppingCapabilityRequested();
188 4784 : if (use_substep && !_models[model_number]->substeppingCapabilityEnabled())
189 : {
190 2 : 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 2 : mooseError(error_message.str());
194 0 : }
195 : }
196 4748 : }
197 :
198 : void
199 23806720 : ComputeMultipleInelasticStressBase::computeQpStress()
200 : {
201 23806720 : if (_damage_model)
202 : {
203 8032 : _undamaged_stress_old = _stress_old[_qp];
204 : _damage_model->setQp(_qp);
205 8032 : _damage_model->computeUndamagedOldStress(_undamaged_stress_old);
206 : }
207 23806720 : computeQpStressIntermediateConfiguration();
208 :
209 23806604 : if (_damage_model)
210 : {
211 8032 : _damage_model->setQp(_qp);
212 8032 : _damage_model->updateDamage();
213 8032 : _damage_model->updateStressForDamage(_stress[_qp]);
214 8032 : _damage_model->finiteStrainRotation(_rotation_increment[_qp]);
215 8032 : _damage_model->updateJacobianMultForDamage(_Jacobian_mult[_qp]);
216 :
217 8032 : const Real damage_timestep_limit = _damage_model->computeTimeStepLimit();
218 8032 : if (_material_timestep_limit[_qp] > damage_timestep_limit)
219 7272 : _material_timestep_limit[_qp] = damage_timestep_limit;
220 : }
221 :
222 23806604 : if (_perform_finite_strain_rotations)
223 22047764 : finiteStrainRotation();
224 23806604 : }
225 :
226 : void
227 23864288 : ComputeMultipleInelasticStressBase::computeQpStressIntermediateConfiguration()
228 : {
229 23864288 : RankTwoTensor elastic_strain_increment;
230 23864288 : RankTwoTensor combined_inelastic_strain_increment;
231 :
232 23864288 : if (_num_models == 0)
233 : {
234 60040 : _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 60040 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
239 60040 : _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 60040 : if (_fe_problem.currentlyComputingJacobian())
248 10320 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
249 :
250 60040 : _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
251 : }
252 : else
253 : {
254 23804248 : if (_num_models == 1 || _cycle_models)
255 21627904 : updateQpStateSingleModel((_t_step - 1) % _num_models,
256 : elastic_strain_increment,
257 : combined_inelastic_strain_increment);
258 : else
259 2176344 : updateQpState(elastic_strain_increment, combined_inelastic_strain_increment);
260 :
261 23804132 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + elastic_strain_increment;
262 23804132 : _inelastic_strain[_qp] = _inelastic_strain_old[_qp] + combined_inelastic_strain_increment;
263 : }
264 23864172 : }
265 :
266 : void
267 23308480 : ComputeMultipleInelasticStressBase::finiteStrainRotation(const bool force_elasticity_rotation)
268 : {
269 23308480 : _elastic_strain[_qp] =
270 23308480 : _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
271 23308480 : _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
272 23308480 : _inelastic_strain[_qp] =
273 23308480 : _rotation_increment[_qp] * _inelastic_strain[_qp] * _rotation_increment[_qp].transpose();
274 :
275 23308480 : if (force_elasticity_rotation ||
276 22105332 : !(_is_elasticity_tensor_guaranteed_isotropic &&
277 22096692 : (_all_models_isotropic ||
278 1489148 : _tangent_calculation_method == TangentCalculationMethod::ELASTIC || _num_models == 0)))
279 1241292 : _Jacobian_mult[_qp].rotate(_rotation_increment[_qp]);
280 23308480 : }
281 :
282 : void
283 219968 : ComputeMultipleInelasticStressBase::computeQpJacobianMult()
284 : {
285 219968 : if (_tangent_calculation_method == TangentCalculationMethod::ELASTIC)
286 212448 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
287 7520 : else if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
288 : {
289 7488 : RankFourTensor A = _identity_symmetric_four;
290 22464 : for (const auto i_rmm : make_range(_num_models))
291 14976 : A += _consistent_tangent_operator[i_rmm];
292 : mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
293 7488 : _Jacobian_mult[_qp] = A.invSymm() * _elasticity_tensor[_qp];
294 : }
295 : else
296 : {
297 32 : const RankFourTensor E_inv = _elasticity_tensor[_qp].invSymm();
298 32 : _Jacobian_mult[_qp] = _consistent_tangent_operator[0];
299 64 : for (const auto i_rmm : make_range(1u, _num_models))
300 32 : _Jacobian_mult[_qp] = _consistent_tangent_operator[i_rmm] * E_inv * _Jacobian_mult[_qp];
301 : }
302 219968 : }
303 :
304 : void
305 21627904 : ComputeMultipleInelasticStressBase::updateQpStateSingleModel(
306 : unsigned model_number,
307 : RankTwoTensor & elastic_strain_increment,
308 : RankTwoTensor & combined_inelastic_strain_increment)
309 : {
310 43273856 : for (auto model : _models)
311 21645952 : model->setQp(_qp);
312 :
313 21627904 : 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 21627904 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
318 21622144 : _stress[_qp] = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
319 : else
320 : {
321 5760 : if (_damage_model)
322 0 : _stress[_qp] = _undamaged_stress_old + _elasticity_tensor[_qp] * elastic_strain_increment;
323 : else
324 5760 : _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
325 : }
326 :
327 21627904 : computeAdmissibleState(model_number,
328 : elastic_strain_increment,
329 : combined_inelastic_strain_increment,
330 : _consistent_tangent_operator[0]);
331 :
332 21627788 : if (_fe_problem.currentlyComputingJacobian())
333 : {
334 1167644 : if (_tangent_calculation_method == TangentCalculationMethod::ELASTIC)
335 663532 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
336 504112 : else if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
337 : {
338 128272 : RankFourTensor A = _identity_symmetric_four + _consistent_tangent_operator[0];
339 : mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
340 128272 : _Jacobian_mult[_qp] = A.invSymm() * _elasticity_tensor[_qp];
341 : }
342 : else
343 375840 : _Jacobian_mult[_qp] = _consistent_tangent_operator[0];
344 : }
345 :
346 21627788 : _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 43273624 : for (const auto i_rmm : make_range(_num_models))
351 21645836 : if (i_rmm != model_number)
352 18048 : _models[i_rmm]->propagateQpStatefulProperties();
353 21627788 : }
354 :
355 : void
356 37800910 : 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 37800910 : _models[model_number]->resetIncrementalMaterialProperties();
364 :
365 37800910 : const bool jac = _fe_problem.currentlyComputingJacobian();
366 37800910 : if (_damage_model)
367 16064 : _models[model_number]->updateState(elastic_strain_increment,
368 : inelastic_strain_increment,
369 8032 : _rotation_increment[_qp],
370 8032 : _stress[_qp],
371 8032 : _undamaged_stress_old,
372 8032 : _elasticity_tensor[_qp],
373 8032 : _elastic_strain_old[_qp],
374 8032 : (jac && _tangent_computation_flag[model_number]),
375 : consistent_tangent_operator);
376 37792878 : else if (_models[model_number]->substeppingCapabilityEnabled() &&
377 30522 : (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations))
378 61044 : _models[model_number]->updateStateSubstep(elastic_strain_increment,
379 : inelastic_strain_increment,
380 30522 : _rotation_increment[_qp],
381 30522 : _stress[_qp],
382 30522 : _stress_old[_qp],
383 30522 : _elasticity_tensor[_qp],
384 30522 : _elastic_strain_old[_qp],
385 30522 : (jac && _tangent_computation_flag[model_number]),
386 : consistent_tangent_operator);
387 : else
388 75524712 : _models[model_number]->updateState(elastic_strain_increment,
389 : inelastic_strain_increment,
390 37762356 : _rotation_increment[_qp],
391 37762356 : _stress[_qp],
392 37762356 : _stress_old[_qp],
393 37762356 : _elasticity_tensor[_qp],
394 37762356 : _elastic_strain_old[_qp],
395 37762356 : (jac && _tangent_computation_flag[model_number]),
396 : consistent_tangent_operator);
397 :
398 37800794 : if (jac && !_tangent_computation_flag[model_number])
399 : {
400 2583916 : if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
401 0 : consistent_tangent_operator.zero();
402 : else
403 2583916 : consistent_tangent_operator = _elasticity_tensor[_qp];
404 : }
405 37800794 : }
|
