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 5930 : ComputeMultipleInelasticStressBase::validParams()
19 : {
20 5930 : InputParameters params = ComputeFiniteStrainElasticStress::validParams();
21 5930 : 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 11860 : params.addParam<unsigned int>("max_iterations",
25 11860 : 30,
26 : "Maximum number of the stress update "
27 : "iterations over the stress change after all "
28 : "update materials are called");
29 11860 : params.addParam<Real>("relative_tolerance",
30 11860 : 1e-5,
31 : "Relative convergence tolerance for the stress "
32 : "update iterations over the stress change "
33 : "after all update materials are called");
34 11860 : params.addParam<Real>("absolute_tolerance",
35 11860 : 1e-5,
36 : "Absolute convergence tolerance for the stress "
37 : "update iterations over the stress change "
38 : "after all update materials are called");
39 11860 : params.addParam<bool>(
40 : "internal_solve_full_iteration_history",
41 11860 : false,
42 : "Set to true to output stress update iteration information over the stress change");
43 11860 : params.addParam<bool>("perform_finite_strain_rotations",
44 11860 : 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 11860 : MooseEnum tangent_operator("elastic nonlinear", "nonlinear");
51 11860 : 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 11860 : 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 11860 : params.addParam<bool>(
64 11860 : "cycle_models", false, "At timestep N use only inelastic model N % num_models.");
65 11860 : params.addParam<MaterialName>("damage_model", "Name of the damage model");
66 :
67 5930 : return params;
68 5930 : }
69 :
70 4442 : ComputeMultipleInelasticStressBase::ComputeMultipleInelasticStressBase(
71 4442 : const InputParameters & parameters)
72 : : ComputeFiniteStrainElasticStress(parameters),
73 4442 : _max_iterations(parameters.get<unsigned int>("max_iterations")),
74 4442 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
75 4442 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
76 8884 : _internal_solve_full_iteration_history(getParam<bool>("internal_solve_full_iteration_history")),
77 8884 : _perform_finite_strain_rotations(getParam<bool>("perform_finite_strain_rotations")),
78 8884 : _elastic_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "elastic_strain")),
79 8884 : _strain_increment(getMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
80 4442 : _inelastic_strain(declareProperty<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
81 4442 : _inelastic_strain_old(
82 4442 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_inelastic_strain")),
83 8884 : _tangent_operator_type(getParam<MooseEnum>("tangent_operator").getEnum<TangentOperatorEnum>()),
84 4442 : _tangent_calculation_method(TangentCalculationMethod::ELASTIC),
85 8884 : _cycle_models(getParam<bool>("cycle_models")),
86 4442 : _material_timestep_limit(declareProperty<Real>(_base_name + "material_timestep_limit")),
87 4442 : _identity_symmetric_four(RankFourTensor::initIdentitySymmetricFour),
88 8884 : _all_models_isotropic(true)
89 : {
90 4442 : }
91 :
92 : void
93 166928 : ComputeMultipleInelasticStressBase::initQpStatefulProperties()
94 : {
95 166928 : ComputeStressBase::initQpStatefulProperties();
96 166928 : _inelastic_strain[_qp].zero();
97 166928 : }
98 :
99 : void
100 4248 : ComputeMultipleInelasticStressBase::initialSetup()
101 : {
102 4248 : _damage_model = isParamValid("damage_model")
103 4284 : ? dynamic_cast<DamageBaseTempl<false> *>(&getMaterial("damage_model"))
104 : : nullptr;
105 :
106 4248 : _is_elasticity_tensor_guaranteed_isotropic =
107 4248 : hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC);
108 :
109 4248 : std::vector<MaterialName> models = getInelasticModelNames();
110 :
111 4248 : _num_models = models.size();
112 4248 : _tangent_computation_flag.resize(_num_models, false);
113 4248 : _consistent_tangent_operator.resize(_num_models);
114 :
115 4248 : _inelastic_weights = isParamValid("combined_inelastic_strain_weights")
116 8532 : ? getParam<std::vector<Real>>("combined_inelastic_strain_weights")
117 8496 : : std::vector<Real>(_num_models, 1.0);
118 :
119 4248 : 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 8516 : for (const auto i : make_range(_num_models))
128 : {
129 4270 : StressUpdateBase * rrr = dynamic_cast<StressUpdateBase *>(&getMaterialByName(models[i]));
130 :
131 4270 : if (rrr)
132 : {
133 4270 : _models.push_back(rrr);
134 4270 : _all_models_isotropic = _all_models_isotropic && rrr->isIsotropic();
135 4270 : 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 4246 : if (_tangent_operator_type != TangentOperatorEnum::elastic)
153 : {
154 : bool full_present = false;
155 : bool partial_present = false;
156 7214 : for (const auto i : make_range(_num_models))
157 : {
158 3546 : if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::FULL)
159 : {
160 : full_present = true;
161 : _tangent_computation_flag[i] = true;
162 2568 : _tangent_calculation_method = TangentCalculationMethod::FULL;
163 : }
164 978 : else if (_models[i]->getTangentCalculationMethod() == TangentCalculationMethod::PARTIAL)
165 : {
166 : partial_present = true;
167 : _tangent_computation_flag[i] = true;
168 750 : _tangent_calculation_method = TangentCalculationMethod::PARTIAL;
169 : }
170 : }
171 3668 : 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 12738 : 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 8512 : for (const auto model_number : index_range(_models))
186 : {
187 4268 : const bool use_substep = _models[model_number]->substeppingCapabilityRequested();
188 4268 : 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 4 : << _models[model_number]->name() << " does not implement substepping yet.";
193 2 : mooseError(error_message.str());
194 0 : }
195 : }
196 4244 : }
197 :
198 : void
199 22032756 : ComputeMultipleInelasticStressBase::computeQpStress()
200 : {
201 22032756 : if (_damage_model)
202 : {
203 6304 : _undamaged_stress_old = _stress_old[_qp];
204 : _damage_model->setQp(_qp);
205 6304 : _damage_model->computeUndamagedOldStress(_undamaged_stress_old);
206 : }
207 22032756 : computeQpStressIntermediateConfiguration();
208 :
209 22032640 : if (_damage_model)
210 : {
211 6304 : _damage_model->setQp(_qp);
212 6304 : _damage_model->updateDamage();
213 6304 : _damage_model->updateStressForDamage(_stress[_qp]);
214 6304 : _damage_model->finiteStrainRotation(_rotation_increment[_qp]);
215 6304 : _damage_model->updateJacobianMultForDamage(_Jacobian_mult[_qp]);
216 :
217 6304 : const Real damage_timestep_limit = _damage_model->computeTimeStepLimit();
218 6304 : if (_material_timestep_limit[_qp] > damage_timestep_limit)
219 5696 : _material_timestep_limit[_qp] = damage_timestep_limit;
220 : }
221 :
222 22032640 : if (_perform_finite_strain_rotations)
223 20503776 : finiteStrainRotation();
224 22032640 : }
225 :
226 : void
227 22078772 : ComputeMultipleInelasticStressBase::computeQpStressIntermediateConfiguration()
228 : {
229 22078772 : RankTwoTensor elastic_strain_increment;
230 22078772 : RankTwoTensor combined_inelastic_strain_increment;
231 :
232 22078772 : if (_num_models == 0)
233 : {
234 197696 : _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 197696 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
239 197696 : _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 197696 : if (_fe_problem.currentlyComputingJacobian())
248 45696 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
249 :
250 197696 : _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
251 : }
252 : else
253 : {
254 21881076 : if (_num_models == 1 || _cycle_models)
255 20115972 : updateQpStateSingleModel((_t_step - 1) % _num_models,
256 : elastic_strain_increment,
257 : combined_inelastic_strain_increment);
258 : else
259 1765104 : updateQpState(elastic_strain_increment, combined_inelastic_strain_increment);
260 :
261 21880960 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + elastic_strain_increment;
262 21880960 : _inelastic_strain[_qp] = _inelastic_strain_old[_qp] + combined_inelastic_strain_increment;
263 : }
264 22078656 : }
265 :
266 : void
267 21510576 : ComputeMultipleInelasticStressBase::finiteStrainRotation(const bool force_elasticity_rotation)
268 : {
269 21510576 : _elastic_strain[_qp] =
270 21510576 : _rotation_increment[_qp] * _elastic_strain[_qp] * _rotation_increment[_qp].transpose();
271 21510576 : _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
272 21510576 : _inelastic_strain[_qp] =
273 21510576 : _rotation_increment[_qp] * _inelastic_strain[_qp] * _rotation_increment[_qp].transpose();
274 :
275 21510576 : if (force_elasticity_rotation ||
276 20549792 : !(_is_elasticity_tensor_guaranteed_isotropic &&
277 20541152 : (_all_models_isotropic ||
278 1177680 : _tangent_calculation_method == TangentCalculationMethod::ELASTIC || _num_models == 0)))
279 998928 : _Jacobian_mult[_qp].rotate(_rotation_increment[_qp]);
280 21510576 : }
281 :
282 : void
283 176672 : ComputeMultipleInelasticStressBase::computeQpJacobianMult()
284 : {
285 176672 : if (_tangent_calculation_method == TangentCalculationMethod::ELASTIC)
286 170880 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
287 5792 : else if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
288 : {
289 5760 : RankFourTensor A = _identity_symmetric_four;
290 17280 : for (const auto i_rmm : make_range(_num_models))
291 11520 : A += _consistent_tangent_operator[i_rmm];
292 : mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
293 5760 : _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 176672 : }
303 :
304 : void
305 20115972 : ComputeMultipleInelasticStressBase::updateQpStateSingleModel(
306 : unsigned model_number,
307 : RankTwoTensor & elastic_strain_increment,
308 : RankTwoTensor & combined_inelastic_strain_increment)
309 : {
310 40245816 : for (auto model : _models)
311 20129844 : model->setQp(_qp);
312 :
313 20115972 : 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 20115972 : if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
318 20110212 : _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 20115972 : computeAdmissibleState(model_number,
328 : elastic_strain_increment,
329 : combined_inelastic_strain_increment,
330 : _consistent_tangent_operator[0]);
331 :
332 20115856 : if (_fe_problem.currentlyComputingJacobian())
333 : {
334 928144 : if (_tangent_calculation_method == TangentCalculationMethod::ELASTIC)
335 528784 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
336 399360 : else if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
337 : {
338 103696 : RankFourTensor A = _identity_symmetric_four + _consistent_tangent_operator[0];
339 : mooseAssert(A.isSymmetric(), "Tangent operator isn't symmetric");
340 103696 : _Jacobian_mult[_qp] = A.invSymm() * _elasticity_tensor[_qp];
341 : }
342 : else
343 295664 : _Jacobian_mult[_qp] = _consistent_tangent_operator[0];
344 : }
345 :
346 20115856 : _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 40245584 : for (const auto i_rmm : make_range(_num_models))
351 20129728 : if (i_rmm != model_number)
352 13872 : _models[i_rmm]->propagateQpStatefulProperties();
353 20115856 : }
354 :
355 : void
356 33063340 : 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 33063340 : _models[model_number]->resetIncrementalMaterialProperties();
364 :
365 33063340 : const bool jac = _fe_problem.currentlyComputingJacobian();
366 33063340 : if (_damage_model)
367 12608 : _models[model_number]->updateState(elastic_strain_increment,
368 : inelastic_strain_increment,
369 6304 : _rotation_increment[_qp],
370 6304 : _stress[_qp],
371 6304 : _undamaged_stress_old,
372 6304 : _elasticity_tensor[_qp],
373 6304 : _elastic_strain_old[_qp],
374 6304 : (jac && _tangent_computation_flag[model_number]),
375 : consistent_tangent_operator);
376 33057036 : else if (_models[model_number]->substeppingCapabilityEnabled() &&
377 26602 : (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations))
378 53204 : _models[model_number]->updateStateSubstep(elastic_strain_increment,
379 : inelastic_strain_increment,
380 26602 : _rotation_increment[_qp],
381 26602 : _stress[_qp],
382 26602 : _stress_old[_qp],
383 26602 : _elasticity_tensor[_qp],
384 26602 : _elastic_strain_old[_qp],
385 26602 : (jac && _tangent_computation_flag[model_number]),
386 : consistent_tangent_operator);
387 : else
388 66060868 : _models[model_number]->updateState(elastic_strain_increment,
389 : inelastic_strain_increment,
390 33030434 : _rotation_increment[_qp],
391 33030434 : _stress[_qp],
392 33030434 : _stress_old[_qp],
393 33030434 : _elasticity_tensor[_qp],
394 33030434 : _elastic_strain_old[_qp],
395 33030434 : (jac && _tangent_computation_flag[model_number]),
396 : consistent_tangent_operator);
397 :
398 33063224 : if (jac && !_tangent_computation_flag[model_number])
399 : {
400 2063920 : if (_tangent_calculation_method == TangentCalculationMethod::PARTIAL)
401 0 : consistent_tangent_operator.zero();
402 : else
403 2063920 : consistent_tangent_operator = _elasticity_tensor[_qp];
404 : }
405 33063224 : }
|
