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 156 : SolidMechanicsPlasticMeanCapTC::validParams()
22 : {
23 156 : InputParameters params = SolidMechanicsPlasticModel::validParams();
24 468 : params.addRangeCheckedParam<unsigned>("max_iterations",
25 312 : 10,
26 : "max_iterations>0",
27 : "Maximum iterations for custom MeanCapTC return map");
28 312 : params.addParam<bool>(
29 312 : "use_custom_returnMap", true, "Whether to use the custom MeanCapTC returnMap algorithm.");
30 312 : params.addParam<bool>("use_custom_cto",
31 312 : true,
32 : "Whether to use the custom consistent tangent operator computations.");
33 312 : 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 312 : 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 156 : params.addClassDescription(
46 : "Associative mean-cap tensile and compressive plasticity with hardening/softening");
47 :
48 156 : return params;
49 0 : }
50 :
51 78 : SolidMechanicsPlasticMeanCapTC::SolidMechanicsPlasticMeanCapTC(const InputParameters & parameters)
52 : : SolidMechanicsPlasticModel(parameters),
53 78 : _max_iters(getParam<unsigned>("max_iterations")),
54 156 : _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
55 156 : _use_custom_cto(getParam<bool>("use_custom_cto")),
56 78 : _strength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
57 156 : _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 78 : 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 78 : }
65 :
66 : Real
67 61472 : SolidMechanicsPlasticMeanCapTC::yieldFunction(const RankTwoTensor & stress, Real intnl) const
68 : {
69 61472 : const Real tr = stress.trace();
70 61472 : const Real t_str = tensile_strength(intnl);
71 61472 : if (tr >= t_str)
72 23936 : return tr - t_str;
73 :
74 37536 : const Real c_str = compressive_strength(intnl);
75 37536 : if (tr <= c_str)
76 33384 : 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 4152 : return (c_str - t_str) / libMesh::pi * std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str));
82 : }
83 :
84 : RankTwoTensor
85 24192 : SolidMechanicsPlasticMeanCapTC::dyieldFunction_dstress(const RankTwoTensor & stress,
86 : Real intnl) const
87 : {
88 24192 : return df_dsig(stress, intnl);
89 : }
90 :
91 : Real
92 24192 : SolidMechanicsPlasticMeanCapTC::dyieldFunction_dintnl(const RankTwoTensor & stress,
93 : Real intnl) const
94 : {
95 24192 : const Real tr = stress.trace();
96 24192 : const Real t_str = tensile_strength(intnl);
97 24192 : if (tr >= t_str)
98 8544 : return -dtensile_strength(intnl);
99 :
100 15648 : const Real c_str = compressive_strength(intnl);
101 15648 : if (tr <= c_str)
102 11656 : return dcompressive_strength(intnl);
103 :
104 3992 : const Real dt = dtensile_strength(intnl);
105 3992 : const Real dc = dcompressive_strength(intnl);
106 3992 : return (dc - dt) / libMesh::pi * std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) +
107 3992 : 1.0 / (t_str - c_str) * std::cos(libMesh::pi * (tr - c_str) / (t_str - c_str)) *
108 3992 : ((tr - c_str) * dt - (tr - t_str) * dc);
109 : }
110 :
111 : RankTwoTensor
112 74688 : SolidMechanicsPlasticMeanCapTC::df_dsig(const RankTwoTensor & stress, Real intnl) const
113 : {
114 74688 : const Real tr = stress.trace();
115 74688 : const Real t_str = tensile_strength(intnl);
116 74688 : if (tr >= t_str)
117 27648 : return stress.dtrace();
118 :
119 47040 : const Real c_str = compressive_strength(intnl);
120 47040 : if (tr <= c_str)
121 34968 : return -stress.dtrace();
122 :
123 12072 : return -std::cos(libMesh::pi * (tr - c_str) / (t_str - c_str)) * stress.dtrace();
124 : }
125 :
126 : RankTwoTensor
127 50496 : SolidMechanicsPlasticMeanCapTC::flowPotential(const RankTwoTensor & stress, Real intnl) const
128 : {
129 50496 : return df_dsig(stress, intnl);
130 : }
131 :
132 : RankFourTensor
133 24192 : SolidMechanicsPlasticMeanCapTC::dflowPotential_dstress(const RankTwoTensor & stress,
134 : Real intnl) const
135 : {
136 24192 : const Real tr = stress.trace();
137 24192 : const Real t_str = tensile_strength(intnl);
138 24192 : if (tr >= t_str)
139 8544 : return RankFourTensor();
140 :
141 15648 : const Real c_str = compressive_strength(intnl);
142 15648 : if (tr <= c_str)
143 11656 : return RankFourTensor();
144 :
145 3992 : return libMesh::pi / (t_str - c_str) * std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) *
146 7984 : stress.dtrace().outerProduct(stress.dtrace());
147 : }
148 :
149 : RankTwoTensor
150 24192 : SolidMechanicsPlasticMeanCapTC::dflowPotential_dintnl(const RankTwoTensor & stress,
151 : Real intnl) const
152 : {
153 24192 : const Real tr = stress.trace();
154 24192 : const Real t_str = tensile_strength(intnl);
155 24192 : if (tr >= t_str)
156 8544 : return RankTwoTensor();
157 :
158 15648 : const Real c_str = compressive_strength(intnl);
159 15648 : if (tr <= c_str)
160 11656 : return RankTwoTensor();
161 :
162 3992 : const Real dt = dtensile_strength(intnl);
163 3992 : const Real dc = dcompressive_strength(intnl);
164 3992 : return std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) * stress.dtrace() * libMesh::pi /
165 3992 : Utility::pow<2>(t_str - c_str) * ((tr - t_str) * dc - (tr - c_str) * dt);
166 : }
167 :
168 : Real
169 50496 : SolidMechanicsPlasticMeanCapTC::hardPotential(const RankTwoTensor & stress, Real intnl) const
170 : {
171 : // This is the key for this whole class!
172 50496 : const Real tr = stress.trace();
173 50496 : const Real t_str = tensile_strength(intnl);
174 :
175 50496 : 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 31392 : const Real c_str = compressive_strength(intnl);
180 31392 : 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 8080 : return std::cos(libMesh::pi * (tr - c_str) /
185 8080 : (t_str - c_str)); // this interpolates C2 smoothly between 1 and -1
186 : }
187 :
188 : RankTwoTensor
189 13152 : SolidMechanicsPlasticMeanCapTC::dhardPotential_dstress(const RankTwoTensor & stress,
190 : Real intnl) const
191 : {
192 13152 : const Real tr = stress.trace();
193 13152 : const Real t_str = tensile_strength(intnl);
194 13152 : if (tr >= t_str)
195 6144 : return RankTwoTensor();
196 :
197 7008 : const Real c_str = compressive_strength(intnl);
198 7008 : if (tr <= c_str)
199 6912 : return RankTwoTensor();
200 :
201 96 : return -std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) * libMesh::pi / (t_str - c_str) *
202 96 : stress.dtrace();
203 : }
204 :
205 : Real
206 13152 : SolidMechanicsPlasticMeanCapTC::dhardPotential_dintnl(const RankTwoTensor & stress,
207 : Real intnl) const
208 : {
209 13152 : const Real tr = stress.trace();
210 13152 : const Real t_str = tensile_strength(intnl);
211 13152 : if (tr >= t_str)
212 : return 0.0;
213 :
214 7008 : const Real c_str = compressive_strength(intnl);
215 7008 : if (tr <= c_str)
216 : return 0.0;
217 :
218 96 : const Real dt = dtensile_strength(intnl);
219 96 : const Real dc = dcompressive_strength(intnl);
220 96 : return -std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) * libMesh::pi /
221 96 : Utility::pow<2>(t_str - c_str) * ((tr - t_str) * dc - (tr - c_str) * dt);
222 : }
223 :
224 : Real
225 336736 : SolidMechanicsPlasticMeanCapTC::tensile_strength(const Real internal_param) const
226 : {
227 336736 : return _strength.value(internal_param);
228 : }
229 :
230 : Real
231 35856 : SolidMechanicsPlasticMeanCapTC::dtensile_strength(const Real internal_param) const
232 : {
233 35856 : return _strength.derivative(internal_param);
234 : }
235 :
236 : Real
237 199648 : SolidMechanicsPlasticMeanCapTC::compressive_strength(const Real internal_param) const
238 : {
239 199648 : return _c_strength.value(internal_param);
240 : }
241 :
242 : Real
243 43448 : SolidMechanicsPlasticMeanCapTC::dcompressive_strength(const Real internal_param) const
244 : {
245 43448 : return _c_strength.derivative(internal_param);
246 : }
247 :
248 : void
249 11040 : 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 11040 : if (f[0] <= _f_tol)
259 : {
260 0 : returned_stress = stress;
261 0 : return;
262 : }
263 :
264 11040 : const Real tr = stress.trace();
265 11040 : const Real t_str = tensile_strength(intnl);
266 : Real str;
267 : Real dirn;
268 11040 : if (tr >= t_str)
269 : {
270 : str = t_str;
271 : dirn = 1.0;
272 : }
273 : else
274 : {
275 6912 : str = compressive_strength(intnl);
276 : dirn = -1.0;
277 : }
278 :
279 11040 : RankTwoTensor n; // flow direction
280 44160 : for (unsigned i = 0; i < 3; ++i)
281 132480 : for (unsigned j = 0; j < 3; ++j)
282 397440 : for (unsigned k = 0; k < 3; ++k)
283 298080 : 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 11040 : Real gamma = (stress.trace() - str) / n.trace();
290 :
291 44160 : for (unsigned i = 0; i < 3; ++i)
292 132480 : for (unsigned j = 0; j < 3; ++j)
293 99360 : returned_stress(i, j) = stress(i, j) - gamma * n(i, j);
294 :
295 : act[0] = true;
296 : }
297 :
298 : bool
299 35168 : 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 35168 : if (!(_use_custom_returnMap))
311 22080 : 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 13088 : yf.resize(1);
323 :
324 13088 : Real yf_orig = yieldFunction(trial_stress, intnl_old);
325 :
326 13088 : yf[0] = yf_orig;
327 :
328 13088 : 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 13024 : 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 13024 : const bool tensile_failure = (trial_stress.trace() >= tensile_strength(intnl_old));
353 13024 : const Real dirn = (tensile_failure ? 1.0 : -1.0);
354 :
355 13024 : RankTwoTensor n; // flow direction, which is E_ijkl * r
356 52096 : for (unsigned i = 0; i < 3; ++i)
357 156288 : for (unsigned j = 0; j < 3; ++j)
358 468864 : for (unsigned k = 0; k < 3; ++k)
359 351648 : n(i, j) += dirn * E_ijkl(i, j, k, k);
360 13024 : 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 13024 : Real trial_trace = trial_stress.trace();
369 : Real residual;
370 : Real jac;
371 13024 : dpm[0] = 0;
372 : unsigned int iter = 0;
373 : do
374 : {
375 29920 : if (tensile_failure)
376 : {
377 14112 : residual = trial_trace - tensile_strength(intnl_old + dpm[0]) - dpm[0] * n_trace;
378 14112 : jac = -dtensile_strength(intnl_old + dpm[0]) - n_trace;
379 : }
380 : else
381 : {
382 15808 : residual = trial_trace - compressive_strength(intnl_old - dpm[0]) - dpm[0] * n_trace;
383 15808 : jac = -dcompressive_strength(intnl_old - dpm[0]) - n_trace;
384 : }
385 29920 : dpm[0] += -residual / jac;
386 29920 : if (iter > _max_iters) // not converging
387 : return false;
388 29920 : iter++;
389 29920 : } while (residual * residual > _f_tol * _f_tol);
390 :
391 : // set the returned values
392 13024 : yf[0] = 0;
393 13024 : returned_intnl = intnl_old + dirn * dpm[0];
394 13024 : returned_stress = trial_stress - dpm[0] * n;
395 13024 : delta_dp = dpm[0] * dirn * returned_stress.dtrace();
396 :
397 13024 : return true;
398 : }
399 :
400 : RankFourTensor
401 13024 : 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 13024 : 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 13024 : if (trial_stress.trace() >= tensile_strength(intnl_old))
416 : {
417 5120 : df_dq = -dtensile_strength(intnl);
418 : alpha = 1.0;
419 : }
420 : else
421 : {
422 7904 : df_dq = dcompressive_strength(intnl);
423 : alpha = -1.0;
424 : }
425 :
426 13024 : RankTwoTensor elas;
427 52096 : for (unsigned int i = 0; i < 3; ++i)
428 156288 : for (unsigned int j = 0; j < 3; ++j)
429 468864 : for (unsigned int k = 0; k < 3; ++k)
430 351648 : elas(i, j) += E_ijkl(i, j, k, k);
431 :
432 13024 : const Real hw = -df_dq + alpha * elas.trace();
433 :
434 26048 : 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 13024 : SolidMechanicsPlasticMeanCapTC::useCustomCTO() const
445 : {
446 13024 : return _use_custom_cto;
447 : }
448 :
449 : std::string
450 0 : SolidMechanicsPlasticMeanCapTC::modelName() const
451 : {
452 0 : return "MeanCapTC";
453 : }
|
