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 "TensorMechanicsPlasticMeanCapTC.h"
11 : #include "RankFourTensor.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("TensorMechanicsApp", TensorMechanicsPlasticMeanCapTC);
15 :
16 : InputParameters
17 78 : TensorMechanicsPlasticMeanCapTC::validParams()
18 : {
19 78 : InputParameters params = TensorMechanicsPlasticModel::validParams();
20 234 : params.addRangeCheckedParam<unsigned>("max_iterations",
21 156 : 10,
22 : "max_iterations>0",
23 : "Maximum iterations for custom MeanCapTC return map");
24 156 : params.addParam<bool>(
25 156 : "use_custom_returnMap", true, "Whether to use the custom MeanCapTC returnMap algorithm.");
26 156 : params.addParam<bool>("use_custom_cto",
27 156 : true,
28 : "Whether to use the custom consistent tangent operator computations.");
29 156 : params.addRequiredParam<UserObjectName>("tensile_strength",
30 : "A TensorMechanicsHardening UserObject that defines "
31 : "hardening of the mean-cap tensile strength (it will "
32 : "typically be positive). Yield function = trace(stress) "
33 : "- tensile_strength for trace(stress)>tensile_strength.");
34 156 : params.addRequiredParam<UserObjectName>(
35 : "compressive_strength",
36 : "A TensorMechanicsHardening UserObject that defines hardening of the "
37 : "mean-cap compressive strength. This should always be less than "
38 : "tensile_strength (it will typically be negative). Yield function = "
39 : "- (trace(stress) - compressive_strength) for "
40 : "trace(stress)<compressive_strength.");
41 78 : params.addClassDescription(
42 : "Associative mean-cap tensile and compressive plasticity with hardening/softening");
43 :
44 78 : return params;
45 0 : }
46 :
47 39 : TensorMechanicsPlasticMeanCapTC::TensorMechanicsPlasticMeanCapTC(const InputParameters & parameters)
48 : : TensorMechanicsPlasticModel(parameters),
49 39 : _max_iters(getParam<unsigned>("max_iterations")),
50 78 : _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
51 78 : _use_custom_cto(getParam<bool>("use_custom_cto")),
52 39 : _strength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
53 78 : _c_strength(getUserObject<TensorMechanicsHardeningModel>("compressive_strength"))
54 : {
55 : // cannot check the following for all values of the internal parameter, but this will catch most
56 : // errors
57 39 : if (_strength.value(0) <= _c_strength.value(0))
58 0 : mooseError("MeanCapTC: tensile strength (which is usually positive) must not be less than "
59 : "compressive strength (which is usually negative)");
60 39 : }
61 :
62 : Real
63 36400 : TensorMechanicsPlasticMeanCapTC::yieldFunction(const RankTwoTensor & stress, Real intnl) const
64 : {
65 36400 : const Real tr = stress.trace();
66 36400 : const Real t_str = tensile_strength(intnl);
67 36400 : if (tr >= t_str)
68 11984 : return tr - t_str;
69 :
70 24416 : const Real c_str = compressive_strength(intnl);
71 24416 : if (tr <= c_str)
72 16708 : return -(tr - c_str);
73 : // the following is zero at tr = t_str, and at tr = c_str
74 : // it also has derivative = 1 at tr = t_str, and derivative = -1 at tr = c_str
75 : // it also has second derivative = 0, at these points.
76 : // This makes the complete yield function C2 continuous.
77 7708 : return (c_str - t_str) / libMesh::pi * std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str));
78 : }
79 :
80 : RankTwoTensor
81 12096 : TensorMechanicsPlasticMeanCapTC::dyieldFunction_dstress(const RankTwoTensor & stress,
82 : Real intnl) const
83 : {
84 12096 : return df_dsig(stress, intnl);
85 : }
86 :
87 : Real
88 12096 : TensorMechanicsPlasticMeanCapTC::dyieldFunction_dintnl(const RankTwoTensor & stress,
89 : Real intnl) const
90 : {
91 12096 : const Real tr = stress.trace();
92 12096 : const Real t_str = tensile_strength(intnl);
93 12096 : if (tr >= t_str)
94 4272 : return -dtensile_strength(intnl);
95 :
96 7824 : const Real c_str = compressive_strength(intnl);
97 7824 : if (tr <= c_str)
98 5828 : return dcompressive_strength(intnl);
99 :
100 1996 : const Real dt = dtensile_strength(intnl);
101 1996 : const Real dc = dcompressive_strength(intnl);
102 1996 : return (dc - dt) / libMesh::pi * std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) +
103 1996 : 1.0 / (t_str - c_str) * std::cos(libMesh::pi * (tr - c_str) / (t_str - c_str)) *
104 1996 : ((tr - c_str) * dt - (tr - t_str) * dc);
105 : }
106 :
107 : RankTwoTensor
108 37344 : TensorMechanicsPlasticMeanCapTC::df_dsig(const RankTwoTensor & stress, Real intnl) const
109 : {
110 37344 : const Real tr = stress.trace();
111 37344 : const Real t_str = tensile_strength(intnl);
112 37344 : if (tr >= t_str)
113 13824 : return stress.dtrace();
114 :
115 23520 : const Real c_str = compressive_strength(intnl);
116 23520 : if (tr <= c_str)
117 17484 : return -stress.dtrace();
118 :
119 6036 : return -std::cos(libMesh::pi * (tr - c_str) / (t_str - c_str)) * stress.dtrace();
120 : }
121 :
122 : RankTwoTensor
123 25248 : TensorMechanicsPlasticMeanCapTC::flowPotential(const RankTwoTensor & stress, Real intnl) const
124 : {
125 25248 : return df_dsig(stress, intnl);
126 : }
127 :
128 : RankFourTensor
129 12096 : TensorMechanicsPlasticMeanCapTC::dflowPotential_dstress(const RankTwoTensor & stress,
130 : Real intnl) const
131 : {
132 12096 : const Real tr = stress.trace();
133 12096 : const Real t_str = tensile_strength(intnl);
134 12096 : if (tr >= t_str)
135 4272 : return RankFourTensor();
136 :
137 7824 : const Real c_str = compressive_strength(intnl);
138 7824 : if (tr <= c_str)
139 5828 : return RankFourTensor();
140 :
141 1996 : return libMesh::pi / (t_str - c_str) * std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) *
142 3992 : stress.dtrace().outerProduct(stress.dtrace());
143 : }
144 :
145 : RankTwoTensor
146 12096 : TensorMechanicsPlasticMeanCapTC::dflowPotential_dintnl(const RankTwoTensor & stress,
147 : Real intnl) const
148 : {
149 12096 : const Real tr = stress.trace();
150 12096 : const Real t_str = tensile_strength(intnl);
151 12096 : if (tr >= t_str)
152 4272 : return RankTwoTensor();
153 :
154 7824 : const Real c_str = compressive_strength(intnl);
155 7824 : if (tr <= c_str)
156 5828 : return RankTwoTensor();
157 :
158 1996 : const Real dt = dtensile_strength(intnl);
159 1996 : const Real dc = dcompressive_strength(intnl);
160 1996 : return std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) * stress.dtrace() * libMesh::pi /
161 1996 : Utility::pow<2>(t_str - c_str) * ((tr - t_str) * dc - (tr - c_str) * dt);
162 : }
163 :
164 : Real
165 25248 : TensorMechanicsPlasticMeanCapTC::hardPotential(const RankTwoTensor & stress, Real intnl) const
166 : {
167 : // This is the key for this whole class!
168 25248 : const Real tr = stress.trace();
169 25248 : const Real t_str = tensile_strength(intnl);
170 :
171 25248 : if (tr >= t_str)
172 : return -1.0; // this will serve to *increase* the internal parameter (so internal parameter will
173 : // be a measure of volumetric strain)
174 :
175 15696 : const Real c_str = compressive_strength(intnl);
176 15696 : if (tr <= c_str)
177 : return 1.0; // this will serve to *decrease* the internal parameter (so internal parameter will
178 : // be a measure of volumetric strain)
179 :
180 4040 : return std::cos(libMesh::pi * (tr - c_str) /
181 4040 : (t_str - c_str)); // this interpolates C2 smoothly between 1 and -1
182 : }
183 :
184 : RankTwoTensor
185 6576 : TensorMechanicsPlasticMeanCapTC::dhardPotential_dstress(const RankTwoTensor & stress,
186 : Real intnl) const
187 : {
188 6576 : const Real tr = stress.trace();
189 6576 : const Real t_str = tensile_strength(intnl);
190 6576 : if (tr >= t_str)
191 3072 : return RankTwoTensor();
192 :
193 3504 : const Real c_str = compressive_strength(intnl);
194 3504 : if (tr <= c_str)
195 3456 : return RankTwoTensor();
196 :
197 48 : return -std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) * libMesh::pi / (t_str - c_str) *
198 48 : stress.dtrace();
199 : }
200 :
201 : Real
202 6576 : TensorMechanicsPlasticMeanCapTC::dhardPotential_dintnl(const RankTwoTensor & stress,
203 : Real intnl) const
204 : {
205 6576 : const Real tr = stress.trace();
206 6576 : const Real t_str = tensile_strength(intnl);
207 6576 : if (tr >= t_str)
208 : return 0.0;
209 :
210 3504 : const Real c_str = compressive_strength(intnl);
211 3504 : if (tr <= c_str)
212 : return 0.0;
213 :
214 48 : const Real dt = dtensile_strength(intnl);
215 48 : const Real dc = dcompressive_strength(intnl);
216 48 : return -std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) * libMesh::pi /
217 48 : Utility::pow<2>(t_str - c_str) * ((tr - t_str) * dc - (tr - c_str) * dt);
218 : }
219 :
220 : Real
221 174160 : TensorMechanicsPlasticMeanCapTC::tensile_strength(const Real internal_param) const
222 : {
223 174160 : return _strength.value(internal_param);
224 : }
225 :
226 : Real
227 18008 : TensorMechanicsPlasticMeanCapTC::dtensile_strength(const Real internal_param) const
228 : {
229 18008 : return _strength.derivative(internal_param);
230 : }
231 :
232 : Real
233 105504 : TensorMechanicsPlasticMeanCapTC::compressive_strength(const Real internal_param) const
234 : {
235 105504 : return _c_strength.value(internal_param);
236 : }
237 :
238 : Real
239 21772 : TensorMechanicsPlasticMeanCapTC::dcompressive_strength(const Real internal_param) const
240 : {
241 21772 : return _c_strength.derivative(internal_param);
242 : }
243 :
244 : void
245 5520 : TensorMechanicsPlasticMeanCapTC::activeConstraints(const std::vector<Real> & f,
246 : const RankTwoTensor & stress,
247 : Real intnl,
248 : const RankFourTensor & Eijkl,
249 : std::vector<bool> & act,
250 : RankTwoTensor & returned_stress) const
251 : {
252 : act.assign(1, false);
253 :
254 5520 : if (f[0] <= _f_tol)
255 : {
256 0 : returned_stress = stress;
257 0 : return;
258 : }
259 :
260 5520 : const Real tr = stress.trace();
261 5520 : const Real t_str = tensile_strength(intnl);
262 : Real str;
263 : Real dirn;
264 5520 : if (tr >= t_str)
265 : {
266 : str = t_str;
267 : dirn = 1.0;
268 : }
269 : else
270 : {
271 3456 : str = compressive_strength(intnl);
272 : dirn = -1.0;
273 : }
274 :
275 5520 : RankTwoTensor n; // flow direction
276 22080 : for (unsigned i = 0; i < 3; ++i)
277 66240 : for (unsigned j = 0; j < 3; ++j)
278 198720 : for (unsigned k = 0; k < 3; ++k)
279 149040 : n(i, j) += dirn * Eijkl(i, j, k, k);
280 :
281 : // returned_stress = stress - gamma*n
282 : // and taking the trace of this and using
283 : // Tr(returned_stress) = str, gives
284 : // gamma = (Tr(stress) - str)/Tr(n)
285 5520 : Real gamma = (stress.trace() - str) / n.trace();
286 :
287 22080 : for (unsigned i = 0; i < 3; ++i)
288 66240 : for (unsigned j = 0; j < 3; ++j)
289 49680 : returned_stress(i, j) = stress(i, j) - gamma * n(i, j);
290 :
291 : act[0] = true;
292 : }
293 :
294 : bool
295 23248 : TensorMechanicsPlasticMeanCapTC::returnMap(const RankTwoTensor & trial_stress,
296 : Real intnl_old,
297 : const RankFourTensor & E_ijkl,
298 : Real ep_plastic_tolerance,
299 : RankTwoTensor & returned_stress,
300 : Real & returned_intnl,
301 : std::vector<Real> & dpm,
302 : RankTwoTensor & delta_dp,
303 : std::vector<Real> & yf,
304 : bool & trial_stress_inadmissible) const
305 : {
306 23248 : if (!(_use_custom_returnMap))
307 13824 : return TensorMechanicsPlasticModel::returnMap(trial_stress,
308 : intnl_old,
309 : E_ijkl,
310 : ep_plastic_tolerance,
311 : returned_stress,
312 : returned_intnl,
313 : dpm,
314 : delta_dp,
315 : yf,
316 : trial_stress_inadmissible);
317 :
318 9424 : yf.resize(1);
319 :
320 9424 : Real yf_orig = yieldFunction(trial_stress, intnl_old);
321 :
322 9424 : yf[0] = yf_orig;
323 :
324 9424 : if (yf_orig < _f_tol)
325 : {
326 : // the trial_stress is admissible
327 2880 : trial_stress_inadmissible = false;
328 2880 : return true;
329 : }
330 :
331 6544 : trial_stress_inadmissible = true;
332 :
333 : // In the following we want to solve
334 : // trial_stress - stress = E_ijkl * dpm * r ...... (1)
335 : // and either
336 : // stress.trace() = tensile_strength(intnl) ...... (2a)
337 : // intnl = intnl_old + dpm ...... (3a)
338 : // or
339 : // stress.trace() = compressive_strength(intnl) ... (2b)
340 : // intnl = intnl_old - dpm ...... (3b)
341 : // The former (2a and 3a) are chosen if
342 : // trial_stress.trace() > tensile_strength(intnl_old)
343 : // while the latter (2b and 3b) are chosen if
344 : // trial_stress.trace() < compressive_strength(intnl_old)
345 : // The variables we want to solve for are stress, dpm
346 : // and intnl. We do this using a Newton approach, starting
347 : // with stress=trial_stress and intnl=intnl_old and dpm=0
348 6544 : const bool tensile_failure = (trial_stress.trace() >= tensile_strength(intnl_old));
349 6544 : const Real dirn = (tensile_failure ? 1.0 : -1.0);
350 :
351 6544 : RankTwoTensor n; // flow direction, which is E_ijkl * r
352 26176 : for (unsigned i = 0; i < 3; ++i)
353 78528 : for (unsigned j = 0; j < 3; ++j)
354 235584 : for (unsigned k = 0; k < 3; ++k)
355 176688 : n(i, j) += dirn * E_ijkl(i, j, k, k);
356 6544 : const Real n_trace = n.trace();
357 :
358 : // Perform a Newton-Raphson to find dpm when
359 : // residual = trial_stress.trace() - tensile_strength(intnl) - dpm * n.trace() [for
360 : // tensile_failure=true]
361 : // or
362 : // residual = trial_stress.trace() - compressive_strength(intnl) - dpm * n.trace() [for
363 : // tensile_failure=false]
364 6544 : Real trial_trace = trial_stress.trace();
365 : Real residual;
366 : Real jac;
367 6544 : dpm[0] = 0;
368 : unsigned int iter = 0;
369 : do
370 : {
371 15056 : if (tensile_failure)
372 : {
373 7120 : residual = trial_trace - tensile_strength(intnl_old + dpm[0]) - dpm[0] * n_trace;
374 7120 : jac = -dtensile_strength(intnl_old + dpm[0]) - n_trace;
375 : }
376 : else
377 : {
378 7936 : residual = trial_trace - compressive_strength(intnl_old - dpm[0]) - dpm[0] * n_trace;
379 7936 : jac = -dcompressive_strength(intnl_old - dpm[0]) - n_trace;
380 : }
381 15056 : dpm[0] += -residual / jac;
382 15056 : if (iter > _max_iters) // not converging
383 : return false;
384 15056 : iter++;
385 15056 : } while (residual * residual > _f_tol * _f_tol);
386 :
387 : // set the returned values
388 6544 : yf[0] = 0;
389 6544 : returned_intnl = intnl_old + dirn * dpm[0];
390 6544 : returned_stress = trial_stress - dpm[0] * n;
391 6544 : delta_dp = dpm[0] * dirn * returned_stress.dtrace();
392 :
393 6544 : return true;
394 : }
395 :
396 : RankFourTensor
397 6544 : TensorMechanicsPlasticMeanCapTC::consistentTangentOperator(
398 : const RankTwoTensor & trial_stress,
399 : Real intnl_old,
400 : const RankTwoTensor & stress,
401 : Real intnl,
402 : const RankFourTensor & E_ijkl,
403 : const std::vector<Real> & cumulative_pm) const
404 : {
405 6544 : if (!_use_custom_cto)
406 0 : return TensorMechanicsPlasticModel::consistentTangentOperator(
407 : trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
408 :
409 : Real df_dq;
410 : Real alpha;
411 6544 : if (trial_stress.trace() >= tensile_strength(intnl_old))
412 : {
413 2576 : df_dq = -dtensile_strength(intnl);
414 : alpha = 1.0;
415 : }
416 : else
417 : {
418 3968 : df_dq = dcompressive_strength(intnl);
419 : alpha = -1.0;
420 : }
421 :
422 6544 : RankTwoTensor elas;
423 26176 : for (unsigned int i = 0; i < 3; ++i)
424 78528 : for (unsigned int j = 0; j < 3; ++j)
425 235584 : for (unsigned int k = 0; k < 3; ++k)
426 176688 : elas(i, j) += E_ijkl(i, j, k, k);
427 :
428 6544 : const Real hw = -df_dq + alpha * elas.trace();
429 :
430 13088 : return E_ijkl - alpha / hw * elas.outerProduct(elas);
431 : }
432 :
433 : bool
434 0 : TensorMechanicsPlasticMeanCapTC::useCustomReturnMap() const
435 : {
436 0 : return _use_custom_returnMap;
437 : }
438 :
439 : bool
440 6544 : TensorMechanicsPlasticMeanCapTC::useCustomCTO() const
441 : {
442 6544 : return _use_custom_cto;
443 : }
444 :
445 : std::string
446 0 : TensorMechanicsPlasticMeanCapTC::modelName() const
447 : {
448 0 : return "MeanCapTC";
449 : }
|
