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 "HillCreepStressUpdate.h"
11 : #include "ElasticityTensorTools.h"
12 : #include "Function.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", ADHillCreepStressUpdate);
15 : registerMooseObject("SolidMechanicsApp", HillCreepStressUpdate);
16 :
17 : template <bool is_ad>
18 : InputParameters
19 224 : HillCreepStressUpdateTempl<is_ad>::validParams()
20 : {
21 224 : InputParameters params = AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::validParams();
22 224 : params.addClassDescription(
23 : "This class uses the stress update material in a generalized radial return anisotropic power "
24 : "law creep "
25 : "model. This class can be used in conjunction with other creep and plasticity materials for "
26 : "more complex simulations.");
27 :
28 : // Linear strain hardening parameters
29 448 : params.addCoupledVar("temperature", "Coupled temperature");
30 448 : params.addRequiredParam<Real>("coefficient", "Leading coefficient in power-law equation");
31 448 : params.addRequiredParam<Real>("n_exponent", "Exponent on effective stress in power-law equation");
32 448 : params.addParam<Real>("m_exponent", 0.0, "Exponent on time in power-law equation");
33 448 : params.addRequiredParam<Real>("activation_energy", "Activation energy");
34 448 : params.addParam<Real>("gas_constant", 8.3143, "Universal gas constant");
35 448 : params.addParam<Real>("start_time", 0.0, "Start time (if not zero)");
36 448 : params.addParam<bool>("anisotropic_elasticity", false, "Enable using anisotropic elasticity");
37 448 : params.addParam<FunctionName>(
38 : "creep_prefactor", "Optional function to use as a scalar prefactor on the creep strain.");
39 224 : return params;
40 0 : }
41 :
42 : template <bool is_ad>
43 168 : HillCreepStressUpdateTempl<is_ad>::HillCreepStressUpdateTempl(const InputParameters & parameters)
44 : : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>(parameters),
45 168 : _has_temp(this->isParamValid("temperature")),
46 168 : _temperature(this->_has_temp ? this->coupledValue("temperature") : this->_zero),
47 336 : _coefficient(this->template getParam<Real>("coefficient")),
48 336 : _n_exponent(this->template getParam<Real>("n_exponent")),
49 336 : _m_exponent(this->template getParam<Real>("m_exponent")),
50 336 : _activation_energy(this->template getParam<Real>("activation_energy")),
51 336 : _gas_constant(this->template getParam<Real>("gas_constant")),
52 336 : _start_time(this->template getParam<Real>("start_time")),
53 168 : _exponential(1.0),
54 168 : _exp_time(1.0),
55 336 : _hill_constants(this->template getMaterialPropertyByName<std::vector<Real>>(this->_base_name +
56 : "hill_constants")),
57 336 : _hill_tensor(this->_use_transformation
58 450 : ? &this->template getMaterialPropertyByName<DenseMatrix<Real>>(
59 : this->_base_name + "hill_tensor")
60 : : nullptr),
61 168 : _C(6, 6),
62 168 : _elasticity_tensor_name(this->_base_name + "elasticity_tensor"),
63 168 : _elasticity_tensor(
64 168 : this->template getGenericMaterialProperty<RankFourTensor, is_ad>(_elasticity_tensor_name)),
65 336 : _anisotropic_elasticity(this->template getParam<bool>("anisotropic_elasticity")),
66 168 : _prefactor_function(
67 504 : this->isParamValid("creep_prefactor") ? &this->getFunction("creep_prefactor") : nullptr)
68 : {
69 168 : if (_start_time < this->_app.getStartTime() && (std::trunc(_m_exponent) != _m_exponent))
70 0 : this->paramError("start_time",
71 : "Start time must be equal to or greater than the Executioner start_time if a "
72 : "non-integer m_exponent is used");
73 168 : }
74 :
75 : template <bool is_ad>
76 : void
77 1368664 : HillCreepStressUpdateTempl<is_ad>::computeStressInitialize(
78 : const GenericDenseVector<is_ad> & /*stress_dev*/,
79 : const GenericDenseVector<is_ad> & /*stress*/,
80 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
81 : {
82 1368664 : if (_has_temp)
83 0 : _exponential = std::exp(-_activation_energy / (_gas_constant * _temperature[_qp]));
84 :
85 1942448 : _two_shear_modulus = 2.0 * ElasticityTensorTools::getIsotropicShearModulus(elasticity_tensor);
86 :
87 1368664 : _exp_time = std::pow(_t - _start_time, _m_exponent);
88 1368664 : }
89 :
90 : template <bool is_ad>
91 : GenericReal<is_ad>
92 1344568 : HillCreepStressUpdateTempl<is_ad>::initialGuess(const GenericDenseVector<is_ad> & /*stress_dev*/)
93 : {
94 1344568 : return 0.0;
95 : }
96 :
97 : template <bool is_ad>
98 : GenericReal<is_ad>
99 6902820 : HillCreepStressUpdateTempl<is_ad>::computeResidual(
100 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
101 : const GenericDenseVector<is_ad> & stress_new,
102 : const GenericReal<is_ad> & delta_gamma)
103 : {
104 : using std::pow;
105 :
106 : GenericReal<is_ad> qsigma_changed;
107 6902820 : GenericRankTwoTensor<is_ad> stress_changed;
108 :
109 : // Get the qsigma_changed and stress_changed (not used)
110 : // delta_gamma (scalar plastic strain) will change the stress (and qsigma as well)
111 6902820 : computeQsigmaChanged(qsigma_changed, stress_new, delta_gamma, stress_changed);
112 :
113 5995080 : GenericReal<is_ad> creep_rate =
114 6902820 : _coefficient * pow(qsigma_changed, _n_exponent) * _exponential * _exp_time;
115 :
116 6902820 : if (_prefactor_function)
117 0 : creep_rate *= _prefactor_function->value(_t, _q_point[_qp]);
118 :
119 6902820 : this->_inelastic_strain_rate[_qp] = MetaPhysicL::raw_value(creep_rate);
120 : // Return iteration difference between creep strain and inelastic strain multiplier
121 9900360 : return creep_rate * _dt - delta_gamma;
122 : }
123 :
124 : template <bool is_ad>
125 : Real
126 6902820 : HillCreepStressUpdateTempl<is_ad>::computeReferenceResidual(
127 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
128 : const GenericDenseVector<is_ad> & /*stress_new*/,
129 : const GenericReal<is_ad> & /*residual*/,
130 : const GenericReal<is_ad> & /*scalar_effective_inelastic_strain*/)
131 : {
132 6902820 : return 1.0;
133 : }
134 :
135 : template <bool is_ad>
136 : GenericReal<is_ad>
137 5558252 : HillCreepStressUpdateTempl<is_ad>::computeDerivative(
138 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
139 : const GenericDenseVector<is_ad> & stress_new,
140 : const GenericReal<is_ad> & delta_gamma)
141 : {
142 : using std::pow;
143 :
144 : GenericReal<is_ad> qsigma_changed;
145 5558252 : GenericRankTwoTensor<is_ad> stress_changed;
146 :
147 : // Get the qsigma_changed and stress_changed
148 : // delta_gamma (scalar plastic strain) will change the stress (and qsigma as well)
149 5558252 : computeQsigmaChanged(qsigma_changed, stress_new, delta_gamma, stress_changed);
150 :
151 5558252 : ElasticityTensorTools::toMooseVoigtNotation<is_ad>(_C, _elasticity_tensor[_qp]);
152 : const unsigned int dimension = _C.n();
153 5558252 : GenericDenseMatrix<is_ad> d_stress_d_inelasticStrainIncrement(dimension, dimension);
154 :
155 38907764 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
156 233446584 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
157 : {
158 200097072 : if (index_j < 3)
159 100048536 : d_stress_d_inelasticStrainIncrement(index_i, index_j) =
160 100048536 : -1.0 * MetaPhysicL::raw_value(_C(index_i, index_j));
161 : else
162 100048536 : d_stress_d_inelasticStrainIncrement(index_i, index_j) =
163 100048536 : -2.0 * MetaPhysicL::raw_value(_C(index_i, index_j));
164 : }
165 :
166 : // Hill constants, some constraints apply
167 5558252 : const Real & F = _hill_constants[_qp][0];
168 : const Real & G = _hill_constants[_qp][1];
169 : const Real & H = _hill_constants[_qp][2];
170 : const Real & L = _hill_constants[_qp][3];
171 : const Real & M = _hill_constants[_qp][4];
172 : const Real & N = _hill_constants[_qp][5];
173 :
174 5558252 : GenericDenseVector<is_ad> d_qsigma_d_inelasticStrainIncrement(6);
175 38907764 : for (unsigned int index_k = 0; index_k < 6; index_k++)
176 : {
177 14641032 : d_qsigma_d_inelasticStrainIncrement(index_k) =
178 18708480 : F * (stress_changed(1, 1) - stress_changed(2, 2)) *
179 18708480 : (d_stress_d_inelasticStrainIncrement(1, index_k) -
180 18708480 : d_stress_d_inelasticStrainIncrement(2, index_k)) +
181 18708480 : G * (stress_changed(2, 2) - stress_changed(0, 0)) *
182 18708480 : (d_stress_d_inelasticStrainIncrement(2, index_k) -
183 18708480 : d_stress_d_inelasticStrainIncrement(0, index_k)) +
184 18708480 : H * (stress_changed(0, 0) - stress_changed(1, 1)) *
185 18708480 : (d_stress_d_inelasticStrainIncrement(0, index_k) -
186 18708480 : d_stress_d_inelasticStrainIncrement(1, index_k)) +
187 33349512 : 2.0 * L * stress_changed(1, 2) * d_stress_d_inelasticStrainIncrement(4, index_k) +
188 33349512 : 2.0 * M * stress_changed(2, 0) * d_stress_d_inelasticStrainIncrement(5, index_k) +
189 33349512 : 2.0 * N * stress_changed(0, 1) * d_stress_d_inelasticStrainIncrement(3, index_k);
190 33349512 : d_qsigma_d_inelasticStrainIncrement(index_k) /= qsigma_changed;
191 : }
192 :
193 5558252 : GenericDenseVector<is_ad> d_qsigma_d_sigma(6);
194 :
195 5558252 : d_qsigma_d_sigma(0) = (H * (stress_changed(0, 0) - stress_changed(1, 1)) -
196 3118080 : G * (stress_changed(2, 2) - stress_changed(0, 0))) /
197 : qsigma_changed;
198 5558252 : d_qsigma_d_sigma(1) = (F * (stress_changed(1, 1) - stress_changed(2, 2)) -
199 3118080 : H * (stress_changed(0, 0) - stress_changed(1, 1))) /
200 : qsigma_changed;
201 5558252 : d_qsigma_d_sigma(2) = (G * (stress_changed(2, 2) - stress_changed(0, 0)) -
202 3118080 : F * (stress_changed(1, 1) - stress_changed(2, 2))) /
203 : qsigma_changed;
204 7998424 : d_qsigma_d_sigma(3) = 2.0 * N * stress_changed(0, 1) / qsigma_changed;
205 7998424 : d_qsigma_d_sigma(4) = 2.0 * L * stress_changed(1, 2) / qsigma_changed;
206 7998424 : d_qsigma_d_sigma(5) = 2.0 * M * stress_changed(0, 2) / qsigma_changed;
207 :
208 2440172 : GenericReal<is_ad> d_qsigma_d_deltaGamma =
209 3118080 : d_qsigma_d_inelasticStrainIncrement.dot(d_qsigma_d_sigma);
210 :
211 5558252 : GenericReal<is_ad> creep_rate_derivative = _coefficient * d_qsigma_d_deltaGamma * _n_exponent *
212 7998424 : pow(qsigma_changed, _n_exponent - 1.0) * _exponential *
213 5558252 : _exp_time;
214 :
215 5558252 : if (_prefactor_function)
216 0 : creep_rate_derivative *= _prefactor_function->value(_t, _q_point[_qp]);
217 :
218 13556676 : return (creep_rate_derivative * _dt - 1.0);
219 5558252 : }
220 :
221 : template <bool is_ad>
222 : void
223 1344568 : HillCreepStressUpdateTempl<is_ad>::computeStrainFinalize(
224 : GenericRankTwoTensor<is_ad> & inelasticStrainIncrement,
225 : const GenericRankTwoTensor<is_ad> & stress,
226 : const GenericDenseVector<is_ad> & stress_dev,
227 : const GenericReal<is_ad> & delta_gamma)
228 : {
229 : using std::sqrt;
230 :
231 : GenericReal<is_ad> qsigma_square;
232 1344568 : if (!this->_use_transformation)
233 : {
234 : // Hill constants, some constraints apply
235 393600 : const Real & F = _hill_constants[_qp][0];
236 : const Real & G = _hill_constants[_qp][1];
237 : const Real & H = _hill_constants[_qp][2];
238 : const Real & L = _hill_constants[_qp][3];
239 : const Real & M = _hill_constants[_qp][4];
240 : const Real & N = _hill_constants[_qp][5];
241 :
242 : // Equivalent deviatoric stress function.
243 393600 : qsigma_square = F * (stress(1, 1) - stress(2, 2)) * (stress(1, 1) - stress(2, 2));
244 393600 : qsigma_square += G * (stress(2, 2) - stress(0, 0)) * (stress(2, 2) - stress(0, 0));
245 393600 : qsigma_square += H * (stress(0, 0) - stress(1, 1)) * (stress(0, 0) - stress(1, 1));
246 393600 : qsigma_square += 2 * L * stress(1, 2) * stress(1, 2);
247 393600 : qsigma_square += 2 * M * stress(0, 2) * stress(0, 2);
248 393600 : qsigma_square += 2 * N * stress(0, 1) * stress(0, 1);
249 : }
250 : else
251 : {
252 950968 : GenericDenseVector<is_ad> Ms(6);
253 950968 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_dev);
254 950968 : qsigma_square = Ms.dot(stress_dev);
255 : }
256 :
257 1344568 : if (qsigma_square == 0)
258 : {
259 0 : inelasticStrainIncrement.zero();
260 :
261 0 : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
262 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
263 :
264 0 : this->_effective_inelastic_strain[_qp] = this->_effective_inelastic_strain_old[_qp];
265 :
266 0 : return;
267 : }
268 :
269 : // Use Hill-type flow rule to compute the time step inelastic increment.
270 1344568 : GenericReal<is_ad> prefactor = delta_gamma / sqrt(qsigma_square);
271 :
272 1344568 : if (!this->_use_transformation)
273 : {
274 : // Hill constants, some constraints apply
275 393600 : const Real & F = _hill_constants[_qp][0];
276 : const Real & G = _hill_constants[_qp][1];
277 : const Real & H = _hill_constants[_qp][2];
278 : const Real & L = _hill_constants[_qp][3];
279 : const Real & M = _hill_constants[_qp][4];
280 : const Real & N = _hill_constants[_qp][5];
281 :
282 393600 : inelasticStrainIncrement(0, 0) =
283 393600 : prefactor * (H * (stress(0, 0) - stress(1, 1)) - G * (stress(2, 2) - stress(0, 0)));
284 393600 : inelasticStrainIncrement(1, 1) =
285 393600 : prefactor * (F * (stress(1, 1) - stress(2, 2)) - H * (stress(0, 0) - stress(1, 1)));
286 393600 : inelasticStrainIncrement(2, 2) =
287 393600 : prefactor * (G * (stress(2, 2) - stress(0, 0)) - F * (stress(1, 1) - stress(2, 2)));
288 :
289 393600 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
290 393600 : prefactor * 2.0 * N * stress(0, 1);
291 393600 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
292 393600 : prefactor * 2.0 * M * stress(0, 2);
293 393600 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
294 393600 : prefactor * 2.0 * L * stress(1, 2);
295 : }
296 : else
297 : {
298 950968 : GenericDenseVector<is_ad> inelastic_strain_increment(6);
299 950968 : GenericDenseVector<is_ad> Ms(6);
300 950968 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_dev);
301 :
302 6656776 : for (unsigned int i = 0; i < 6; i++)
303 5705808 : inelastic_strain_increment(i) = Ms(i) * prefactor;
304 :
305 950968 : inelasticStrainIncrement(0, 0) = inelastic_strain_increment(0);
306 950968 : inelasticStrainIncrement(1, 1) = inelastic_strain_increment(1);
307 950968 : inelasticStrainIncrement(2, 2) = inelastic_strain_increment(2);
308 :
309 950968 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) = inelastic_strain_increment(3);
310 950968 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) = inelastic_strain_increment(4);
311 950968 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) = inelastic_strain_increment(5);
312 : }
313 :
314 1344568 : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
315 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
316 :
317 1901936 : this->_effective_inelastic_strain[_qp] = this->_effective_inelastic_strain_old[_qp] + delta_gamma;
318 : }
319 :
320 : template <bool is_ad>
321 : void
322 1368664 : HillCreepStressUpdateTempl<is_ad>::computeStressFinalize(
323 : const GenericRankTwoTensor<is_ad> & creepStrainIncrement,
324 : const GenericReal<is_ad> & /*delta_gamma*/,
325 : GenericRankTwoTensor<is_ad> & stress_new,
326 : const GenericDenseVector<is_ad> & /*stress_dev*/,
327 : const GenericRankTwoTensor<is_ad> & stress_old,
328 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
329 : {
330 : using std::abs;
331 :
332 : // Class only valid for isotropic elasticity (for now)
333 1368664 : stress_new -= elasticity_tensor * creepStrainIncrement;
334 :
335 : // Compute the maximum time step allowed due to creep strain numerical integration error
336 1942448 : Real stress_dif = MetaPhysicL::raw_value(stress_new - stress_old).L2norm();
337 :
338 : // Get a representative value of the elasticity tensor
339 1368664 : Real elasticity_value =
340 : 1.0 / 3.0 *
341 1368664 : MetaPhysicL::raw_value((elasticity_tensor(0, 0, 0, 0) + elasticity_tensor(1, 1, 1, 1) +
342 : elasticity_tensor(2, 2, 2, 2)));
343 :
344 1368664 : if (abs(stress_dif) > TOLERANCE * TOLERANCE)
345 1179172 : this->_max_integration_error_time_step =
346 1179172 : _dt / (stress_dif / elasticity_value / this->_max_integration_error);
347 : else
348 189492 : this->_max_integration_error_time_step = std::numeric_limits<Real>::max();
349 1368664 : }
350 :
351 : template <bool is_ad>
352 : void
353 12461072 : HillCreepStressUpdateTempl<is_ad>::computeQsigmaChanged(
354 : GenericReal<is_ad> & qsigma_changed,
355 : const GenericDenseVector<is_ad> & stress_new,
356 : const GenericReal<is_ad> & delta_gamma,
357 : GenericRankTwoTensor<is_ad> & stress_changed)
358 : {
359 : using std::sqrt;
360 :
361 : GenericReal<is_ad> qsigma_square;
362 :
363 : // Hill constants, some constraints apply
364 12461072 : const Real & F = _hill_constants[_qp][0];
365 : const Real & G = _hill_constants[_qp][1];
366 : const Real & H = _hill_constants[_qp][2];
367 : const Real & L = _hill_constants[_qp][3];
368 : const Real & M = _hill_constants[_qp][4];
369 : const Real & N = _hill_constants[_qp][5];
370 :
371 12461072 : if (!this->_use_transformation)
372 : {
373 3511680 : qsigma_square = F * (stress_new(1) - stress_new(2)) * (stress_new(1) - stress_new(2));
374 3511680 : qsigma_square += G * (stress_new(2) - stress_new(0)) * (stress_new(2) - stress_new(0));
375 3511680 : qsigma_square += H * (stress_new(0) - stress_new(1)) * (stress_new(0) - stress_new(1));
376 3511680 : qsigma_square += 2 * L * stress_new(4) * stress_new(4);
377 3511680 : qsigma_square += 2 * M * stress_new(5) * stress_new(5);
378 3511680 : qsigma_square += 2 * N * stress_new(3) * stress_new(3);
379 : }
380 : else
381 : {
382 8949392 : GenericDenseVector<is_ad> Ms(6);
383 8949392 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_new);
384 8949392 : qsigma_square = Ms.dot(stress_new);
385 : }
386 :
387 12461072 : GenericReal<is_ad> qsigma = sqrt(qsigma_square);
388 :
389 12461072 : if (!_anisotropic_elasticity)
390 17856240 : qsigma_changed = qsigma - 1.5 * _two_shear_modulus * delta_gamma;
391 : else
392 : {
393 21272 : GenericRankTwoTensor<is_ad> stress;
394 21272 : stress(0, 0) = stress_new(0);
395 21272 : stress(1, 1) = stress_new(1);
396 21272 : stress(2, 2) = stress_new(2);
397 21272 : stress(0, 1) = stress(1, 0) = stress_new(3);
398 21272 : stress(1, 2) = stress(2, 1) = stress_new(4);
399 21272 : stress(0, 2) = stress(2, 0) = stress_new(5);
400 :
401 21272 : GenericDenseVector<is_ad> b(6);
402 21272 : b(0) = H * (stress(0, 0) - stress(1, 1)) - G * (stress(2, 2) - stress(0, 0));
403 21272 : b(1) = F * (stress(1, 1) - stress(2, 2)) - H * (stress(0, 0) - stress(1, 1));
404 21272 : b(2) = G * (stress(2, 2) - stress(0, 0)) - F * (stress(1, 1) - stress(2, 2));
405 42544 : b(3) = 2.0 * N * stress(0, 1);
406 42544 : b(4) = 2.0 * L * stress(1, 2);
407 42544 : b(5) = 2.0 * M * stress(0, 2);
408 :
409 21272 : GenericRankTwoTensor<is_ad> inelasticStrainIncrement;
410 21272 : inelasticStrainIncrement(0, 0) = delta_gamma * b(0) / qsigma;
411 21272 : inelasticStrainIncrement(1, 1) = delta_gamma * b(1) / qsigma;
412 21272 : inelasticStrainIncrement(2, 2) = delta_gamma * b(2) / qsigma;
413 21272 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) = delta_gamma * b(3) / qsigma;
414 21272 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) = delta_gamma * b(4) / qsigma;
415 21272 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) = delta_gamma * b(5) / qsigma;
416 :
417 42544 : stress_changed = stress - _elasticity_tensor[_qp] * inelasticStrainIncrement;
418 :
419 : GenericReal<is_ad> qsigma_square_changed;
420 21272 : qsigma_square_changed = F * (stress_changed(1, 1) - stress_changed(2, 2)) *
421 : (stress_changed(1, 1) - stress_changed(2, 2));
422 21272 : qsigma_square_changed += G * (stress_changed(2, 2) - stress_changed(0, 0)) *
423 : (stress_changed(2, 2) - stress_changed(0, 0));
424 21272 : qsigma_square_changed += H * (stress_changed(0, 0) - stress_changed(1, 1)) *
425 : (stress_changed(0, 0) - stress_changed(1, 1));
426 42544 : qsigma_square_changed += 2 * L * stress_changed(1, 2) * stress_changed(1, 2);
427 42544 : qsigma_square_changed += 2 * M * stress_changed(0, 2) * stress_changed(0, 2);
428 42544 : qsigma_square_changed += 2 * N * stress_changed(0, 1) * stress_changed(0, 1);
429 21272 : qsigma_changed = sqrt(qsigma_square_changed);
430 : }
431 12461072 : }
432 :
433 : template class HillCreepStressUpdateTempl<false>;
434 : template class HillCreepStressUpdateTempl<true>;
|
