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 4722 : RadialReturnStressUpdateTempl<is_ad>::validParams()
19 : {
20 4722 : InputParameters params = StressUpdateBaseTempl<is_ad>::validParams();
21 4722 : 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 4722 : params += SingleVariableReturnMappingSolutionTempl<is_ad>::validParams();
26 9444 : params.addParam<Real>("max_inelastic_increment",
27 9444 : 1e-4,
28 : "The maximum inelastic strain increment allowed in a time step");
29 9444 : params.addRequiredParam<std::string>(
30 : "effective_inelastic_strain_name",
31 : "Name of the material property that stores the effective inelastic strain");
32 9444 : params.addParam<bool>("use_substep", false, "Whether to use substepping");
33 :
34 9444 : MooseEnum substeppingType("NONE ERROR_BASED INCREMENT_BASED", "NONE");
35 9444 : substeppingType.addDocumentation("NONE", "Do not use substepping");
36 9444 : 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 9444 : 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 9444 : params.addParam<MooseEnum>(
45 : "use_substepping", substeppingType, "Whether and how to use substepping");
46 9444 : params.addParam<bool>(
47 : "adaptive_substepping",
48 9444 : 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 14166 : params.addDeprecatedParam<bool>(
53 9444 : "use_substep", false, "Whether to use substepping", "Use `use_substepping` instead");
54 9444 : params.addParam<Real>("substep_strain_tolerance",
55 9444 : 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 9444 : params.addParam<bool>("apply_strain", true, "Flag to apply strain. Used for testing.");
60 9444 : params.addParamNamesToGroup(
61 : "effective_inelastic_strain_name substep_strain_tolerance apply_strain", "Advanced");
62 9444 : params.addParam<bool>("use_substep_integration_error",
63 9444 : 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 9444 : params.addParam<unsigned>("maximum_number_substeps",
67 9444 : 25,
68 : "The maximum number of substeps allowed before cutting the time step.");
69 4722 : return params;
70 4722 : }
71 :
72 : template <bool is_ad>
73 3544 : RadialReturnStressUpdateTempl<is_ad>::RadialReturnStressUpdateTempl(
74 : const InputParameters & parameters)
75 : : StressUpdateBaseTempl<is_ad>(parameters),
76 : SingleVariableReturnMappingSolutionTempl<is_ad>(parameters),
77 4920 : _effective_inelastic_strain(this->template declareGenericProperty<Real, is_ad>(
78 3544 : this->_base_name +
79 : this->template getParam<std::string>("effective_inelastic_strain_name"))),
80 12008 : _effective_inelastic_strain_old(this->template getMaterialPropertyOld<Real>(
81 : this->_base_name +
82 : this->template getParam<std::string>("effective_inelastic_strain_name"))),
83 7088 : _max_inelastic_increment(this->template getParam<Real>("max_inelastic_increment")),
84 7088 : _substep_tolerance(this->template getParam<Real>("substep_strain_tolerance")),
85 3544 : _identity_two(RankTwoTensor::initIdentity),
86 3544 : _identity_symmetric_four(RankFourTensor::initIdentitySymmetricFour),
87 3544 : _deviatoric_projection_four(_identity_symmetric_four -
88 3544 : _identity_two.outerProduct(_identity_two) / 3.0),
89 7088 : _apply_strain(this->template getParam<bool>("apply_strain")),
90 3544 : _use_substepping(
91 3544 : this->template getParam<MooseEnum>("use_substepping").template getEnum<SubsteppingType>()),
92 7088 : _adaptive_substepping(this->template getParam<bool>("adaptive_substepping")),
93 10632 : _maximum_number_substeps(this->template getParam<unsigned>("maximum_number_substeps"))
94 : {
95 7088 : 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 7088 : if (this->_pars.isParamSetByUser("maximum_number_substeps") &&
111 60 : _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 3544 : 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 3544 : }
123 :
124 : template <bool is_ad>
125 : void
126 25552 : RadialReturnStressUpdateTempl<is_ad>::initQpStatefulProperties()
127 : {
128 69088 : _effective_inelastic_strain[_qp] = 0.0;
129 25552 : }
130 :
131 : template <bool is_ad>
132 : void
133 27318506 : 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 83255658 : _three_shear_modulus = 3.0 * ElasticityTensorTools::getIsotropicShearModulus(elasticity_tensor);
139 27318506 : }
140 :
141 : template <bool is_ad>
142 : void
143 517968 : RadialReturnStressUpdateTempl<is_ad>::propagateQpStatefulPropertiesRadialReturn()
144 : {
145 517968 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
146 517968 : }
147 :
148 : template <bool is_ad>
149 : int
150 47908 : RadialReturnStressUpdateTempl<is_ad>::calculateNumberSubsteps(
151 : const GenericRankTwoTensor<is_ad> & strain_increment)
152 : {
153 : // compute an effective elastic strain measure
154 21306 : const GenericReal<is_ad> contracted_elastic_strain =
155 26602 : strain_increment.doubleContraction(strain_increment);
156 : const Real effective_elastic_strain =
157 47908 : std::sqrt(3.0 / 2.0 * MetaPhysicL::raw_value(contracted_elastic_strain));
158 :
159 47908 : if (MooseUtils::absoluteFuzzyEqual(effective_elastic_strain, 0.0))
160 : return 1;
161 :
162 33380 : if (_use_substepping == SubsteppingType::INCREMENT_BASED)
163 : {
164 22350 : const Real ratio = effective_elastic_strain / _max_inelastic_increment;
165 :
166 22350 : if (ratio > _substep_tolerance)
167 22344 : return std::ceil(ratio / _substep_tolerance);
168 : return 1;
169 : }
170 :
171 11030 : if (_use_substepping == SubsteppingType::ERROR_BASED)
172 : {
173 11030 : const Real accurate_time_step_ratio = _substep_tolerance / effective_elastic_strain;
174 :
175 11030 : if (accurate_time_step_ratio < 1.0)
176 4292 : 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 589552 : RadialReturnStressUpdateTempl<false>::computeTangentOperator(Real effective_trial_stress,
195 : const RankTwoTensor & stress_new,
196 : RankFourTensor & tangent_operator)
197 : {
198 589552 : if (getTangentCalculationMethod() == TangentCalculationMethod::PARTIAL)
199 : {
200 508992 : if (MooseUtils::absoluteFuzzyEqual(_effective_inelastic_strain_increment, 0.0))
201 65320 : 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 443672 : const RankTwoTensor deviatoric_stress = stress_new.deviatoric();
214 443672 : const Real norm_dev_stress_squared = deviatoric_stress.doubleContraction(deviatoric_stress);
215 443672 : if (MooseUtils::absoluteFuzzyEqual(norm_dev_stress_squared, 0.0))
216 : {
217 0 : tangent_operator.zero();
218 : return;
219 : }
220 443672 : const Real norm_dev_stress = std::sqrt(norm_dev_stress_squared);
221 :
222 443672 : 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 443672 : computeStressDerivative(effective_trial_stress, _effective_inelastic_strain_increment);
226 443672 : const Real scalar_one = _three_shear_modulus * _effective_inelastic_strain_increment /
227 443672 : std::sqrt(1.5) / norm_dev_stress;
228 :
229 443672 : tangent_operator = scalar_one * _deviatoric_projection_four +
230 887344 : (_three_shear_modulus * deriv - scalar_one) * flow_direction_dyad;
231 : }
232 : }
233 : }
234 :
235 : template <bool is_ad>
236 : void
237 58043632 : 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 :
249 : // compute the deviatoric trial stress and trial strain from the current intermediate
250 : // configuration
251 58043632 : GenericRankTwoTensor<is_ad> deviatoric_trial_stress = stress_new.deviatoric();
252 :
253 : // compute the effective trial stress
254 27318506 : GenericReal<is_ad> dev_trial_stress_squared =
255 30725126 : deviatoric_trial_stress.doubleContraction(deviatoric_trial_stress);
256 57241184 : GenericReal<is_ad> effective_trial_stress = MetaPhysicL::raw_value(dev_trial_stress_squared)
257 58043632 : ? std::sqrt(3.0 / 2.0 * dev_trial_stress_squared)
258 27318506 : : 0.0;
259 :
260 58043632 : computeStressInitialize(effective_trial_stress, elasticity_tensor);
261 :
262 : mooseAssert(
263 : _three_shear_modulus != 0.0,
264 : "Shear modulus is zero. Ensure that the base class computeStressInitialize() is called.");
265 :
266 : // Use Newton iteration to determine the scalar effective inelastic strain increment
267 58043628 : _effective_inelastic_strain_increment = 0.0;
268 58043628 : if (!MooseUtils::absoluteFuzzyEqual(effective_trial_stress, 0.0))
269 : {
270 56366236 : this->returnMappingSolve(
271 56366236 : effective_trial_stress, _effective_inelastic_strain_increment, this->_console);
272 56366208 : if (_effective_inelastic_strain_increment != 0.0)
273 53454856 : inelastic_strain_increment =
274 : deviatoric_trial_stress *
275 79235006 : (1.5 * _effective_inelastic_strain_increment / effective_trial_stress);
276 : else
277 735890 : inelastic_strain_increment.zero();
278 : }
279 : else
280 802448 : inelastic_strain_increment.zero();
281 :
282 58043600 : if (_apply_strain)
283 : {
284 58031504 : strain_increment -= inelastic_strain_increment;
285 58031504 : updateEffectiveInelasticStrain(_effective_inelastic_strain_increment);
286 :
287 : // Use the old elastic strain here because we require tensors used by this class
288 : // to be isotropic and this method natively allows for changing in time
289 : // elasticity tensors
290 85347912 : stress_new = elasticity_tensor * (strain_increment + elastic_strain_old);
291 : }
292 :
293 58043600 : computeStressFinalize(inelastic_strain_increment);
294 :
295 : if constexpr (!is_ad)
296 : {
297 30725112 : if (compute_full_tangent_operator)
298 181024 : computeTangentOperator(effective_trial_stress, stress_new, tangent_operator);
299 : }
300 : else
301 : {
302 : libmesh_ignore(compute_full_tangent_operator);
303 : libmesh_ignore(tangent_operator);
304 : }
305 58043600 : }
306 :
307 : template <bool is_ad>
308 : void
309 47916 : RadialReturnStressUpdateTempl<is_ad>::updateStateSubstepInternal(
310 : GenericRankTwoTensor<is_ad> & strain_increment,
311 : GenericRankTwoTensor<is_ad> & inelastic_strain_increment,
312 : const GenericRankTwoTensor<is_ad> & rotation_increment,
313 : GenericRankTwoTensor<is_ad> & stress_new,
314 : const RankTwoTensor & stress_old,
315 : const GenericRankFourTensor<is_ad> & elasticity_tensor,
316 : const RankTwoTensor & elastic_strain_old,
317 : unsigned int total_number_substeps,
318 : bool compute_full_tangent_operator,
319 : RankFourTensor & tangent_operator)
320 : {
321 : // if only one substep is needed, then call the original update state method
322 47916 : if (total_number_substeps == 1)
323 : {
324 21272 : updateState(strain_increment,
325 : inelastic_strain_increment,
326 : rotation_increment,
327 : stress_new,
328 : stress_old,
329 : elasticity_tensor,
330 : elastic_strain_old,
331 : compute_full_tangent_operator,
332 : tangent_operator);
333 :
334 21264 : this->storeIncrementalMaterialProperties(total_number_substeps);
335 21264 : return;
336 : }
337 :
338 26644 : if (total_number_substeps > _maximum_number_substeps)
339 76 : mooseException("The number of substeps computed exceeds 'maximum_number_substeps'.");
340 :
341 : // cut the original timestep
342 26568 : _dt = _dt_original / total_number_substeps;
343 :
344 : // initialize the inputs
345 26568 : const GenericRankTwoTensor<is_ad> strain_increment_per_step =
346 : strain_increment / total_number_substeps;
347 26568 : GenericRankTwoTensor<is_ad> sub_stress_new = elasticity_tensor * elastic_strain_old;
348 8900 : GenericRankTwoTensor<is_ad> sub_elastic_strain_old = elastic_strain_old;
349 :
350 : // clear the original inputs
351 26568 : MathUtils::mooseSetToZero(strain_increment);
352 26568 : MathUtils::mooseSetToZero(inelastic_strain_increment);
353 26568 : MathUtils::mooseSetToZero(stress_new);
354 :
355 26568 : GenericReal<is_ad> sub_effective_inelastic_strain_increment = 0.0;
356 8900 : GenericRankTwoTensor<is_ad> sub_inelastic_strain_increment = inelastic_strain_increment;
357 :
358 651272 : for (unsigned int step = 0; step < total_number_substeps; ++step)
359 : {
360 : // set up input for this substep
361 247412 : GenericRankTwoTensor<is_ad> sub_strain_increment = strain_increment_per_step;
362 624712 : sub_stress_new += elasticity_tensor * sub_strain_increment;
363 :
364 : Real effective_sub_stress_new;
365 : if constexpr (!is_ad)
366 : {
367 : // compute effective_sub_stress_new
368 247412 : const RankTwoTensor deviatoric_sub_stress_new = sub_stress_new.deviatoric();
369 : const Real dev_sub_stress_new_squared =
370 247412 : deviatoric_sub_stress_new.doubleContraction(deviatoric_sub_stress_new);
371 247412 : effective_sub_stress_new = std::sqrt(3.0 / 2.0 * dev_sub_stress_new_squared);
372 : }
373 : else
374 : libmesh_ignore(effective_sub_stress_new);
375 :
376 : // update stress and strain based on the strain increment
377 624712 : updateState(sub_strain_increment,
378 : sub_inelastic_strain_increment,
379 : rotation_increment, // not used in updateState
380 : sub_stress_new,
381 : stress_old, // not used in updateState
382 : elasticity_tensor,
383 : elastic_strain_old,
384 : false);
385 : // do not compute tangent until the end of this substep (or not at all for is_ad == true)
386 :
387 : // update strain and stress
388 624704 : strain_increment += sub_strain_increment;
389 624704 : inelastic_strain_increment += sub_inelastic_strain_increment;
390 624704 : sub_elastic_strain_old += sub_strain_increment;
391 624708 : sub_stress_new = elasticity_tensor * sub_elastic_strain_old;
392 :
393 : // accumulate scalar_effective_inelastic_strain
394 624704 : sub_effective_inelastic_strain_increment += _effective_inelastic_strain_increment;
395 :
396 : if constexpr (!is_ad)
397 247408 : computeTangentOperator(effective_sub_stress_new, sub_stress_new, tangent_operator);
398 :
399 : // store incremental material properties for this step
400 624704 : this->storeIncrementalMaterialProperties(total_number_substeps);
401 : }
402 :
403 : // update stress
404 8896 : stress_new = sub_stress_new;
405 :
406 : // update effective inelastic strain
407 26560 : updateEffectiveInelasticStrain(sub_effective_inelastic_strain_increment);
408 : }
409 :
410 : template <bool is_ad>
411 : void
412 47908 : RadialReturnStressUpdateTempl<is_ad>::updateStateSubstep(
413 : GenericRankTwoTensor<is_ad> & strain_increment,
414 : GenericRankTwoTensor<is_ad> & inelastic_strain_increment,
415 : const GenericRankTwoTensor<is_ad> & rotation_increment,
416 : GenericRankTwoTensor<is_ad> & stress_new,
417 : const RankTwoTensor & stress_old,
418 : const GenericRankFourTensor<is_ad> & elasticity_tensor,
419 : const RankTwoTensor & elastic_strain_old,
420 : bool compute_full_tangent_operator,
421 : RankFourTensor & tangent_operator)
422 : {
423 47908 : unsigned int num_substeps = calculateNumberSubsteps(strain_increment);
424 47908 : _dt_original = _dt;
425 : while (true)
426 : {
427 : try
428 : {
429 47916 : updateStateSubstepInternal(strain_increment,
430 : inelastic_strain_increment,
431 : rotation_increment,
432 : stress_new,
433 : stress_old,
434 : elasticity_tensor,
435 : elastic_strain_old,
436 : num_substeps,
437 : compute_full_tangent_operator,
438 : tangent_operator);
439 : }
440 176 : catch (MooseException & e)
441 : {
442 : // if we are not using adaptive substepping we just rethrow the exception
443 92 : if (!_adaptive_substepping)
444 76 : throw e;
445 :
446 : // otherwise we double the number of substeps and try again
447 16 : num_substeps *= 2;
448 16 : if (num_substeps <= _maximum_number_substeps)
449 : continue;
450 :
451 : // too meany substeps, break out of the loop
452 : break;
453 : }
454 :
455 : // updateStateSubstepInternal was successful (didn't throw)
456 47824 : _dt = _dt_original;
457 47824 : return;
458 : }
459 :
460 : // recover the original timestep
461 8 : _dt = _dt_original;
462 :
463 8 : mooseException("Adaptive substepping failed. Maximum number of substeps exceeded.");
464 : }
465 :
466 : template <bool is_ad>
467 : Real
468 291524968 : RadialReturnStressUpdateTempl<is_ad>::computeReferenceResidual(
469 : const GenericReal<is_ad> & effective_trial_stress,
470 : const GenericReal<is_ad> & scalar_effective_inelastic_strain)
471 : {
472 291524968 : return MetaPhysicL::raw_value(effective_trial_stress / _three_shear_modulus) -
473 291524968 : MetaPhysicL::raw_value(scalar_effective_inelastic_strain);
474 : }
475 :
476 : template <bool is_ad>
477 : GenericReal<is_ad>
478 56112772 : RadialReturnStressUpdateTempl<is_ad>::maximumPermissibleValue(
479 : const GenericReal<is_ad> & effective_trial_stress) const
480 : {
481 56112772 : return effective_trial_stress / _three_shear_modulus;
482 : }
483 :
484 : template <bool is_ad>
485 : Real
486 29877152 : RadialReturnStressUpdateTempl<is_ad>::computeTimeStepLimit()
487 : {
488 : const Real scalar_inelastic_strain_incr =
489 29877152 : std::abs(MetaPhysicL::raw_value(_effective_inelastic_strain[_qp]) -
490 29877152 : _effective_inelastic_strain_old[_qp]);
491 29877152 : if (!scalar_inelastic_strain_incr)
492 : return std::numeric_limits<Real>::max();
493 :
494 25092984 : return _dt * _max_inelastic_increment / scalar_inelastic_strain_incr;
495 : }
496 :
497 : template <bool is_ad>
498 : void
499 25364 : RadialReturnStressUpdateTempl<is_ad>::outputIterationSummary(std::stringstream * iter_output,
500 : const unsigned int total_it)
501 : {
502 25364 : if (iter_output)
503 : {
504 25364 : *iter_output << "At element " << _current_elem->id() << " _qp=" << _qp << " Coordinates "
505 25364 : << _q_point[_qp] << " block=" << _current_elem->subdomain_id() << '\n';
506 : }
507 25364 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationSummary(iter_output, total_it);
508 25364 : }
509 :
510 : template class RadialReturnStressUpdateTempl<false>;
511 : template class RadialReturnStressUpdateTempl<true>;
|
