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