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 "SolidMechanicsPlasticTensileMulti.h"
11 : #include "RankFourTensor.h"
12 :
13 : // Following is for perturbing eigvenvalues. This looks really bodgy, but works quite well!
14 : #include "MooseRandom.h"
15 :
16 : registerMooseObject("SolidMechanicsApp", SolidMechanicsPlasticTensileMulti);
17 : registerMooseObjectRenamed("SolidMechanicsApp",
18 : TensorMechanicsPlasticTensileMulti,
19 : "01/01/2025 00:00",
20 : SolidMechanicsPlasticTensileMulti);
21 :
22 : InputParameters
23 234 : SolidMechanicsPlasticTensileMulti::validParams()
24 : {
25 234 : InputParameters params = SolidMechanicsPlasticModel::validParams();
26 234 : params.addClassDescription("Associative tensile plasticity with hardening/softening");
27 468 : params.addRequiredParam<UserObjectName>(
28 : "tensile_strength",
29 : "A SolidMechanicsHardening UserObject that defines hardening of the tensile strength");
30 468 : params.addParam<Real>("shift",
31 : "Yield surface is shifted by this amount to avoid problems with "
32 : "defining derivatives when eigenvalues are equal. If this is "
33 : "larger than f_tol, a warning will be issued. Default = f_tol.");
34 468 : params.addParam<unsigned int>("max_iterations",
35 468 : 10,
36 : "Maximum number of Newton-Raphson iterations "
37 : "allowed in the custom return-map algorithm. "
38 : " For highly nonlinear hardening this may "
39 : "need to be higher than 10.");
40 468 : params.addParam<bool>("use_custom_returnMap",
41 468 : true,
42 : "Whether to use the custom returnMap "
43 : "algorithm. Set to true if you are using "
44 : "isotropic elasticity.");
45 468 : params.addParam<bool>("use_custom_cto",
46 468 : true,
47 : "Whether to use the custom consistent tangent "
48 : "operator computations. Set to true if you are "
49 : "using isotropic elasticity.");
50 234 : return params;
51 0 : }
52 :
53 117 : SolidMechanicsPlasticTensileMulti::SolidMechanicsPlasticTensileMulti(
54 117 : const InputParameters & parameters)
55 : : SolidMechanicsPlasticModel(parameters),
56 117 : _strength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
57 234 : _max_iters(getParam<unsigned int>("max_iterations")),
58 337 : _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
59 234 : _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
60 351 : _use_custom_cto(getParam<bool>("use_custom_cto"))
61 : {
62 117 : if (_shift < 0)
63 0 : mooseError("Value of 'shift' in SolidMechanicsPlasticTensileMulti must not be negative\n");
64 117 : if (_shift > _f_tol)
65 : _console << "WARNING: value of 'shift' in SolidMechanicsPlasticTensileMulti is probably set "
66 0 : "too high"
67 0 : << std::endl;
68 : if (LIBMESH_DIM != 3)
69 : mooseError("SolidMechanicsPlasticTensileMulti is only defined for LIBMESH_DIM=3");
70 : MooseRandom::seed(0);
71 117 : }
72 :
73 : unsigned int
74 4175768 : SolidMechanicsPlasticTensileMulti::numberSurfaces() const
75 : {
76 4175768 : return 3;
77 : }
78 :
79 : void
80 165288 : SolidMechanicsPlasticTensileMulti::yieldFunctionV(const RankTwoTensor & stress,
81 : Real intnl,
82 : std::vector<Real> & f) const
83 : {
84 : std::vector<Real> eigvals;
85 165288 : stress.symmetricEigenvalues(eigvals);
86 165288 : const Real str = tensile_strength(intnl);
87 :
88 165288 : f.resize(3);
89 165288 : f[0] = eigvals[0] + _shift - str;
90 165288 : f[1] = eigvals[1] - str;
91 165288 : f[2] = eigvals[2] - _shift - str;
92 165288 : }
93 :
94 : void
95 118064 : SolidMechanicsPlasticTensileMulti::dyieldFunction_dstressV(
96 : const RankTwoTensor & stress, Real /*intnl*/, std::vector<RankTwoTensor> & df_dstress) const
97 : {
98 : std::vector<Real> eigvals;
99 118064 : stress.dsymmetricEigenvalues(eigvals, df_dstress);
100 :
101 118064 : if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
102 : {
103 : Real small_perturbation;
104 0 : RankTwoTensor shifted_stress = stress;
105 0 : while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
106 : {
107 0 : for (unsigned i = 0; i < 3; ++i)
108 0 : for (unsigned j = 0; j <= i; ++j)
109 : {
110 0 : small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
111 0 : shifted_stress(i, j) += small_perturbation;
112 0 : shifted_stress(j, i) += small_perturbation;
113 : }
114 0 : shifted_stress.dsymmetricEigenvalues(eigvals, df_dstress);
115 : }
116 : }
117 118064 : }
118 :
119 : void
120 36088 : SolidMechanicsPlasticTensileMulti::dyieldFunction_dintnlV(const RankTwoTensor & /*stress*/,
121 : Real intnl,
122 : std::vector<Real> & df_dintnl) const
123 : {
124 36088 : df_dintnl.assign(3, -dtensile_strength(intnl));
125 36088 : }
126 :
127 : void
128 78344 : SolidMechanicsPlasticTensileMulti::flowPotentialV(const RankTwoTensor & stress,
129 : Real intnl,
130 : std::vector<RankTwoTensor> & r) const
131 : {
132 : // This plasticity is associative so
133 78344 : dyieldFunction_dstressV(stress, intnl, r);
134 78344 : }
135 :
136 : void
137 36088 : SolidMechanicsPlasticTensileMulti::dflowPotential_dstressV(
138 : const RankTwoTensor & stress, Real /*intnl*/, std::vector<RankFourTensor> & dr_dstress) const
139 : {
140 36088 : stress.d2symmetricEigenvalues(dr_dstress);
141 36088 : }
142 :
143 : void
144 36088 : SolidMechanicsPlasticTensileMulti::dflowPotential_dintnlV(
145 : const RankTwoTensor & /*stress*/, Real /*intnl*/, std::vector<RankTwoTensor> & dr_dintnl) const
146 : {
147 36088 : dr_dintnl.assign(3, RankTwoTensor());
148 36088 : }
149 :
150 : Real
151 353544 : SolidMechanicsPlasticTensileMulti::tensile_strength(const Real internal_param) const
152 : {
153 353544 : return _strength.value(internal_param);
154 : }
155 :
156 : Real
157 89928 : SolidMechanicsPlasticTensileMulti::dtensile_strength(const Real internal_param) const
158 : {
159 89928 : return _strength.derivative(internal_param);
160 : }
161 :
162 : void
163 24496 : SolidMechanicsPlasticTensileMulti::activeConstraints(const std::vector<Real> & f,
164 : const RankTwoTensor & stress,
165 : Real intnl,
166 : const RankFourTensor & Eijkl,
167 : std::vector<bool> & act,
168 : RankTwoTensor & returned_stress) const
169 : {
170 : act.assign(3, false);
171 :
172 24496 : if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol)
173 : {
174 1368 : returned_stress = stress;
175 1368 : return;
176 : }
177 :
178 : Real returned_intnl;
179 23128 : std::vector<Real> dpm(3);
180 23128 : RankTwoTensor delta_dp;
181 23128 : std::vector<Real> yf(3);
182 : bool trial_stress_inadmissible;
183 23128 : doReturnMap(stress,
184 : intnl,
185 : Eijkl,
186 : 0.0,
187 : returned_stress,
188 : returned_intnl,
189 : dpm,
190 : delta_dp,
191 : yf,
192 : trial_stress_inadmissible);
193 :
194 92512 : for (unsigned i = 0; i < 3; ++i)
195 69384 : act[i] = (dpm[i] > 0);
196 23128 : }
197 :
198 : Real
199 0 : SolidMechanicsPlasticTensileMulti::dot(const std::vector<Real> & a,
200 : const std::vector<Real> & b) const
201 : {
202 0 : return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
203 : }
204 :
205 : Real
206 0 : SolidMechanicsPlasticTensileMulti::triple(const std::vector<Real> & a,
207 : const std::vector<Real> & b,
208 : const std::vector<Real> & c) const
209 : {
210 0 : return a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) +
211 0 : a[2] * (b[0] * c[1] - b[1] * c[0]);
212 : }
213 :
214 : std::string
215 51 : SolidMechanicsPlasticTensileMulti::modelName() const
216 : {
217 51 : return "TensileMulti";
218 : }
219 :
220 : bool
221 60144 : SolidMechanicsPlasticTensileMulti::returnMap(const RankTwoTensor & trial_stress,
222 : Real intnl_old,
223 : const RankFourTensor & E_ijkl,
224 : Real ep_plastic_tolerance,
225 : RankTwoTensor & returned_stress,
226 : Real & returned_intnl,
227 : std::vector<Real> & dpm,
228 : RankTwoTensor & delta_dp,
229 : std::vector<Real> & yf,
230 : bool & trial_stress_inadmissible) const
231 : {
232 60144 : if (!_use_custom_returnMap)
233 46608 : return SolidMechanicsPlasticModel::returnMap(trial_stress,
234 : intnl_old,
235 : E_ijkl,
236 : ep_plastic_tolerance,
237 : returned_stress,
238 : returned_intnl,
239 : dpm,
240 : delta_dp,
241 : yf,
242 46608 : trial_stress_inadmissible);
243 :
244 13536 : return doReturnMap(trial_stress,
245 : intnl_old,
246 : E_ijkl,
247 : ep_plastic_tolerance,
248 : returned_stress,
249 : returned_intnl,
250 : dpm,
251 : delta_dp,
252 : yf,
253 13536 : trial_stress_inadmissible);
254 : }
255 :
256 : bool
257 36664 : SolidMechanicsPlasticTensileMulti::doReturnMap(const RankTwoTensor & trial_stress,
258 : Real intnl_old,
259 : const RankFourTensor & E_ijkl,
260 : Real ep_plastic_tolerance,
261 : RankTwoTensor & returned_stress,
262 : Real & returned_intnl,
263 : std::vector<Real> & dpm,
264 : RankTwoTensor & delta_dp,
265 : std::vector<Real> & yf,
266 : bool & trial_stress_inadmissible) const
267 : {
268 : mooseAssert(dpm.size() == 3,
269 : "SolidMechanicsPlasticTensileMulti size of dpm should be 3 but it is " << dpm.size());
270 :
271 : std::vector<Real> eigvals;
272 36664 : RankTwoTensor eigvecs;
273 36664 : trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
274 36664 : eigvals[0] += _shift;
275 36664 : eigvals[2] -= _shift;
276 :
277 36664 : Real str = tensile_strength(intnl_old);
278 :
279 36664 : yf.resize(3);
280 36664 : yf[0] = eigvals[0] - str;
281 36664 : yf[1] = eigvals[1] - str;
282 36664 : yf[2] = eigvals[2] - str;
283 :
284 36664 : if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol)
285 : {
286 : // purely elastic (trial_stress, intnl_old)
287 376 : trial_stress_inadmissible = false;
288 376 : return true;
289 : }
290 :
291 36288 : trial_stress_inadmissible = true;
292 : delta_dp.zero();
293 : returned_stress.zero();
294 :
295 : // In the following i often assume that E_ijkl is
296 : // for an isotropic situation. This reduces FLOPS
297 : // substantially which is important since the returnMap
298 : // is potentially the most compute-intensive function
299 : // of a simulation.
300 : // In many comments i write the general expression, and
301 : // i hope that might guide future coders if they are
302 : // generalising to a non-istropic E_ijkl
303 :
304 : // n[alpha] = E_ijkl*r[alpha]_kl expressed in principal stress space
305 : // (alpha = 0, 1, 2, corresponding to the three surfaces)
306 : // Note that in principal stress space, the flow
307 : // directions are, expressed in 'vector' form,
308 : // r[0] = (1,0,0), r[1] = (0,1,0), r[2] = (0,0,1).
309 : // Similar for _n:
310 : // so _n[0] = E_ij00*r[0], _n[1] = E_ij11*r[1], _n[2] = E_ij22*r[2]
311 : // In the following I assume that the E_ijkl is
312 : // for an isotropic situation.
313 : // In the anisotropic situation, we couldn't express
314 : // the flow directions as vectors in the same principal
315 : // stress space as the stress: they'd be full rank-2 tensors
316 36288 : std::vector<RealVectorValue> n(3);
317 36288 : n[0](0) = E_ijkl(0, 0, 0, 0);
318 36288 : n[0](1) = E_ijkl(1, 1, 0, 0);
319 36288 : n[0](2) = E_ijkl(2, 2, 0, 0);
320 36288 : n[1](0) = E_ijkl(0, 0, 1, 1);
321 36288 : n[1](1) = E_ijkl(1, 1, 1, 1);
322 36288 : n[1](2) = E_ijkl(2, 2, 1, 1);
323 36288 : n[2](0) = E_ijkl(0, 0, 2, 2);
324 36288 : n[2](1) = E_ijkl(1, 1, 2, 2);
325 36288 : n[2](2) = E_ijkl(2, 2, 2, 2);
326 :
327 : // With non-zero Poisson's ratio and hardening
328 : // it is not computationally cheap to know whether
329 : // the trial stress will return to the tip, edge,
330 : // or plane. The following is correct for zero
331 : // Poisson's ratio and no hardening, and at least
332 : // gives a not-completely-stupid guess in the
333 : // more general case.
334 : // trial_order[0] = type of return to try first
335 : // trial_order[1] = type of return to try second
336 : // trial_order[2] = type of return to try third
337 : const unsigned int number_of_return_paths = 3;
338 36288 : std::vector<int> trial_order(number_of_return_paths);
339 36288 : if (yf[0] > _f_tol) // all the yield functions are positive, since eigvals are ordered eigvals[0]
340 : // <= eigvals[1] <= eigvals[2]
341 : {
342 7344 : trial_order[0] = tip;
343 7344 : trial_order[1] = edge;
344 7344 : trial_order[2] = plane;
345 : }
346 28944 : else if (yf[1] > _f_tol) // two yield functions are positive
347 : {
348 12448 : trial_order[0] = edge;
349 12448 : trial_order[1] = tip;
350 12448 : trial_order[2] = plane;
351 : }
352 : else
353 : {
354 16496 : trial_order[0] = plane;
355 16496 : trial_order[1] = edge;
356 16496 : trial_order[2] = tip;
357 : }
358 :
359 : unsigned trial;
360 : bool nr_converged = false;
361 53568 : for (trial = 0; trial < number_of_return_paths; ++trial)
362 : {
363 53568 : switch (trial_order[trial])
364 : {
365 12384 : case tip:
366 12384 : nr_converged = returnTip(eigvals, n, dpm, returned_stress, intnl_old, 0);
367 : break;
368 17072 : case edge:
369 17072 : nr_converged = returnEdge(eigvals, n, dpm, returned_stress, intnl_old, 0);
370 : break;
371 24112 : case plane:
372 24112 : nr_converged = returnPlane(eigvals, n, dpm, returned_stress, intnl_old, 0);
373 : break;
374 : }
375 :
376 53568 : str = tensile_strength(intnl_old + dpm[0] + dpm[1] + dpm[2]);
377 :
378 53568 : if (nr_converged && KuhnTuckerOK(returned_stress, dpm, str, ep_plastic_tolerance))
379 : break;
380 : }
381 :
382 36288 : if (trial == number_of_return_paths)
383 : {
384 : Moose::err << "Trial stress = \n";
385 0 : trial_stress.print(Moose::err);
386 : Moose::err << "Internal parameter = " << intnl_old << std::endl;
387 0 : mooseError("SolidMechanicsPlasticTensileMulti: FAILURE! You probably need to implement a "
388 : "line search\n");
389 : // failure - must place yield function values at trial stress into yf
390 : str = tensile_strength(intnl_old);
391 : yf[0] = eigvals[0] - str;
392 : yf[1] = eigvals[1] - str;
393 : yf[2] = eigvals[2] - str;
394 : return false;
395 : }
396 :
397 : // success
398 :
399 36288 : returned_intnl = intnl_old;
400 145152 : for (unsigned i = 0; i < 3; ++i)
401 : {
402 108864 : yf[i] = returned_stress(i, i) - str;
403 108864 : delta_dp(i, i) = dpm[i];
404 108864 : returned_intnl += dpm[i];
405 : }
406 36288 : returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
407 36288 : delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
408 : return true;
409 36664 : }
410 :
411 : bool
412 12384 : SolidMechanicsPlasticTensileMulti::returnTip(const std::vector<Real> & eigvals,
413 : const std::vector<RealVectorValue> & n,
414 : std::vector<Real> & dpm,
415 : RankTwoTensor & returned_stress,
416 : Real intnl_old,
417 : Real initial_guess) const
418 : {
419 : // The returned point is defined by f0=f1=f2=0.
420 : // that is, returned_stress = diag(str, str, str), where
421 : // str = tensile_strength(intnl),
422 : // where intnl = intnl_old + dpm[0] + dpm[1] + dpm[2]
423 : // The 3 plastic multipliers, dpm, are defiend by the normality condition
424 : // eigvals - str = dpm[0]*n[0] + dpm[1]*n[1] + dpm[2]*n[2]
425 : // (Kuhn-Tucker demands that all dpm are non-negative, but we leave
426 : // that checking for later.)
427 : // This is a vector equation with solution (A):
428 : // dpm[0] = triple(eigvals - str, n[1], n[2])/trip;
429 : // dpm[1] = triple(eigvals - str, n[2], n[0])/trip;
430 : // dpm[2] = triple(eigvals - str, n[0], n[1])/trip;
431 : // where trip = triple(n[0], n[1], n[2]).
432 : // By adding the three components of that solution together
433 : // we can get an equation for x = dpm[0] + dpm[1] + dpm[2],
434 : // and then our Newton-Raphson only involves one variable (x).
435 : // In the following, i specialise to the isotropic situation.
436 :
437 : Real x = initial_guess;
438 12384 : const Real denom = (n[0](0) - n[0](1)) * (n[0](0) + 2 * n[0](1));
439 12384 : Real str = tensile_strength(intnl_old + x);
440 :
441 12384 : if (_strength.modelName().compare("Constant") != 0)
442 : {
443 : // Finding x is expensive. Therefore
444 : // although x!=0 for Constant Hardening, solution (A)
445 : // demonstrates that we don't
446 : // actually need to know x to find the dpm for
447 : // Constant Hardening.
448 : //
449 : // However, for nontrivial Hardening, the following
450 : // is necessary
451 2736 : const Real eig = eigvals[0] + eigvals[1] + eigvals[2];
452 2736 : const Real bul = (n[0](0) + 2 * n[0](1));
453 :
454 : // and finally, the equation we want to solve is:
455 : // bul*x - eig + 3*str = 0
456 : // where str=tensile_strength(intnl_old + x)
457 : // and x = dpm[0] + dpm[1] + dpm[2]
458 : // (Note this has units of stress, so using _f_tol as a convergence check is reasonable.)
459 : // Use Netwon-Raphson with initial guess x = 0
460 2736 : Real residual = bul * x - eig + 3 * str;
461 : Real jacobian;
462 : unsigned int iter = 0;
463 : do
464 : {
465 8096 : jacobian = bul + 3 * dtensile_strength(intnl_old + x);
466 8096 : x += -residual / jacobian;
467 8096 : if (iter > _max_iters) // not converging
468 : return false;
469 8096 : str = tensile_strength(intnl_old + x);
470 8096 : residual = bul * x - eig + 3 * str;
471 8096 : iter++;
472 8096 : } while (residual * residual > _f_tol * _f_tol);
473 : }
474 :
475 : // The following is the solution (A) written above
476 : // (dpm[0] = triple(eigvals - str, n[1], n[2])/trip, etc)
477 : // in the isotropic situation
478 12384 : dpm[0] = (n[0](0) * (eigvals[0] - str) + n[0](1) * (eigvals[0] - eigvals[1] - eigvals[2] + str)) /
479 : denom;
480 12384 : dpm[1] = (n[0](0) * (eigvals[1] - str) + n[0](1) * (eigvals[1] - eigvals[2] - eigvals[0] + str)) /
481 : denom;
482 12384 : dpm[2] = (n[0](0) * (eigvals[2] - str) + n[0](1) * (eigvals[2] - eigvals[0] - eigvals[1] + str)) /
483 : denom;
484 12384 : returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = str;
485 12384 : return true;
486 : }
487 :
488 : bool
489 17072 : SolidMechanicsPlasticTensileMulti::returnEdge(const std::vector<Real> & eigvals,
490 : const std::vector<RealVectorValue> & n,
491 : std::vector<Real> & dpm,
492 : RankTwoTensor & returned_stress,
493 : Real intnl_old,
494 : Real initial_guess) const
495 : {
496 : // work out the point to which we would return, "a". It is defined by
497 : // f1 = 0 = f2, and the normality condition:
498 : // (eigvals - a).(n1 x n2) = 0,
499 : // where eigvals is the starting position
500 : // (it is a vector in principal stress space).
501 : // To get f1=0=f2, we need a = (a0, str, str), and a0 is found
502 : // by expanding the normality condition to yield:
503 : // a0 = (-str*n1xn2[1] - str*n1xn2[2] + edotn1xn2)/n1xn2[0];
504 : // where edotn1xn2 = eigvals.(n1 x n2)
505 : //
506 : // We need to find the plastic multipliers, dpm, defined by
507 : // eigvals - a = dpm[1]*n1 + dpm[2]*n2
508 : // For the isotropic case, and defining eminusa = eigvals - a,
509 : // the solution is easy:
510 : // dpm[0] = 0;
511 : // dpm[1] = (eminusa[1] - eminusa[0])/(n[1][1] - n[1][0]);
512 : // dpm[2] = (eminusa[2] - eminusa[0])/(n[2][2] - n[2][0]);
513 : //
514 : // Now specialise to the isotropic case. Define
515 : // x = dpm[1] + dpm[2] = (eigvals[1] + eigvals[2] - 2*str)/(n[0][0] + n[0][1])
516 : // Notice that the RHS is a function of x, so we solve using
517 : // Newton-Raphson starting with x=initial_guess
518 : Real x = initial_guess;
519 17072 : const Real denom = n[0](0) + n[0](1);
520 17072 : Real str = tensile_strength(intnl_old + x);
521 :
522 17072 : if (_strength.modelName().compare("Constant") != 0)
523 : {
524 : // Finding x is expensive. Therefore
525 : // although x!=0 for Constant Hardening, solution
526 : // for dpm above demonstrates that we don't
527 : // actually need to know x to find the dpm for
528 : // Constant Hardening.
529 : //
530 : // However, for nontrivial Hardening, the following
531 : // is necessary
532 2592 : const Real eig = eigvals[1] + eigvals[2];
533 2592 : Real residual = denom * x - eig + 2 * str;
534 : Real jacobian;
535 : unsigned int iter = 0;
536 : do
537 : {
538 7520 : jacobian = denom + 2 * dtensile_strength(intnl_old + x);
539 7520 : x += -residual / jacobian;
540 7520 : if (iter > _max_iters) // not converging
541 : return false;
542 7520 : str = tensile_strength(intnl_old + x);
543 7520 : residual = denom * x - eig + 2 * str;
544 7520 : iter++;
545 7520 : } while (residual * residual > _f_tol * _f_tol);
546 : }
547 :
548 17072 : dpm[0] = 0;
549 17072 : dpm[1] = ((eigvals[1] * n[0](0) - eigvals[2] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
550 17072 : dpm[2] = ((eigvals[2] * n[0](0) - eigvals[1] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
551 :
552 17072 : returned_stress(0, 0) = eigvals[0] - n[0](1) * (dpm[1] + dpm[2]);
553 17072 : returned_stress(1, 1) = returned_stress(2, 2) = str;
554 17072 : return true;
555 : }
556 :
557 : bool
558 24112 : SolidMechanicsPlasticTensileMulti::returnPlane(const std::vector<Real> & eigvals,
559 : const std::vector<RealVectorValue> & n,
560 : std::vector<Real> & dpm,
561 : RankTwoTensor & returned_stress,
562 : Real intnl_old,
563 : Real initial_guess) const
564 : {
565 : // the returned point, "a", is defined by f2=0 and
566 : // a = p - dpm[2]*n2.
567 : // This is a vector equation in
568 : // principal stress space, and dpm[2] is the third
569 : // plasticity multiplier (dpm[0]=0=dpm[1] for return
570 : // to the plane) and "p" is the starting
571 : // position (p=eigvals).
572 : // (Kuhn-Tucker demands that dpm[2]>=0, but we leave checking
573 : // that condition for later.)
574 : // Since f2=0, we must have a[2]=tensile_strength,
575 : // so we can just look at the [2] component of the
576 : // equation, which yields
577 : // n[2][2]*dpm[2] - eigvals[2] + str = 0
578 : // For hardening, str=tensile_strength(intnl_old+dpm[2]),
579 : // and we want to solve for dpm[2].
580 : // Use Newton-Raphson with initial guess dpm[2] = initial_guess
581 24112 : dpm[2] = initial_guess;
582 24112 : Real residual = n[2](2) * dpm[2] - eigvals[2] + tensile_strength(intnl_old + dpm[2]);
583 : Real jacobian;
584 : unsigned int iter = 0;
585 : do
586 : {
587 28840 : jacobian = n[2](2) + dtensile_strength(intnl_old + dpm[2]);
588 28840 : dpm[2] += -residual / jacobian;
589 28840 : if (iter > _max_iters) // not converging
590 : return false;
591 28840 : residual = n[2](2) * dpm[2] - eigvals[2] + tensile_strength(intnl_old + dpm[2]);
592 28840 : iter++;
593 28840 : } while (residual * residual > _f_tol * _f_tol);
594 :
595 24112 : dpm[0] = 0;
596 24112 : dpm[1] = 0;
597 24112 : returned_stress(0, 0) = eigvals[0] - dpm[2] * n[2](0);
598 24112 : returned_stress(1, 1) = eigvals[1] - dpm[2] * n[2](1);
599 24112 : returned_stress(2, 2) = eigvals[2] - dpm[2] * n[2](2);
600 24112 : return true;
601 : }
602 :
603 : bool
604 53568 : SolidMechanicsPlasticTensileMulti::KuhnTuckerOK(const RankTwoTensor & returned_diagonal_stress,
605 : const std::vector<Real> & dpm,
606 : Real str,
607 : Real ep_plastic_tolerance) const
608 : {
609 170048 : for (unsigned i = 0; i < 3; ++i)
610 133760 : if (!SolidMechanicsPlasticModel::KuhnTuckerSingleSurface(
611 133760 : returned_diagonal_stress(i, i) - str, dpm[i], ep_plastic_tolerance))
612 : return false;
613 : return true;
614 : }
615 :
616 : RankFourTensor
617 9384 : SolidMechanicsPlasticTensileMulti::consistentTangentOperator(
618 : const RankTwoTensor & trial_stress,
619 : Real intnl_old,
620 : const RankTwoTensor & stress,
621 : Real intnl,
622 : const RankFourTensor & E_ijkl,
623 : const std::vector<Real> & cumulative_pm) const
624 : {
625 9384 : if (!_use_custom_cto)
626 0 : return SolidMechanicsPlasticModel::consistentTangentOperator(
627 : trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
628 :
629 : mooseAssert(cumulative_pm.size() == 3,
630 : "SolidMechanicsPlasticTensileMulti size of cumulative_pm should be 3 but it is "
631 : << cumulative_pm.size());
632 :
633 9384 : if (cumulative_pm[2] <= 0) // All cumulative_pm are non-positive, so this is admissible
634 0 : return E_ijkl;
635 :
636 : // Need the eigenvalues at the returned configuration
637 : std::vector<Real> eigvals;
638 9384 : stress.symmetricEigenvalues(eigvals);
639 :
640 : // need to rotate to and from principal stress space
641 : // using the eigenvectors of the trial configuration
642 : // (not the returned configuration).
643 : std::vector<Real> trial_eigvals;
644 9384 : RankTwoTensor trial_eigvecs;
645 9384 : trial_stress.symmetricEigenvaluesEigenvectors(trial_eigvals, trial_eigvecs);
646 :
647 : // The returnMap will have returned to the Tip, Edge or
648 : // Plane. The consistentTangentOperator describes the
649 : // change in stress for an arbitrary change in applied
650 : // strain. I assume that the change in strain will not
651 : // change the type of return (Tip remains Tip, Edge remains
652 : // Edge, Plane remains Plane).
653 : // I assume isotropic elasticity.
654 : //
655 : // The consistent tangent operator is a little different
656 : // than cases where no rotation to principal stress space
657 : // is made during the returnMap. Let S_ij be the stress
658 : // in original coordinates, and s_ij be the stress in the
659 : // principal stress coordinates, so that
660 : // s_ij = diag(eigvals[0], eigvals[1], eigvals[2])
661 : // We want dS_ij under an arbitrary change in strain (ep->ep+dep)
662 : // dS = S(ep+dep) - S(ep)
663 : // = R(ep+dep) s(ep+dep) R(ep+dep)^T - R(ep) s(ep) R(ep)^T
664 : // Here R = the rotation to principal-stress space, ie
665 : // R_ij = eigvecs[i][j] = i^th component of j^th eigenvector
666 : // Expanding to first order in dep,
667 : // dS = R(ep) (s(ep+dep) - s(ep)) R(ep)^T
668 : // + dR/dep s(ep) R^T + R(ep) s(ep) dR^T/dep
669 : // The first line is all that is usually calculated in the
670 : // consistent tangent operator calculation, and is called
671 : // cto below.
672 : // The second line involves changes in the eigenvectors, and
673 : // is called sec below.
674 :
675 9384 : RankFourTensor cto;
676 9384 : const Real hard = dtensile_strength(intnl);
677 9384 : const Real la = E_ijkl(0, 0, 1, 1);
678 9384 : const Real mu = 0.5 * (E_ijkl(0, 0, 0, 0) - la);
679 :
680 9384 : if (cumulative_pm[1] <= 0)
681 : {
682 : // only cumulative_pm[2] is positive, so this is return to the Plane
683 3112 : const Real denom = hard + la + 2 * mu;
684 3112 : const Real al = la * la / denom;
685 3112 : const Real be = la * (la + 2 * mu) / denom;
686 3112 : const Real ga = hard * (la + 2 * mu) / denom;
687 3112 : std::vector<Real> comps(9);
688 3112 : comps[0] = comps[4] = la + 2 * mu - al;
689 3112 : comps[1] = comps[3] = la - al;
690 3112 : comps[2] = comps[5] = comps[6] = comps[7] = la - be;
691 3112 : comps[8] = ga;
692 3112 : cto.fillFromInputVector(comps, RankFourTensor::principal);
693 3112 : }
694 6272 : else if (cumulative_pm[0] <= 0)
695 : {
696 : // both cumulative_pm[2] and cumulative_pm[1] are positive, so Edge
697 4080 : const Real denom = 2 * hard + 2 * la + 2 * mu;
698 4080 : const Real al = hard * 2 * la / denom;
699 4080 : const Real be = hard * (2 * la + 2 * mu) / denom;
700 4080 : std::vector<Real> comps(9);
701 4080 : comps[0] = la + 2 * mu - 2 * la * la / denom;
702 4080 : comps[1] = comps[2] = al;
703 4080 : comps[3] = comps[6] = al;
704 4080 : comps[4] = comps[5] = comps[7] = comps[8] = be;
705 4080 : cto.fillFromInputVector(comps, RankFourTensor::principal);
706 4080 : }
707 : else
708 : {
709 : // all cumulative_pm are positive, so Tip
710 2192 : const Real denom = 3 * hard + 3 * la + 2 * mu;
711 2192 : std::vector<Real> comps(2);
712 2192 : comps[0] = hard * (3 * la + 2 * mu) / denom;
713 2192 : comps[1] = 0;
714 2192 : cto.fillFromInputVector(comps, RankFourTensor::symmetric_isotropic);
715 2192 : }
716 :
717 9384 : cto.rotate(trial_eigvecs);
718 :
719 : // drdsig = change in eigenvectors under a small stress change
720 : // drdsig(i,j,m,n) = dR(i,j)/dS_mn
721 : // The formula below is fairly easily derived:
722 : // S R = R s, so taking the variation
723 : // dS R + S dR = dR s + R ds, and multiplying by R^T
724 : // R^T dS R + R^T S dR = R^T dR s + ds .... (eqn 1)
725 : // I demand that RR^T = 1 = R^T R, and also that
726 : // (R+dR)(R+dR)^T = 1 = (R+dT)^T (R+dR), which means
727 : // that dR = R*c, for some antisymmetric c, so Eqn1 reads
728 : // R^T dS R + s c = c s + ds
729 : // Grabbing the components of this gives ds/dS (already
730 : // in RankTwoTensor), and c, which is:
731 : // dR_ik/dS_mn = drdsig(i, k, m, n) = trial_eigvecs(m, b)*trial_eigvecs(n, k)*trial_eigvecs(i,
732 : // b)/(trial_eigvals[k] - trial_eigvals[b]);
733 : // (sum over b!=k).
734 :
735 9384 : RankFourTensor drdsig;
736 37536 : for (unsigned k = 0; k < 3; ++k)
737 112608 : for (unsigned b = 0; b < 3; ++b)
738 : {
739 84456 : if (b == k)
740 28152 : continue;
741 225216 : for (unsigned m = 0; m < 3; ++m)
742 675648 : for (unsigned n = 0; n < 3; ++n)
743 2026944 : for (unsigned i = 0; i < 3; ++i)
744 1520208 : drdsig(i, k, m, n) += trial_eigvecs(m, b) * trial_eigvecs(n, k) * trial_eigvecs(i, b) /
745 1520208 : (trial_eigvals[k] - trial_eigvals[b]);
746 : }
747 :
748 : // With diagla = diag(eigvals[0], eigvals[1], digvals[2])
749 : // The following implements
750 : // ans(i, j, a, b) += (drdsig(i, k, m, n)*trial_eigvecs(j, l)*diagla(k, l) + trial_eigvecs(i,
751 : // k)*drdsig(j, l, m, n)*diagla(k, l))*E_ijkl(m, n, a, b);
752 : // (sum over k, l, m and n)
753 :
754 9384 : RankFourTensor ans;
755 37536 : for (unsigned i = 0; i < 3; ++i)
756 112608 : for (unsigned j = 0; j < 3; ++j)
757 337824 : for (unsigned a = 0; a < 3; ++a)
758 1013472 : for (unsigned k = 0; k < 3; ++k)
759 3040416 : for (unsigned m = 0; m < 3; ++m)
760 2280312 : ans(i, j, a, a) += (drdsig(i, k, m, m) * trial_eigvecs(j, k) +
761 2280312 : trial_eigvecs(i, k) * drdsig(j, k, m, m)) *
762 2280312 : eigvals[k] * la; // E_ijkl(m, n, a, b) = la*(m==n)*(a==b);
763 :
764 37536 : for (unsigned i = 0; i < 3; ++i)
765 112608 : for (unsigned j = 0; j < 3; ++j)
766 337824 : for (unsigned a = 0; a < 3; ++a)
767 1013472 : for (unsigned b = 0; b < 3; ++b)
768 3040416 : for (unsigned k = 0; k < 3; ++k)
769 : {
770 2280312 : ans(i, j, a, b) += (drdsig(i, k, a, b) * trial_eigvecs(j, k) +
771 2280312 : trial_eigvecs(i, k) * drdsig(j, k, a, b)) *
772 2280312 : eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==a)*(n==b)
773 2280312 : ans(i, j, a, b) += (drdsig(i, k, b, a) * trial_eigvecs(j, k) +
774 2280312 : trial_eigvecs(i, k) * drdsig(j, k, b, a)) *
775 2280312 : eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==b)*(n==a)
776 : }
777 :
778 9384 : return cto + ans;
779 9384 : }
780 :
781 : bool
782 0 : SolidMechanicsPlasticTensileMulti::useCustomReturnMap() const
783 : {
784 0 : return _use_custom_returnMap;
785 : }
786 :
787 : bool
788 9384 : SolidMechanicsPlasticTensileMulti::useCustomCTO() const
789 : {
790 9384 : return _use_custom_cto;
791 : }
|
