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