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