21 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.");
28 "The maximum inelastic strain increment allowed in a time step");
30 "effective_inelastic_strain_name",
31 "Name of the material property that stores the effective inelastic strain");
32 params.
addParam<
bool>(
"use_substep",
false,
"Whether to use substepping");
34 MooseEnum substeppingType(
"NONE ERROR_BASED INCREMENT_BASED",
"NONE");
35 substeppingType.addDocumentation(
"NONE",
"Do not use substepping");
36 substeppingType.addDocumentation(
38 "Use substepping with a substep size that will yield, at most, the creep numerical " 39 "integration error given by substep_strain_tolerance.");
40 substeppingType.addDocumentation(
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.");
45 "use_substepping", substeppingType,
"Whether and how to use substepping");
47 "adaptive_substepping",
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. ");
53 "use_substep",
false,
"Whether to use substepping",
"Use `use_substepping` instead");
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 params.
addParam<
bool>(
"apply_strain",
true,
"Flag to apply strain. Used for testing.");
61 "effective_inelastic_strain_name substep_strain_tolerance apply_strain",
"Advanced");
62 params.
addParam<
bool>(
"use_substep_integration_error",
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 params.
addParam<
unsigned>(
"maximum_number_substeps",
68 "The maximum number of substeps allowed before cutting the time step.");
77 _effective_inelastic_strain(this->template declareGenericProperty<
Real, is_ad>(
79 this->template getParam<
std::string>(
"effective_inelastic_strain_name"))),
80 _effective_inelastic_strain_old(this->template getMaterialPropertyOld<
Real>(
82 this->template getParam<
std::string>(
"effective_inelastic_strain_name"))),
83 _max_inelastic_increment(this->template getParam<
Real>(
"max_inelastic_increment")),
84 _substep_tolerance(this->template getParam<
Real>(
"substep_strain_tolerance")),
86 _identity_symmetric_four(
RankFourTensor::initIdentitySymmetricFour),
87 _deviatoric_projection_four(_identity_symmetric_four -
88 _identity_two.outerProduct(_identity_two) / 3.0),
89 _apply_strain(this->template getParam<bool>(
"apply_strain")),
92 _adaptive_substepping(this->template getParam<bool>(
"adaptive_substepping")),
93 _maximum_number_substeps(this->template getParam<unsigned>(
"maximum_number_substeps"))
99 "use_substep",
"Remove this parameter and just keep `use_substepping` in the input");
101 if (parameters.
get<
bool>(
"use_substep"))
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");
119 "adaptive_substepping",
120 "The parameter adaptive_substepping can only be used when the substepping option " 121 "(use_substepping) is not set to NONE");
124 template <
bool is_ad>
128 _effective_inelastic_strain[_qp] = 0.0;
131 template <
bool is_ad>
141 template <
bool is_ad>
145 _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
148 template <
bool is_ad>
155 strain_increment.doubleContraction(strain_increment);
156 const Real effective_elastic_strain =
162 if (_use_substepping == SubsteppingType::INCREMENT_BASED)
164 const Real ratio = effective_elastic_strain / _max_inelastic_increment;
166 if (ratio > _substep_tolerance)
167 return std::ceil(ratio / _substep_tolerance);
171 if (_use_substepping == SubsteppingType::ERROR_BASED)
173 const Real accurate_time_step_ratio = _substep_tolerance / effective_elastic_strain;
175 if (accurate_time_step_ratio < 1.0)
176 return std::ceil(1.0 / accurate_time_step_ratio);
180 mooseError(
"calculateNumberSubsteps should not have been called. Notify a developer.");
183 template <
bool is_ad>
189 mooseError(
"computeTangentOperator called: no tangent computation is needed for AD");
201 tangent_operator.
zero();
211 mooseAssert(_three_shear_modulus != 0.0,
"Shear modulus is zero");
217 tangent_operator.
zero();
220 const Real norm_dev_stress =
std::sqrt(norm_dev_stress_squared);
222 const RankTwoTensor flow_direction = deviatoric_stress / norm_dev_stress;
225 computeStressDerivative(effective_trial_stress, _effective_inelastic_strain_increment);
226 const Real scalar_one = _three_shear_modulus * _effective_inelastic_strain_increment /
229 tangent_operator = scalar_one * _deviatoric_projection_four +
230 (_three_shear_modulus *
deriv - scalar_one) * flow_direction_dyad;
235 template <
bool is_ad>
245 bool compute_full_tangent_operator,
255 deviatoric_trial_stress.doubleContraction(deviatoric_trial_stress);
257 ?
std::sqrt(3.0 / 2.0 * dev_trial_stress_squared)
263 _three_shear_modulus != 0.0,
264 "Shear modulus is zero. Ensure that the base class computeStressInitialize() is called.");
267 _effective_inelastic_strain_increment = 0.0;
270 this->returnMappingSolve(
271 effective_trial_stress, _effective_inelastic_strain_increment, this->_console);
272 if (_effective_inelastic_strain_increment != 0.0)
273 inelastic_strain_increment =
274 deviatoric_trial_stress *
275 (1.5 * _effective_inelastic_strain_increment / effective_trial_stress);
277 inelastic_strain_increment.zero();
280 inelastic_strain_increment.zero();
284 strain_increment -= inelastic_strain_increment;
285 updateEffectiveInelasticStrain(_effective_inelastic_strain_increment);
293 computeStressFinalize(inelastic_strain_increment);
295 if constexpr (!is_ad)
297 if (compute_full_tangent_operator)
298 computeTangentOperator(effective_trial_stress, stress_new, tangent_operator);
307 template <
bool is_ad>
317 unsigned int total_number_substeps,
318 bool compute_full_tangent_operator,
322 if (total_number_substeps == 1)
324 updateState(strain_increment,
325 inelastic_strain_increment,
331 compute_full_tangent_operator,
334 this->storeIncrementalMaterialProperties(total_number_substeps);
338 if (total_number_substeps > _maximum_number_substeps)
339 mooseException(
"The number of substeps computed exceeds 'maximum_number_substeps'.");
342 _dt = _dt_original / total_number_substeps;
346 strain_increment / total_number_substeps;
358 for (
unsigned int step = 0; step < total_number_substeps; ++step)
364 Real effective_sub_stress_new;
365 if constexpr (!is_ad)
369 const Real dev_sub_stress_new_squared =
371 effective_sub_stress_new =
std::sqrt(3.0 / 2.0 * dev_sub_stress_new_squared);
377 updateState(sub_strain_increment,
378 sub_inelastic_strain_increment,
388 strain_increment += sub_strain_increment;
389 inelastic_strain_increment += sub_inelastic_strain_increment;
390 sub_elastic_strain_old += sub_strain_increment;
394 sub_effective_inelastic_strain_increment += _effective_inelastic_strain_increment;
396 if constexpr (!is_ad)
397 computeTangentOperator(effective_sub_stress_new, sub_stress_new, tangent_operator);
400 this->storeIncrementalMaterialProperties(total_number_substeps);
404 stress_new = sub_stress_new;
407 updateEffectiveInelasticStrain(sub_effective_inelastic_strain_increment);
410 template <
bool is_ad>
420 bool compute_full_tangent_operator,
423 unsigned int num_substeps = calculateNumberSubsteps(strain_increment);
429 updateStateSubstepInternal(strain_increment,
430 inelastic_strain_increment,
437 compute_full_tangent_operator,
443 if (!_adaptive_substepping)
448 if (num_substeps <= _maximum_number_substeps)
463 mooseException(
"Adaptive substepping failed. Maximum number of substeps exceeded.");
466 template <
bool is_ad>
476 template <
bool is_ad>
481 return effective_trial_stress / _three_shear_modulus;
484 template <
bool is_ad>
488 const Real scalar_inelastic_strain_incr =
490 _effective_inelastic_strain_old[_qp]);
491 if (!scalar_inelastic_strain_incr)
492 return std::numeric_limits<Real>::max();
494 return _dt * _max_inelastic_increment / scalar_inelastic_strain_incr;
497 template <
bool is_ad>
500 const unsigned int total_it)
504 *iter_output <<
"At element " << _current_elem->id() <<
" _qp=" << _qp <<
" Coordinates " 505 << _q_point[_qp] <<
" block=" << _current_elem->subdomain_id() <<
'\n';
RankFourTensorTempl< Real > outerProduct(const RankTwoTensorTempl< Real > &b) const
SubsteppingType _use_substepping
Whether user has requested the use of substepping technique to improve convergence [make const later]...
typename Moose::GenericType< RankFourTensor, is_ad > GenericRankFourTensor
void mooseSetToZero(T &v)
void outputIterationSummary(std::stringstream *iter_output, const unsigned int total_it) override
Output summary information for the convergence history of the model.
virtual void updateState(GenericRankTwoTensor< is_ad > &strain_increment, GenericRankTwoTensor< is_ad > &inelastic_strain_increment, const GenericRankTwoTensor< is_ad > &rotation_increment, GenericRankTwoTensor< is_ad > &stress_new, const RankTwoTensor &stress_old, const GenericRankFourTensor< is_ad > &elasticity_tensor, const RankTwoTensor &elastic_strain_old, bool compute_full_tangent_operator=false, RankFourTensor &tangent_operator=StressUpdateBaseTempl< is_ad >::_identityTensor) override
A radial return (J2) mapping method is performed with return mapping iterations.
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
static InputParameters validParams()
void mooseError(Args &&... args)
virtual void computeStressInitialize(const GenericReal< is_ad > &effective_trial_stress, const GenericRankFourTensor< is_ad > &elasticity_tensor)
Perform any necessary initialization before return mapping iterations.
RadialReturnStressUpdate computes the radial return stress increment for an isotropic elastic-viscopl...
virtual int calculateNumberSubsteps(const GenericRankTwoTensor< is_ad > &strain_increment) override
If substepping is enabled, calculate the number of substeps as a function of the elastic strain incre...
static InputParameters validParams()
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
virtual void updateStateSubstepInternal(GenericRankTwoTensor< is_ad > &, GenericRankTwoTensor< is_ad > &, const GenericRankTwoTensor< is_ad > &, GenericRankTwoTensor< is_ad > &, const RankTwoTensor &, const GenericRankFourTensor< is_ad > &, const RankTwoTensor &, unsigned int total_number_substeps, bool compute_full_tangent_operator=false, RankFourTensor &tangent_operator=StressUpdateBaseTempl< is_ad >::_identityTensor)
static InputParameters validParams()
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
RankTwoTensorTempl< Real > deviatoric() const
void libmesh_ignore(const Args &...)
virtual void initQpStatefulProperties() override
Real deriv(unsigned n, unsigned alpha, unsigned beta, Real x)
virtual void outputIterationSummary(std::stringstream *iter_output, const unsigned int total_it)
Output summary information for the convergence history of the model.
Base class that provides capability for Newton return mapping iterations on a single variable...
Real doubleContraction(const RankTwoTensorTempl< Real > &a) const
virtual Real computeTimeStepLimit() override
Compute the limiting value of the time step for this material.
void paramError(const std::string ¶m, Args... args) const
StressUpdateBase is a material that is not called by MOOSE because of the compute=false flag set in t...
void computeTangentOperator(Real effective_trial_stress, const RankTwoTensor &stress_new, RankFourTensor &tangent_operator)
Calculate the tangent_operator.
virtual GenericReal< is_ad > maximumPermissibleValue(const GenericReal< is_ad > &effective_trial_stress) const override
Compute the maximum permissible value of the scalar.
RadialReturnStressUpdateTempl(const InputParameters ¶meters)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void propagateQpStatefulPropertiesRadialReturn()
Propagate the properties pertaining to this intermediate class.
const InputParameters & _pars
typename Moose::GenericType< RankTwoTensor, is_ad > GenericRankTwoTensor
const InputParameters & parameters() const
typename Moose::GenericType< Real, is_ad > GenericReal
virtual void updateStateSubstep(GenericRankTwoTensor< is_ad > &, GenericRankTwoTensor< is_ad > &, const GenericRankTwoTensor< is_ad > &, GenericRankTwoTensor< is_ad > &, const RankTwoTensor &, const GenericRankFourTensor< is_ad > &, const RankTwoTensor &, bool compute_full_tangent_operator=false, RankFourTensor &tangent_operator=StressUpdateBaseTempl< is_ad >::_identityTensor) override
Similar to the updateState function, this method updates the strain and stress for one substep...
const bool _adaptive_substepping
Use adaptive substepping, cutting substep sizes until convergence is achieved.
virtual Real computeReferenceResidual(const GenericReal< is_ad > &effective_trial_stress, const GenericReal< is_ad > &scalar_effective_inelastic_strain) override
Compute a reference quantity to be used for checking relative convergence.