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 "GeneralizedRadialReturnStressUpdate.h"
11 :
12 : #include "MooseMesh.h"
13 : #include "MooseTypes.h"
14 : #include "ElasticityTensorTools.h"
15 : #include "libmesh/ignore_warnings.h"
16 : #include "Eigen/Dense"
17 : #include "Eigen/Eigenvalues"
18 : #include "libmesh/restore_warnings.h"
19 :
20 : template <bool is_ad>
21 : InputParameters
22 424 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::validParams()
23 : {
24 424 : InputParameters params = StressUpdateBaseTempl<is_ad>::validParams();
25 424 : params.addClassDescription("Calculates the effective inelastic strain increment required to "
26 : "return the isotropic stress state to a J2 yield surface. This class "
27 : "is intended to be a parent class for classes with specific "
28 : "constitutive models.");
29 424 : params += GeneralizedReturnMappingSolutionTempl<is_ad>::validParams();
30 848 : params.addParam<Real>("max_inelastic_increment",
31 848 : 1e-4,
32 : "The maximum inelastic strain increment allowed in a time step");
33 848 : params.addParam<Real>("max_integration_error",
34 848 : 5e-4,
35 : "The maximum inelastic strain increment integration error allowed");
36 848 : params.addRequiredParam<std::string>(
37 : "effective_inelastic_strain_name",
38 : "Name of the material property that stores the effective inelastic strain");
39 848 : params.addRequiredParam<std::string>(
40 : "inelastic_strain_rate_name",
41 : "Name of the material property that stores the inelastic strain rate");
42 848 : params.addParam<bool>(
43 : "use_transformation",
44 848 : true,
45 : "Whether to employ updated Hill's tensor due to rigid body or large "
46 : "deformation kinematic rotations. If an initial rigid body rotation is provided by the user "
47 : "in increments of 90 degrees (e.g. 90, 180, 270), this option can be set to false, in which "
48 : "case the Hill's coefficients are extracted from the transformed Hill's tensor.");
49 :
50 424 : return params;
51 0 : }
52 : template <bool is_ad>
53 318 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::GeneralizedRadialReturnStressUpdateTempl(
54 : const InputParameters & parameters)
55 : : StressUpdateBaseTempl<is_ad>(parameters),
56 : GeneralizedReturnMappingSolutionTempl<is_ad>(parameters),
57 588 : _effective_inelastic_strain(this->template declareGenericProperty<Real, is_ad>(
58 318 : this->_base_name +
59 : this->template getParam<std::string>("effective_inelastic_strain_name"))),
60 954 : _effective_inelastic_strain_old(this->template getMaterialPropertyOld<Real>(
61 : this->_base_name +
62 : this->template getParam<std::string>("effective_inelastic_strain_name"))),
63 636 : _inelastic_strain_rate(this->template declareProperty<Real>(
64 : this->_base_name + this->template getParam<std::string>("inelastic_strain_rate_name"))),
65 954 : _inelastic_strain_rate_old(this->template getMaterialPropertyOld<Real>(
66 : this->_base_name + this->template getParam<std::string>("inelastic_strain_rate_name"))),
67 636 : _max_inelastic_increment(this->template getParam<Real>("max_inelastic_increment")),
68 636 : _max_integration_error(this->template getParam<Real>("max_integration_error")),
69 318 : _max_integration_error_time_step(std::numeric_limits<Real>::max()),
70 954 : _use_transformation(this->template getParam<bool>("use_transformation"))
71 : {
72 318 : }
73 :
74 : template <bool is_ad>
75 : void
76 24096 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::initQpStatefulProperties()
77 : {
78 24096 : _effective_inelastic_strain[_qp] = 0.0;
79 24096 : _inelastic_strain_rate[_qp] = 0.0;
80 24096 : }
81 :
82 : template <bool is_ad>
83 : void
84 0 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::propagateQpStatefulPropertiesRadialReturn()
85 : {
86 0 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
87 0 : _inelastic_strain_rate[_qp] = _inelastic_strain_rate_old[_qp];
88 0 : }
89 :
90 : template <bool is_ad>
91 : void
92 1934240 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::updateState(
93 : GenericRankTwoTensor<is_ad> & elastic_strain_increment,
94 : GenericRankTwoTensor<is_ad> & inelastic_strain_increment,
95 : const GenericRankTwoTensor<is_ad> & /*rotation_increment*/,
96 : GenericRankTwoTensor<is_ad> & stress_new,
97 : const RankTwoTensor & stress_old,
98 : const GenericRankFourTensor<is_ad> & elasticity_tensor,
99 : const RankTwoTensor & /*elastic_strain_old*/,
100 : bool /*compute_full_tangent_operator = false*/,
101 : RankFourTensor & /*tangent_operator = StressUpdateBaseTempl<is_ad>::_identityTensor*/)
102 : {
103 : // Prepare initial trial stress for generalized return mapping
104 1934240 : GenericRankTwoTensor<is_ad> deviatoric_trial_stress = stress_new.deviatoric();
105 :
106 1934240 : GenericDenseVector<is_ad> stress_new_vector(6);
107 1934240 : stress_new_vector(0) = stress_new(0, 0);
108 1934240 : stress_new_vector(1) = stress_new(1, 1);
109 1934240 : stress_new_vector(2) = stress_new(2, 2);
110 1934240 : stress_new_vector(3) = stress_new(0, 1);
111 1934240 : stress_new_vector(4) = stress_new(1, 2);
112 1934240 : stress_new_vector(5) = stress_new(0, 2);
113 :
114 1934240 : GenericDenseVector<is_ad> stress_dev(6);
115 1934240 : stress_dev(0) = deviatoric_trial_stress(0, 0);
116 1934240 : stress_dev(1) = deviatoric_trial_stress(1, 1);
117 1934240 : stress_dev(2) = deviatoric_trial_stress(2, 2);
118 1934240 : stress_dev(3) = deviatoric_trial_stress(0, 1);
119 1934240 : stress_dev(4) = deviatoric_trial_stress(1, 2);
120 1934240 : stress_dev(5) = deviatoric_trial_stress(0, 2);
121 :
122 1934240 : computeStressInitialize(stress_dev, stress_new_vector, elasticity_tensor);
123 :
124 : // Use Newton iteration to determine a plastic multiplier variable
125 1934240 : GenericReal<is_ad> delta_gamma = 0.0;
126 :
127 : // Use Newton iteration to determine the scalar effective inelastic strain increment
128 3868480 : if (!MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(stress_dev).l2_norm(), 0.0))
129 : {
130 1888576 : this->returnMappingSolve(stress_dev, stress_new_vector, delta_gamma, this->_console);
131 :
132 1888576 : computeStrainFinalize(inelastic_strain_increment, stress_new, stress_dev, delta_gamma);
133 : }
134 : else
135 35168 : inelastic_strain_increment.zero();
136 :
137 1934240 : elastic_strain_increment -= inelastic_strain_increment;
138 :
139 1934240 : computeStressFinalize(inelastic_strain_increment,
140 : delta_gamma,
141 : stress_new,
142 : stress_dev,
143 : stress_old,
144 : elasticity_tensor);
145 1934240 : }
146 :
147 : template <bool is_ad>
148 : Real
149 0 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::computeReferenceResidual(
150 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
151 : const GenericDenseVector<is_ad> & /*stress_new*/,
152 : const GenericReal<is_ad> & /*residual*/,
153 : const GenericReal<is_ad> & /*scalar_effective_inelastic_strain*/)
154 : {
155 0 : mooseError("GeneralizedRadialReturnStressUpdate::computeReferenceResidual must be implemented "
156 : "by child classes");
157 :
158 : return 0.0;
159 : }
160 :
161 : template <bool is_ad>
162 : GenericReal<is_ad>
163 1888576 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::maximumPermissibleValue(
164 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/) const
165 : {
166 1888576 : return std::numeric_limits<Real>::max();
167 : }
168 :
169 : template <bool is_ad>
170 : Real
171 1889360 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::computeTimeStepLimit()
172 : {
173 :
174 : // Add a new criterion including numerical integration error
175 1889360 : Real scalar_inelastic_strain_incr = MetaPhysicL::raw_value(_effective_inelastic_strain[_qp]) -
176 1889360 : _effective_inelastic_strain_old[_qp];
177 :
178 1889360 : if (MooseUtils::absoluteFuzzyEqual(scalar_inelastic_strain_incr, 0.0))
179 : return std::numeric_limits<Real>::max();
180 :
181 3486784 : return std::min(_dt * _max_inelastic_increment / scalar_inelastic_strain_incr,
182 1797200 : computeIntegrationErrorTimeStep());
183 : }
184 :
185 : template <bool is_ad>
186 : void
187 0 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::outputIterationSummary(
188 : std::stringstream * iter_output, const unsigned int total_it)
189 : {
190 0 : if (iter_output)
191 : {
192 0 : *iter_output << "At element " << _current_elem->id() << " _qp=" << _qp << " Coordinates "
193 0 : << _q_point[_qp] << " block=" << _current_elem->subdomain_id() << '\n';
194 : }
195 0 : GeneralizedReturnMappingSolutionTempl<is_ad>::outputIterationSummary(iter_output, total_it);
196 0 : }
197 :
198 : template <bool is_ad>
199 : bool
200 119488 : GeneralizedRadialReturnStressUpdateTempl<is_ad>::isBlockDiagonal(const AnisotropyMatrixReal & A)
201 : {
202 119488 : AnisotropyMatrixRealBlock A_bottom_left(A.block<3, 3>(0, 3));
203 119488 : AnisotropyMatrixRealBlock A_top_right(A.block<3, 3>(3, 0));
204 :
205 135800 : for (unsigned int index_i = 0; index_i < 3; index_i++)
206 186842 : for (unsigned int index_j = 0; index_j < 3; index_j++)
207 170530 : if (A_bottom_left(index_i, index_j) != 0 || A_top_right(index_i, index_j) != 0)
208 : return false;
209 :
210 : return true;
211 : }
212 :
213 : template class GeneralizedRadialReturnStressUpdateTempl<false>;
214 : template class GeneralizedRadialReturnStressUpdateTempl<true>;
|
