Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
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 "HillCreepStressUpdate.h"
11 : #include "ElasticityTensorTools.h"
12 :
13 : registerMooseObject("TensorMechanicsApp", ADHillCreepStressUpdate);
14 : registerMooseObject("TensorMechanicsApp", HillCreepStressUpdate);
15 :
16 : template <bool is_ad>
17 : InputParameters
18 128 : HillCreepStressUpdateTempl<is_ad>::validParams()
19 : {
20 128 : InputParameters params = AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::validParams();
21 128 : params.addClassDescription(
22 : "This class uses the stress update material in a generalized radial return anisotropic power "
23 : "law creep "
24 : "model. This class can be used in conjunction with other creep and plasticity materials for "
25 : "more complex simulations.");
26 :
27 : // Linear strain hardening parameters
28 256 : params.addCoupledVar("temperature", "Coupled temperature");
29 256 : params.addRequiredParam<Real>("coefficient", "Leading coefficient in power-law equation");
30 256 : params.addRequiredParam<Real>("n_exponent", "Exponent on effective stress in power-law equation");
31 256 : params.addParam<Real>("m_exponent", 0.0, "Exponent on time in power-law equation");
32 256 : params.addRequiredParam<Real>("activation_energy", "Activation energy");
33 256 : params.addParam<Real>("gas_constant", 8.3143, "Universal gas constant");
34 256 : params.addParam<Real>("start_time", 0.0, "Start time (if not zero)");
35 :
36 128 : return params;
37 0 : }
38 :
39 : template <bool is_ad>
40 96 : HillCreepStressUpdateTempl<is_ad>::HillCreepStressUpdateTempl(const InputParameters & parameters)
41 : : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>(parameters),
42 96 : _has_temp(this->isParamValid("temperature")),
43 96 : _temperature(this->_has_temp ? this->coupledValue("temperature") : this->_zero),
44 192 : _coefficient(this->template getParam<Real>("coefficient")),
45 192 : _n_exponent(this->template getParam<Real>("n_exponent")),
46 192 : _m_exponent(this->template getParam<Real>("m_exponent")),
47 192 : _activation_energy(this->template getParam<Real>("activation_energy")),
48 192 : _gas_constant(this->template getParam<Real>("gas_constant")),
49 192 : _start_time(this->template getParam<Real>("start_time")),
50 96 : _exponential(1.0),
51 96 : _exp_time(1.0),
52 192 : _hill_constants(this->template getMaterialPropertyByName<std::vector<Real>>(this->_base_name +
53 : "hill_constants")),
54 192 : _hill_tensor(this->_use_transformation
55 180 : ? &this->template getMaterialPropertyByName<DenseMatrix<Real>>(
56 : this->_base_name + "hill_tensor")
57 : : nullptr),
58 168 : _qsigma(0.0)
59 : {
60 96 : if (_start_time < this->_app.getStartTime() && (std::trunc(_m_exponent) != _m_exponent))
61 0 : this->template paramError(
62 : "start_time",
63 : "Start time must be equal to or greater than the Executioner start_time if a "
64 : "non-integer m_exponent is used");
65 96 : }
66 :
67 : template <bool is_ad>
68 : void
69 985520 : HillCreepStressUpdateTempl<is_ad>::computeStressInitialize(
70 : const GenericDenseVector<is_ad> & /*stress_dev*/,
71 : const GenericDenseVector<is_ad> & /*stress*/,
72 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
73 : {
74 985520 : if (_has_temp)
75 0 : _exponential = std::exp(-_activation_energy / (_gas_constant * _temperature[_qp]));
76 :
77 1417312 : _two_shear_modulus = 2.0 * ElasticityTensorTools::getIsotropicShearModulus(elasticity_tensor);
78 :
79 985520 : _exp_time = std::pow(_t - _start_time, _m_exponent);
80 985520 : }
81 :
82 : template <bool is_ad>
83 : GenericReal<is_ad>
84 960968 : HillCreepStressUpdateTempl<is_ad>::initialGuess(const GenericDenseVector<is_ad> & /*stress_dev*/)
85 : {
86 960968 : return 0.0;
87 : }
88 :
89 : template <bool is_ad>
90 : GenericReal<is_ad>
91 2853716 : HillCreepStressUpdateTempl<is_ad>::computeResidual(
92 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
93 : const GenericDenseVector<is_ad> & stress_new,
94 : const GenericReal<is_ad> & delta_gamma)
95 : {
96 : GenericReal<is_ad> qsigma_square;
97 2853716 : if (!this->_use_transformation)
98 : {
99 : // Hill constants, some constraints apply
100 813536 : const Real & F = _hill_constants[_qp][0];
101 : const Real & G = _hill_constants[_qp][1];
102 : const Real & H = _hill_constants[_qp][2];
103 : const Real & L = _hill_constants[_qp][3];
104 : const Real & M = _hill_constants[_qp][4];
105 : const Real & N = _hill_constants[_qp][5];
106 :
107 813536 : qsigma_square = F * (stress_new(1) - stress_new(2)) * (stress_new(1) - stress_new(2));
108 813536 : qsigma_square += G * (stress_new(2) - stress_new(0)) * (stress_new(2) - stress_new(0));
109 813536 : qsigma_square += H * (stress_new(0) - stress_new(1)) * (stress_new(0) - stress_new(1));
110 813536 : qsigma_square += 2 * L * stress_new(4) * stress_new(4);
111 813536 : qsigma_square += 2 * M * stress_new(5) * stress_new(5);
112 813536 : qsigma_square += 2 * N * stress_new(3) * stress_new(3);
113 : }
114 : else
115 : {
116 2040180 : GenericDenseVector<is_ad> Ms(6);
117 2040180 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_new);
118 2040180 : qsigma_square = Ms.dot(stress_new);
119 : }
120 :
121 2853716 : qsigma_square = std::sqrt(qsigma_square);
122 4080360 : qsigma_square -= 1.5 * _two_shear_modulus * delta_gamma;
123 :
124 1226644 : const GenericReal<is_ad> creep_rate =
125 2853716 : _coefficient * std::pow(qsigma_square, _n_exponent) * _exponential * _exp_time;
126 :
127 2853716 : this->_inelastic_strain_rate[_qp] = MetaPhysicL::raw_value(creep_rate);
128 : // Return iteration difference between creep strain and inelastic strain multiplier
129 4080360 : return creep_rate * _dt - delta_gamma;
130 : }
131 :
132 : template <bool is_ad>
133 : Real
134 2853716 : HillCreepStressUpdateTempl<is_ad>::computeReferenceResidual(
135 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
136 : const GenericDenseVector<is_ad> & /*stress_new*/,
137 : const GenericReal<is_ad> & /*residual*/,
138 : const GenericReal<is_ad> & /*scalar_effective_inelastic_strain*/)
139 : {
140 2853716 : return 1.0;
141 : }
142 :
143 : template <bool is_ad>
144 : GenericReal<is_ad>
145 1892748 : HillCreepStressUpdateTempl<is_ad>::computeDerivative(
146 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
147 : const GenericDenseVector<is_ad> & stress_new,
148 : const GenericReal<is_ad> & delta_gamma)
149 : {
150 : GenericReal<is_ad> qsigma_square;
151 1892748 : if (!this->_use_transformation)
152 : {
153 : // Hill constants, some constraints apply
154 540608 : const Real & F = _hill_constants[_qp][0];
155 : const Real & G = _hill_constants[_qp][1];
156 : const Real & H = _hill_constants[_qp][2];
157 : const Real & L = _hill_constants[_qp][3];
158 : const Real & M = _hill_constants[_qp][4];
159 : const Real & N = _hill_constants[_qp][5];
160 :
161 : // Equivalent deviatoric stress function.
162 540608 : qsigma_square = F * (stress_new(1) - stress_new(2)) * (stress_new(1) - stress_new(2));
163 540608 : qsigma_square += G * (stress_new(2) - stress_new(0)) * (stress_new(2) - stress_new(0));
164 540608 : qsigma_square += H * (stress_new(0) - stress_new(1)) * (stress_new(0) - stress_new(1));
165 540608 : qsigma_square += 2 * L * stress_new(4) * stress_new(4);
166 540608 : qsigma_square += 2 * M * stress_new(5) * stress_new(5);
167 540608 : qsigma_square += 2 * N * stress_new(3) * stress_new(3);
168 : }
169 : else
170 : {
171 1352140 : GenericDenseVector<is_ad> Ms(6);
172 1352140 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_new);
173 1352140 : qsigma_square = Ms.dot(stress_new);
174 : }
175 :
176 1892748 : qsigma_square = std::sqrt(qsigma_square);
177 2704280 : qsigma_square -= 1.5 * _two_shear_modulus * delta_gamma;
178 :
179 1892748 : _qsigma = qsigma_square;
180 :
181 2704280 : const GenericReal<is_ad> creep_rate_derivative =
182 1892748 : -_coefficient * 1.5 * _two_shear_modulus * _n_exponent *
183 1892748 : std::pow(qsigma_square, _n_exponent - 1.0) * _exponential * _exp_time;
184 2704280 : return (creep_rate_derivative * _dt - 1.0);
185 : }
186 :
187 : template <bool is_ad>
188 : void
189 960968 : HillCreepStressUpdateTempl<is_ad>::computeStrainFinalize(
190 : GenericRankTwoTensor<is_ad> & inelasticStrainIncrement,
191 : const GenericRankTwoTensor<is_ad> & stress,
192 : const GenericDenseVector<is_ad> & stress_dev,
193 : const GenericReal<is_ad> & delta_gamma)
194 : {
195 : GenericReal<is_ad> qsigma_square;
196 960968 : if (!this->_use_transformation)
197 : {
198 : // Hill constants, some constraints apply
199 272928 : const Real & F = _hill_constants[_qp][0];
200 : const Real & G = _hill_constants[_qp][1];
201 : const Real & H = _hill_constants[_qp][2];
202 : const Real & L = _hill_constants[_qp][3];
203 : const Real & M = _hill_constants[_qp][4];
204 : const Real & N = _hill_constants[_qp][5];
205 :
206 : // Equivalent deviatoric stress function.
207 272928 : qsigma_square = F * (stress(1, 1) - stress(2, 2)) * (stress(1, 1) - stress(2, 2));
208 272928 : qsigma_square += G * (stress(2, 2) - stress(0, 0)) * (stress(2, 2) - stress(0, 0));
209 272928 : qsigma_square += H * (stress(0, 0) - stress(1, 1)) * (stress(0, 0) - stress(1, 1));
210 272928 : qsigma_square += 2 * L * stress(1, 2) * stress(1, 2);
211 272928 : qsigma_square += 2 * M * stress(0, 2) * stress(0, 2);
212 272928 : qsigma_square += 2 * N * stress(0, 1) * stress(0, 1);
213 : }
214 : else
215 : {
216 688040 : GenericDenseVector<is_ad> Ms(6);
217 688040 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_dev);
218 688040 : qsigma_square = Ms.dot(stress_dev);
219 : }
220 :
221 960968 : if (qsigma_square == 0)
222 : {
223 0 : inelasticStrainIncrement.zero();
224 :
225 0 : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
226 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
227 :
228 0 : this->_effective_inelastic_strain[_qp] = this->_effective_inelastic_strain_old[_qp];
229 :
230 0 : return;
231 : }
232 :
233 : // Use Hill-type flow rule to compute the time step inelastic increment.
234 960968 : GenericReal<is_ad> prefactor = delta_gamma / std::sqrt(qsigma_square);
235 :
236 960968 : if (!this->_use_transformation)
237 : {
238 : // Hill constants, some constraints apply
239 272928 : const Real & F = _hill_constants[_qp][0];
240 : const Real & G = _hill_constants[_qp][1];
241 : const Real & H = _hill_constants[_qp][2];
242 : const Real & L = _hill_constants[_qp][3];
243 : const Real & M = _hill_constants[_qp][4];
244 : const Real & N = _hill_constants[_qp][5];
245 :
246 272928 : inelasticStrainIncrement(0, 0) =
247 272928 : prefactor * (H * (stress(0, 0) - stress(1, 1)) - G * (stress(2, 2) - stress(0, 0)));
248 272928 : inelasticStrainIncrement(1, 1) =
249 272928 : prefactor * (F * (stress(1, 1) - stress(2, 2)) - H * (stress(0, 0) - stress(1, 1)));
250 272928 : inelasticStrainIncrement(2, 2) =
251 272928 : prefactor * (G * (stress(2, 2) - stress(0, 0)) - F * (stress(1, 1) - stress(2, 2)));
252 :
253 272928 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
254 272928 : prefactor * 2.0 * N * stress(0, 1);
255 272928 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
256 272928 : prefactor * 2.0 * M * stress(0, 2);
257 272928 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
258 272928 : prefactor * 2.0 * L * stress(1, 2);
259 : }
260 : else
261 : {
262 688040 : GenericDenseVector<is_ad> inelastic_strain_increment(6);
263 688040 : GenericDenseVector<is_ad> Ms(6);
264 688040 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_dev);
265 :
266 4816280 : for (unsigned int i = 0; i < 6; i++)
267 4128240 : inelastic_strain_increment(i) = Ms(i) * prefactor;
268 :
269 688040 : inelasticStrainIncrement(0, 0) = inelastic_strain_increment(0);
270 688040 : inelasticStrainIncrement(1, 1) = inelastic_strain_increment(1);
271 688040 : inelasticStrainIncrement(2, 2) = inelastic_strain_increment(2);
272 :
273 688040 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) = inelastic_strain_increment(3);
274 688040 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) = inelastic_strain_increment(4);
275 688040 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) = inelastic_strain_increment(5);
276 : }
277 :
278 960968 : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
279 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
280 :
281 1376080 : this->_effective_inelastic_strain[_qp] = this->_effective_inelastic_strain_old[_qp] + delta_gamma;
282 : }
283 :
284 : template <bool is_ad>
285 : void
286 985520 : HillCreepStressUpdateTempl<is_ad>::computeStressFinalize(
287 : const GenericRankTwoTensor<is_ad> & creepStrainIncrement,
288 : const GenericReal<is_ad> & /*delta_gamma*/,
289 : GenericRankTwoTensor<is_ad> & stress_new,
290 : const GenericDenseVector<is_ad> & /*stress_dev*/,
291 : const GenericRankTwoTensor<is_ad> & stress_old,
292 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
293 : {
294 : // Class only valid for isotropic elasticity (for now)
295 985520 : stress_new -= elasticity_tensor * creepStrainIncrement;
296 :
297 : // Compute the maximum time step allowed due to creep strain numerical integration error
298 1417312 : Real stress_dif = MetaPhysicL::raw_value(stress_new - stress_old).L2norm();
299 :
300 : // Get a representative value of the elasticity tensor
301 985520 : Real elasticity_value =
302 : 1.0 / 3.0 *
303 985520 : MetaPhysicL::raw_value((elasticity_tensor(0, 0, 0, 0) + elasticity_tensor(1, 1, 1, 1) +
304 : elasticity_tensor(2, 2, 2, 2)));
305 :
306 985520 : if (std::abs(stress_dif) > TOLERANCE * TOLERANCE)
307 853600 : this->_max_integration_error_time_step =
308 853600 : _dt / (stress_dif / elasticity_value / this->_max_integration_error);
309 : else
310 131920 : this->_max_integration_error_time_step = std::numeric_limits<Real>::max();
311 985520 : }
312 :
313 : template class HillCreepStressUpdateTempl<false>;
314 : template class HillCreepStressUpdateTempl<true>;
|
