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 "TensorMechanicsPlasticMohrCoulombMulti.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 : #include "libmesh/utility.h"
16 :
17 : registerMooseObject("TensorMechanicsApp", TensorMechanicsPlasticMohrCoulombMulti);
18 :
19 : InputParameters
20 68 : TensorMechanicsPlasticMohrCoulombMulti::validParams()
21 : {
22 68 : InputParameters params = TensorMechanicsPlasticModel::validParams();
23 68 : params.addClassDescription("Non-associative Mohr-Coulomb plasticity with hardening/softening");
24 136 : params.addRequiredParam<UserObjectName>(
25 : "cohesion", "A TensorMechanicsHardening UserObject that defines hardening of the cohesion");
26 136 : params.addRequiredParam<UserObjectName>("friction_angle",
27 : "A TensorMechanicsHardening UserObject "
28 : "that defines hardening of the "
29 : "friction angle (in radians)");
30 136 : params.addRequiredParam<UserObjectName>("dilation_angle",
31 : "A TensorMechanicsHardening UserObject "
32 : "that defines hardening of the "
33 : "dilation angle (in radians)");
34 136 : params.addParam<unsigned int>("max_iterations",
35 136 : 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 136 : params.addParam<Real>("shift",
41 : "Yield surface is shifted by this amount to avoid problems with "
42 : "defining derivatives when eigenvalues are equal. If this is "
43 : "larger than f_tol, a warning will be issued. This may be set "
44 : "very small when using the custom returnMap. Default = f_tol.");
45 136 : params.addParam<bool>("use_custom_returnMap",
46 136 : true,
47 : "Use a custom return-map algorithm for this "
48 : "plasticity model, which may speed up "
49 : "computations considerably. Set to true "
50 : "only for isotropic elasticity with no "
51 : "hardening of the dilation angle. In this "
52 : "case you may set 'shift' very small.");
53 :
54 68 : return params;
55 0 : }
56 :
57 34 : TensorMechanicsPlasticMohrCoulombMulti::TensorMechanicsPlasticMohrCoulombMulti(
58 34 : const InputParameters & parameters)
59 : : TensorMechanicsPlasticModel(parameters),
60 34 : _cohesion(getUserObject<TensorMechanicsHardeningModel>("cohesion")),
61 34 : _phi(getUserObject<TensorMechanicsHardeningModel>("friction_angle")),
62 34 : _psi(getUserObject<TensorMechanicsHardeningModel>("dilation_angle")),
63 68 : _max_iters(getParam<unsigned int>("max_iterations")),
64 84 : _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
65 102 : _use_custom_returnMap(getParam<bool>("use_custom_returnMap"))
66 : {
67 34 : if (_shift < 0)
68 0 : mooseError("Value of 'shift' in TensorMechanicsPlasticMohrCoulombMulti must not be negative\n");
69 34 : if (_shift > _f_tol)
70 : _console << "WARNING: value of 'shift' in TensorMechanicsPlasticMohrCoulombMulti is probably "
71 0 : "set too high"
72 0 : << std::endl;
73 : if (LIBMESH_DIM != 3)
74 : mooseError("TensorMechanicsPlasticMohrCoulombMulti is only defined for LIBMESH_DIM=3");
75 : MooseRandom::seed(0);
76 34 : }
77 :
78 : unsigned int
79 2528288 : TensorMechanicsPlasticMohrCoulombMulti::numberSurfaces() const
80 : {
81 2528288 : return 6;
82 : }
83 :
84 : void
85 68408 : TensorMechanicsPlasticMohrCoulombMulti::yieldFunctionV(const RankTwoTensor & stress,
86 : Real intnl,
87 : std::vector<Real> & f) const
88 : {
89 : std::vector<Real> eigvals;
90 68408 : stress.symmetricEigenvalues(eigvals);
91 68408 : eigvals[0] += _shift;
92 68408 : eigvals[2] -= _shift;
93 :
94 68408 : const Real sinphi = std::sin(phi(intnl));
95 68408 : const Real cosphi = std::cos(phi(intnl));
96 68408 : const Real cohcos = cohesion(intnl) * cosphi;
97 :
98 68408 : yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, f);
99 68408 : }
100 :
101 : void
102 94088 : TensorMechanicsPlasticMohrCoulombMulti::yieldFunctionEigvals(
103 : Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector<Real> & f) const
104 : {
105 : // Naively it seems a shame to have 6 yield functions active instead of just
106 : // 3. But 3 won't do. Eg, think of a loading with eigvals[0]=eigvals[1]=eigvals[2]
107 : // Then to return to the yield surface would require 2 positive plastic multipliers
108 : // and one negative one. Boo hoo.
109 :
110 94088 : f.resize(6);
111 94088 : f[0] = 0.5 * (e0 - e1) + 0.5 * (e0 + e1) * sinphi - cohcos;
112 94088 : f[1] = 0.5 * (e1 - e0) + 0.5 * (e0 + e1) * sinphi - cohcos;
113 94088 : f[2] = 0.5 * (e0 - e2) + 0.5 * (e0 + e2) * sinphi - cohcos;
114 94088 : f[3] = 0.5 * (e2 - e0) + 0.5 * (e0 + e2) * sinphi - cohcos;
115 94088 : f[4] = 0.5 * (e1 - e2) + 0.5 * (e1 + e2) * sinphi - cohcos;
116 94088 : f[5] = 0.5 * (e2 - e1) + 0.5 * (e1 + e2) * sinphi - cohcos;
117 94088 : }
118 :
119 : void
120 7632 : TensorMechanicsPlasticMohrCoulombMulti::perturbStress(const RankTwoTensor & stress,
121 : std::vector<Real> & eigvals,
122 : std::vector<RankTwoTensor> & deigvals) const
123 : {
124 : Real small_perturbation;
125 7632 : RankTwoTensor shifted_stress = stress;
126 20827 : while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
127 : {
128 52780 : for (unsigned i = 0; i < 3; ++i)
129 118755 : for (unsigned j = 0; j <= i; ++j)
130 : {
131 79170 : small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
132 79170 : shifted_stress(i, j) += small_perturbation;
133 79170 : shifted_stress(j, i) += small_perturbation;
134 : }
135 13195 : shifted_stress.dsymmetricEigenvalues(eigvals, deigvals);
136 : }
137 7632 : }
138 :
139 : void
140 55568 : TensorMechanicsPlasticMohrCoulombMulti::df_dsig(const RankTwoTensor & stress,
141 : Real sin_angle,
142 : std::vector<RankTwoTensor> & df) const
143 : {
144 : std::vector<Real> eigvals;
145 : std::vector<RankTwoTensor> deigvals;
146 55568 : stress.dsymmetricEigenvalues(eigvals, deigvals);
147 :
148 55568 : if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
149 5088 : perturbStress(stress, eigvals, deigvals);
150 :
151 55568 : df.resize(6);
152 55568 : df[0] = 0.5 * (deigvals[0] - deigvals[1]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
153 55568 : df[1] = 0.5 * (deigvals[1] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
154 55568 : df[2] = 0.5 * (deigvals[0] - deigvals[2]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
155 55568 : df[3] = 0.5 * (deigvals[2] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
156 55568 : df[4] = 0.5 * (deigvals[1] - deigvals[2]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
157 55568 : df[5] = 0.5 * (deigvals[2] - deigvals[1]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
158 55568 : }
159 :
160 : void
161 20168 : TensorMechanicsPlasticMohrCoulombMulti::dyieldFunction_dstressV(
162 : const RankTwoTensor & stress, Real intnl, std::vector<RankTwoTensor> & df_dstress) const
163 : {
164 20168 : const Real sinphi = std::sin(phi(intnl));
165 20168 : df_dsig(stress, sinphi, df_dstress);
166 20168 : }
167 :
168 : void
169 18936 : TensorMechanicsPlasticMohrCoulombMulti::dyieldFunction_dintnlV(const RankTwoTensor & stress,
170 : Real intnl,
171 : std::vector<Real> & df_dintnl) const
172 : {
173 : std::vector<Real> eigvals;
174 18936 : stress.symmetricEigenvalues(eigvals);
175 18936 : eigvals[0] += _shift;
176 18936 : eigvals[2] -= _shift;
177 :
178 18936 : const Real sin_angle = std::sin(phi(intnl));
179 18936 : const Real cos_angle = std::cos(phi(intnl));
180 18936 : const Real dsin_angle = cos_angle * dphi(intnl);
181 18936 : const Real dcos_angle = -sin_angle * dphi(intnl);
182 18936 : const Real dcohcos = dcohesion(intnl) * cos_angle + cohesion(intnl) * dcos_angle;
183 :
184 18936 : df_dintnl.resize(6);
185 18936 : df_dintnl[0] = df_dintnl[1] = 0.5 * (eigvals[0] + eigvals[1]) * dsin_angle - dcohcos;
186 18936 : df_dintnl[2] = df_dintnl[3] = 0.5 * (eigvals[0] + eigvals[2]) * dsin_angle - dcohcos;
187 18936 : df_dintnl[4] = df_dintnl[5] = 0.5 * (eigvals[1] + eigvals[2]) * dsin_angle - dcohcos;
188 18936 : }
189 :
190 : void
191 35400 : TensorMechanicsPlasticMohrCoulombMulti::flowPotentialV(const RankTwoTensor & stress,
192 : Real intnl,
193 : std::vector<RankTwoTensor> & r) const
194 : {
195 35400 : const Real sinpsi = std::sin(psi(intnl));
196 35400 : df_dsig(stress, sinpsi, r);
197 35400 : }
198 :
199 : void
200 18936 : TensorMechanicsPlasticMohrCoulombMulti::dflowPotential_dstressV(
201 : const RankTwoTensor & stress, Real intnl, std::vector<RankFourTensor> & dr_dstress) const
202 : {
203 : std::vector<RankFourTensor> d2eigvals;
204 18936 : stress.d2symmetricEigenvalues(d2eigvals);
205 :
206 18936 : const Real sinpsi = std::sin(psi(intnl));
207 :
208 18936 : dr_dstress.resize(6);
209 : dr_dstress[0] =
210 56808 : 0.5 * (d2eigvals[0] - d2eigvals[1]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
211 : dr_dstress[1] =
212 56808 : 0.5 * (d2eigvals[1] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
213 : dr_dstress[2] =
214 56808 : 0.5 * (d2eigvals[0] - d2eigvals[2]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
215 : dr_dstress[3] =
216 56808 : 0.5 * (d2eigvals[2] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
217 : dr_dstress[4] =
218 56808 : 0.5 * (d2eigvals[1] - d2eigvals[2]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
219 : dr_dstress[5] =
220 56808 : 0.5 * (d2eigvals[2] - d2eigvals[1]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
221 18936 : }
222 :
223 : void
224 18936 : TensorMechanicsPlasticMohrCoulombMulti::dflowPotential_dintnlV(
225 : const RankTwoTensor & stress, Real intnl, std::vector<RankTwoTensor> & dr_dintnl) const
226 : {
227 18936 : const Real cos_angle = std::cos(psi(intnl));
228 18936 : const Real dsin_angle = cos_angle * dpsi(intnl);
229 :
230 : std::vector<Real> eigvals;
231 : std::vector<RankTwoTensor> deigvals;
232 18936 : stress.dsymmetricEigenvalues(eigvals, deigvals);
233 :
234 18936 : if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
235 2544 : perturbStress(stress, eigvals, deigvals);
236 :
237 18936 : dr_dintnl.resize(6);
238 18936 : dr_dintnl[0] = dr_dintnl[1] = 0.5 * (deigvals[0] + deigvals[1]) * dsin_angle;
239 18936 : dr_dintnl[2] = dr_dintnl[3] = 0.5 * (deigvals[0] + deigvals[2]) * dsin_angle;
240 18936 : dr_dintnl[4] = dr_dintnl[5] = 0.5 * (deigvals[1] + deigvals[2]) * dsin_angle;
241 18936 : }
242 :
243 : void
244 6176 : TensorMechanicsPlasticMohrCoulombMulti::activeConstraints(const std::vector<Real> & f,
245 : const RankTwoTensor & stress,
246 : Real intnl,
247 : const RankFourTensor & Eijkl,
248 : std::vector<bool> & act,
249 : RankTwoTensor & returned_stress) const
250 : {
251 : act.assign(6, false);
252 :
253 6176 : if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol && f[3] <= _f_tol && f[4] <= _f_tol &&
254 0 : f[5] <= _f_tol)
255 : {
256 0 : returned_stress = stress;
257 0 : return;
258 : }
259 :
260 : Real returned_intnl;
261 6176 : std::vector<Real> dpm(6);
262 6176 : RankTwoTensor delta_dp;
263 6176 : std::vector<Real> yf(6);
264 : bool trial_stress_inadmissible;
265 6176 : doReturnMap(stress,
266 : intnl,
267 : Eijkl,
268 : 0.0,
269 : returned_stress,
270 : returned_intnl,
271 : dpm,
272 : delta_dp,
273 : yf,
274 : trial_stress_inadmissible);
275 :
276 43232 : for (unsigned i = 0; i < 6; ++i)
277 37056 : act[i] = (dpm[i] > 0);
278 : }
279 :
280 : Real
281 147720 : TensorMechanicsPlasticMohrCoulombMulti::cohesion(const Real internal_param) const
282 : {
283 147720 : return _cohesion.value(internal_param);
284 : }
285 :
286 : Real
287 47272 : TensorMechanicsPlasticMohrCoulombMulti::dcohesion(const Real internal_param) const
288 : {
289 47272 : return _cohesion.derivative(internal_param);
290 : }
291 :
292 : Real
293 315608 : TensorMechanicsPlasticMohrCoulombMulti::phi(const Real internal_param) const
294 : {
295 315608 : return _phi.value(internal_param);
296 : }
297 :
298 : Real
299 66208 : TensorMechanicsPlasticMohrCoulombMulti::dphi(const Real internal_param) const
300 : {
301 66208 : return _phi.derivative(internal_param);
302 : }
303 :
304 : Real
305 83976 : TensorMechanicsPlasticMohrCoulombMulti::psi(const Real internal_param) const
306 : {
307 83976 : return _psi.value(internal_param);
308 : }
309 :
310 : Real
311 18936 : TensorMechanicsPlasticMohrCoulombMulti::dpsi(const Real internal_param) const
312 : {
313 18936 : return _psi.derivative(internal_param);
314 : }
315 :
316 : std::string
317 24 : TensorMechanicsPlasticMohrCoulombMulti::modelName() const
318 : {
319 24 : return "MohrCoulombMulti";
320 : }
321 :
322 : bool
323 13152 : TensorMechanicsPlasticMohrCoulombMulti::returnMap(const RankTwoTensor & trial_stress,
324 : Real intnl_old,
325 : const RankFourTensor & E_ijkl,
326 : Real ep_plastic_tolerance,
327 : RankTwoTensor & returned_stress,
328 : Real & returned_intnl,
329 : std::vector<Real> & dpm,
330 : RankTwoTensor & delta_dp,
331 : std::vector<Real> & yf,
332 : bool & trial_stress_inadmissible) const
333 : {
334 13152 : if (!_use_custom_returnMap)
335 5376 : return TensorMechanicsPlasticModel::returnMap(trial_stress,
336 : intnl_old,
337 : E_ijkl,
338 : ep_plastic_tolerance,
339 : returned_stress,
340 : returned_intnl,
341 : dpm,
342 : delta_dp,
343 : yf,
344 5376 : trial_stress_inadmissible);
345 :
346 7776 : return doReturnMap(trial_stress,
347 : intnl_old,
348 : E_ijkl,
349 : ep_plastic_tolerance,
350 : returned_stress,
351 : returned_intnl,
352 : dpm,
353 : delta_dp,
354 : yf,
355 7776 : trial_stress_inadmissible);
356 : }
357 :
358 : bool
359 13952 : TensorMechanicsPlasticMohrCoulombMulti::doReturnMap(const RankTwoTensor & trial_stress,
360 : Real intnl_old,
361 : const RankFourTensor & E_ijkl,
362 : Real ep_plastic_tolerance,
363 : RankTwoTensor & returned_stress,
364 : Real & returned_intnl,
365 : std::vector<Real> & dpm,
366 : RankTwoTensor & delta_dp,
367 : std::vector<Real> & yf,
368 : bool & trial_stress_inadmissible) const
369 : {
370 : mooseAssert(dpm.size() == 6,
371 : "TensorMechanicsPlasticMohrCoulombMulti size of dpm should be 6 but it is "
372 : << dpm.size());
373 :
374 : std::vector<Real> eigvals;
375 13952 : RankTwoTensor eigvecs;
376 13952 : trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
377 13952 : eigvals[0] += _shift;
378 13952 : eigvals[2] -= _shift;
379 :
380 13952 : Real sinphi = std::sin(phi(intnl_old));
381 13952 : Real cosphi = std::cos(phi(intnl_old));
382 13952 : Real coh = cohesion(intnl_old);
383 13952 : Real cohcos = coh * cosphi;
384 :
385 13952 : yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
386 :
387 13952 : if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol && yf[3] <= _f_tol && yf[4] <= _f_tol &&
388 3248 : yf[5] <= _f_tol)
389 : {
390 : // purely elastic (trial_stress, intnl_old)
391 3248 : trial_stress_inadmissible = false;
392 3248 : return true;
393 : }
394 :
395 10704 : trial_stress_inadmissible = true;
396 : delta_dp.zero();
397 10704 : returned_stress = RankTwoTensor();
398 :
399 : // these are the normals to the 6 yield surfaces, which are const because of the assumption of no
400 : // psi hardening
401 10704 : std::vector<RealVectorValue> norm(6);
402 10704 : const Real sinpsi = std::sin(psi(intnl_old));
403 10704 : const Real oneminus = 0.5 * (1 - sinpsi);
404 10704 : const Real oneplus = 0.5 * (1 + sinpsi);
405 10704 : norm[0](0) = oneplus;
406 10704 : norm[0](1) = -oneminus;
407 10704 : norm[0](2) = 0;
408 10704 : norm[1](0) = -oneminus;
409 10704 : norm[1](1) = oneplus;
410 10704 : norm[1](2) = 0;
411 10704 : norm[2](0) = oneplus;
412 10704 : norm[2](1) = 0;
413 10704 : norm[2](2) = -oneminus;
414 10704 : norm[3](0) = -oneminus;
415 10704 : norm[3](1) = 0;
416 10704 : norm[3](2) = oneplus;
417 10704 : norm[4](0) = 0;
418 10704 : norm[4](1) = oneplus;
419 10704 : norm[4](2) = -oneminus;
420 10704 : norm[5](0) = 0;
421 10704 : norm[5](1) = -oneminus;
422 10704 : norm[5](2) = oneplus;
423 :
424 : // the flow directions are these norm multiplied by Eijkl.
425 : // I call the flow directions "n".
426 : // In the following I assume that the Eijkl is
427 : // for an isotropic situation. Then I don't have to
428 : // rotate to the principal-stress frame, and i don't
429 : // have to worry about strange off-diagonal things
430 10704 : std::vector<RealVectorValue> n(6);
431 74928 : for (unsigned ys = 0; ys < 6; ++ys)
432 256896 : for (unsigned i = 0; i < 3; ++i)
433 770688 : for (unsigned j = 0; j < 3; ++j)
434 578016 : n[ys](i) += E_ijkl(i, i, j, j) * norm[ys](j);
435 10704 : const Real mag_E = E_ijkl(0, 0, 0, 0);
436 :
437 : // With non-zero Poisson's ratio and hardening
438 : // it is not computationally cheap to know whether
439 : // the trial stress will return to the tip, edge,
440 : // or plane. The following at least
441 : // gives a not-completely-stupid guess
442 : // trial_order[0] = type of return to try first
443 : // trial_order[1] = type of return to try second
444 : // trial_order[2] = type of return to try third
445 : // trial_order[3] = type of return to try fourth
446 : // trial_order[4] = type of return to try fifth
447 : // In the following the "binary" stuff indicates the
448 : // deactive (0) and active (1) surfaces, eg
449 : // 110100 means that surfaces 0, 1 and 3 are active
450 : // and 2, 4 and 5 are deactive
451 : const unsigned int number_of_return_paths = 5;
452 10704 : std::vector<int> trial_order(number_of_return_paths);
453 10704 : if (yf[1] > _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
454 : {
455 6048 : trial_order[0] = tip110100;
456 6048 : trial_order[1] = edge010100;
457 6048 : trial_order[2] = plane000100;
458 6048 : trial_order[3] = edge000101;
459 6048 : trial_order[4] = tip010101;
460 : }
461 4656 : else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
462 : {
463 3456 : trial_order[0] = edge000101;
464 3456 : trial_order[1] = plane000100;
465 3456 : trial_order[2] = tip110100;
466 3456 : trial_order[3] = tip010101;
467 3456 : trial_order[4] = edge010100;
468 : }
469 1200 : else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] <= _f_tol)
470 : {
471 800 : trial_order[0] = plane000100;
472 800 : trial_order[1] = edge000101;
473 800 : trial_order[2] = edge010100;
474 800 : trial_order[3] = tip110100;
475 800 : trial_order[4] = tip010101;
476 : }
477 : else
478 : {
479 400 : trial_order[0] = edge010100;
480 400 : trial_order[1] = plane000100;
481 400 : trial_order[2] = edge000101;
482 400 : trial_order[3] = tip110100;
483 400 : trial_order[4] = tip010101;
484 : }
485 :
486 : unsigned trial;
487 10704 : bool nr_converged = false;
488 : bool kt_success = false;
489 10704 : std::vector<RealVectorValue> ntip(3);
490 10704 : std::vector<Real> dpmtip(3);
491 :
492 18088 : for (trial = 0; trial < number_of_return_paths; ++trial)
493 : {
494 18088 : switch (trial_order[trial])
495 : {
496 : case tip110100:
497 24992 : for (unsigned int i = 0; i < 3; ++i)
498 : {
499 18744 : ntip[0](i) = n[0](i);
500 18744 : ntip[1](i) = n[1](i);
501 18744 : ntip[2](i) = n[3](i);
502 : }
503 6248 : kt_success = returnTip(eigvals,
504 : ntip,
505 : dpmtip,
506 : returned_stress,
507 : intnl_old,
508 : sinphi,
509 : cohcos,
510 : 0,
511 : nr_converged,
512 : ep_plastic_tolerance,
513 : yf);
514 6248 : if (nr_converged && kt_success)
515 : {
516 3000 : dpm[0] = dpmtip[0];
517 3000 : dpm[1] = dpmtip[1];
518 3000 : dpm[3] = dpmtip[2];
519 3000 : dpm[2] = dpm[4] = dpm[5] = 0;
520 : }
521 : break;
522 :
523 : case tip010101:
524 448 : for (unsigned int i = 0; i < 3; ++i)
525 : {
526 336 : ntip[0](i) = n[1](i);
527 336 : ntip[1](i) = n[3](i);
528 336 : ntip[2](i) = n[5](i);
529 : }
530 112 : kt_success = returnTip(eigvals,
531 : ntip,
532 : dpmtip,
533 : returned_stress,
534 : intnl_old,
535 : sinphi,
536 : cohcos,
537 : 0,
538 : nr_converged,
539 : ep_plastic_tolerance,
540 : yf);
541 112 : if (nr_converged && kt_success)
542 : {
543 112 : dpm[1] = dpmtip[0];
544 112 : dpm[3] = dpmtip[1];
545 112 : dpm[5] = dpmtip[2];
546 112 : dpm[0] = dpm[2] = dpm[4] = 0;
547 : }
548 : break;
549 :
550 3560 : case edge000101:
551 3560 : kt_success = returnEdge000101(eigvals,
552 : n,
553 : dpm,
554 : returned_stress,
555 : intnl_old,
556 : sinphi,
557 : cohcos,
558 : 0,
559 : mag_E,
560 : nr_converged,
561 : ep_plastic_tolerance,
562 : yf);
563 : break;
564 :
565 3592 : case edge010100:
566 3592 : kt_success = returnEdge010100(eigvals,
567 : n,
568 : dpm,
569 : returned_stress,
570 : intnl_old,
571 : sinphi,
572 : cohcos,
573 : 0,
574 : mag_E,
575 : nr_converged,
576 : ep_plastic_tolerance,
577 : yf);
578 : break;
579 :
580 4576 : case plane000100:
581 4576 : kt_success = returnPlane(eigvals,
582 : n,
583 : dpm,
584 : returned_stress,
585 : intnl_old,
586 : sinphi,
587 : cohcos,
588 : 0,
589 : nr_converged,
590 : ep_plastic_tolerance,
591 : yf);
592 : break;
593 : }
594 :
595 18088 : if (nr_converged && kt_success)
596 : break;
597 : }
598 :
599 10704 : if (trial == number_of_return_paths)
600 : {
601 0 : sinphi = std::sin(phi(intnl_old));
602 0 : cosphi = std::cos(phi(intnl_old));
603 0 : coh = cohesion(intnl_old);
604 0 : cohcos = coh * cosphi;
605 0 : yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
606 : Moose::err << "Trial stress = \n";
607 0 : trial_stress.print(Moose::err);
608 : Moose::err << "which has eigenvalues = " << eigvals[0] << " " << eigvals[1] << " " << eigvals[2]
609 : << "\n";
610 : Moose::err << "and yield functions = " << yf[0] << " " << yf[1] << " " << yf[2] << " " << yf[3]
611 : << " " << yf[4] << " " << yf[5] << "\n";
612 : Moose::err << "Internal parameter = " << intnl_old << std::endl;
613 0 : mooseError("TensorMechanicsPlasticMohrCoulombMulti: FAILURE! You probably need to implement a "
614 : "line search if your hardening is too severe, or you need to tune your tolerances "
615 : "(eg, yield_function_tolerance should be a little smaller than (young "
616 : "modulus)*ep_plastic_tolerance).\n");
617 : return false;
618 : }
619 :
620 : // success
621 :
622 10704 : returned_intnl = intnl_old;
623 74928 : for (unsigned i = 0; i < 6; ++i)
624 64224 : returned_intnl += dpm[i];
625 74928 : for (unsigned i = 0; i < 6; ++i)
626 256896 : for (unsigned j = 0; j < 3; ++j)
627 192672 : delta_dp(j, j) += dpm[i] * norm[i](j);
628 10704 : returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
629 10704 : delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
630 : return true;
631 : }
632 :
633 : bool
634 6360 : TensorMechanicsPlasticMohrCoulombMulti::returnTip(const std::vector<Real> & eigvals,
635 : const std::vector<RealVectorValue> & n,
636 : std::vector<Real> & dpm,
637 : RankTwoTensor & returned_stress,
638 : Real intnl_old,
639 : Real & sinphi,
640 : Real & cohcos,
641 : Real initial_guess,
642 : bool & nr_converged,
643 : Real ep_plastic_tolerance,
644 : std::vector<Real> & yf) const
645 : {
646 : // This returns to the Mohr-Coulomb tip using the THREE directions
647 : // given in n, and yields the THREE dpm values. Note that you
648 : // must supply THREE suitable n vectors out of the total of SIX
649 : // flow directions, and then interpret the THREE dpm values appropriately.
650 : //
651 : // Eg1. You supply the flow directions n[0], n[1] and n[3] as
652 : // the "n" vectors. This is return-to-the-tip via 110100.
653 : // Then the three returned dpm values will be dpm[0], dpm[1] and dpm[3].
654 :
655 : // Eg2. You supply the flow directions n[1], n[3] and n[5] as
656 : // the "n" vectors. This is return-to-the-tip via 010101.
657 : // Then the three returned dpm values will be dpm[1], dpm[3] and dpm[5].
658 :
659 : // The returned point is defined by the three yield functions (corresonding
660 : // to the three supplied flow directions) all being zero.
661 : // that is, returned_stress = diag(cohcot, cohcot, cohcot), where
662 : // cohcot = cohesion*cosphi/sinphi
663 : // where intnl = intnl_old + dpm[0] + dpm[1] + dpm[2]
664 : // The 3 plastic multipliers, dpm, are defiend by the normality condition
665 : // eigvals - cohcot = dpm[0]*n[0] + dpm[1]*n[1] + dpm[2]*n[2]
666 : // (Kuhn-Tucker demands that all dpm are non-negative, but we leave
667 : // that checking for the end.)
668 : // (Remember that these "n" vectors and "dpm" values must be interpreted
669 : // differently depending on what you pass into this function.)
670 : // This is a vector equation with solution (A):
671 : // dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip;
672 : // dpm[1] = triple(eigvals - cohcot, n[2], n[0])/trip;
673 : // dpm[2] = triple(eigvals - cohcot, n[0], n[1])/trip;
674 : // where trip = triple(n[0], n[1], n[2]).
675 : // By adding the three components of that solution together
676 : // we can get an equation for x = dpm[0] + dpm[1] + dpm[2],
677 : // and then our Newton-Raphson only involves one variable (x).
678 : // In the following, i specialise to the isotropic situation.
679 :
680 : mooseAssert(n.size() == 3,
681 : "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
682 : "supplied with n of size 3, whereas yours is "
683 : << n.size());
684 : mooseAssert(dpm.size() == 3,
685 : "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
686 : "supplied with dpm of size 3, whereas yours is "
687 : << dpm.size());
688 : mooseAssert(yf.size() == 6,
689 : "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
690 : "supplied with yf of size 6, whereas yours is "
691 : << yf.size());
692 :
693 : Real x = initial_guess;
694 6360 : const Real trip = triple_product(n[0], n[1], n[2]);
695 6360 : sinphi = std::sin(phi(intnl_old + x));
696 6360 : Real cosphi = std::cos(phi(intnl_old + x));
697 6360 : Real coh = cohesion(intnl_old + x);
698 6360 : cohcos = coh * cosphi;
699 6360 : Real cohcot = cohcos / sinphi;
700 :
701 10184 : if (_cohesion.modelName().compare("Constant") != 0 || _phi.modelName().compare("Constant") != 0)
702 : {
703 : // Finding x is expensive. Therefore
704 : // although x!=0 for Constant Hardening, solution (A)
705 : // demonstrates that we don't
706 : // actually need to know x to find the dpm for
707 : // Constant Hardening.
708 : //
709 : // However, for nontrivial Hardening, the following
710 : // is necessary
711 : // cohcot_coeff = [1,1,1].(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
712 : Real cohcot_coeff =
713 3016 : (n[0](0) * (n[1](1) - n[1](2) - n[2](1)) + (n[1](2) - n[1](1)) * n[2](0) +
714 3016 : (n[1](0) - n[1](2)) * n[2](1) + n[0](2) * (n[1](0) - n[1](1) - n[2](0) + n[2](1)) +
715 3016 : n[0](1) * (n[1](2) - n[1](0) + n[2](0) - n[2](2)) +
716 3016 : (n[0](0) - n[1](0) + n[1](1)) * n[2](2)) /
717 3016 : trip;
718 : // eig_term = eigvals.(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
719 3016 : Real eig_term = eigvals[0] *
720 3016 : (-n[0](2) * n[1](1) + n[0](1) * n[1](2) + n[0](2) * n[2](1) -
721 3016 : n[1](2) * n[2](1) - n[0](1) * n[2](2) + n[1](1) * n[2](2)) /
722 3016 : trip;
723 3016 : eig_term += eigvals[1] *
724 3016 : (n[0](2) * n[1](0) - n[0](0) * n[1](2) - n[0](2) * n[2](0) + n[1](2) * n[2](0) +
725 3016 : n[0](0) * n[2](2) - n[1](0) * n[2](2)) /
726 : trip;
727 3016 : eig_term += eigvals[2] *
728 3016 : (n[0](0) * n[1](1) - n[1](1) * n[2](0) + n[0](1) * n[2](0) - n[0](1) * n[1](0) -
729 3016 : n[0](0) * n[2](1) + n[1](0) * n[2](1)) /
730 : trip;
731 : // and finally, the equation we want to solve is:
732 : // x - eig_term + cohcot*cohcot_coeff = 0
733 : // but i divide by cohcot_coeff so the result has the units of
734 : // stress, so using _f_tol as a convergence check is reasonable
735 3016 : eig_term /= cohcot_coeff;
736 3016 : Real residual = x / cohcot_coeff - eig_term + cohcot;
737 : Real jacobian;
738 : Real deriv_phi;
739 : Real deriv_coh;
740 : unsigned int iter = 0;
741 : do
742 : {
743 6104 : deriv_phi = dphi(intnl_old + x);
744 6104 : deriv_coh = dcohesion(intnl_old + x);
745 6104 : jacobian = 1.0 / cohcot_coeff + deriv_coh * cosphi / sinphi -
746 6104 : coh * deriv_phi / Utility::pow<2>(sinphi);
747 6104 : x += -residual / jacobian;
748 :
749 6104 : if (iter > _max_iters) // not converging
750 : {
751 0 : nr_converged = false;
752 0 : return false;
753 : }
754 :
755 6104 : sinphi = std::sin(phi(intnl_old + x));
756 6104 : cosphi = std::cos(phi(intnl_old + x));
757 6104 : coh = cohesion(intnl_old + x);
758 6104 : cohcos = coh * cosphi;
759 6104 : cohcot = cohcos / sinphi;
760 6104 : residual = x / cohcot_coeff - eig_term + cohcot;
761 6104 : iter++;
762 6104 : } while (residual * residual > _f_tol * _f_tol / 100);
763 : }
764 :
765 : // so the NR process converged, but we must
766 : // calculate the individual dpm values and
767 : // check Kuhn-Tucker
768 6360 : nr_converged = true;
769 6360 : if (x < -3 * ep_plastic_tolerance)
770 : // obviously at least one of the dpm are < -ep_plastic_tolerance. No point in proceeding. This
771 : // is a potential weak-point: if the user has set _f_tol quite large, and ep_plastic_tolerance
772 : // quite small, the above NR process will quickly converge, but the solution may be wrong and
773 : // violate Kuhn-Tucker.
774 : return false;
775 :
776 : // The following is the solution (A) written above
777 : // (dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip, etc)
778 : // in the isotropic situation
779 : RealVectorValue v;
780 5664 : v(0) = eigvals[0] - cohcot;
781 5664 : v(1) = eigvals[1] - cohcot;
782 5664 : v(2) = eigvals[2] - cohcot;
783 5664 : dpm[0] = triple_product(v, n[1], n[2]) / trip;
784 5664 : dpm[1] = triple_product(v, n[2], n[0]) / trip;
785 5664 : dpm[2] = triple_product(v, n[0], n[1]) / trip;
786 :
787 5664 : if (dpm[0] < -ep_plastic_tolerance || dpm[1] < -ep_plastic_tolerance ||
788 : dpm[2] < -ep_plastic_tolerance)
789 : // Kuhn-Tucker failure. No point in proceeding
790 : return false;
791 :
792 : // Kuhn-Tucker has succeeded: just need returned_stress and yf values
793 : // I do not use the dpm to calculate returned_stress, because that
794 : // might add error (and computational effort), simply:
795 3112 : returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = cohcot;
796 : // So by construction the yield functions are all zero
797 3112 : yf[0] = yf[1] = yf[2] = yf[3] = yf[4] = yf[5] = 0;
798 3112 : return true;
799 : }
800 :
801 : bool
802 4576 : TensorMechanicsPlasticMohrCoulombMulti::returnPlane(const std::vector<Real> & eigvals,
803 : const std::vector<RealVectorValue> & n,
804 : std::vector<Real> & dpm,
805 : RankTwoTensor & returned_stress,
806 : Real intnl_old,
807 : Real & sinphi,
808 : Real & cohcos,
809 : Real initial_guess,
810 : bool & nr_converged,
811 : Real ep_plastic_tolerance,
812 : std::vector<Real> & yf) const
813 : {
814 : // This returns to the Mohr-Coulomb plane using n[3] (ie 000100)
815 : //
816 : // The returned point is defined by the f[3]=0 and
817 : // a = eigvals - dpm[3]*n[3]
818 : // where "a" is the returned point and dpm[3] is the plastic multiplier.
819 : // This equation is a vector equation in principal stress space.
820 : // (Kuhn-Tucker also demands that dpm[3]>=0, but we leave checking
821 : // that condition for the end.)
822 : // Since f[3]=0, we must have
823 : // a[2]*(1+sinphi) + a[0]*(-1+sinphi) - 2*coh*cosphi = 0
824 : // which gives dpm[3] as the solution of
825 : // alpha*dpm[3] + eigvals[2] - eigvals[0] + beta*sinphi - 2*coh*cosphi = 0
826 : // with alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0))*sinphi
827 : // beta = eigvals[2] + eigvals[0]
828 :
829 : mooseAssert(n.size() == 6,
830 : "TensorMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
831 : "supplied with n of size 6, whereas yours is "
832 : << n.size());
833 : mooseAssert(dpm.size() == 6,
834 : "TensorMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
835 : "supplied with dpm of size 6, whereas yours is "
836 : << dpm.size());
837 : mooseAssert(yf.size() == 6,
838 : "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
839 : "supplied with yf of size 6, whereas yours is "
840 : << yf.size());
841 :
842 4576 : dpm[3] = initial_guess;
843 4576 : sinphi = std::sin(phi(intnl_old + dpm[3]));
844 4576 : Real cosphi = std::cos(phi(intnl_old + dpm[3]));
845 4576 : Real coh = cohesion(intnl_old + dpm[3]);
846 4576 : cohcos = coh * cosphi;
847 :
848 4576 : Real alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
849 : Real deriv_phi;
850 : Real dalpha;
851 4576 : const Real beta = eigvals[2] + eigvals[0];
852 : Real deriv_coh;
853 :
854 : Real residual =
855 4576 : alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos; // this is 2*yf[3]
856 : Real jacobian;
857 :
858 : const Real f_tol2 = Utility::pow<2>(_f_tol);
859 : unsigned int iter = 0;
860 : do
861 : {
862 7832 : deriv_phi = dphi(intnl_old + dpm[3]);
863 7832 : dalpha = -(n[3](2) + n[3](0)) * cosphi * deriv_phi;
864 7832 : deriv_coh = dcohesion(intnl_old + dpm[3]);
865 7832 : jacobian = alpha + dalpha * dpm[3] + beta * cosphi * deriv_phi - 2.0 * deriv_coh * cosphi +
866 7832 : 2.0 * coh * sinphi * deriv_phi;
867 :
868 7832 : dpm[3] -= residual / jacobian;
869 7832 : if (iter > _max_iters) // not converging
870 : {
871 0 : nr_converged = false;
872 0 : return false;
873 : }
874 :
875 7832 : sinphi = std::sin(phi(intnl_old + dpm[3]));
876 7832 : cosphi = std::cos(phi(intnl_old + dpm[3]));
877 7832 : coh = cohesion(intnl_old + dpm[3]);
878 7832 : cohcos = coh * cosphi;
879 7832 : alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
880 7832 : residual = alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos;
881 7832 : iter++;
882 7832 : } while (residual * residual > f_tol2);
883 :
884 : // so the NR process converged, but we must
885 : // check Kuhn-Tucker
886 4576 : nr_converged = true;
887 4576 : if (dpm[3] < -ep_plastic_tolerance)
888 : // Kuhn-Tucker failure
889 : return false;
890 :
891 18304 : for (unsigned i = 0; i < 3; ++i)
892 13728 : returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i);
893 4576 : yieldFunctionEigvals(
894 : returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
895 :
896 : // by construction abs(yf[3]) = abs(residual/2) < _f_tol/2
897 4576 : if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
898 : // Kuhn-Tucker failure
899 : return false;
900 :
901 : // success!
902 4272 : dpm[0] = dpm[1] = dpm[2] = dpm[4] = dpm[5] = 0;
903 4272 : return true;
904 : }
905 :
906 : bool
907 3560 : TensorMechanicsPlasticMohrCoulombMulti::returnEdge000101(const std::vector<Real> & eigvals,
908 : const std::vector<RealVectorValue> & n,
909 : std::vector<Real> & dpm,
910 : RankTwoTensor & returned_stress,
911 : Real intnl_old,
912 : Real & sinphi,
913 : Real & cohcos,
914 : Real initial_guess,
915 : Real mag_E,
916 : bool & nr_converged,
917 : Real ep_plastic_tolerance,
918 : std::vector<Real> & yf) const
919 : {
920 : // This returns to the Mohr-Coulomb edge
921 : // with 000101 being active. This means that
922 : // f3=f5=0. Denoting the returned stress by "a"
923 : // (in principal stress space), this means that
924 : // a0=a1 = (2Ccosphi + a2(1+sinphi))/(sinphi-1)
925 : //
926 : // Also, a = eigvals - dpm3*n3 - dpm5*n5,
927 : // where the dpm are the plastic multipliers
928 : // and the n are the flow directions.
929 : //
930 : // Hence we have 5 equations and 5 unknowns (a,
931 : // dpm3 and dpm5). Eliminating all unknowns
932 : // except for x = dpm3+dpm5 gives the residual below
933 : // (I used mathematica)
934 :
935 : Real x = initial_guess;
936 3560 : sinphi = std::sin(phi(intnl_old + x));
937 3560 : Real cosphi = std::cos(phi(intnl_old + x));
938 3560 : Real coh = cohesion(intnl_old + x);
939 3560 : cohcos = coh * cosphi;
940 :
941 : const Real numer_const =
942 3560 : -n[3](2) * eigvals[0] - n[5](1) * eigvals[0] + n[5](2) * eigvals[0] + n[3](2) * eigvals[1] +
943 3560 : n[5](0) * eigvals[1] - n[5](2) * eigvals[1] - n[5](0) * eigvals[2] + n[5](1) * eigvals[2] +
944 3560 : n[3](0) * (-eigvals[1] + eigvals[2]) - n[3](1) * (-eigvals[0] + eigvals[2]);
945 3560 : const Real numer_coeff1 = 2 * (-n[3](0) + n[3](1) + n[5](0) - n[5](1));
946 : const Real numer_coeff2 =
947 3560 : n[5](2) * (eigvals[0] - eigvals[1]) + n[3](2) * (-eigvals[0] + eigvals[1]) +
948 3560 : n[5](1) * (eigvals[0] + eigvals[2]) + (n[3](0) - n[5](0)) * (eigvals[1] + eigvals[2]) -
949 3560 : n[3](1) * (eigvals[0] + eigvals[2]);
950 3560 : Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
951 3560 : const Real denom_const = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (-n[5](0) + n[5](2)) +
952 3560 : n[3](0) * (-n[5](1) + n[5](2));
953 3560 : const Real denom_coeff = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (n[5](0) + n[5](2)) +
954 3560 : n[3](0) * (n[5](1) + n[5](2));
955 3560 : Real denom = denom_const + denom_coeff * sinphi;
956 3560 : Real residual = -x + numer / denom;
957 :
958 : Real deriv_phi;
959 : Real deriv_coh;
960 : Real jacobian;
961 3560 : const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
962 : unsigned int iter = 0;
963 : do
964 : {
965 : do
966 : {
967 6240 : deriv_phi = dphi(intnl_old + dpm[3]);
968 6240 : deriv_coh = dcohesion(intnl_old + dpm[3]);
969 6240 : jacobian = -1 +
970 6240 : (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
971 6240 : numer_coeff2 * cosphi * deriv_phi) /
972 : denom -
973 6240 : numer * denom_coeff * cosphi * deriv_phi / denom / denom;
974 6240 : x -= residual / jacobian;
975 6240 : if (iter > _max_iters) // not converging
976 : {
977 0 : nr_converged = false;
978 0 : return false;
979 : }
980 :
981 6240 : sinphi = std::sin(phi(intnl_old + x));
982 6240 : cosphi = std::cos(phi(intnl_old + x));
983 6240 : coh = cohesion(intnl_old + x);
984 6240 : cohcos = coh * cosphi;
985 6240 : numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
986 6240 : denom = denom_const + denom_coeff * sinphi;
987 6240 : residual = -x + numer / denom;
988 6240 : iter++;
989 6240 : } while (residual * residual > tol);
990 :
991 : // now must ensure that yf[3] and yf[5] are both "zero"
992 : const Real dpm3minusdpm5 =
993 3560 : (2.0 * (eigvals[0] - eigvals[1]) + x * (n[3](1) - n[3](0) + n[5](1) - n[5](0))) /
994 3560 : (n[3](0) - n[3](1) + n[5](1) - n[5](0));
995 3560 : dpm[3] = (x + dpm3minusdpm5) / 2.0;
996 3560 : dpm[5] = (x - dpm3minusdpm5) / 2.0;
997 :
998 14240 : for (unsigned i = 0; i < 3; ++i)
999 10680 : returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i) - dpm[5] * n[5](i);
1000 3560 : yieldFunctionEigvals(
1001 : returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
1002 3560 : } while (yf[3] * yf[3] > _f_tol * _f_tol && yf[5] * yf[5] > _f_tol * _f_tol);
1003 :
1004 : // so the NR process converged, but we must
1005 : // check Kuhn-Tucker
1006 3560 : nr_converged = true;
1007 :
1008 3560 : if (dpm[3] < -ep_plastic_tolerance || dpm[5] < -ep_plastic_tolerance)
1009 : // Kuhn-Tucker failure. This is a potential weak-point: if the user has set _f_tol quite
1010 : // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
1011 : // the solution may be wrong and violate Kuhn-Tucker.
1012 : return false;
1013 :
1014 1888 : if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol)
1015 : // Kuhn-Tucker failure
1016 : return false;
1017 :
1018 : // success
1019 1632 : dpm[0] = dpm[1] = dpm[2] = dpm[4] = 0;
1020 1632 : return true;
1021 : }
1022 :
1023 : bool
1024 3592 : TensorMechanicsPlasticMohrCoulombMulti::returnEdge010100(const std::vector<Real> & eigvals,
1025 : const std::vector<RealVectorValue> & n,
1026 : std::vector<Real> & dpm,
1027 : RankTwoTensor & returned_stress,
1028 : Real intnl_old,
1029 : Real & sinphi,
1030 : Real & cohcos,
1031 : Real initial_guess,
1032 : Real mag_E,
1033 : bool & nr_converged,
1034 : Real ep_plastic_tolerance,
1035 : std::vector<Real> & yf) const
1036 : {
1037 : // This returns to the Mohr-Coulomb edge
1038 : // with 010100 being active. This means that
1039 : // f1=f3=0. Denoting the returned stress by "a"
1040 : // (in principal stress space), this means that
1041 : // a1=a2 = (2Ccosphi + a0(1-sinphi))/(sinphi+1)
1042 : //
1043 : // Also, a = eigvals - dpm1*n1 - dpm3*n3,
1044 : // where the dpm are the plastic multipliers
1045 : // and the n are the flow directions.
1046 : //
1047 : // Hence we have 5 equations and 5 unknowns (a,
1048 : // dpm1 and dpm3). Eliminating all unknowns
1049 : // except for x = dpm1+dpm3 gives the residual below
1050 : // (I used mathematica)
1051 :
1052 : Real x = initial_guess;
1053 3592 : sinphi = std::sin(phi(intnl_old + x));
1054 3592 : Real cosphi = std::cos(phi(intnl_old + x));
1055 3592 : Real coh = cohesion(intnl_old + x);
1056 3592 : cohcos = coh * cosphi;
1057 :
1058 3592 : const Real numer_const = -n[1](2) * eigvals[0] - n[3](1) * eigvals[0] + n[3](2) * eigvals[0] -
1059 3592 : n[1](0) * eigvals[1] + n[1](2) * eigvals[1] + n[3](0) * eigvals[1] -
1060 3592 : n[3](2) * eigvals[1] + n[1](0) * eigvals[2] - n[3](0) * eigvals[2] +
1061 3592 : n[3](1) * eigvals[2] - n[1](1) * (-eigvals[0] + eigvals[2]);
1062 3592 : const Real numer_coeff1 = 2 * (n[1](1) - n[1](2) - n[3](1) + n[3](2));
1063 : const Real numer_coeff2 =
1064 3592 : n[3](2) * (-eigvals[0] - eigvals[1]) + n[1](2) * (eigvals[0] + eigvals[1]) +
1065 3592 : n[3](1) * (eigvals[0] + eigvals[2]) + (n[1](0) - n[3](0)) * (eigvals[1] - eigvals[2]) -
1066 3592 : n[1](1) * (eigvals[0] + eigvals[2]);
1067 3592 : Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
1068 3592 : const Real denom_const = -n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (-n[3](0) + n[3](1)) +
1069 3592 : n[1](1) * (-n[3](2) + n[3](0));
1070 : const Real denom_coeff =
1071 3592 : n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (n[3](0) + n[3](1)) - n[1](1) * (n[3](0) + n[3](2));
1072 3592 : Real denom = denom_const + denom_coeff * sinphi;
1073 3592 : Real residual = -x + numer / denom;
1074 :
1075 : Real deriv_phi;
1076 : Real deriv_coh;
1077 : Real jacobian;
1078 3592 : const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
1079 : unsigned int iter = 0;
1080 : do
1081 : {
1082 : do
1083 : {
1084 8160 : deriv_phi = dphi(intnl_old + dpm[3]);
1085 8160 : deriv_coh = dcohesion(intnl_old + dpm[3]);
1086 8160 : jacobian = -1 +
1087 8160 : (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
1088 8160 : numer_coeff2 * cosphi * deriv_phi) /
1089 : denom -
1090 8160 : numer * denom_coeff * cosphi * deriv_phi / denom / denom;
1091 8160 : x -= residual / jacobian;
1092 8160 : if (iter > _max_iters) // not converging
1093 : {
1094 0 : nr_converged = false;
1095 0 : return false;
1096 : }
1097 :
1098 8160 : sinphi = std::sin(phi(intnl_old + x));
1099 8160 : cosphi = std::cos(phi(intnl_old + x));
1100 8160 : coh = cohesion(intnl_old + x);
1101 8160 : cohcos = coh * cosphi;
1102 8160 : numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
1103 8160 : denom = denom_const + denom_coeff * sinphi;
1104 8160 : residual = -x + numer / denom;
1105 8160 : iter++;
1106 8160 : } while (residual * residual > tol);
1107 :
1108 : // now must ensure that yf[1] and yf[3] are both "zero"
1109 : Real dpm1minusdpm3 =
1110 3592 : (2 * (eigvals[1] - eigvals[2]) + x * (n[1](2) - n[1](1) + n[3](2) - n[3](1))) /
1111 3592 : (n[1](1) - n[1](2) + n[3](2) - n[3](1));
1112 3592 : dpm[1] = (x + dpm1minusdpm3) / 2.0;
1113 3592 : dpm[3] = (x - dpm1minusdpm3) / 2.0;
1114 :
1115 14368 : for (unsigned i = 0; i < 3; ++i)
1116 10776 : returned_stress(i, i) = eigvals[i] - dpm[1] * n[1](i) - dpm[3] * n[3](i);
1117 3592 : yieldFunctionEigvals(
1118 : returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
1119 3592 : } while (yf[1] * yf[1] > _f_tol * _f_tol && yf[3] * yf[3] > _f_tol * _f_tol);
1120 :
1121 : // so the NR process converged, but we must
1122 : // check Kuhn-Tucker
1123 3592 : nr_converged = true;
1124 :
1125 3592 : if (dpm[1] < -ep_plastic_tolerance || dpm[3] < -ep_plastic_tolerance)
1126 : // Kuhn-Tucker failure. This is a potential weak-point: if the user has set _f_tol quite
1127 : // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
1128 : // the solution may be wrong and violate Kuhn-Tucker
1129 : return false;
1130 :
1131 1688 : if (yf[0] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
1132 : // Kuhn-Tucker failure
1133 : return false;
1134 :
1135 : // success
1136 1688 : dpm[0] = dpm[2] = dpm[4] = dpm[5] = 0;
1137 1688 : return true;
1138 : }
1139 :
1140 : bool
1141 0 : TensorMechanicsPlasticMohrCoulombMulti::KuhnTuckerOK(const std::vector<Real> & yf,
1142 : const std::vector<Real> & dpm,
1143 : Real ep_plastic_tolerance) const
1144 : {
1145 : mooseAssert(
1146 : yf.size() == 6,
1147 : "TensorMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires yf to be size 6, but your is "
1148 : << yf.size());
1149 : mooseAssert(
1150 : dpm.size() == 6,
1151 : "TensorMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires dpm to be size 6, but your is "
1152 : << dpm.size());
1153 0 : for (unsigned i = 0; i < 6; ++i)
1154 0 : if (!TensorMechanicsPlasticModel::KuhnTuckerSingleSurface(yf[i], dpm[i], ep_plastic_tolerance))
1155 : return false;
1156 : return true;
1157 : }
1158 :
1159 : bool
1160 0 : TensorMechanicsPlasticMohrCoulombMulti::useCustomReturnMap() const
1161 : {
1162 0 : return _use_custom_returnMap;
1163 : }
|
