Line data Source code
1 :
2 : //* This file is part of the MOOSE framework
3 : //* https://www.mooseframework.org
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 2178 : RadialReturnStressUpdateTempl<is_ad>::validParams()
19 : {
20 2178 : InputParameters params = StressUpdateBaseTempl<is_ad>::validParams();
21 2178 : 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 2178 : params += SingleVariableReturnMappingSolutionTempl<is_ad>::validParams();
26 4356 : params.addParam<Real>("max_inelastic_increment",
27 4356 : 1e-4,
28 : "The maximum inelastic strain increment allowed in a time step");
29 4356 : params.addRequiredParam<std::string>(
30 : "effective_inelastic_strain_name",
31 : "Name of the material property that stores the effective inelastic strain");
32 4356 : params.addParam<bool>("use_substep", false, "Whether to use substepping");
33 :
34 4356 : MooseEnum substeppingType("NONE ERROR_BASED INCREMENT_BASED", "NONE");
35 4356 : substeppingType.addDocumentation("NONE", "Do not use substepping");
36 4356 : 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 4356 : 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 4356 : params.addParam<MooseEnum>(
45 : "use_substepping", substeppingType, "Whether and how to use substepping");
46 4356 : params.addParam<bool>(
47 : "adaptive_substepping",
48 4356 : 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 6534 : params.addDeprecatedParam<bool>(
53 4356 : "use_substep", false, "Whether to use substepping", "Use `use_substepping` instead");
54 4356 : params.addParam<Real>("substep_strain_tolerance",
55 4356 : 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 4356 : params.addParam<bool>("apply_strain", true, "Flag to apply strain. Used for testing.");
60 4356 : params.addParamNamesToGroup(
61 : "effective_inelastic_strain_name substep_strain_tolerance apply_strain", "Advanced");
62 4356 : params.addParam<bool>("use_substep_integration_error",
63 4356 : 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 4356 : params.addParam<unsigned>("maximum_number_substeps",
67 4356 : 25,
68 : "The maximum number of substeps allowed before cutting the time step.");
69 2178 : return params;
70 2178 : }
71 :
72 : template <bool is_ad>
73 1634 : RadialReturnStressUpdateTempl<is_ad>::RadialReturnStressUpdateTempl(
74 : const InputParameters & parameters)
75 : : StressUpdateBaseTempl<is_ad>(parameters),
76 : SingleVariableReturnMappingSolutionTempl<is_ad>(parameters),
77 2322 : _effective_inelastic_strain(this->template declareGenericProperty<Real, is_ad>(
78 1634 : this->_base_name +
79 : this->template getParam<std::string>("effective_inelastic_strain_name"))),
80 5590 : _effective_inelastic_strain_old(this->template getMaterialPropertyOld<Real>(
81 : this->_base_name +
82 : this->template getParam<std::string>("effective_inelastic_strain_name"))),
83 3268 : _max_inelastic_increment(this->template getParam<Real>("max_inelastic_increment")),
84 3268 : _substep_tolerance(this->template getParam<Real>("substep_strain_tolerance")),
85 1634 : _identity_two(RankTwoTensor::initIdentity),
86 1634 : _identity_symmetric_four(RankFourTensor::initIdentitySymmetricFour),
87 1634 : _deviatoric_projection_four(_identity_symmetric_four -
88 1634 : _identity_two.outerProduct(_identity_two) / 3.0),
89 3268 : _apply_strain(this->template getParam<bool>("apply_strain")),
90 1634 : _use_substepping(
91 1634 : this->template getParam<MooseEnum>("use_substepping").template getEnum<SubsteppingType>()),
92 3268 : _adaptive_substepping(this->template getParam<bool>("adaptive_substepping")),
93 4902 : _maximum_number_substeps(this->template getParam<unsigned>("maximum_number_substeps"))
94 : {
95 3268 : if (this->_pars.isParamSetByUser("use_substep"))
96 : {
97 6 : if (this->_pars.isParamSetByUser("use_substepping"))
98 0 : this->template paramError(
99 : "use_substep", "Remove this parameter and just keep `use_substepping` in the input");
100 :
101 3 : if (parameters.get<bool>("use_substep"))
102 : {
103 3 : if (parameters.get<bool>("use_substep_integration_error"))
104 0 : _use_substepping = SubsteppingType::ERROR_BASED;
105 : else
106 3 : _use_substepping = SubsteppingType::INCREMENT_BASED;
107 : }
108 : }
109 :
110 3268 : if (this->_pars.isParamSetByUser("maximum_number_substeps") &&
111 30 : _use_substepping == SubsteppingType::NONE)
112 0 : this->template 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 1634 : if (_adaptive_substepping && _use_substepping == SubsteppingType::NONE)
118 0 : this->template 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 1634 : }
123 :
124 : template <bool is_ad>
125 : void
126 12776 : RadialReturnStressUpdateTempl<is_ad>::initQpStatefulProperties()
127 : {
128 32816 : _effective_inelastic_strain[_qp] = 0.0;
129 12776 : }
130 :
131 : template <bool is_ad>
132 : void
133 13718532 : 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 40918796 : _three_shear_modulus = 3.0 * ElasticityTensorTools::getIsotropicShearModulus(elasticity_tensor);
139 13718532 : }
140 :
141 : template <bool is_ad>
142 : void
143 260784 : RadialReturnStressUpdateTempl<is_ad>::propagateQpStatefulPropertiesRadialReturn()
144 : {
145 260784 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
146 260784 : }
147 :
148 : template <bool is_ad>
149 : int
150 29618 : RadialReturnStressUpdateTempl<is_ad>::calculateNumberSubsteps(
151 : const GenericRankTwoTensor<is_ad> & strain_increment)
152 : {
153 : // compute an effective elastic strain measure
154 13085 : const GenericReal<is_ad> contracted_elastic_strain =
155 16533 : strain_increment.doubleContraction(strain_increment);
156 : const Real effective_elastic_strain =
157 29618 : std::sqrt(3.0 / 2.0 * MetaPhysicL::raw_value(contracted_elastic_strain));
158 :
159 29618 : if (MooseUtils::absoluteFuzzyEqual(effective_elastic_strain, 0.0))
160 : return 1;
161 :
162 16690 : if (_use_substepping == SubsteppingType::INCREMENT_BASED)
163 : {
164 11175 : const Real ratio = effective_elastic_strain / _max_inelastic_increment;
165 :
166 11175 : if (ratio > _substep_tolerance)
167 11172 : return std::ceil(ratio / _substep_tolerance);
168 : return 1;
169 : }
170 :
171 5515 : if (_use_substepping == SubsteppingType::ERROR_BASED)
172 : {
173 5515 : const Real accurate_time_step_ratio = _substep_tolerance / effective_elastic_strain;
174 :
175 5515 : if (accurate_time_step_ratio < 1.0)
176 2146 : 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 297252 : RadialReturnStressUpdateTempl<false>::computeTangentOperator(Real effective_trial_stress,
195 : const RankTwoTensor & stress_new,
196 : RankFourTensor & tangent_operator)
197 : {
198 297252 : if (getTangentCalculationMethod() == TangentCalculationMethod::PARTIAL)
199 : {
200 255918 : if (MooseUtils::absoluteFuzzyEqual(_effective_inelastic_strain_increment, 0.0))
201 33666 : 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 222252 : const RankTwoTensor deviatoric_stress = stress_new.deviatoric();
214 222252 : const Real norm_dev_stress_squared = deviatoric_stress.doubleContraction(deviatoric_stress);
215 222252 : if (MooseUtils::absoluteFuzzyEqual(norm_dev_stress_squared, 0.0))
216 : {
217 0 : tangent_operator.zero();
218 : return;
219 : }
220 222252 : const Real norm_dev_stress = std::sqrt(norm_dev_stress_squared);
221 :
222 222252 : 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 222252 : computeStressDerivative(effective_trial_stress, _effective_inelastic_strain_increment);
226 222252 : const Real scalar_one = _three_shear_modulus * _effective_inelastic_strain_increment /
227 222252 : std::sqrt(1.5) / norm_dev_stress;
228 :
229 222252 : tangent_operator = scalar_one * _deviatoric_projection_four +
230 444504 : (_three_shear_modulus * deriv - scalar_one) * flow_direction_dyad;
231 : }
232 : }
233 : }
234 :
235 : template <bool is_ad>
236 : void
237 27158930 : 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 27158930 : GenericRankTwoTensor<is_ad> deviatoric_trial_stress = stress_new.deviatoric();
252 :
253 : // compute the effective trial stress
254 13718532 : GenericReal<is_ad> dev_trial_stress_squared =
255 13440398 : deviatoric_trial_stress.doubleContraction(deviatoric_trial_stress);
256 26750822 : GenericReal<is_ad> effective_trial_stress = MetaPhysicL::raw_value(dev_trial_stress_squared)
257 27158930 : ? std::sqrt(3.0 / 2.0 * dev_trial_stress_squared)
258 13718532 : : 0.0;
259 :
260 27158930 : 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 27158928 : _effective_inelastic_strain_increment = 0.0;
268 27158928 : if (!MooseUtils::absoluteFuzzyEqual(effective_trial_stress, 0.0))
269 : {
270 26306692 : this->returnMappingSolve(
271 26306692 : effective_trial_stress, _effective_inelastic_strain_increment, this->_console);
272 26306678 : if (_effective_inelastic_strain_increment != 0.0)
273 37772736 : inelastic_strain_increment =
274 : deviatoric_trial_stress *
275 38293437 : (1.5 * _effective_inelastic_strain_increment / effective_trial_stress);
276 : else
277 372369 : inelastic_strain_increment.zero();
278 : }
279 : else
280 408108 : inelastic_strain_increment.zero();
281 :
282 27158914 : if (_apply_strain)
283 : {
284 27151202 : strain_increment -= inelastic_strain_increment;
285 27151202 : 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 27151202 : stress_new = elasticity_tensor * (strain_increment + elastic_strain_old);
291 : }
292 :
293 27158914 : computeStressFinalize(inelastic_strain_increment);
294 :
295 : if constexpr (!is_ad)
296 : {
297 13440391 : if (compute_full_tangent_operator)
298 90880 : 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 27158914 : }
306 :
307 : template <bool is_ad>
308 : void
309 29622 : 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 29622 : if (total_number_substeps == 1)
323 : {
324 16300 : 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 16296 : this->storeIncrementalMaterialProperties(total_number_substeps);
335 16296 : return;
336 : }
337 :
338 13322 : if (total_number_substeps > _maximum_number_substeps)
339 38 : mooseException("The number of substeps computed exceeds the maximum_number_substeps. The "
340 : "system time step will be cut.");
341 :
342 : // cut the original timestep
343 13284 : _dt = _dt_original / total_number_substeps;
344 :
345 : // initialize the inputs
346 13284 : const GenericRankTwoTensor<is_ad> strain_increment_per_step =
347 : strain_increment / total_number_substeps;
348 13284 : GenericRankTwoTensor<is_ad> sub_stress_new = elasticity_tensor * elastic_strain_old;
349 4450 : GenericRankTwoTensor<is_ad> sub_elastic_strain_old = elastic_strain_old;
350 :
351 : // clear the original inputs
352 13284 : MathUtils::mooseSetToZero(strain_increment);
353 13284 : MathUtils::mooseSetToZero(inelastic_strain_increment);
354 13284 : MathUtils::mooseSetToZero(stress_new);
355 :
356 13284 : GenericReal<is_ad> sub_effective_inelastic_strain_increment = 0.0;
357 4450 : GenericRankTwoTensor<is_ad> sub_inelastic_strain_increment = inelastic_strain_increment;
358 :
359 325636 : for (unsigned int step = 0; step < total_number_substeps; ++step)
360 : {
361 : // set up input for this substep
362 123706 : GenericRankTwoTensor<is_ad> sub_strain_increment = strain_increment_per_step;
363 312356 : sub_stress_new += elasticity_tensor * sub_strain_increment;
364 :
365 : Real effective_sub_stress_new;
366 : if constexpr (!is_ad)
367 : {
368 : // compute effective_sub_stress_new
369 123706 : const RankTwoTensor deviatoric_sub_stress_new = sub_stress_new.deviatoric();
370 : const Real dev_sub_stress_new_squared =
371 123706 : deviatoric_sub_stress_new.doubleContraction(deviatoric_sub_stress_new);
372 123706 : effective_sub_stress_new = std::sqrt(3.0 / 2.0 * dev_sub_stress_new_squared);
373 : }
374 : else
375 : libmesh_ignore(effective_sub_stress_new);
376 :
377 : // update stress and strain based on the strain increment
378 312356 : updateState(sub_strain_increment,
379 : sub_inelastic_strain_increment,
380 : rotation_increment, // not used in updateState
381 : sub_stress_new,
382 : stress_old, // not used in updateState
383 : elasticity_tensor,
384 : elastic_strain_old,
385 : false);
386 : // do not compute tangent until the end of this substep (or not at all for is_ad == true)
387 :
388 : // update strain and stress
389 312352 : strain_increment += sub_strain_increment;
390 312352 : inelastic_strain_increment += sub_inelastic_strain_increment;
391 312352 : sub_elastic_strain_old += sub_strain_increment;
392 312354 : sub_stress_new = elasticity_tensor * sub_elastic_strain_old;
393 :
394 : // accumulate scalar_effective_inelastic_strain
395 312352 : sub_effective_inelastic_strain_increment += _effective_inelastic_strain_increment;
396 :
397 : if constexpr (!is_ad)
398 123704 : computeTangentOperator(effective_sub_stress_new, sub_stress_new, tangent_operator);
399 :
400 : // store incremental material properties for this step
401 312352 : this->storeIncrementalMaterialProperties(total_number_substeps);
402 : }
403 :
404 : // update stress
405 13280 : stress_new = sub_stress_new;
406 :
407 : // update effective inelastic strain
408 13280 : updateEffectiveInelasticStrain(sub_effective_inelastic_strain_increment);
409 : }
410 :
411 : template <bool is_ad>
412 : void
413 29618 : RadialReturnStressUpdateTempl<is_ad>::updateStateSubstep(
414 : GenericRankTwoTensor<is_ad> & strain_increment,
415 : GenericRankTwoTensor<is_ad> & inelastic_strain_increment,
416 : const GenericRankTwoTensor<is_ad> & rotation_increment,
417 : GenericRankTwoTensor<is_ad> & stress_new,
418 : const RankTwoTensor & stress_old,
419 : const GenericRankFourTensor<is_ad> & elasticity_tensor,
420 : const RankTwoTensor & elastic_strain_old,
421 : bool compute_full_tangent_operator,
422 : RankFourTensor & tangent_operator)
423 : {
424 29618 : unsigned int num_substeps = calculateNumberSubsteps(strain_increment);
425 29618 : _dt_original = _dt;
426 : while (true)
427 : {
428 : try
429 : {
430 29622 : updateStateSubstepInternal(strain_increment,
431 : inelastic_strain_increment,
432 : rotation_increment,
433 : stress_new,
434 : stress_old,
435 : elasticity_tensor,
436 : elastic_strain_old,
437 : num_substeps,
438 : compute_full_tangent_operator,
439 : tangent_operator);
440 : }
441 88 : catch (MooseException & e)
442 : {
443 : // if we are not using adaptive substepping we just rethrow the exception
444 46 : if (!_adaptive_substepping)
445 38 : throw e;
446 :
447 : // otherwise we double the number of substeps and try again
448 8 : num_substeps *= 2;
449 8 : if (num_substeps <= _maximum_number_substeps)
450 : continue;
451 :
452 : // too meany substeps, break out of the loop
453 : break;
454 : }
455 :
456 : // updateStateSubstepInternal was successful (didn't throw)
457 29576 : _dt = _dt_original;
458 29576 : return;
459 : }
460 :
461 : // recover the original timestep
462 4 : _dt = _dt_original;
463 :
464 4 : mooseException("Adaptive substepping failed. Maximum number of substeps exceeded.");
465 : }
466 :
467 : template <bool is_ad>
468 : Real
469 132524329 : RadialReturnStressUpdateTempl<is_ad>::computeReferenceResidual(
470 : const GenericReal<is_ad> & effective_trial_stress,
471 : const GenericReal<is_ad> & scalar_effective_inelastic_strain)
472 : {
473 132524329 : return MetaPhysicL::raw_value(effective_trial_stress / _three_shear_modulus) -
474 132524329 : MetaPhysicL::raw_value(scalar_effective_inelastic_strain);
475 : }
476 :
477 : template <bool is_ad>
478 : GenericReal<is_ad>
479 26162144 : RadialReturnStressUpdateTempl<is_ad>::maximumPermissibleValue(
480 : const GenericReal<is_ad> & effective_trial_stress) const
481 : {
482 26162144 : return effective_trial_stress / _three_shear_modulus;
483 : }
484 :
485 : template <bool is_ad>
486 : Real
487 13026140 : RadialReturnStressUpdateTempl<is_ad>::computeTimeStepLimit()
488 : {
489 : const Real scalar_inelastic_strain_incr =
490 13026140 : std::abs(MetaPhysicL::raw_value(_effective_inelastic_strain[_qp]) -
491 13026140 : _effective_inelastic_strain_old[_qp]);
492 13026140 : if (!scalar_inelastic_strain_incr)
493 : return std::numeric_limits<Real>::max();
494 :
495 11125191 : return _dt * _max_inelastic_increment / scalar_inelastic_strain_incr;
496 : }
497 :
498 : template <bool is_ad>
499 : void
500 12682 : RadialReturnStressUpdateTempl<is_ad>::outputIterationSummary(std::stringstream * iter_output,
501 : const unsigned int total_it)
502 : {
503 12682 : if (iter_output)
504 : {
505 12682 : *iter_output << "At element " << _current_elem->id() << " _qp=" << _qp << " Coordinates "
506 12682 : << _q_point[_qp] << " block=" << _current_elem->subdomain_id() << '\n';
507 : }
508 12682 : SingleVariableReturnMappingSolutionTempl<is_ad>::outputIterationSummary(iter_output, total_it);
509 12682 : }
510 :
511 : template class RadialReturnStressUpdateTempl<false>;
512 : template class RadialReturnStressUpdateTempl<true>;
|
