Line data Source code
1 :
2 : //* This file is part of the MOOSE framework
3 : //* https://mooseframework.inl.gov
4 : //*
5 : //* All rights reserved, see COPYRIGHT for full restrictions
6 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
7 : //*
8 : //* Licensed under LGPL 2.1, please see LICENSE for details
9 : //* https://www.gnu.org/licenses/lgpl-2.1.html
10 :
11 : #include "RadialReturnStressUpdate.h"
12 :
13 : #include "MooseMesh.h"
14 : #include "ElasticityTensorTools.h"
15 :
16 : template <bool is_ad>
17 : InputParameters
18 5390 : RadialReturnStressUpdateTempl<is_ad>::validParams()
19 : {
20 5390 : InputParameters params = StressUpdateBaseTempl<is_ad>::validParams();
21 5390 : params.addClassDescription("Calculates the effective inelastic strain increment required to "
22 : "return the isotropic stress state to a J2 yield surface. This class "
23 : "is intended to be a parent class for classes with specific "
24 : "constitutive models.");
25 5390 : params += SingleVariableReturnMappingSolutionTempl<is_ad>::validParams();
26 10780 : params.addParam<Real>("max_inelastic_increment",
27 10780 : 1e-4,
28 : "The maximum inelastic strain increment allowed in a time step");
29 10780 : params.addRequiredParam<std::string>(
30 : "effective_inelastic_strain_name",
31 : "Name of the material property that stores the effective inelastic strain");
32 10780 : params.addParam<bool>("use_substep", false, "Whether to use substepping");
33 :
34 10780 : MooseEnum substeppingType("NONE ERROR_BASED INCREMENT_BASED", "NONE");
35 10780 : substeppingType.addDocumentation("NONE", "Do not use substepping");
36 10780 : substeppingType.addDocumentation(
37 : "ERROR_BASED",
38 : "Use substepping with a substep size that will yield, at most, the creep numerical "
39 : "integration error given by substep_strain_tolerance.");
40 10780 : substeppingType.addDocumentation(
41 : "INCREMENT_BASED",
42 : "Use substepping with a substep size that will keep each inelastic strain increment below "
43 : "the maximum inelastic strain increment allowed in a time step.");
44 10780 : params.addParam<MooseEnum>(
45 : "use_substepping", substeppingType, "Whether and how to use substepping");
46 10780 : params.addParam<bool>(
47 : "adaptive_substepping",
48 10780 : false,
49 : "Use adaptive substepping, where the number of substeps is successively doubled until the "
50 : "return mapping model successfully converges or the maximum number of substeps is reached. ");
51 :
52 16170 : params.addDeprecatedParam<bool>(
53 10780 : "use_substep", false, "Whether to use substepping", "Use `use_substepping` instead");
54 10780 : params.addParam<Real>("substep_strain_tolerance",
55 10780 : 0.1,
56 : "Maximum ratio of the initial elastic strain increment at start of the "
57 : "return mapping solve to the maximum inelastic strain allowable in a "
58 : "single substep. Reduce this value to increase the number of substeps");
59 10780 : params.addParam<bool>("apply_strain", true, "Flag to apply strain. Used for testing.");
60 10780 : params.addParamNamesToGroup(
61 : "effective_inelastic_strain_name substep_strain_tolerance apply_strain", "Advanced");
62 10780 : params.addParam<bool>("use_substep_integration_error",
63 10780 : false,
64 : "If true, it establishes a substep size that will yield, at most,"
65 : "the creep numerical integration error given by substep_strain_tolerance.");
66 10780 : params.addParam<unsigned>("maximum_number_substeps",
67 10780 : 25,
68 : "The maximum number of substeps allowed before cutting the time step.");
69 5390 : return params;
70 5390 : }
71 :
72 : template <bool is_ad>
73 4045 : RadialReturnStressUpdateTempl<is_ad>::RadialReturnStressUpdateTempl(
74 : const InputParameters & parameters)
75 : : StressUpdateBaseTempl<is_ad>(parameters),
76 : SingleVariableReturnMappingSolutionTempl<is_ad>(parameters),
77 5538 : _effective_inelastic_strain(this->template declareGenericProperty<Real, is_ad>(
78 4045 : this->_base_name +
79 : this->template getParam<std::string>("effective_inelastic_strain_name"))),
80 13628 : _effective_inelastic_strain_old(this->template getMaterialPropertyOld<Real>(
81 : this->_base_name +
82 : this->template getParam<std::string>("effective_inelastic_strain_name"))),
83 8090 : _max_inelastic_increment(this->template getParam<Real>("max_inelastic_increment")),
84 8090 : _substep_tolerance(this->template getParam<Real>("substep_strain_tolerance")),
85 4045 : _identity_two(RankTwoTensor::initIdentity),
86 4045 : _identity_symmetric_four(RankFourTensor::initIdentitySymmetricFour),
87 4045 : _deviatoric_projection_four(_identity_symmetric_four -
88 4045 : _identity_two.outerProduct(_identity_two) / 3.0),
89 8090 : _apply_strain(this->template getParam<bool>("apply_strain")),
90 4045 : _use_substepping(
91 4045 : this->template getParam<MooseEnum>("use_substepping").template getEnum<SubsteppingType>()),
92 8090 : _adaptive_substepping(this->template getParam<bool>("adaptive_substepping")),
93 12135 : _maximum_number_substeps(this->template getParam<unsigned>("maximum_number_substeps"))
94 : {
95 8090 : if (this->_pars.isParamSetByUser("use_substep"))
96 : {
97 12 : if (this->_pars.isParamSetByUser("use_substepping"))
98 0 : this->paramError("use_substep",
99 : "Remove this parameter and just keep `use_substepping` in the input");
100 :
101 6 : if (parameters.get<bool>("use_substep"))
102 : {
103 6 : if (parameters.get<bool>("use_substep_integration_error"))
104 0 : _use_substepping = SubsteppingType::ERROR_BASED;
105 : else
106 6 : _use_substepping = SubsteppingType::INCREMENT_BASED;
107 : }
108 : }
109 :
110 8090 : if (this->_pars.isParamSetByUser("maximum_number_substeps") &&
111 66 : _use_substepping == SubsteppingType::NONE)
112 0 : this->paramError(
113 : "maximum_number_substeps",
114 : "The parameter maximum_number_substeps can only be used when the substepping option "
115 : "(use_substepping) is not set to NONE");
116 :
117 4045 : if (_adaptive_substepping && _use_substepping == SubsteppingType::NONE)
118 0 : this->paramError(
119 : "adaptive_substepping",
120 : "The parameter adaptive_substepping can only be used when the substepping option "
121 : "(use_substepping) is not set to NONE");
122 4045 : }
123 :
124 : template <bool is_ad>
125 : void
126 26376 : RadialReturnStressUpdateTempl<is_ad>::initQpStatefulProperties()
127 : {
128 75232 : _effective_inelastic_strain[_qp] = 0.0;
129 26376 : }
130 :
131 : template <bool is_ad>
132 : void
133 29766204 : RadialReturnStressUpdateTempl<is_ad>::computeStressInitialize(
134 : const GenericReal<is_ad> & /*effective_trial_stress*/,
135 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
136 : {
137 : // Set the value of 3 * shear modulus for use as a reference residual value
138 86597218 : _three_shear_modulus = 3.0 * ElasticityTensorTools::getIsotropicShearModulus(elasticity_tensor);
139 29766204 : }
140 :
141 : template <bool is_ad>
142 : void
143 549504 : RadialReturnStressUpdateTempl<is_ad>::propagateQpStatefulPropertiesRadialReturn()
144 : {
145 549504 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
146 549504 : }
147 :
148 : template <bool is_ad>
149 : int
150 26238 : RadialReturnStressUpdateTempl<is_ad>::calculateNumberSubsteps(
151 : const GenericRankTwoTensor<is_ad> & strain_increment)
152 : {
153 : // compute an effective elastic strain measure
154 4740 : const GenericReal<is_ad> contracted_elastic_strain =
155 21498 : strain_increment.doubleContraction(strain_increment);
156 : const Real effective_elastic_strain =
157 26238 : std::sqrt(3.0 / 2.0 * MetaPhysicL::raw_value(contracted_elastic_strain));
158 :
159 26238 : if (MooseUtils::absoluteFuzzyEqual(effective_elastic_strain, 0.0))
160 : return 1;
161 :
162 18166 : if (_use_substepping == SubsteppingType::INCREMENT_BASED)
163 : {
164 10008 : const Real ratio = effective_elastic_strain / _max_inelastic_increment;
165 :
166 10008 : if (ratio > _substep_tolerance)
167 10002 : return std::ceil(ratio / _substep_tolerance);
168 : return 1;
169 : }
170 :
171 8158 : if (_use_substepping == SubsteppingType::ERROR_BASED)
172 : {
173 8158 : const Real accurate_time_step_ratio = _substep_tolerance / effective_elastic_strain;
174 :
175 8158 : if (accurate_time_step_ratio < 1.0)
176 2764 : return std::ceil(1.0 / accurate_time_step_ratio);
177 : return 1;
178 : }
179 :
180 0 : mooseError("calculateNumberSubsteps should not have been called. Notify a developer.");
181 : }
182 :
183 : template <bool is_ad>
184 : void
185 0 : RadialReturnStressUpdateTempl<is_ad>::computeTangentOperator(Real /*effective_trial_stress*/,
186 : const RankTwoTensor & /*stress_new*/,
187 : RankFourTensor & /*tangent_operator*/)
188 : {
189 0 : mooseError("computeTangentOperator called: no tangent computation is needed for AD");
190 : }
191 :
192 : template <>
193 : void
194 678956 : RadialReturnStressUpdateTempl<false>::computeTangentOperator(Real effective_trial_stress,
195 : const RankTwoTensor & stress_new,
196 : RankFourTensor & tangent_operator)
197 : {
198 678956 : if (getTangentCalculationMethod() == TangentCalculationMethod::PARTIAL)
199 : {
200 603284 : if (MooseUtils::absoluteFuzzyEqual(_effective_inelastic_strain_increment, 0.0))
201 75614 : tangent_operator.zero();
202 : else
203 : {
204 : // mu = _three_shear_modulus / 3.0;
205 : // norm_dev_stress = ||s_n+1||
206 : // effective_trial_stress = von mises trial stress = std::sqrt(3.0 / 2.0) * ||s_n+1^trial||
207 : // scalar_effective_inelastic_strain = Delta epsilon^cr_n+1
208 : // deriv = derivative of scalar_effective_inelastic_strain w.r.t. von mises stress
209 : // deriv = std::sqrt(3.0 / 2.0) partial Delta epsilon^cr_n+1n over partial ||s_n+1^trial||
210 :
211 : mooseAssert(_three_shear_modulus != 0.0, "Shear modulus is zero");
212 :
213 527670 : const RankTwoTensor deviatoric_stress = stress_new.deviatoric();
214 527670 : const Real norm_dev_stress_squared = deviatoric_stress.doubleContraction(deviatoric_stress);
215 527670 : if (MooseUtils::absoluteFuzzyEqual(norm_dev_stress_squared, 0.0))
216 : {
217 0 : tangent_operator.zero();
218 : return;
219 : }
220 527670 : const Real norm_dev_stress = std::sqrt(norm_dev_stress_squared);
221 :
222 527670 : const RankTwoTensor flow_direction = deviatoric_stress / norm_dev_stress;
223 : const RankFourTensor flow_direction_dyad = flow_direction.outerProduct(flow_direction);
224 : const Real deriv =
225 527670 : computeStressDerivative(effective_trial_stress, _effective_inelastic_strain_increment);
226 527670 : const Real scalar_one = _three_shear_modulus * _effective_inelastic_strain_increment /
227 527670 : std::sqrt(1.5) / norm_dev_stress;
228 :
229 527670 : tangent_operator = scalar_one * _deviatoric_projection_four +
230 1055340 : (_three_shear_modulus * deriv - scalar_one) * flow_direction_dyad;
231 : }
232 : }
233 : }
234 :
235 : template <bool is_ad>
236 : void
237 58555054 : RadialReturnStressUpdateTempl<is_ad>::updateState(
238 : GenericRankTwoTensor<is_ad> & strain_increment,
239 : GenericRankTwoTensor<is_ad> & inelastic_strain_increment,
240 : const GenericRankTwoTensor<is_ad> & /*rotation_increment*/,
241 : GenericRankTwoTensor<is_ad> & stress_new,
242 : const RankTwoTensor & /*stress_old*/,
243 : const GenericRankFourTensor<is_ad> & elasticity_tensor,
244 : const RankTwoTensor & elastic_strain_old,
245 : bool compute_full_tangent_operator,
246 : RankFourTensor & tangent_operator)
247 : {
248 : using std::sqrt;
249 :
250 : // compute the deviatoric trial stress and trial strain from the current intermediate
251 : // configuration
252 58555054 : GenericRankTwoTensor<is_ad> deviatoric_trial_stress = stress_new.deviatoric();
253 :
254 : // compute the effective trial stress
255 29766204 : GenericReal<is_ad> dev_trial_stress_squared =
256 28788850 : deviatoric_trial_stress.doubleContraction(deviatoric_trial_stress);
257 57781970 : GenericReal<is_ad> effective_trial_stress = MetaPhysicL::raw_value(dev_trial_stress_squared)
258 58555054 : ? sqrt(3.0 / 2.0 * dev_trial_stress_squared)
259 29766204 : : 0.0;
260 :
261 58555054 : computeStressInitialize(effective_trial_stress, elasticity_tensor);
262 :
263 : mooseAssert(
264 : _three_shear_modulus != 0.0,
265 : "Shear modulus is zero. Ensure that the base class computeStressInitialize() is called.");
266 :
267 : // Use Newton iteration to determine the scalar effective inelastic strain increment
268 58555050 : _effective_inelastic_strain_increment = 0.0;
269 58555050 : if (!MooseUtils::absoluteFuzzyEqual(effective_trial_stress, 0.0))
270 : {
271 56889686 : this->returnMappingSolve(
272 56889686 : effective_trial_stress, _effective_inelastic_strain_increment, this->_console);
273 56889658 : if (_effective_inelastic_strain_increment != 0.0)
274 54197785 : inelastic_strain_increment =
275 : deviatoric_trial_stress *
276 82442780 : (1.5 * _effective_inelastic_strain_increment / effective_trial_stress);
277 : else
278 748107 : inelastic_strain_increment.zero();
279 : }
280 : else
281 773084 : inelastic_strain_increment.zero();
282 :
283 58555022 : if (_apply_strain)
284 : {
285 58540702 : strain_increment -= inelastic_strain_increment;
286 58540702 : updateEffectiveInelasticStrain(_effective_inelastic_strain_increment);
287 :
288 : // Use the old elastic strain here because we require tensors used by this class
289 : // to be isotropic and this method natively allows for changing in time
290 : // elasticity tensors
291 88304736 : stress_new = elasticity_tensor * (strain_increment + elastic_strain_old);
292 : }
293 :
294 58555022 : computeStressFinalize(inelastic_strain_increment);
295 :
296 : if constexpr (!is_ad)
297 : {
298 28788836 : if (compute_full_tangent_operator)
299 203888 : computeTangentOperator(effective_trial_stress, stress_new, tangent_operator);
300 : }
301 : else
302 : {
303 : libmesh_ignore(compute_full_tangent_operator);
304 : libmesh_ignore(tangent_operator);
305 : }
306 58555022 : }
307 :
308 : template <bool is_ad>
309 : void
310 26246 : RadialReturnStressUpdateTempl<is_ad>::updateStateSubstepInternal(
311 : GenericRankTwoTensor<is_ad> & strain_increment,
312 : GenericRankTwoTensor<is_ad> & inelastic_strain_increment,
313 : const GenericRankTwoTensor<is_ad> & rotation_increment,
314 : GenericRankTwoTensor<is_ad> & stress_new,
315 : const RankTwoTensor & stress_old,
316 : const GenericRankFourTensor<is_ad> & elasticity_tensor,
317 : const RankTwoTensor & elastic_strain_old,
318 : unsigned int total_number_substeps,
319 : bool compute_full_tangent_operator,
320 : RankFourTensor & tangent_operator)
321 : {
322 : using std::sqrt;
323 :
324 : // if only one substep is needed, then call the original update state method
325 26246 : if (total_number_substeps == 1)
326 : {
327 13472 : updateState(strain_increment,
328 : inelastic_strain_increment,
329 : rotation_increment,
330 : stress_new,
331 : stress_old,
332 : elasticity_tensor,
333 : elastic_strain_old,
334 : compute_full_tangent_operator,
335 : tangent_operator);
336 :
337 13464 : this->storeIncrementalMaterialProperties(total_number_substeps);
338 13464 : return;
339 : }
340 :
341 12774 : if (total_number_substeps > _maximum_number_substeps)
342 22 : mooseException("The number of substeps computed exceeds 'maximum_number_substeps'.");
343 :
344 : // cut the original timestep
345 12752 : _dt = _dt_original / total_number_substeps;
346 :
347 : // initialize the inputs
348 12752 : const GenericRankTwoTensor<is_ad> strain_increment_per_step =
349 : strain_increment / total_number_substeps;
350 12752 : GenericRankTwoTensor<is_ad> sub_stress_new = elasticity_tensor * elastic_strain_old;
351 8524 : GenericRankTwoTensor<is_ad> sub_elastic_strain_old = elastic_strain_old;
352 :
353 : // clear the original inputs
354 12752 : MathUtils::mooseSetToZero(strain_increment);
355 12752 : MathUtils::mooseSetToZero(inelastic_strain_increment);
356 12752 : MathUtils::mooseSetToZero(stress_new);
357 :
358 12752 : GenericReal<is_ad> sub_effective_inelastic_strain_increment = 0.0;
359 8524 : GenericRankTwoTensor<is_ad> sub_inelastic_strain_increment = inelastic_strain_increment;
360 :
361 666236 : for (unsigned int step = 0; step < total_number_substeps; ++step)
362 : {
363 : // set up input for this substep
364 323728 : GenericRankTwoTensor<is_ad> sub_strain_increment = strain_increment_per_step;
365 653492 : sub_stress_new += elasticity_tensor * sub_strain_increment;
366 :
367 : Real effective_sub_stress_new;
368 : if constexpr (!is_ad)
369 : {
370 : // compute effective_sub_stress_new
371 323728 : const RankTwoTensor deviatoric_sub_stress_new = sub_stress_new.deviatoric();
372 : const Real dev_sub_stress_new_squared =
373 323728 : deviatoric_sub_stress_new.doubleContraction(deviatoric_sub_stress_new);
374 323728 : effective_sub_stress_new = sqrt(3.0 / 2.0 * dev_sub_stress_new_squared);
375 : }
376 : else
377 : libmesh_ignore(effective_sub_stress_new);
378 :
379 : // update stress and strain based on the strain increment
380 653492 : updateState(sub_strain_increment,
381 : sub_inelastic_strain_increment,
382 : rotation_increment, // not used in updateState
383 : sub_stress_new,
384 : stress_old, // not used in updateState
385 : elasticity_tensor,
386 : elastic_strain_old,
387 : false);
388 : // do not compute tangent until the end of this substep (or not at all for is_ad == true)
389 :
390 : // update strain and stress
391 653484 : strain_increment += sub_strain_increment;
392 653484 : inelastic_strain_increment += sub_inelastic_strain_increment;
393 653484 : sub_elastic_strain_old += sub_strain_increment;
394 653488 : sub_stress_new = elasticity_tensor * sub_elastic_strain_old;
395 :
396 : // accumulate scalar_effective_inelastic_strain
397 653484 : sub_effective_inelastic_strain_increment += _effective_inelastic_strain_increment;
398 :
399 : if constexpr (!is_ad)
400 323724 : computeTangentOperator(effective_sub_stress_new, sub_stress_new, tangent_operator);
401 :
402 : // store incremental material properties for this step
403 653484 : this->storeIncrementalMaterialProperties(total_number_substeps);
404 : }
405 :
406 : // update stress
407 8520 : stress_new = sub_stress_new;
408 :
409 : // update effective inelastic strain
410 12744 : updateEffectiveInelasticStrain(sub_effective_inelastic_strain_increment);
411 : }
412 :
413 : template <bool is_ad>
414 : void
415 26238 : RadialReturnStressUpdateTempl<is_ad>::updateStateSubstep(
416 : GenericRankTwoTensor<is_ad> & strain_increment,
417 : GenericRankTwoTensor<is_ad> & inelastic_strain_increment,
418 : const GenericRankTwoTensor<is_ad> & rotation_increment,
419 : GenericRankTwoTensor<is_ad> & stress_new,
420 : const RankTwoTensor & stress_old,
421 : const GenericRankFourTensor<is_ad> & elasticity_tensor,
422 : const RankTwoTensor & elastic_strain_old,
423 : bool compute_full_tangent_operator,
424 : RankFourTensor & tangent_operator)
425 : {
426 26238 : unsigned int num_substeps = calculateNumberSubsteps(strain_increment);
427 26238 : _dt_original = _dt;
428 : while (true)
429 : {
430 : try
431 : {
432 26246 : updateStateSubstepInternal(strain_increment,
433 : inelastic_strain_increment,
434 : rotation_increment,
435 : stress_new,
436 : stress_old,
437 : elasticity_tensor,
438 : elastic_strain_old,
439 : num_substeps,
440 : compute_full_tangent_operator,
441 : tangent_operator);
442 : }
443 68 : catch (MooseException & e)
444 : {
445 : // if we are not using adaptive substepping we just rethrow the exception
446 38 : if (!_adaptive_substepping)
447 22 : throw e;
448 :
449 : // otherwise we double the number of substeps and try again
450 16 : num_substeps *= 2;
451 16 : if (num_substeps <= _maximum_number_substeps)
452 : continue;
453 :
454 : // too meany substeps, break out of the loop
455 : break;
456 : }
457 :
458 : // updateStateSubstepInternal was successful (didn't throw)
459 26208 : _dt = _dt_original;
460 26208 : return;
461 : }
462 :
463 : // recover the original timestep
464 8 : _dt = _dt_original;
465 :
466 8 : mooseException("Adaptive substepping failed. Maximum number of substeps exceeded.");
467 : }
468 :
469 : template <bool is_ad>
470 : Real
471 297632479 : RadialReturnStressUpdateTempl<is_ad>::computeReferenceResidual(
472 : const GenericReal<is_ad> & effective_trial_stress,
473 : const GenericReal<is_ad> & scalar_effective_inelastic_strain)
474 : {
475 297632479 : return MetaPhysicL::raw_value(effective_trial_stress / _three_shear_modulus) -
476 297632479 : MetaPhysicL::raw_value(scalar_effective_inelastic_strain);
477 : }
478 :
479 : template <bool is_ad>
480 : GenericReal<is_ad>
481 56697310 : RadialReturnStressUpdateTempl<is_ad>::maximumPermissibleValue(
482 : const GenericReal<is_ad> & effective_trial_stress) const
483 : {
484 56697310 : return effective_trial_stress / _three_shear_modulus;
485 : }
486 :
487 : template <bool is_ad>
488 : Real
489 29247740 : RadialReturnStressUpdateTempl<is_ad>::computeTimeStepLimit()
490 : {
491 : const Real scalar_inelastic_strain_incr =
492 29247740 : std::abs(MetaPhysicL::raw_value(_effective_inelastic_strain[_qp]) -
493 29247740 : _effective_inelastic_strain_old[_qp]);
494 29247740 : if (!scalar_inelastic_strain_incr)
495 : return std::numeric_limits<Real>::max();
496 :
497 24680287 : return _dt * _max_inelastic_increment / scalar_inelastic_strain_incr;
498 : }
499 :
500 : template <bool is_ad>
501 : void
502 25364 : RadialReturnStressUpdateTempl<is_ad>::outputIterationSummary(std::stringstream * iter_output,
503 : const unsigned int total_it)
504 : {
505 25364 : if (iter_output)
506 : {
507 25364 : *iter_output << "At element " << _current_elem->id() << " _qp=" << _qp << " Coordinates "
508 25364 : << _q_point[_qp] << " block=" << _current_elem->subdomain_id() << '\n';
509 : }
510 25364 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationSummary(iter_output, total_it);
511 25364 : }
512 :
513 : template class RadialReturnStressUpdateTempl<false>;
514 : template class RadialReturnStressUpdateTempl<true>;
|
