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 2310908 : 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 2310908 : if (_has_temp)
83 0 : _exponential = std::exp(-_activation_energy / (_gas_constant * _temperature[_qp]));
84 :
85 3279480 : _two_shear_modulus = 2.0 * ElasticityTensorTools::getIsotropicShearModulus(elasticity_tensor);
86 :
87 2310908 : _exp_time = std::pow(_t - _start_time, _m_exponent);
88 2310908 : }
89 :
90 : template <bool is_ad>
91 : GenericReal<is_ad>
92 2269588 : HillCreepStressUpdateTempl<is_ad>::initialGuess(const GenericDenseVector<is_ad> & /*stress_dev*/)
93 : {
94 2269588 : return 0.0;
95 : }
96 :
97 : template <bool is_ad>
98 : GenericReal<is_ad>
99 11629065 : 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 : GenericReal<is_ad> qsigma_changed;
105 11629065 : GenericRankTwoTensor<is_ad> stress_changed;
106 :
107 : // Get the qsigma_changed and stress_changed (not used)
108 : // delta_gamma (scalar plastic strain) will change the stress (and qsigma as well)
109 11629065 : computeQsigmaChanged(qsigma_changed, stress_new, delta_gamma, stress_changed);
110 :
111 10070930 : GenericReal<is_ad> creep_rate =
112 11629065 : _coefficient * std::pow(qsigma_changed, _n_exponent) * _exponential * _exp_time;
113 :
114 11629065 : if (_prefactor_function)
115 0 : creep_rate *= _prefactor_function->value(_t, _q_point[_qp]);
116 :
117 11629065 : this->_inelastic_strain_rate[_qp] = MetaPhysicL::raw_value(creep_rate);
118 : // Return iteration difference between creep strain and inelastic strain multiplier
119 16664530 : return creep_rate * _dt - delta_gamma;
120 : }
121 :
122 : template <bool is_ad>
123 : Real
124 11629065 : HillCreepStressUpdateTempl<is_ad>::computeReferenceResidual(
125 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
126 : const GenericDenseVector<is_ad> & /*stress_new*/,
127 : const GenericReal<is_ad> & /*residual*/,
128 : const GenericReal<is_ad> & /*scalar_effective_inelastic_strain*/)
129 : {
130 11629065 : return 1.0;
131 : }
132 :
133 : template <bool is_ad>
134 : GenericReal<is_ad>
135 9359477 : HillCreepStressUpdateTempl<is_ad>::computeDerivative(
136 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
137 : const GenericDenseVector<is_ad> & stress_new,
138 : const GenericReal<is_ad> & delta_gamma)
139 : {
140 : GenericReal<is_ad> qsigma_changed;
141 9359477 : GenericRankTwoTensor<is_ad> stress_changed;
142 :
143 : // Get the qsigma_changed and stress_changed
144 : // delta_gamma (scalar plastic strain) will change the stress (and qsigma as well)
145 9359477 : computeQsigmaChanged(qsigma_changed, stress_new, delta_gamma, stress_changed);
146 :
147 9359477 : ElasticityTensorTools::toMooseVoigtNotation<is_ad>(_C, _elasticity_tensor[_qp]);
148 : const unsigned int dimension = _C.n();
149 9359477 : GenericDenseMatrix<is_ad> d_stress_d_inelasticStrainIncrement(dimension, dimension);
150 :
151 65516339 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
152 393098034 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
153 : {
154 336941172 : if (index_j < 3)
155 168470586 : d_stress_d_inelasticStrainIncrement(index_i, index_j) =
156 168470586 : -1.0 * MetaPhysicL::raw_value(_C(index_i, index_j));
157 : else
158 168470586 : d_stress_d_inelasticStrainIncrement(index_i, index_j) =
159 168470586 : -2.0 * MetaPhysicL::raw_value(_C(index_i, index_j));
160 : }
161 :
162 : // Hill constants, some constraints apply
163 9359477 : const Real & F = _hill_constants[_qp][0];
164 : const Real & G = _hill_constants[_qp][1];
165 : const Real & H = _hill_constants[_qp][2];
166 : const Real & L = _hill_constants[_qp][3];
167 : const Real & M = _hill_constants[_qp][4];
168 : const Real & N = _hill_constants[_qp][5];
169 :
170 9359477 : GenericDenseVector<is_ad> d_qsigma_d_inelasticStrainIncrement(6);
171 65516339 : for (unsigned int index_k = 0; index_k < 6; index_k++)
172 : {
173 24570558 : d_qsigma_d_inelasticStrainIncrement(index_k) =
174 31586304 : F * (stress_changed(1, 1) - stress_changed(2, 2)) *
175 31586304 : (d_stress_d_inelasticStrainIncrement(1, index_k) -
176 31586304 : d_stress_d_inelasticStrainIncrement(2, index_k)) +
177 31586304 : G * (stress_changed(2, 2) - stress_changed(0, 0)) *
178 31586304 : (d_stress_d_inelasticStrainIncrement(2, index_k) -
179 31586304 : d_stress_d_inelasticStrainIncrement(0, index_k)) +
180 31586304 : H * (stress_changed(0, 0) - stress_changed(1, 1)) *
181 31586304 : (d_stress_d_inelasticStrainIncrement(0, index_k) -
182 31586304 : d_stress_d_inelasticStrainIncrement(1, index_k)) +
183 56156862 : 2.0 * L * stress_changed(1, 2) * d_stress_d_inelasticStrainIncrement(4, index_k) +
184 56156862 : 2.0 * M * stress_changed(2, 0) * d_stress_d_inelasticStrainIncrement(5, index_k) +
185 56156862 : 2.0 * N * stress_changed(0, 1) * d_stress_d_inelasticStrainIncrement(3, index_k);
186 56156862 : d_qsigma_d_inelasticStrainIncrement(index_k) /= qsigma_changed;
187 : }
188 :
189 9359477 : GenericDenseVector<is_ad> d_qsigma_d_sigma(6);
190 :
191 9359477 : d_qsigma_d_sigma(0) = (H * (stress_changed(0, 0) - stress_changed(1, 1)) -
192 5264384 : G * (stress_changed(2, 2) - stress_changed(0, 0))) /
193 : qsigma_changed;
194 9359477 : d_qsigma_d_sigma(1) = (F * (stress_changed(1, 1) - stress_changed(2, 2)) -
195 5264384 : H * (stress_changed(0, 0) - stress_changed(1, 1))) /
196 : qsigma_changed;
197 9359477 : d_qsigma_d_sigma(2) = (G * (stress_changed(2, 2) - stress_changed(0, 0)) -
198 5264384 : F * (stress_changed(1, 1) - stress_changed(2, 2))) /
199 : qsigma_changed;
200 13454570 : d_qsigma_d_sigma(3) = 2.0 * N * stress_changed(0, 1) / qsigma_changed;
201 13454570 : d_qsigma_d_sigma(4) = 2.0 * L * stress_changed(1, 2) / qsigma_changed;
202 13454570 : d_qsigma_d_sigma(5) = 2.0 * M * stress_changed(0, 2) / qsigma_changed;
203 :
204 4095093 : GenericReal<is_ad> d_qsigma_d_deltaGamma =
205 5264384 : d_qsigma_d_inelasticStrainIncrement.dot(d_qsigma_d_sigma);
206 :
207 9359477 : GenericReal<is_ad> creep_rate_derivative = _coefficient * d_qsigma_d_deltaGamma * _n_exponent *
208 13454570 : std::pow(qsigma_changed, _n_exponent - 1.0) *
209 9359477 : _exponential * _exp_time;
210 :
211 9359477 : if (_prefactor_function)
212 0 : creep_rate_derivative *= _prefactor_function->value(_t, _q_point[_qp]);
213 :
214 22814047 : return (creep_rate_derivative * _dt - 1.0);
215 9359477 : }
216 :
217 : template <bool is_ad>
218 : void
219 2269588 : HillCreepStressUpdateTempl<is_ad>::computeStrainFinalize(
220 : GenericRankTwoTensor<is_ad> & inelasticStrainIncrement,
221 : const GenericRankTwoTensor<is_ad> & stress,
222 : const GenericDenseVector<is_ad> & stress_dev,
223 : const GenericReal<is_ad> & delta_gamma)
224 : {
225 : GenericReal<is_ad> qsigma_square;
226 2269588 : if (!this->_use_transformation)
227 : {
228 : // Hill constants, some constraints apply
229 664608 : const Real & F = _hill_constants[_qp][0];
230 : const Real & G = _hill_constants[_qp][1];
231 : const Real & H = _hill_constants[_qp][2];
232 : const Real & L = _hill_constants[_qp][3];
233 : const Real & M = _hill_constants[_qp][4];
234 : const Real & N = _hill_constants[_qp][5];
235 :
236 : // Equivalent deviatoric stress function.
237 664608 : qsigma_square = F * (stress(1, 1) - stress(2, 2)) * (stress(1, 1) - stress(2, 2));
238 664608 : qsigma_square += G * (stress(2, 2) - stress(0, 0)) * (stress(2, 2) - stress(0, 0));
239 664608 : qsigma_square += H * (stress(0, 0) - stress(1, 1)) * (stress(0, 0) - stress(1, 1));
240 664608 : qsigma_square += 2 * L * stress(1, 2) * stress(1, 2);
241 664608 : qsigma_square += 2 * M * stress(0, 2) * stress(0, 2);
242 664608 : qsigma_square += 2 * N * stress(0, 1) * stress(0, 1);
243 : }
244 : else
245 : {
246 1604980 : GenericDenseVector<is_ad> Ms(6);
247 1604980 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_dev);
248 1604980 : qsigma_square = Ms.dot(stress_dev);
249 : }
250 :
251 2269588 : if (qsigma_square == 0)
252 : {
253 0 : inelasticStrainIncrement.zero();
254 :
255 0 : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
256 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
257 :
258 0 : this->_effective_inelastic_strain[_qp] = this->_effective_inelastic_strain_old[_qp];
259 :
260 0 : return;
261 : }
262 :
263 : // Use Hill-type flow rule to compute the time step inelastic increment.
264 2269588 : GenericReal<is_ad> prefactor = delta_gamma / std::sqrt(qsigma_square);
265 :
266 2269588 : if (!this->_use_transformation)
267 : {
268 : // Hill constants, some constraints apply
269 664608 : const Real & F = _hill_constants[_qp][0];
270 : const Real & G = _hill_constants[_qp][1];
271 : const Real & H = _hill_constants[_qp][2];
272 : const Real & L = _hill_constants[_qp][3];
273 : const Real & M = _hill_constants[_qp][4];
274 : const Real & N = _hill_constants[_qp][5];
275 :
276 664608 : inelasticStrainIncrement(0, 0) =
277 664608 : prefactor * (H * (stress(0, 0) - stress(1, 1)) - G * (stress(2, 2) - stress(0, 0)));
278 664608 : inelasticStrainIncrement(1, 1) =
279 664608 : prefactor * (F * (stress(1, 1) - stress(2, 2)) - H * (stress(0, 0) - stress(1, 1)));
280 664608 : inelasticStrainIncrement(2, 2) =
281 664608 : prefactor * (G * (stress(2, 2) - stress(0, 0)) - F * (stress(1, 1) - stress(2, 2)));
282 :
283 664608 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
284 664608 : prefactor * 2.0 * N * stress(0, 1);
285 664608 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
286 664608 : prefactor * 2.0 * M * stress(0, 2);
287 664608 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
288 664608 : prefactor * 2.0 * L * stress(1, 2);
289 : }
290 : else
291 : {
292 1604980 : GenericDenseVector<is_ad> inelastic_strain_increment(6);
293 1604980 : GenericDenseVector<is_ad> Ms(6);
294 1604980 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_dev);
295 :
296 11234860 : for (unsigned int i = 0; i < 6; i++)
297 9629880 : inelastic_strain_increment(i) = Ms(i) * prefactor;
298 :
299 1604980 : inelasticStrainIncrement(0, 0) = inelastic_strain_increment(0);
300 1604980 : inelasticStrainIncrement(1, 1) = inelastic_strain_increment(1);
301 1604980 : inelasticStrainIncrement(2, 2) = inelastic_strain_increment(2);
302 :
303 1604980 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) = inelastic_strain_increment(3);
304 1604980 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) = inelastic_strain_increment(4);
305 1604980 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) = inelastic_strain_increment(5);
306 : }
307 :
308 2269588 : AnisotropicReturnCreepStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
309 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
310 :
311 3209960 : this->_effective_inelastic_strain[_qp] = this->_effective_inelastic_strain_old[_qp] + delta_gamma;
312 : }
313 :
314 : template <bool is_ad>
315 : void
316 2310908 : HillCreepStressUpdateTempl<is_ad>::computeStressFinalize(
317 : const GenericRankTwoTensor<is_ad> & creepStrainIncrement,
318 : const GenericReal<is_ad> & /*delta_gamma*/,
319 : GenericRankTwoTensor<is_ad> & stress_new,
320 : const GenericDenseVector<is_ad> & /*stress_dev*/,
321 : const GenericRankTwoTensor<is_ad> & stress_old,
322 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
323 : {
324 : // Class only valid for isotropic elasticity (for now)
325 2310908 : stress_new -= elasticity_tensor * creepStrainIncrement;
326 :
327 : // Compute the maximum time step allowed due to creep strain numerical integration error
328 3279480 : Real stress_dif = MetaPhysicL::raw_value(stress_new - stress_old).L2norm();
329 :
330 : // Get a representative value of the elasticity tensor
331 2310908 : Real elasticity_value =
332 : 1.0 / 3.0 *
333 2310908 : MetaPhysicL::raw_value((elasticity_tensor(0, 0, 0, 0) + elasticity_tensor(1, 1, 1, 1) +
334 : elasticity_tensor(2, 2, 2, 2)));
335 :
336 2310908 : if (std::abs(stress_dif) > TOLERANCE * TOLERANCE)
337 1989532 : this->_max_integration_error_time_step =
338 1989532 : _dt / (stress_dif / elasticity_value / this->_max_integration_error);
339 : else
340 321376 : this->_max_integration_error_time_step = std::numeric_limits<Real>::max();
341 2310908 : }
342 :
343 : template <bool is_ad>
344 : void
345 20988542 : HillCreepStressUpdateTempl<is_ad>::computeQsigmaChanged(
346 : GenericReal<is_ad> & qsigma_changed,
347 : const GenericDenseVector<is_ad> & stress_new,
348 : const GenericReal<is_ad> & delta_gamma,
349 : GenericRankTwoTensor<is_ad> & stress_changed)
350 : {
351 : GenericReal<is_ad> qsigma_square;
352 :
353 : // Hill constants, some constraints apply
354 20988542 : const Real & F = _hill_constants[_qp][0];
355 : const Real & G = _hill_constants[_qp][1];
356 : const Real & H = _hill_constants[_qp][2];
357 : const Real & L = _hill_constants[_qp][3];
358 : const Real & M = _hill_constants[_qp][4];
359 : const Real & N = _hill_constants[_qp][5];
360 :
361 20988542 : if (!this->_use_transformation)
362 : {
363 5928992 : qsigma_square = F * (stress_new(1) - stress_new(2)) * (stress_new(1) - stress_new(2));
364 5928992 : qsigma_square += G * (stress_new(2) - stress_new(0)) * (stress_new(2) - stress_new(0));
365 5928992 : qsigma_square += H * (stress_new(0) - stress_new(1)) * (stress_new(0) - stress_new(1));
366 5928992 : qsigma_square += 2 * L * stress_new(4) * stress_new(4);
367 5928992 : qsigma_square += 2 * M * stress_new(5) * stress_new(5);
368 5928992 : qsigma_square += 2 * N * stress_new(3) * stress_new(3);
369 : }
370 : else
371 : {
372 15059550 : GenericDenseVector<is_ad> Ms(6);
373 15059550 : (*_hill_tensor)[_qp].vector_mult(Ms, stress_new);
374 15059550 : qsigma_square = Ms.dot(stress_new);
375 : }
376 :
377 20988542 : GenericReal<is_ad> qsigma = std::sqrt(qsigma_square);
378 :
379 20988542 : if (!_anisotropic_elasticity)
380 30032620 : qsigma_changed = qsigma - 1.5 * _two_shear_modulus * delta_gamma;
381 : else
382 : {
383 43240 : GenericRankTwoTensor<is_ad> stress;
384 43240 : stress(0, 0) = stress_new(0);
385 43240 : stress(1, 1) = stress_new(1);
386 43240 : stress(2, 2) = stress_new(2);
387 43240 : stress(0, 1) = stress(1, 0) = stress_new(3);
388 43240 : stress(1, 2) = stress(2, 1) = stress_new(4);
389 43240 : stress(0, 2) = stress(2, 0) = stress_new(5);
390 :
391 43240 : GenericDenseVector<is_ad> b(6);
392 43240 : b(0) = H * (stress(0, 0) - stress(1, 1)) - G * (stress(2, 2) - stress(0, 0));
393 43240 : b(1) = F * (stress(1, 1) - stress(2, 2)) - H * (stress(0, 0) - stress(1, 1));
394 43240 : b(2) = G * (stress(2, 2) - stress(0, 0)) - F * (stress(1, 1) - stress(2, 2));
395 86480 : b(3) = 2.0 * N * stress(0, 1);
396 86480 : b(4) = 2.0 * L * stress(1, 2);
397 86480 : b(5) = 2.0 * M * stress(0, 2);
398 :
399 43240 : GenericRankTwoTensor<is_ad> inelasticStrainIncrement;
400 43240 : inelasticStrainIncrement(0, 0) = delta_gamma * b(0) / qsigma;
401 43240 : inelasticStrainIncrement(1, 1) = delta_gamma * b(1) / qsigma;
402 43240 : inelasticStrainIncrement(2, 2) = delta_gamma * b(2) / qsigma;
403 43240 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) = delta_gamma * b(3) / qsigma;
404 43240 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) = delta_gamma * b(4) / qsigma;
405 43240 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) = delta_gamma * b(5) / qsigma;
406 :
407 86480 : stress_changed = stress - _elasticity_tensor[_qp] * inelasticStrainIncrement;
408 :
409 : GenericReal<is_ad> qsigma_square_changed;
410 43240 : qsigma_square_changed = F * (stress_changed(1, 1) - stress_changed(2, 2)) *
411 : (stress_changed(1, 1) - stress_changed(2, 2));
412 43240 : qsigma_square_changed += G * (stress_changed(2, 2) - stress_changed(0, 0)) *
413 : (stress_changed(2, 2) - stress_changed(0, 0));
414 43240 : qsigma_square_changed += H * (stress_changed(0, 0) - stress_changed(1, 1)) *
415 : (stress_changed(0, 0) - stress_changed(1, 1));
416 86480 : qsigma_square_changed += 2 * L * stress_changed(1, 2) * stress_changed(1, 2);
417 86480 : qsigma_square_changed += 2 * M * stress_changed(0, 2) * stress_changed(0, 2);
418 86480 : qsigma_square_changed += 2 * N * stress_changed(0, 1) * stress_changed(0, 1);
419 43240 : qsigma_changed = std::sqrt(qsigma_square_changed);
420 : }
421 20988542 : }
422 :
423 : template class HillCreepStressUpdateTempl<false>;
424 : template class HillCreepStressUpdateTempl<true>;
|
