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 "CappedDruckerPragerStressUpdate.h"
11 :
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("TensorMechanicsApp", CappedDruckerPragerStressUpdate);
15 :
16 : InputParameters
17 300 : CappedDruckerPragerStressUpdate::validParams()
18 : {
19 300 : InputParameters params = TwoParameterPlasticityStressUpdate::validParams();
20 300 : params.addClassDescription("Capped Drucker-Prager plasticity stress calculator");
21 600 : params.addRequiredParam<UserObjectName>(
22 : "DP_model",
23 : "A TensorMechanicsPlasticDruckerPrager UserObject that defines the "
24 : "Drucker-Prager parameters (cohesion, friction angle and dilation angle)");
25 600 : params.addRequiredParam<UserObjectName>(
26 : "tensile_strength",
27 : "A TensorMechanicsHardening 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 600 : params.addRequiredParam<UserObjectName>(
31 : "compressive_strength",
32 : "A TensorMechanicsHardening UserObject that defines hardening of the "
33 : "compressive strength. In physical situations this is positive.");
34 600 : 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 600 : params.addParam<bool>("perfect_guess",
40 600 : 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 600 : params.addParam<bool>("small_dilation",
46 600 : 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 300 : return params;
52 0 : }
53 :
54 225 : CappedDruckerPragerStressUpdate::CappedDruckerPragerStressUpdate(const InputParameters & parameters)
55 : : TwoParameterPlasticityStressUpdate(parameters, 3, 2),
56 225 : _dp(getUserObject<TensorMechanicsPlasticDruckerPrager>("DP_model")),
57 225 : _tstrength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
58 225 : _cstrength(getUserObject<TensorMechanicsHardeningModel>("compressive_strength")),
59 450 : _small_smoother2(Utility::pow<2>(getParam<Real>("tip_smoother"))),
60 450 : _perfect_guess(getParam<bool>("perfect_guess")),
61 225 : _stress_return_type(StressReturnType::nothing_special),
62 450 : _small_dilation(getParam<bool>("small_dilation")),
63 225 : _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 225 : 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 225 : }
71 :
72 : void
73 55792 : CappedDruckerPragerStressUpdate::initializeReturnProcess()
74 : {
75 55792 : _stress_return_type = StressReturnType::nothing_special;
76 55792 : }
77 :
78 : void
79 36736 : CappedDruckerPragerStressUpdate::finalizeReturnProcess(const RankTwoTensor & /*rotation_increment*/)
80 : {
81 36736 : _stress_return_type = StressReturnType::nothing_special;
82 36736 : }
83 :
84 : void
85 36736 : 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 36736 : if (yf[2] >= 0)
94 1752 : _stress_return_type = StressReturnType::no_tension;
95 34984 : else if (_small_dilation && yf[1] >= 0)
96 7488 : _stress_return_type = StressReturnType::no_compression;
97 :
98 36736 : _in_q_trial = q_trial;
99 36736 : }
100 :
101 : void
102 81568 : CappedDruckerPragerStressUpdate::computePQ(const RankTwoTensor & stress, Real & p, Real & q) const
103 : {
104 81568 : p = stress.trace();
105 81568 : q = std::sqrt(stress.secondInvariant());
106 81568 : }
107 :
108 : void
109 33056 : CappedDruckerPragerStressUpdate::setEppEqq(const RankFourTensor & Eijkl,
110 : Real & Epp,
111 : Real & Eqq) const
112 : {
113 33056 : Epp = Eijkl.sum3x3();
114 33056 : Eqq = Eijkl(0, 1, 0, 1);
115 33056 : }
116 :
117 : void
118 100552 : 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 100552 : _dp.onlyB(intnl[0], _dp.dilation, tanpsi);
127 : Real dtanpsi;
128 100552 : _dp.donlyB(intnl[0], _dp.dilation, dtanpsi);
129 100552 : dintnl[0][0] = 0.0;
130 100552 : dintnl[0][1] = -1.0 / _Eqq;
131 100552 : dintnl[1][0] = -1.0 / _Epp;
132 100552 : dintnl[1][1] = tanpsi / _Eqq - (q_trial - q) * dtanpsi * dintnl[0][1] / _Eqq;
133 100552 : }
134 :
135 : void
136 92528 : 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 92528 : _dp.bothAB(intnl[0], aaa, bbb);
144 92528 : yf[0] = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * bbb - aaa;
145 :
146 92528 : if (_stress_return_type == StressReturnType::no_tension)
147 0 : yf[1] = std::numeric_limits<Real>::lowest();
148 : else
149 92528 : yf[1] = p - _tstrength.value(intnl[1]);
150 :
151 92528 : if (_stress_return_type == StressReturnType::no_compression)
152 0 : yf[2] = std::numeric_limits<Real>::lowest();
153 : else
154 92528 : yf[2] = -p - _cstrength.value(intnl[1]);
155 92528 : }
156 :
157 : void
158 137192 : 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 137192 : _dp.bothAB(intnl[0], aaa, bbb);
166 : Real daaa;
167 : Real dbbb;
168 137192 : _dp.dbothAB(intnl[0], daaa, dbbb);
169 : Real tanpsi;
170 137192 : _dp.onlyB(intnl[0], _dp.dilation, tanpsi);
171 : Real dtanpsi;
172 137192 : _dp.donlyB(intnl[0], _dp.dilation, dtanpsi);
173 :
174 : // yield function values
175 137192 : all_q[0].f = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * bbb - aaa;
176 137192 : if (_stress_return_type == StressReturnType::no_tension)
177 10432 : all_q[1].f = std::numeric_limits<Real>::lowest();
178 : else
179 126760 : all_q[1].f = p - _tstrength.value(intnl[1]);
180 137192 : if (_stress_return_type == StressReturnType::no_compression)
181 21120 : all_q[2].f = std::numeric_limits<Real>::lowest();
182 : else
183 116072 : all_q[2].f = -p - _cstrength.value(intnl[1]);
184 :
185 : // d(yield Function)/d(p, q)
186 : // derivatives wrt p
187 137192 : all_q[0].df[0] = bbb;
188 137192 : all_q[1].df[0] = 1.0;
189 137192 : all_q[2].df[0] = -1.0;
190 :
191 : // derivatives wrt q
192 137192 : if (_small_smoother2 == 0.0)
193 0 : all_q[0].df[1] = 1.0;
194 : else
195 137192 : all_q[0].df[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
196 137192 : all_q[1].df[1] = 0.0;
197 137192 : all_q[2].df[1] = 0.0;
198 :
199 : // d(yield Function)/d(intnl)
200 : // derivatives wrt intnl[0] (shear plastic strain)
201 137192 : all_q[0].df_di[0] = p * dbbb - daaa;
202 137192 : all_q[1].df_di[0] = 0.0;
203 137192 : all_q[2].df_di[0] = 0.0;
204 : // derivatives wrt intnl[q] (tensile plastic strain)
205 137192 : all_q[0].df_di[1] = 0.0;
206 137192 : all_q[1].df_di[1] = -_tstrength.derivative(intnl[1]);
207 137192 : all_q[2].df_di[1] = -_cstrength.derivative(intnl[1]);
208 :
209 : // d(flowPotential)/d(p, q)
210 : // derivatives wrt p
211 137192 : all_q[0].dg[0] = tanpsi;
212 137192 : all_q[1].dg[0] = 1.0;
213 137192 : all_q[2].dg[0] = -1.0;
214 : // derivatives wrt q
215 137192 : if (_small_smoother2 == 0.0)
216 0 : all_q[0].dg[1] = 1.0;
217 : else
218 137192 : all_q[0].dg[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
219 137192 : all_q[1].dg[1] = 0.0;
220 137192 : all_q[2].dg[1] = 0.0;
221 :
222 : // d2(flowPotential)/d(p, q)/d(intnl)
223 : // d(dg/dp)/dintnl[0]
224 137192 : all_q[0].d2g_di[0][0] = dtanpsi;
225 137192 : all_q[1].d2g_di[0][0] = 0.0;
226 137192 : all_q[2].d2g_di[0][0] = 0.0;
227 : // d(dg/dp)/dintnl[1]
228 137192 : all_q[0].d2g_di[0][1] = 0.0;
229 137192 : all_q[1].d2g_di[0][1] = 0.0;
230 137192 : all_q[2].d2g_di[0][1] = 0.0;
231 : // d(dg/dq)/dintnl[0]
232 137192 : all_q[0].d2g_di[1][0] = 0.0;
233 137192 : all_q[1].d2g_di[1][0] = 0.0;
234 137192 : all_q[2].d2g_di[1][0] = 0.0;
235 : // d(dg/dq)/dintnl[1]
236 137192 : all_q[0].d2g_di[1][1] = 0.0;
237 137192 : all_q[1].d2g_di[1][1] = 0.0;
238 137192 : 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 137192 : all_q[0].d2g[0][0] = 0.0;
243 137192 : all_q[1].d2g[0][0] = 0.0;
244 137192 : all_q[2].d2g[0][0] = 0.0;
245 : // d(dg/dp)/dq
246 137192 : all_q[0].d2g[0][1] = 0.0;
247 137192 : all_q[1].d2g[0][1] = 0.0;
248 137192 : all_q[2].d2g[0][1] = 0.0;
249 : // d(dg/dq)/dp
250 137192 : all_q[0].d2g[1][0] = 0.0;
251 137192 : all_q[1].d2g[1][0] = 0.0;
252 137192 : all_q[2].d2g[1][0] = 0.0;
253 : // d(dg/dq)/dq
254 137192 : if (_small_smoother2 == 0.0)
255 0 : all_q[0].d2g[1][1] = 0.0;
256 : else
257 137192 : all_q[0].d2g[1][1] = _small_smoother2 / std::pow(Utility::pow<2>(q) + _small_smoother2, 1.5);
258 137192 : all_q[1].d2g[1][1] = 0.0;
259 137192 : all_q[2].d2g[1][1] = 0.0;
260 137192 : }
261 :
262 : void
263 36736 : 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 36736 : 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 36736 : _dp.bothAB(intnl[0], coh, tanphi);
282 : Real tanpsi;
283 36736 : _dp.onlyB(intnl_old[0], _dp.dilation, tanpsi);
284 36736 : const Real tens = _tstrength.value(intnl_old[1]);
285 36736 : const Real comp = _cstrength.value(intnl_old[1]);
286 36736 : const Real q_at_T = coh - tens * tanphi;
287 36736 : const Real q_at_C = coh + comp * tanphi;
288 :
289 36736 : if ((p_trial >= tens) && (q_trial <= q_at_T))
290 : {
291 : // pure tensile failure
292 4128 : p = tens;
293 4128 : q = q_trial;
294 4128 : gaE = p_trial - p;
295 : }
296 32608 : else if ((p_trial <= -comp) && (q_trial <= q_at_C))
297 : {
298 : // pure compressive failure
299 48 : p = -comp;
300 48 : q = q_trial;
301 48 : gaE = p - p_trial;
302 : }
303 : else
304 : {
305 : // shear failure or a mixture
306 : // Calculate ga assuming a pure shear failure
307 32560 : const Real ga = (q_trial + p_trial * tanphi - coh) / (_Eqq + _Epp * tanphi * tanpsi);
308 32560 : 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 64 : p = p_trial;
313 64 : q = q_trial;
314 64 : gaE = 0.0;
315 : }
316 : else
317 : {
318 32496 : q = q_trial - _Eqq * ga;
319 32496 : 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 32496 : else if (q <= q_at_T)
328 : {
329 : // pure shear is incorrect: mixture of tension and shear is correct
330 2864 : q = q_at_T;
331 2864 : p = tens;
332 2864 : gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
333 : }
334 29632 : else if (q >= q_at_C)
335 : {
336 : // pure shear is incorrect: mixture of compression and shear is correct
337 1936 : q = q_at_C;
338 1936 : p = -comp;
339 1936 : if (p - p_trial < _Epp * tanpsi * (q_trial - q) / _Eqq)
340 : // trial point is sitting almost directly above corner
341 320 : gaE = (q_trial - q) * _Epp / _Eqq;
342 : else
343 : // trial point is sitting to the left of the corner
344 1616 : gaE = (p - p_trial) / tanpsi;
345 : }
346 : else
347 : {
348 : // pure shear was correct
349 27696 : p = p_trial - _Epp * ga * tanpsi;
350 27696 : gaE = ga * _Epp;
351 : }
352 : }
353 : }
354 : }
355 36736 : setIntnlValues(p_trial, q_trial, p, q, intnl_old, intnl);
356 36736 : }
357 :
358 : void
359 173928 : 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 173928 : intnl[0] = intnl_old[0] + (q_trial - q) / _Eqq;
367 : Real tanpsi;
368 173928 : _dp.onlyB(intnl[0], _dp.dilation, tanpsi);
369 173928 : intnl[1] = intnl_old[1] + (p_trial - p) / _Epp - (q_trial - q) * tanpsi / _Eqq;
370 173928 : }
371 :
372 : void
373 33056 : 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 33056 : const Real p_trial = stress_trial.trace();
385 99168 : stress = RankTwoTensor(RankTwoTensor::initIdentity) / 3.0 *
386 33056 : (p_ok - (_in_q_trial == 0.0 ? 0.0 : p_trial * q_ok / _in_q_trial));
387 33056 : if (_in_q_trial > 0)
388 32992 : stress += q_ok / _in_q_trial * stress_trial;
389 33056 : }
390 :
391 : RankTwoTensor
392 36736 : CappedDruckerPragerStressUpdate::dpdstress(const RankTwoTensor & stress) const
393 : {
394 36736 : 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 36736 : CappedDruckerPragerStressUpdate::dqdstress(const RankTwoTensor & stress) const
405 : {
406 36736 : const Real j2 = stress.secondInvariant();
407 36736 : if (j2 == 0.0)
408 64 : return RankTwoTensor();
409 36672 : return 0.5 * stress.dsecondInvariant() / std::sqrt(j2);
410 : }
411 :
412 : RankFourTensor
413 96 : CappedDruckerPragerStressUpdate::d2qdstress2(const RankTwoTensor & stress) const
414 : {
415 96 : const Real j2 = stress.secondInvariant();
416 96 : if (j2 == 0.0)
417 0 : return RankFourTensor();
418 :
419 96 : const RankTwoTensor dj2 = stress.dsecondInvariant();
420 96 : return -0.25 * dj2.outerProduct(dj2) / std::pow(j2, 1.5) +
421 192 : 0.5 * stress.d2secondInvariant() / std::sqrt(j2);
422 : }
423 :
424 : void
425 64 : 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 64 : cto = Eijkl;
438 64 : if (!compute_full_tangent_operator)
439 0 : return;
440 :
441 : const RankTwoTensor s_over_q =
442 64 : (q == 0.0 ? RankTwoTensor()
443 256 : : (stress - stress.trace() * RankTwoTensor(RankTwoTensor::initIdentity) / 3.0) / q);
444 :
445 256 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
446 768 : for (unsigned j = 0; j < _tensor_dimensionality; ++j)
447 2304 : for (unsigned k = 0; k < _tensor_dimensionality; ++k)
448 6912 : for (unsigned l = 0; l < _tensor_dimensionality; ++l)
449 : {
450 5184 : cto(i, j, k, l) -=
451 10368 : (i == j) * (1.0 / 3.0) *
452 8640 : (_Epp * (1.0 - _dp_dpt) / 3.0 * (k == l) + s_over_q(k, l) * _Eqq * (-_dp_dqt));
453 5184 : cto(i, j, k, l) -= s_over_q(i, j) * (_Epp * (-_dq_dpt) / 3.0 * (k == l) +
454 5184 : s_over_q(k, l) * _Eqq * (1.0 - _dq_dqt));
455 : }
456 :
457 64 : if (smoothed_q.dg[1] != 0.0)
458 : {
459 96 : const RankFourTensor Tijab = Eijkl * (gaE / _Epp) * smoothed_q.dg[1] * d2qdstress2(stress);
460 48 : RankFourTensor inv = RankFourTensor(RankFourTensor::initIdentityFour) + Tijab;
461 : try
462 : {
463 48 : 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 48 : cto = (cto.transposeMajor() * inv).transposeMajor();
477 : }
478 : }
|
