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 "CappedDruckerPragerStressUpdate.h"
11 :
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", CappedDruckerPragerStressUpdate);
15 :
16 : InputParameters
17 600 : CappedDruckerPragerStressUpdate::validParams()
18 : {
19 600 : InputParameters params = TwoParameterPlasticityStressUpdate::validParams();
20 600 : params.addClassDescription("Capped Drucker-Prager plasticity stress calculator");
21 1200 : params.addRequiredParam<UserObjectName>(
22 : "DP_model",
23 : "A SolidMechanicsPlasticDruckerPrager UserObject that defines the "
24 : "Drucker-Prager parameters (cohesion, friction angle and dilation angle)");
25 1200 : params.addRequiredParam<UserObjectName>(
26 : "tensile_strength",
27 : "A SolidMechanicsHardening UserObject that defines hardening of the "
28 : "tensile strength. In physical situations this is positive (and always "
29 : "must be greater than negative compressive-strength.");
30 1200 : params.addRequiredParam<UserObjectName>(
31 : "compressive_strength",
32 : "A SolidMechanicsHardening UserObject that defines hardening of the "
33 : "compressive strength. In physical situations this is positive.");
34 1200 : params.addRequiredRangeCheckedParam<Real>(
35 : "tip_smoother",
36 : "tip_smoother>=0",
37 : "The cone vertex at J2 = 0 will be smoothed by the given "
38 : "amount. Typical value is 0.1*cohesion");
39 1200 : params.addParam<bool>("perfect_guess",
40 1200 : true,
41 : "Provide a guess to the Newton-Raphson procedure "
42 : "that is the result from perfect plasticity. With "
43 : "severe hardening/softening this may be "
44 : "suboptimal.");
45 1200 : params.addParam<bool>("small_dilation",
46 1200 : true,
47 : "If true, and if the trial stress exceeds the "
48 : "tensile strength, then the user guarantees that "
49 : "the returned stress will be independent of the "
50 : "compressive strength.");
51 600 : return params;
52 0 : }
53 :
54 450 : CappedDruckerPragerStressUpdate::CappedDruckerPragerStressUpdate(const InputParameters & parameters)
55 : : TwoParameterPlasticityStressUpdate(parameters, 3, 2),
56 450 : _dp(getUserObject<SolidMechanicsPlasticDruckerPrager>("DP_model")),
57 450 : _tstrength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
58 450 : _cstrength(getUserObject<SolidMechanicsHardeningModel>("compressive_strength")),
59 900 : _small_smoother2(Utility::pow<2>(getParam<Real>("tip_smoother"))),
60 900 : _perfect_guess(getParam<bool>("perfect_guess")),
61 450 : _stress_return_type(StressReturnType::nothing_special),
62 900 : _small_dilation(getParam<bool>("small_dilation")),
63 450 : _in_q_trial(0.0)
64 : {
65 : // With arbitary UserObjects, it is impossible to check everything,
66 : // but this will catch the common errors
67 450 : if (_tstrength.value(0) + _cstrength.value(0) <= _smoothing_tol)
68 0 : mooseError("CappedDruckerPragerStressUpdate: Tensile strength plus compressive strength must "
69 : "be greater than smoothing_tol");
70 450 : }
71 :
72 : void
73 71680 : CappedDruckerPragerStressUpdate::initializeReturnProcess()
74 : {
75 71680 : _stress_return_type = StressReturnType::nothing_special;
76 71680 : }
77 :
78 : void
79 71328 : CappedDruckerPragerStressUpdate::finalizeReturnProcess(const RankTwoTensor & /*rotation_increment*/)
80 : {
81 71328 : _stress_return_type = StressReturnType::nothing_special;
82 71328 : }
83 :
84 : void
85 71328 : CappedDruckerPragerStressUpdate::preReturnMap(Real /*p_trial*/,
86 : Real q_trial,
87 : const RankTwoTensor & /*stress_trial*/,
88 : const std::vector<Real> & /*intnl_old*/,
89 : const std::vector<Real> & yf,
90 : const RankFourTensor & /*Eijkl*/)
91 : {
92 : // If it's obvious, then simplify the return-type
93 71328 : if (yf[2] >= 0)
94 3504 : _stress_return_type = StressReturnType::no_tension;
95 67824 : else if (_small_dilation && yf[1] >= 0)
96 12960 : _stress_return_type = StressReturnType::no_compression;
97 :
98 71328 : _in_q_trial = q_trial;
99 71328 : }
100 :
101 : void
102 123232 : CappedDruckerPragerStressUpdate::computePQ(const RankTwoTensor & stress, Real & p, Real & q) const
103 : {
104 123232 : p = stress.trace();
105 123232 : q = std::sqrt(stress.secondInvariant());
106 123232 : }
107 :
108 : void
109 64128 : CappedDruckerPragerStressUpdate::setEppEqq(const RankFourTensor & Eijkl,
110 : Real & Epp,
111 : Real & Eqq) const
112 : {
113 64128 : Epp = Eijkl.sum3x3();
114 64128 : Eqq = Eijkl(0, 1, 0, 1);
115 64128 : }
116 :
117 : void
118 200176 : CappedDruckerPragerStressUpdate::setIntnlDerivatives(Real /*p_trial*/,
119 : Real q_trial,
120 : Real /*p*/,
121 : Real q,
122 : const std::vector<Real> & intnl,
123 : std::vector<std::vector<Real>> & dintnl) const
124 : {
125 : Real tanpsi;
126 200176 : _dp.onlyB(intnl[0], _dp.dilation, tanpsi);
127 : Real dtanpsi;
128 200176 : _dp.donlyB(intnl[0], _dp.dilation, dtanpsi);
129 200176 : dintnl[0][0] = 0.0;
130 200176 : dintnl[0][1] = -1.0 / _Eqq;
131 200176 : dintnl[1][0] = -1.0 / _Epp;
132 200176 : dintnl[1][1] = tanpsi / _Eqq - (q_trial - q) * dtanpsi * dintnl[0][1] / _Eqq;
133 200176 : }
134 :
135 : void
136 143008 : CappedDruckerPragerStressUpdate::yieldFunctionValues(Real p,
137 : Real q,
138 : const std::vector<Real> & intnl,
139 : std::vector<Real> & yf) const
140 : {
141 : Real aaa;
142 : Real bbb;
143 143008 : _dp.bothAB(intnl[0], aaa, bbb);
144 143008 : yf[0] = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * bbb - aaa;
145 :
146 143008 : if (_stress_return_type == StressReturnType::no_tension)
147 0 : yf[1] = std::numeric_limits<Real>::lowest();
148 : else
149 143008 : yf[1] = p - _tstrength.value(intnl[1]);
150 :
151 143008 : if (_stress_return_type == StressReturnType::no_compression)
152 0 : yf[2] = std::numeric_limits<Real>::lowest();
153 : else
154 143008 : yf[2] = -p - _cstrength.value(intnl[1]);
155 143008 : }
156 :
157 : void
158 271312 : CappedDruckerPragerStressUpdate::computeAllQ(Real p,
159 : Real q,
160 : const std::vector<Real> & intnl,
161 : std::vector<yieldAndFlow> & all_q) const
162 : {
163 : Real aaa;
164 : Real bbb;
165 271312 : _dp.bothAB(intnl[0], aaa, bbb);
166 : Real daaa;
167 : Real dbbb;
168 271312 : _dp.dbothAB(intnl[0], daaa, dbbb);
169 : Real tanpsi;
170 271312 : _dp.onlyB(intnl[0], _dp.dilation, tanpsi);
171 : Real dtanpsi;
172 271312 : _dp.donlyB(intnl[0], _dp.dilation, dtanpsi);
173 :
174 : // yield function values
175 271312 : all_q[0].f = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * bbb - aaa;
176 271312 : if (_stress_return_type == StressReturnType::no_tension)
177 20864 : all_q[1].f = std::numeric_limits<Real>::lowest();
178 : else
179 250448 : all_q[1].f = p - _tstrength.value(intnl[1]);
180 271312 : if (_stress_return_type == StressReturnType::no_compression)
181 39424 : all_q[2].f = std::numeric_limits<Real>::lowest();
182 : else
183 231888 : all_q[2].f = -p - _cstrength.value(intnl[1]);
184 :
185 : // d(yield Function)/d(p, q)
186 : // derivatives wrt p
187 271312 : all_q[0].df[0] = bbb;
188 271312 : all_q[1].df[0] = 1.0;
189 271312 : all_q[2].df[0] = -1.0;
190 :
191 : // derivatives wrt q
192 271312 : if (_small_smoother2 == 0.0)
193 0 : all_q[0].df[1] = 1.0;
194 : else
195 271312 : all_q[0].df[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
196 271312 : all_q[1].df[1] = 0.0;
197 271312 : all_q[2].df[1] = 0.0;
198 :
199 : // d(yield Function)/d(intnl)
200 : // derivatives wrt intnl[0] (shear plastic strain)
201 271312 : all_q[0].df_di[0] = p * dbbb - daaa;
202 271312 : all_q[1].df_di[0] = 0.0;
203 271312 : all_q[2].df_di[0] = 0.0;
204 : // derivatives wrt intnl[q] (tensile plastic strain)
205 271312 : all_q[0].df_di[1] = 0.0;
206 271312 : all_q[1].df_di[1] = -_tstrength.derivative(intnl[1]);
207 271312 : all_q[2].df_di[1] = -_cstrength.derivative(intnl[1]);
208 :
209 : // d(flowPotential)/d(p, q)
210 : // derivatives wrt p
211 271312 : all_q[0].dg[0] = tanpsi;
212 271312 : all_q[1].dg[0] = 1.0;
213 271312 : all_q[2].dg[0] = -1.0;
214 : // derivatives wrt q
215 271312 : if (_small_smoother2 == 0.0)
216 0 : all_q[0].dg[1] = 1.0;
217 : else
218 271312 : all_q[0].dg[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
219 271312 : all_q[1].dg[1] = 0.0;
220 271312 : all_q[2].dg[1] = 0.0;
221 :
222 : // d2(flowPotential)/d(p, q)/d(intnl)
223 : // d(dg/dp)/dintnl[0]
224 271312 : all_q[0].d2g_di[0][0] = dtanpsi;
225 271312 : all_q[1].d2g_di[0][0] = 0.0;
226 271312 : all_q[2].d2g_di[0][0] = 0.0;
227 : // d(dg/dp)/dintnl[1]
228 271312 : all_q[0].d2g_di[0][1] = 0.0;
229 271312 : all_q[1].d2g_di[0][1] = 0.0;
230 271312 : all_q[2].d2g_di[0][1] = 0.0;
231 : // d(dg/dq)/dintnl[0]
232 271312 : all_q[0].d2g_di[1][0] = 0.0;
233 271312 : all_q[1].d2g_di[1][0] = 0.0;
234 271312 : all_q[2].d2g_di[1][0] = 0.0;
235 : // d(dg/dq)/dintnl[1]
236 271312 : all_q[0].d2g_di[1][1] = 0.0;
237 271312 : all_q[1].d2g_di[1][1] = 0.0;
238 271312 : all_q[2].d2g_di[1][1] = 0.0;
239 :
240 : // d2(flowPotential)/d(p, q)/d(p, q)
241 : // d(dg/dp)/dp
242 271312 : all_q[0].d2g[0][0] = 0.0;
243 271312 : all_q[1].d2g[0][0] = 0.0;
244 271312 : all_q[2].d2g[0][0] = 0.0;
245 : // d(dg/dp)/dq
246 271312 : all_q[0].d2g[0][1] = 0.0;
247 271312 : all_q[1].d2g[0][1] = 0.0;
248 271312 : all_q[2].d2g[0][1] = 0.0;
249 : // d(dg/dq)/dp
250 271312 : all_q[0].d2g[1][0] = 0.0;
251 271312 : all_q[1].d2g[1][0] = 0.0;
252 271312 : all_q[2].d2g[1][0] = 0.0;
253 : // d(dg/dq)/dq
254 271312 : if (_small_smoother2 == 0.0)
255 0 : all_q[0].d2g[1][1] = 0.0;
256 : else
257 271312 : all_q[0].d2g[1][1] = _small_smoother2 / std::pow(Utility::pow<2>(q) + _small_smoother2, 1.5);
258 271312 : all_q[1].d2g[1][1] = 0.0;
259 271312 : all_q[2].d2g[1][1] = 0.0;
260 271312 : }
261 :
262 : void
263 71328 : CappedDruckerPragerStressUpdate::initializeVars(Real p_trial,
264 : Real q_trial,
265 : const std::vector<Real> & intnl_old,
266 : Real & p,
267 : Real & q,
268 : Real & gaE,
269 : std::vector<Real> & intnl) const
270 : {
271 71328 : if (!_perfect_guess)
272 : {
273 0 : p = p_trial;
274 0 : q = q_trial;
275 0 : gaE = 0.0;
276 : }
277 : else
278 : {
279 : Real coh;
280 : Real tanphi;
281 71328 : _dp.bothAB(intnl[0], coh, tanphi);
282 : Real tanpsi;
283 71328 : _dp.onlyB(intnl_old[0], _dp.dilation, tanpsi);
284 71328 : const Real tens = _tstrength.value(intnl_old[1]);
285 71328 : const Real comp = _cstrength.value(intnl_old[1]);
286 71328 : const Real q_at_T = coh - tens * tanphi;
287 71328 : const Real q_at_C = coh + comp * tanphi;
288 :
289 71328 : if ((p_trial >= tens) && (q_trial <= q_at_T))
290 : {
291 : // pure tensile failure
292 6432 : p = tens;
293 6432 : q = q_trial;
294 6432 : gaE = p_trial - p;
295 : }
296 64896 : else if ((p_trial <= -comp) && (q_trial <= q_at_C))
297 : {
298 : // pure compressive failure
299 96 : p = -comp;
300 96 : q = q_trial;
301 96 : gaE = p - p_trial;
302 : }
303 : else
304 : {
305 : // shear failure or a mixture
306 : // Calculate ga assuming a pure shear failure
307 64800 : const Real ga = (q_trial + p_trial * tanphi - coh) / (_Eqq + _Epp * tanphi * tanpsi);
308 64800 : if (ga <= 0 && p_trial <= tens && p_trial >= -comp)
309 : {
310 : // very close to one of the rounded corners: there is no advantage to guessing the
311 : // solution, so:
312 0 : p = p_trial;
313 0 : q = q_trial;
314 0 : gaE = 0.0;
315 : }
316 : else
317 : {
318 64800 : q = q_trial - _Eqq * ga;
319 64800 : if (q <= 0.0 && q_at_T <= 0.0)
320 : {
321 : // user has set tensile strength so large that it is irrelevant: return to the tip of the
322 : // shear cone
323 0 : q = 0.0;
324 0 : p = coh / tanphi;
325 0 : gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
326 : }
327 64800 : else if (q <= q_at_T)
328 : {
329 : // pure shear is incorrect: mixture of tension and shear is correct
330 5600 : q = q_at_T;
331 5600 : p = tens;
332 5600 : gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
333 : }
334 59200 : else if (q >= q_at_C)
335 : {
336 : // pure shear is incorrect: mixture of compression and shear is correct
337 3872 : q = q_at_C;
338 3872 : p = -comp;
339 3872 : if (p - p_trial < _Epp * tanpsi * (q_trial - q) / _Eqq)
340 : // trial point is sitting almost directly above corner
341 640 : gaE = (q_trial - q) * _Epp / _Eqq;
342 : else
343 : // trial point is sitting to the left of the corner
344 3232 : gaE = (p - p_trial) / tanpsi;
345 : }
346 : else
347 : {
348 : // pure shear was correct
349 55328 : p = p_trial - _Epp * ga * tanpsi;
350 55328 : gaE = ga * _Epp;
351 : }
352 : }
353 : }
354 : }
355 71328 : setIntnlValues(p_trial, q_trial, p, q, intnl_old, intnl);
356 71328 : }
357 :
358 : void
359 342640 : CappedDruckerPragerStressUpdate::setIntnlValues(Real p_trial,
360 : Real q_trial,
361 : Real p,
362 : Real q,
363 : const std::vector<Real> & intnl_old,
364 : std::vector<Real> & intnl) const
365 : {
366 342640 : intnl[0] = intnl_old[0] + (q_trial - q) / _Eqq;
367 : Real tanpsi;
368 342640 : _dp.onlyB(intnl[0], _dp.dilation, tanpsi);
369 342640 : intnl[1] = intnl_old[1] + (p_trial - p) / _Epp - (q_trial - q) * tanpsi / _Eqq;
370 342640 : }
371 :
372 : void
373 64128 : CappedDruckerPragerStressUpdate::setStressAfterReturn(const RankTwoTensor & stress_trial,
374 : Real p_ok,
375 : Real q_ok,
376 : Real /*gaE*/,
377 : const std::vector<Real> & /*intnl*/,
378 : const yieldAndFlow & /*smoothed_q*/,
379 : const RankFourTensor & /*Eijkl*/,
380 : RankTwoTensor & stress) const
381 : {
382 : // stress = s_ij + de_ij tr(stress) / 3 = q / q_trial * s_ij^trial + de_ij p / 3 = q / q_trial *
383 : // (stress_ij^trial - de_ij tr(stress^trial) / 3) + de_ij p / 3
384 64128 : const Real p_trial = stress_trial.trace();
385 64128 : stress = RankTwoTensor(RankTwoTensor::initIdentity) / 3.0 *
386 64128 : (p_ok - (_in_q_trial == 0.0 ? 0.0 : p_trial * q_ok / _in_q_trial));
387 64128 : if (_in_q_trial > 0)
388 64128 : stress += q_ok / _in_q_trial * stress_trial;
389 64128 : }
390 :
391 : RankTwoTensor
392 71328 : CappedDruckerPragerStressUpdate::dpdstress(const RankTwoTensor & stress) const
393 : {
394 71328 : return stress.dtrace();
395 : }
396 :
397 : RankFourTensor
398 0 : CappedDruckerPragerStressUpdate::d2pdstress2(const RankTwoTensor & /*stress*/) const
399 : {
400 0 : return RankFourTensor();
401 : }
402 :
403 : RankTwoTensor
404 71328 : CappedDruckerPragerStressUpdate::dqdstress(const RankTwoTensor & stress) const
405 : {
406 71328 : const Real j2 = stress.secondInvariant();
407 71328 : if (j2 == 0.0)
408 0 : return RankTwoTensor();
409 71328 : return 0.5 * stress.dsecondInvariant() / std::sqrt(j2);
410 : }
411 :
412 : RankFourTensor
413 192 : CappedDruckerPragerStressUpdate::d2qdstress2(const RankTwoTensor & stress) const
414 : {
415 192 : const Real j2 = stress.secondInvariant();
416 192 : if (j2 == 0.0)
417 0 : return RankFourTensor();
418 :
419 192 : const RankTwoTensor dj2 = stress.dsecondInvariant();
420 192 : return -0.25 * dj2.outerProduct(dj2) / std::pow(j2, 1.5) +
421 384 : 0.5 * stress.d2secondInvariant() / std::sqrt(j2);
422 : }
423 :
424 : void
425 128 : CappedDruckerPragerStressUpdate::consistentTangentOperator(const RankTwoTensor & /*stress_trial*/,
426 : Real /*p_trial*/,
427 : Real /*q_trial*/,
428 : const RankTwoTensor & stress,
429 : Real /*p*/,
430 : Real q,
431 : Real gaE,
432 : const yieldAndFlow & smoothed_q,
433 : const RankFourTensor & Eijkl,
434 : bool compute_full_tangent_operator,
435 : RankFourTensor & cto) const
436 : {
437 128 : cto = Eijkl;
438 128 : if (!compute_full_tangent_operator)
439 0 : return;
440 :
441 : const RankTwoTensor s_over_q =
442 128 : (q == 0.0 ? RankTwoTensor()
443 512 : : (stress - stress.trace() * RankTwoTensor(RankTwoTensor::initIdentity) / 3.0) / q);
444 :
445 512 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
446 1536 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
447 4608 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
448 13824 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
449 : {
450 10368 : cto(i, j, k, l) -=
451 20736 : (i == j) * (1.0 / 3.0) *
452 17280 : (_Epp * (1.0 - _dp_dpt) / 3.0 * (k == l) + s_over_q(k, l) * _Eqq * (-_dp_dqt));
453 10368 : cto(i, j, k, l) -= s_over_q(i, j) * (_Epp * (-_dq_dpt) / 3.0 * (k == l) +
454 10368 : s_over_q(k, l) * _Eqq * (1.0 - _dq_dqt));
455 : }
456 :
457 128 : if (smoothed_q.dg[1] != 0.0)
458 : {
459 192 : const RankFourTensor Tijab = Eijkl * (gaE / _Epp) * smoothed_q.dg[1] * d2qdstress2(stress);
460 96 : RankFourTensor inv = RankFourTensor(RankFourTensor::initIdentityFour) + Tijab;
461 : try
462 : {
463 96 : inv = inv.transposeMajor().invSymm();
464 : }
465 0 : catch (const MooseException & e)
466 : {
467 : // Cannot form the inverse, so probably at some degenerate place in stress space.
468 : // Just return with the "best estimate" of the cto.
469 0 : mooseWarning("CappedDruckerPragerStressUpdate: Cannot invert 1+T in consistent tangent "
470 : "operator computation at quadpoint ",
471 0 : _qp,
472 : " of element ",
473 0 : _current_elem->id());
474 : return;
475 0 : }
476 96 : cto = (cto.transposeMajor() * inv).transposeMajor();
477 : }
478 : }
|
