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 208 : SolidMechanicsPlasticTensileMulti::validParams()
24 : {
25 208 : InputParameters params = SolidMechanicsPlasticModel::validParams();
26 208 : params.addClassDescription("Associative tensile plasticity with hardening/softening");
27 416 : params.addRequiredParam<UserObjectName>(
28 : "tensile_strength",
29 : "A SolidMechanicsHardening UserObject that defines hardening of the tensile strength");
30 416 : 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 416 : params.addParam<unsigned int>("max_iterations",
35 416 : 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 416 : params.addParam<bool>("use_custom_returnMap",
41 416 : true,
42 : "Whether to use the custom returnMap "
43 : "algorithm. Set to true if you are using "
44 : "isotropic elasticity.");
45 416 : params.addParam<bool>("use_custom_cto",
46 416 : true,
47 : "Whether to use the custom consistent tangent "
48 : "operator computations. Set to true if you are "
49 : "using isotropic elasticity.");
50 208 : return params;
51 0 : }
52 :
53 104 : SolidMechanicsPlasticTensileMulti::SolidMechanicsPlasticTensileMulti(
54 104 : const InputParameters & parameters)
55 : : SolidMechanicsPlasticModel(parameters),
56 104 : _strength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
57 208 : _max_iters(getParam<unsigned int>("max_iterations")),
58 300 : _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
59 208 : _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
60 312 : _use_custom_cto(getParam<bool>("use_custom_cto"))
61 : {
62 104 : if (_shift < 0)
63 0 : mooseError("Value of 'shift' in SolidMechanicsPlasticTensileMulti must not be negative\n");
64 104 : 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 104 : }
72 :
73 : unsigned int
74 2901320 : SolidMechanicsPlasticTensileMulti::numberSurfaces() const
75 : {
76 2901320 : return 3;
77 : }
78 :
79 : void
80 113424 : SolidMechanicsPlasticTensileMulti::yieldFunctionV(const RankTwoTensor & stress,
81 : Real intnl,
82 : std::vector<Real> & f) const
83 : {
84 : std::vector<Real> eigvals;
85 113424 : stress.symmetricEigenvalues(eigvals);
86 113424 : const Real str = tensile_strength(intnl);
87 :
88 113424 : f.resize(3);
89 113424 : f[0] = eigvals[0] + _shift - str;
90 113424 : f[1] = eigvals[1] - str;
91 113424 : f[2] = eigvals[2] - _shift - str;
92 113424 : }
93 :
94 : void
95 79008 : SolidMechanicsPlasticTensileMulti::dyieldFunction_dstressV(
96 : const RankTwoTensor & stress, Real /*intnl*/, std::vector<RankTwoTensor> & df_dstress) const
97 : {
98 : std::vector<Real> eigvals;
99 79008 : stress.dsymmetricEigenvalues(eigvals, df_dstress);
100 :
101 79008 : 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 79008 : }
118 :
119 : void
120 24144 : SolidMechanicsPlasticTensileMulti::dyieldFunction_dintnlV(const RankTwoTensor & /*stress*/,
121 : Real intnl,
122 : std::vector<Real> & df_dintnl) const
123 : {
124 24144 : df_dintnl.assign(3, -dtensile_strength(intnl));
125 24144 : }
126 :
127 : void
128 52400 : SolidMechanicsPlasticTensileMulti::flowPotentialV(const RankTwoTensor & stress,
129 : Real intnl,
130 : std::vector<RankTwoTensor> & r) const
131 : {
132 : // This plasticity is associative so
133 52400 : dyieldFunction_dstressV(stress, intnl, r);
134 52400 : }
135 :
136 : void
137 24144 : SolidMechanicsPlasticTensileMulti::dflowPotential_dstressV(
138 : const RankTwoTensor & stress, Real /*intnl*/, std::vector<RankFourTensor> & dr_dstress) const
139 : {
140 24144 : stress.d2symmetricEigenvalues(dr_dstress);
141 24144 : }
142 :
143 : void
144 24144 : SolidMechanicsPlasticTensileMulti::dflowPotential_dintnlV(
145 : const RankTwoTensor & /*stress*/, Real /*intnl*/, std::vector<RankTwoTensor> & dr_dintnl) const
146 : {
147 24144 : dr_dintnl.assign(3, RankTwoTensor());
148 24144 : }
149 :
150 : Real
151 260336 : SolidMechanicsPlasticTensileMulti::tensile_strength(const Real internal_param) const
152 : {
153 260336 : return _strength.value(internal_param);
154 : }
155 :
156 : Real
157 70064 : SolidMechanicsPlasticTensileMulti::dtensile_strength(const Real internal_param) const
158 : {
159 70064 : return _strength.derivative(internal_param);
160 : }
161 :
162 : void
163 16960 : 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 16960 : if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol)
173 : {
174 912 : returned_stress = stress;
175 912 : return;
176 : }
177 :
178 : Real returned_intnl;
179 16048 : std::vector<Real> dpm(3);
180 16048 : RankTwoTensor delta_dp;
181 16048 : std::vector<Real> yf(3);
182 : bool trial_stress_inadmissible;
183 16048 : 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 64192 : for (unsigned i = 0; i < 3; ++i)
195 48144 : act[i] = (dpm[i] > 0);
196 : }
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 48 : SolidMechanicsPlasticTensileMulti::modelName() const
216 : {
217 48 : return "TensileMulti";
218 : }
219 :
220 : bool
221 43872 : 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 43872 : if (!_use_custom_returnMap)
233 31072 : 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 31072 : trial_stress_inadmissible);
243 :
244 12800 : 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 12800 : trial_stress_inadmissible);
254 : }
255 :
256 : bool
257 28848 : 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 28848 : RankTwoTensor eigvecs;
273 28848 : trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
274 28848 : eigvals[0] += _shift;
275 28848 : eigvals[2] -= _shift;
276 :
277 28848 : Real str = tensile_strength(intnl_old);
278 :
279 28848 : yf.resize(3);
280 28848 : yf[0] = eigvals[0] - str;
281 28848 : yf[1] = eigvals[1] - str;
282 28848 : yf[2] = eigvals[2] - str;
283 :
284 28848 : if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol)
285 : {
286 : // purely elastic (trial_stress, intnl_old)
287 336 : trial_stress_inadmissible = false;
288 336 : return true;
289 : }
290 :
291 28512 : 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 28512 : std::vector<RealVectorValue> n(3);
317 28512 : n[0](0) = E_ijkl(0, 0, 0, 0);
318 28512 : n[0](1) = E_ijkl(1, 1, 0, 0);
319 28512 : n[0](2) = E_ijkl(2, 2, 0, 0);
320 28512 : n[1](0) = E_ijkl(0, 0, 1, 1);
321 28512 : n[1](1) = E_ijkl(1, 1, 1, 1);
322 28512 : n[1](2) = E_ijkl(2, 2, 1, 1);
323 28512 : n[2](0) = E_ijkl(0, 0, 2, 2);
324 28512 : n[2](1) = E_ijkl(1, 1, 2, 2);
325 28512 : 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 28512 : std::vector<int> trial_order(number_of_return_paths);
339 28512 : if (yf[0] > _f_tol) // all the yield functions are positive, since eigvals are ordered eigvals[0]
340 : // <= eigvals[1] <= eigvals[2]
341 : {
342 5696 : trial_order[0] = tip;
343 5696 : trial_order[1] = edge;
344 5696 : trial_order[2] = plane;
345 : }
346 22816 : else if (yf[1] > _f_tol) // two yield functions are positive
347 : {
348 9472 : trial_order[0] = edge;
349 9472 : trial_order[1] = tip;
350 9472 : trial_order[2] = plane;
351 : }
352 : else
353 : {
354 13344 : trial_order[0] = plane;
355 13344 : trial_order[1] = edge;
356 13344 : trial_order[2] = tip;
357 : }
358 :
359 : unsigned trial;
360 : bool nr_converged = false;
361 40416 : for (trial = 0; trial < number_of_return_paths; ++trial)
362 : {
363 40416 : switch (trial_order[trial])
364 : {
365 9056 : case tip:
366 9056 : nr_converged = returnTip(eigvals, n, dpm, returned_stress, intnl_old, 0);
367 : break;
368 12768 : case edge:
369 12768 : nr_converged = returnEdge(eigvals, n, dpm, returned_stress, intnl_old, 0);
370 : break;
371 18592 : case plane:
372 18592 : nr_converged = returnPlane(eigvals, n, dpm, returned_stress, intnl_old, 0);
373 : break;
374 : }
375 :
376 40416 : str = tensile_strength(intnl_old + dpm[0] + dpm[1] + dpm[2]);
377 :
378 40416 : if (nr_converged && KuhnTuckerOK(returned_stress, dpm, str, ep_plastic_tolerance))
379 : break;
380 : }
381 :
382 28512 : 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 28512 : returned_intnl = intnl_old;
400 114048 : for (unsigned i = 0; i < 3; ++i)
401 : {
402 85536 : yf[i] = returned_stress(i, i) - str;
403 85536 : delta_dp(i, i) = dpm[i];
404 85536 : returned_intnl += dpm[i];
405 : }
406 28512 : returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
407 28512 : delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
408 : return true;
409 : }
410 :
411 : bool
412 9056 : 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 9056 : const Real denom = (n[0](0) - n[0](1)) * (n[0](0) + 2 * n[0](1));
439 9056 : Real str = tensile_strength(intnl_old + x);
440 :
441 18112 : 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 2496 : const Real eig = eigvals[0] + eigvals[1] + eigvals[2];
452 2496 : 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 2496 : Real residual = bul * x - eig + 3 * str;
461 : Real jacobian;
462 : unsigned int iter = 0;
463 : do
464 : {
465 7328 : jacobian = bul + 3 * dtensile_strength(intnl_old + x);
466 7328 : x += -residual / jacobian;
467 7328 : if (iter > _max_iters) // not converging
468 : return false;
469 7328 : str = tensile_strength(intnl_old + x);
470 7328 : residual = bul * x - eig + 3 * str;
471 7328 : iter++;
472 7328 : } 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 9056 : dpm[0] = (n[0](0) * (eigvals[0] - str) + n[0](1) * (eigvals[0] - eigvals[1] - eigvals[2] + str)) /
479 : denom;
480 9056 : dpm[1] = (n[0](0) * (eigvals[1] - str) + n[0](1) * (eigvals[1] - eigvals[2] - eigvals[0] + str)) /
481 : denom;
482 9056 : dpm[2] = (n[0](0) * (eigvals[2] - str) + n[0](1) * (eigvals[2] - eigvals[0] - eigvals[1] + str)) /
483 : denom;
484 9056 : returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = str;
485 9056 : return true;
486 : }
487 :
488 : bool
489 12768 : 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 12768 : const Real denom = n[0](0) + n[0](1);
520 12768 : Real str = tensile_strength(intnl_old + x);
521 :
522 25536 : 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 2400 : const Real eig = eigvals[1] + eigvals[2];
533 2400 : Real residual = denom * x - eig + 2 * str;
534 : Real jacobian;
535 : unsigned int iter = 0;
536 : do
537 : {
538 6944 : jacobian = denom + 2 * dtensile_strength(intnl_old + x);
539 6944 : x += -residual / jacobian;
540 6944 : if (iter > _max_iters) // not converging
541 : return false;
542 6944 : str = tensile_strength(intnl_old + x);
543 6944 : residual = denom * x - eig + 2 * str;
544 6944 : iter++;
545 6944 : } while (residual * residual > _f_tol * _f_tol);
546 : }
547 :
548 12768 : dpm[0] = 0;
549 12768 : dpm[1] = ((eigvals[1] * n[0](0) - eigvals[2] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
550 12768 : dpm[2] = ((eigvals[2] * n[0](0) - eigvals[1] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
551 :
552 12768 : returned_stress(0, 0) = eigvals[0] - n[0](1) * (dpm[1] + dpm[2]);
553 12768 : returned_stress(1, 1) = returned_stress(2, 2) = str;
554 12768 : return true;
555 : }
556 :
557 : bool
558 18592 : 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 18592 : dpm[2] = initial_guess;
582 18592 : 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 22960 : jacobian = n[2](2) + dtensile_strength(intnl_old + dpm[2]);
588 22960 : dpm[2] += -residual / jacobian;
589 22960 : if (iter > _max_iters) // not converging
590 : return false;
591 22960 : residual = n[2](2) * dpm[2] - eigvals[2] + tensile_strength(intnl_old + dpm[2]);
592 22960 : iter++;
593 22960 : } while (residual * residual > _f_tol * _f_tol);
594 :
595 18592 : dpm[0] = 0;
596 18592 : dpm[1] = 0;
597 18592 : returned_stress(0, 0) = eigvals[0] - dpm[2] * n[2](0);
598 18592 : returned_stress(1, 1) = eigvals[1] - dpm[2] * n[2](1);
599 18592 : returned_stress(2, 2) = eigvals[2] - dpm[2] * n[2](2);
600 18592 : return true;
601 : }
602 :
603 : bool
604 40416 : SolidMechanicsPlasticTensileMulti::KuhnTuckerOK(const RankTwoTensor & returned_diagonal_stress,
605 : const std::vector<Real> & dpm,
606 : Real str,
607 : Real ep_plastic_tolerance) const
608 : {
609 131200 : for (unsigned i = 0; i < 3; ++i)
610 102688 : if (!SolidMechanicsPlasticModel::KuhnTuckerSingleSurface(
611 102688 : returned_diagonal_stress(i, i) - str, dpm[i], ep_plastic_tolerance))
612 : return false;
613 : return true;
614 : }
615 :
616 : RankFourTensor
617 8688 : 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 8688 : 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 8688 : 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 8688 : 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 8688 : RankTwoTensor trial_eigvecs;
645 8688 : 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 8688 : RankFourTensor cto;
676 8688 : const Real hard = dtensile_strength(intnl);
677 8688 : const Real la = E_ijkl(0, 0, 1, 1);
678 8688 : const Real mu = 0.5 * (E_ijkl(0, 0, 0, 0) - la);
679 :
680 8688 : if (cumulative_pm[1] <= 0)
681 : {
682 : // only cumulative_pm[2] is positive, so this is return to the Plane
683 2704 : const Real denom = hard + la + 2 * mu;
684 2704 : const Real al = la * la / denom;
685 2704 : const Real be = la * (la + 2 * mu) / denom;
686 2704 : const Real ga = hard * (la + 2 * mu) / denom;
687 2704 : std::vector<Real> comps(9);
688 2704 : comps[0] = comps[4] = la + 2 * mu - al;
689 2704 : comps[1] = comps[3] = la - al;
690 2704 : comps[2] = comps[5] = comps[6] = comps[7] = la - be;
691 2704 : comps[8] = ga;
692 2704 : cto.fillFromInputVector(comps, RankFourTensor::principal);
693 : }
694 5984 : else if (cumulative_pm[0] <= 0)
695 : {
696 : // both cumulative_pm[2] and cumulative_pm[1] are positive, so Edge
697 3936 : const Real denom = 2 * hard + 2 * la + 2 * mu;
698 3936 : const Real al = hard * 2 * la / denom;
699 3936 : const Real be = hard * (2 * la + 2 * mu) / denom;
700 3936 : std::vector<Real> comps(9);
701 3936 : comps[0] = la + 2 * mu - 2 * la * la / denom;
702 3936 : comps[1] = comps[2] = al;
703 3936 : comps[3] = comps[6] = al;
704 3936 : comps[4] = comps[5] = comps[7] = comps[8] = be;
705 3936 : cto.fillFromInputVector(comps, RankFourTensor::principal);
706 : }
707 : else
708 : {
709 : // all cumulative_pm are positive, so Tip
710 2048 : const Real denom = 3 * hard + 3 * la + 2 * mu;
711 2048 : std::vector<Real> comps(2);
712 2048 : comps[0] = hard * (3 * la + 2 * mu) / denom;
713 2048 : comps[1] = 0;
714 2048 : cto.fillFromInputVector(comps, RankFourTensor::symmetric_isotropic);
715 : }
716 :
717 8688 : 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 8688 : RankFourTensor drdsig;
736 34752 : for (unsigned k = 0; k < 3; ++k)
737 104256 : for (unsigned b = 0; b < 3; ++b)
738 : {
739 78192 : if (b == k)
740 26064 : continue;
741 208512 : for (unsigned m = 0; m < 3; ++m)
742 625536 : for (unsigned n = 0; n < 3; ++n)
743 1876608 : for (unsigned i = 0; i < 3; ++i)
744 1407456 : drdsig(i, k, m, n) += trial_eigvecs(m, b) * trial_eigvecs(n, k) * trial_eigvecs(i, b) /
745 1407456 : (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 8688 : RankFourTensor ans;
755 34752 : for (unsigned i = 0; i < 3; ++i)
756 104256 : for (unsigned j = 0; j < 3; ++j)
757 312768 : for (unsigned a = 0; a < 3; ++a)
758 938304 : for (unsigned k = 0; k < 3; ++k)
759 2814912 : for (unsigned m = 0; m < 3; ++m)
760 2111184 : ans(i, j, a, a) += (drdsig(i, k, m, m) * trial_eigvecs(j, k) +
761 2111184 : trial_eigvecs(i, k) * drdsig(j, k, m, m)) *
762 2111184 : eigvals[k] * la; // E_ijkl(m, n, a, b) = la*(m==n)*(a==b);
763 :
764 34752 : for (unsigned i = 0; i < 3; ++i)
765 104256 : for (unsigned j = 0; j < 3; ++j)
766 312768 : for (unsigned a = 0; a < 3; ++a)
767 938304 : for (unsigned b = 0; b < 3; ++b)
768 2814912 : for (unsigned k = 0; k < 3; ++k)
769 : {
770 2111184 : ans(i, j, a, b) += (drdsig(i, k, a, b) * trial_eigvecs(j, k) +
771 2111184 : trial_eigvecs(i, k) * drdsig(j, k, a, b)) *
772 2111184 : eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==a)*(n==b)
773 2111184 : ans(i, j, a, b) += (drdsig(i, k, b, a) * trial_eigvecs(j, k) +
774 2111184 : trial_eigvecs(i, k) * drdsig(j, k, b, a)) *
775 2111184 : eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==b)*(n==a)
776 : }
777 :
778 8688 : return cto + ans;
779 : }
780 :
781 : bool
782 0 : SolidMechanicsPlasticTensileMulti::useCustomReturnMap() const
783 : {
784 0 : return _use_custom_returnMap;
785 : }
786 :
787 : bool
788 8688 : SolidMechanicsPlasticTensileMulti::useCustomCTO() const
789 : {
790 8688 : return _use_custom_cto;
791 : }
|
