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 376 : HillCreepStressUpdateTempl<is_ad>::validParams()
20 : {
21 376 : InputParameters params = AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::validParams();
22 376 : 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 752 : params.addCoupledVar("temperature", "Coupled temperature");
30 752 : params.addRequiredParam<Real>("coefficient", "Leading coefficient in power-law equation");
31 752 : params.addRequiredParam<Real>("n_exponent", "Exponent on effective stress in power-law equation");
32 752 : params.addParam<Real>("m_exponent", 0.0, "Exponent on time in power-law equation");
33 752 : params.addRequiredParam<Real>("activation_energy", "Activation energy");
34 752 : params.addParam<Real>("gas_constant", 8.3143, "Universal gas constant");
35 752 : params.addParam<Real>("start_time", 0.0, "Start time (if not zero)");
36 752 : params.addParam<bool>("anisotropic_elasticity", false, "Enable using anisotropic elasticity");
37 752 : params.addParam<FunctionName>(
38 : "creep_prefactor", "Optional function to use as a scalar prefactor on the creep strain.");
39 376 : return params;
40 0 : }
41 :
42 : template <bool is_ad>
43 282 : HillCreepStressUpdateTempl<is_ad>::HillCreepStressUpdateTempl(const InputParameters & parameters)
44 : : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>(parameters),
45 282 : _has_temp(this->isParamValid("temperature")),
46 282 : _temperature(this->_has_temp ? this->coupledValue("temperature") : this->_zero),
47 564 : _coefficient(this->template getParam<Real>("coefficient")),
48 564 : _n_exponent(this->template getParam<Real>("n_exponent")),
49 564 : _m_exponent(this->template getParam<Real>("m_exponent")),
50 564 : _activation_energy(this->template getParam<Real>("activation_energy")),
51 564 : _gas_constant(this->template getParam<Real>("gas_constant")),
52 564 : _start_time(this->template getParam<Real>("start_time")),
53 282 : _exponential(1.0),
54 282 : _exp_time(1.0),
55 564 : _hill_constants(this->template getMaterialPropertyByName<std::vector<Real>>(this->_base_name +
56 : "hill_constants")),
57 564 : _hill_tensor(this->_use_transformation
58 756 : ? &this->template getMaterialPropertyByName<DenseMatrix<Real>>(
59 : this->_base_name + "hill_tensor")
60 : : nullptr),
61 282 : _C(6, 6),
62 282 : _elasticity_tensor_name(this->_base_name + "elasticity_tensor"),
63 282 : _elasticity_tensor(
64 282 : this->template getGenericMaterialProperty<RankFourTensor, is_ad>(_elasticity_tensor_name)),
65 564 : _anisotropic_elasticity(this->template getParam<bool>("anisotropic_elasticity")),
66 282 : _prefactor_function(
67 846 : this->isParamValid("creep_prefactor") ? &this->getFunction("creep_prefactor") : nullptr)
68 : {
69 282 : 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 282 : }
74 :
75 : template <bool is_ad>
76 : void
77 2247224 : 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 2247224 : if (_has_temp)
83 0 : _exponential = std::exp(-_activation_energy / (_gas_constant * _temperature[_qp]));
84 :
85 3185008 : _two_shear_modulus = 2.0 * ElasticityTensorTools::getIsotropicShearModulus(elasticity_tensor);
86 :
87 2247224 : _exp_time = std::pow(_t - _start_time, _m_exponent);
88 2247224 : }
89 :
90 : template <bool is_ad>
91 : GenericReal<is_ad>
92 2207064 : HillCreepStressUpdateTempl<is_ad>::initialGuess(const GenericDenseVector<is_ad> & /*stress_dev*/)
93 : {
94 2207064 : return 0.0;
95 : }
96 :
97 : template <bool is_ad>
98 : GenericReal<is_ad>
99 11326874 : 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 11326874 : 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 11326874 : computeQsigmaChanged(qsigma_changed, stress_new, delta_gamma, stress_changed);
112 :
113 9789748 : GenericReal<is_ad> creep_rate =
114 11326874 : _coefficient * pow(qsigma_changed, _n_exponent) * _exponential * _exp_time;
115 :
116 11326874 : if (_prefactor_function)
117 0 : creep_rate *= _prefactor_function->value(_t, _q_point[_qp]);
118 :
119 11326874 : this->_inelastic_strain_rate[_qp] = MetaPhysicL::raw_value(creep_rate);
120 : // Return iteration difference between creep strain and inelastic strain multiplier
121 16221748 : return creep_rate * _dt - delta_gamma;
122 : }
123 :
124 : template <bool is_ad>
125 : Real
126 11326874 : 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 11326874 : return 1.0;
133 : }
134 :
135 : template <bool is_ad>
136 : GenericReal<is_ad>
137 9119810 : 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 9119810 : 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 9119810 : computeQsigmaChanged(qsigma_changed, stress_new, delta_gamma, stress_changed);
150 :
151 9119810 : ElasticityTensorTools::toMooseVoigtNotation<is_ad>(_C, _elasticity_tensor[_qp]);
152 : const unsigned int dimension = _C.n();
153 9119810 : GenericDenseMatrix<is_ad> d_stress_d_inelasticStrainIncrement(dimension, dimension);
154 :
155 63838670 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
156 383032020 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
157 : {
158 328313160 : if (index_j < 3)
159 164156580 : d_stress_d_inelasticStrainIncrement(index_i, index_j) =
160 164156580 : -1.0 * MetaPhysicL::raw_value(_C(index_i, index_j));
161 : else
162 164156580 : d_stress_d_inelasticStrainIncrement(index_i, index_j) =
163 164156580 : -2.0 * MetaPhysicL::raw_value(_C(index_i, index_j));
164 : }
165 :
166 : // Hill constants, some constraints apply
167 9119810 : 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 9119810 : GenericDenseVector<is_ad> d_qsigma_d_inelasticStrainIncrement(6);
175 63838670 : for (unsigned int index_k = 0; index_k < 6; index_k++)
176 : {
177 23906700 : d_qsigma_d_inelasticStrainIncrement(index_k) =
178 30812160 : F * (stress_changed(1, 1) - stress_changed(2, 2)) *
179 30812160 : (d_stress_d_inelasticStrainIncrement(1, index_k) -
180 30812160 : d_stress_d_inelasticStrainIncrement(2, index_k)) +
181 30812160 : G * (stress_changed(2, 2) - stress_changed(0, 0)) *
182 30812160 : (d_stress_d_inelasticStrainIncrement(2, index_k) -
183 30812160 : d_stress_d_inelasticStrainIncrement(0, index_k)) +
184 30812160 : H * (stress_changed(0, 0) - stress_changed(1, 1)) *
185 30812160 : (d_stress_d_inelasticStrainIncrement(0, index_k) -
186 30812160 : d_stress_d_inelasticStrainIncrement(1, index_k)) +
187 54718860 : 2.0 * L * stress_changed(1, 2) * d_stress_d_inelasticStrainIncrement(4, index_k) +
188 54718860 : 2.0 * M * stress_changed(2, 0) * d_stress_d_inelasticStrainIncrement(5, index_k) +
189 54718860 : 2.0 * N * stress_changed(0, 1) * d_stress_d_inelasticStrainIncrement(3, index_k);
190 54718860 : d_qsigma_d_inelasticStrainIncrement(index_k) /= qsigma_changed;
191 : }
192 :
193 9119810 : GenericDenseVector<is_ad> d_qsigma_d_sigma(6);
194 :
195 9119810 : d_qsigma_d_sigma(0) = (H * (stress_changed(0, 0) - stress_changed(1, 1)) -
196 5135360 : G * (stress_changed(2, 2) - stress_changed(0, 0))) /
197 : qsigma_changed;
198 9119810 : d_qsigma_d_sigma(1) = (F * (stress_changed(1, 1) - stress_changed(2, 2)) -
199 5135360 : H * (stress_changed(0, 0) - stress_changed(1, 1))) /
200 : qsigma_changed;
201 9119810 : d_qsigma_d_sigma(2) = (G * (stress_changed(2, 2) - stress_changed(0, 0)) -
202 5135360 : F * (stress_changed(1, 1) - stress_changed(2, 2))) /
203 : qsigma_changed;
204 13104260 : d_qsigma_d_sigma(3) = 2.0 * N * stress_changed(0, 1) / qsigma_changed;
205 13104260 : d_qsigma_d_sigma(4) = 2.0 * L * stress_changed(1, 2) / qsigma_changed;
206 13104260 : d_qsigma_d_sigma(5) = 2.0 * M * stress_changed(0, 2) / qsigma_changed;
207 :
208 3984450 : GenericReal<is_ad> d_qsigma_d_deltaGamma =
209 5135360 : d_qsigma_d_inelasticStrainIncrement.dot(d_qsigma_d_sigma);
210 :
211 9119810 : GenericReal<is_ad> creep_rate_derivative = _coefficient * d_qsigma_d_deltaGamma * _n_exponent *
212 13104260 : pow(qsigma_changed, _n_exponent - 1.0) * _exponential *
213 9119810 : _exp_time;
214 :
215 9119810 : if (_prefactor_function)
216 0 : creep_rate_derivative *= _prefactor_function->value(_t, _q_point[_qp]);
217 :
218 22224070 : return (creep_rate_derivative * _dt - 1.0);
219 9119810 : }
220 :
221 : template <bool is_ad>
222 : void
223 2207064 : 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 2207064 : if (!this->_use_transformation)
233 : {
234 : // Hill constants, some constraints apply
235 648320 : 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 648320 : qsigma_square = F * (stress(1, 1) - stress(2, 2)) * (stress(1, 1) - stress(2, 2));
244 648320 : qsigma_square += G * (stress(2, 2) - stress(0, 0)) * (stress(2, 2) - stress(0, 0));
245 648320 : qsigma_square += H * (stress(0, 0) - stress(1, 1)) * (stress(0, 0) - stress(1, 1));
246 648320 : qsigma_square += 2 * L * stress(1, 2) * stress(1, 2);
247 648320 : qsigma_square += 2 * M * stress(0, 2) * stress(0, 2);
248 648320 : qsigma_square += 2 * N * stress(0, 1) * stress(0, 1);
249 : }
250 : else
251 : {
252 1558744 : GenericDenseVector<is_ad> Ms(6);
253 1558744 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_dev);
254 1558744 : qsigma_square = Ms.dot(stress_dev);
255 : }
256 :
257 2207064 : 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 2207064 : GenericReal<is_ad> prefactor = delta_gamma / sqrt(qsigma_square);
271 :
272 2207064 : if (!this->_use_transformation)
273 : {
274 : // Hill constants, some constraints apply
275 648320 : 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 648320 : inelasticStrainIncrement(0, 0) =
283 648320 : prefactor * (H * (stress(0, 0) - stress(1, 1)) - G * (stress(2, 2) - stress(0, 0)));
284 648320 : inelasticStrainIncrement(1, 1) =
285 648320 : prefactor * (F * (stress(1, 1) - stress(2, 2)) - H * (stress(0, 0) - stress(1, 1)));
286 648320 : inelasticStrainIncrement(2, 2) =
287 648320 : prefactor * (G * (stress(2, 2) - stress(0, 0)) - F * (stress(1, 1) - stress(2, 2)));
288 :
289 648320 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
290 648320 : prefactor * 2.0 * N * stress(0, 1);
291 648320 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
292 648320 : prefactor * 2.0 * M * stress(0, 2);
293 648320 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
294 648320 : prefactor * 2.0 * L * stress(1, 2);
295 : }
296 : else
297 : {
298 1558744 : GenericDenseVector<is_ad> inelastic_strain_increment(6);
299 1558744 : GenericDenseVector<is_ad> Ms(6);
300 1558744 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_dev);
301 :
302 10911208 : for (unsigned int i = 0; i < 6; i++)
303 9352464 : inelastic_strain_increment(i) = Ms(i) * prefactor;
304 :
305 1558744 : inelasticStrainIncrement(0, 0) = inelastic_strain_increment(0);
306 1558744 : inelasticStrainIncrement(1, 1) = inelastic_strain_increment(1);
307 1558744 : inelasticStrainIncrement(2, 2) = inelastic_strain_increment(2);
308 :
309 1558744 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) = inelastic_strain_increment(3);
310 1558744 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) = inelastic_strain_increment(4);
311 1558744 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) = inelastic_strain_increment(5);
312 : }
313 :
314 2207064 : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
315 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
316 :
317 3117488 : this->_effective_inelastic_strain[_qp] = this->_effective_inelastic_strain_old[_qp] + delta_gamma;
318 : }
319 :
320 : template <bool is_ad>
321 : void
322 2247224 : 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 2247224 : stress_new -= elasticity_tensor * creepStrainIncrement;
334 :
335 : // Compute the maximum time step allowed due to creep strain numerical integration error
336 3185008 : Real stress_dif = MetaPhysicL::raw_value(stress_new - stress_old).L2norm();
337 :
338 : // Get a representative value of the elasticity tensor
339 2247224 : Real elasticity_value =
340 : 1.0 / 3.0 *
341 2247224 : MetaPhysicL::raw_value((elasticity_tensor(0, 0, 0, 0) + elasticity_tensor(1, 1, 1, 1) +
342 : elasticity_tensor(2, 2, 2, 2)));
343 :
344 2247224 : if (abs(stress_dif) > TOLERANCE * TOLERANCE)
345 1931664 : this->_max_integration_error_time_step =
346 1931664 : _dt / (stress_dif / elasticity_value / this->_max_integration_error);
347 : else
348 315560 : this->_max_integration_error_time_step = std::numeric_limits<Real>::max();
349 2247224 : }
350 :
351 : template <bool is_ad>
352 : void
353 20446684 : 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 20446684 : 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 20446684 : if (!this->_use_transformation)
372 : {
373 5783680 : qsigma_square = F * (stress_new(1) - stress_new(2)) * (stress_new(1) - stress_new(2));
374 5783680 : qsigma_square += G * (stress_new(2) - stress_new(0)) * (stress_new(2) - stress_new(0));
375 5783680 : qsigma_square += H * (stress_new(0) - stress_new(1)) * (stress_new(0) - stress_new(1));
376 5783680 : qsigma_square += 2 * L * stress_new(4) * stress_new(4);
377 5783680 : qsigma_square += 2 * M * stress_new(5) * stress_new(5);
378 5783680 : qsigma_square += 2 * N * stress_new(3) * stress_new(3);
379 : }
380 : else
381 : {
382 14663004 : GenericDenseVector<is_ad> Ms(6);
383 14663004 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_new);
384 14663004 : qsigma_square = Ms.dot(stress_new);
385 : }
386 :
387 20446684 : GenericReal<is_ad> qsigma = sqrt(qsigma_square);
388 :
389 20446684 : if (!_anisotropic_elasticity)
390 29256648 : qsigma_changed = qsigma - 1.5 * _two_shear_modulus * delta_gamma;
391 : else
392 : {
393 34680 : GenericRankTwoTensor<is_ad> stress;
394 34680 : stress(0, 0) = stress_new(0);
395 34680 : stress(1, 1) = stress_new(1);
396 34680 : stress(2, 2) = stress_new(2);
397 34680 : stress(0, 1) = stress(1, 0) = stress_new(3);
398 34680 : stress(1, 2) = stress(2, 1) = stress_new(4);
399 34680 : stress(0, 2) = stress(2, 0) = stress_new(5);
400 :
401 34680 : GenericDenseVector<is_ad> b(6);
402 34680 : b(0) = H * (stress(0, 0) - stress(1, 1)) - G * (stress(2, 2) - stress(0, 0));
403 34680 : b(1) = F * (stress(1, 1) - stress(2, 2)) - H * (stress(0, 0) - stress(1, 1));
404 34680 : b(2) = G * (stress(2, 2) - stress(0, 0)) - F * (stress(1, 1) - stress(2, 2));
405 69360 : b(3) = 2.0 * N * stress(0, 1);
406 69360 : b(4) = 2.0 * L * stress(1, 2);
407 69360 : b(5) = 2.0 * M * stress(0, 2);
408 :
409 34680 : GenericRankTwoTensor<is_ad> inelasticStrainIncrement;
410 34680 : inelasticStrainIncrement(0, 0) = delta_gamma * b(0) / qsigma;
411 34680 : inelasticStrainIncrement(1, 1) = delta_gamma * b(1) / qsigma;
412 34680 : inelasticStrainIncrement(2, 2) = delta_gamma * b(2) / qsigma;
413 34680 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) = delta_gamma * b(3) / qsigma;
414 34680 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) = delta_gamma * b(4) / qsigma;
415 34680 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) = delta_gamma * b(5) / qsigma;
416 :
417 69360 : stress_changed = stress - _elasticity_tensor[_qp] * inelasticStrainIncrement;
418 :
419 : GenericReal<is_ad> qsigma_square_changed;
420 34680 : qsigma_square_changed = F * (stress_changed(1, 1) - stress_changed(2, 2)) *
421 : (stress_changed(1, 1) - stress_changed(2, 2));
422 34680 : qsigma_square_changed += G * (stress_changed(2, 2) - stress_changed(0, 0)) *
423 : (stress_changed(2, 2) - stress_changed(0, 0));
424 34680 : qsigma_square_changed += H * (stress_changed(0, 0) - stress_changed(1, 1)) *
425 : (stress_changed(0, 0) - stress_changed(1, 1));
426 69360 : qsigma_square_changed += 2 * L * stress_changed(1, 2) * stress_changed(1, 2);
427 69360 : qsigma_square_changed += 2 * M * stress_changed(0, 2) * stress_changed(0, 2);
428 69360 : qsigma_square_changed += 2 * N * stress_changed(0, 1) * stress_changed(0, 1);
429 34680 : qsigma_changed = sqrt(qsigma_square_changed);
430 : }
431 20446684 : }
432 :
433 : template class HillCreepStressUpdateTempl<false>;
434 : template class HillCreepStressUpdateTempl<true>;
|
