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 "SolidMechanicsPlasticTensile.h"
11 : #include "RankFourTensor.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", SolidMechanicsPlasticTensile);
15 : registerMooseObjectRenamed("SolidMechanicsApp",
16 : TensorMechanicsPlasticTensile,
17 : "01/01/2025 00:00",
18 : SolidMechanicsPlasticTensile);
19 :
20 : InputParameters
21 286 : SolidMechanicsPlasticTensile::validParams()
22 : {
23 286 : InputParameters params = SolidMechanicsPlasticModel::validParams();
24 572 : params.addRequiredParam<UserObjectName>(
25 : "tensile_strength",
26 : "A SolidMechanicsHardening UserObject that defines hardening of the tensile strength");
27 858 : params.addRangeCheckedParam<Real>(
28 : "tensile_edge_smoother",
29 572 : 25.0,
30 : "tensile_edge_smoother>=0 & tensile_edge_smoother<=30",
31 : "Smoothing parameter: the edges of the cone are smoothed by the given amount.");
32 572 : MooseEnum tip_scheme("hyperbolic cap", "hyperbolic");
33 572 : params.addParam<MooseEnum>(
34 : "tip_scheme", tip_scheme, "Scheme by which the pyramid's tip will be smoothed.");
35 572 : params.addRequiredRangeCheckedParam<Real>("tensile_tip_smoother",
36 : "tensile_tip_smoother>=0",
37 : "For the 'hyperbolic' tip_scheme, the pyramid vertex "
38 : "will be smoothed by the given amount. For the 'cap' "
39 : "tip_scheme, additional smoothing will occur. Typical "
40 : "value is 0.1*tensile_strength");
41 572 : params.addParam<Real>(
42 : "cap_start",
43 572 : 0.0,
44 : "For the 'cap' tip_scheme, smoothing is performed in the stress_mean > cap_start region");
45 858 : params.addRangeCheckedParam<Real>("cap_rate",
46 572 : 0.0,
47 : "cap_rate>=0",
48 : "For the 'cap' tip_scheme, this controls how quickly the cap "
49 : "degenerates to a hemisphere: small values mean a slow "
50 : "degeneration to a hemisphere (and zero means the 'cap' will "
51 : "be totally inactive). Typical value is 1/tensile_strength");
52 572 : params.addParam<Real>("tensile_lode_cutoff",
53 : "If the second invariant of stress is less than "
54 : "this amount, the Lode angle is assumed to be zero. "
55 : "This is to guard against precision-loss problems, "
56 : "and this parameter should be set small. Default = "
57 : "0.00001*((yield_Function_tolerance)^2)");
58 286 : params.addClassDescription(
59 : "Associative tensile plasticity with hardening/softening, and tensile_strength = 1");
60 :
61 286 : return params;
62 286 : }
63 :
64 143 : SolidMechanicsPlasticTensile::SolidMechanicsPlasticTensile(const InputParameters & parameters)
65 : : SolidMechanicsPlasticModel(parameters),
66 143 : _strength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
67 286 : _tip_scheme(getParam<MooseEnum>("tip_scheme")),
68 286 : _small_smoother2(Utility::pow<2>(getParam<Real>("tensile_tip_smoother"))),
69 286 : _cap_start(getParam<Real>("cap_start")),
70 286 : _cap_rate(getParam<Real>("cap_rate")),
71 286 : _tt(getParam<Real>("tensile_edge_smoother") * libMesh::pi / 180.0),
72 143 : _sin3tt(std::sin(3.0 * _tt)),
73 143 : _lode_cutoff(parameters.isParamValid("tensile_lode_cutoff")
74 143 : ? getParam<Real>("tensile_lode_cutoff")
75 143 : : 1.0e-5 * Utility::pow<2>(_f_tol))
76 :
77 : {
78 143 : if (_lode_cutoff < 0)
79 0 : mooseError("tensile_lode_cutoff must not be negative");
80 143 : _ccc = (-std::cos(3.0 * _tt) * (std::cos(_tt) - std::sin(_tt) / std::sqrt(3.0)) -
81 143 : 3.0 * _sin3tt * (std::sin(_tt) + std::cos(_tt) / std::sqrt(3.0))) /
82 143 : (18.0 * Utility::pow<3>(std::cos(3.0 * _tt)));
83 143 : _bbb = (std::sin(6.0 * _tt) * (std::cos(_tt) - std::sin(_tt) / std::sqrt(3.0)) -
84 143 : 6.0 * std::cos(6.0 * _tt) * (std::sin(_tt) + std::cos(_tt) / std::sqrt(3.0))) /
85 : (18.0 * Utility::pow<3>(std::cos(3.0 * _tt)));
86 143 : _aaa = -std::sin(_tt) / std::sqrt(3.0) - _bbb * _sin3tt - _ccc * Utility::pow<2>(_sin3tt) +
87 : std::cos(_tt);
88 143 : }
89 :
90 : Real
91 348432 : SolidMechanicsPlasticTensile::yieldFunction(const RankTwoTensor & stress, Real intnl) const
92 : {
93 348432 : Real mean_stress = stress.trace() / 3.0;
94 348432 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
95 348432 : if (sin3Lode <= _sin3tt)
96 : {
97 : // the non-edge-smoothed version
98 : std::vector<Real> eigvals;
99 237432 : stress.symmetricEigenvalues(eigvals);
100 237432 : return mean_stress + std::sqrt(smooth(stress) + Utility::pow<2>(eigvals[2] - mean_stress)) -
101 237432 : tensile_strength(intnl);
102 237432 : }
103 : else
104 : {
105 : // the edge-smoothed version
106 111000 : Real kk = _aaa + _bbb * sin3Lode + _ccc * Utility::pow<2>(sin3Lode);
107 111000 : Real sibar2 = stress.secondInvariant();
108 111000 : return mean_stress + std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk)) -
109 111000 : tensile_strength(intnl);
110 : }
111 : }
112 :
113 : RankTwoTensor
114 559152 : SolidMechanicsPlasticTensile::dyieldFunction_dstress(const RankTwoTensor & stress,
115 : Real /*intnl*/) const
116 : {
117 559152 : Real mean_stress = stress.trace() / 3.0;
118 559152 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
119 559152 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
120 559152 : if (sin3Lode <= _sin3tt)
121 : {
122 : // the non-edge-smoothed version
123 : std::vector<Real> eigvals;
124 : std::vector<RankTwoTensor> deigvals;
125 380160 : stress.dsymmetricEigenvalues(eigvals, deigvals);
126 380160 : Real denom = std::sqrt(smooth(stress) + Utility::pow<2>(eigvals[2] - mean_stress));
127 380160 : return dmean_stress + (0.5 * dsmooth(stress) * dmean_stress +
128 380160 : (eigvals[2] - mean_stress) * (deigvals[2] - dmean_stress)) /
129 380160 : denom;
130 380160 : }
131 : else
132 : {
133 : // the edge-smoothed version
134 178992 : Real kk = _aaa + _bbb * sin3Lode + _ccc * Utility::pow<2>(sin3Lode);
135 178992 : RankTwoTensor dkk = (_bbb + 2.0 * _ccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff);
136 178992 : Real sibar2 = stress.secondInvariant();
137 178992 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
138 178992 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
139 178992 : return dmean_stress + (0.5 * dsmooth(stress) * dmean_stress +
140 178992 : 0.5 * dsibar2 * Utility::pow<2>(kk) + sibar2 * kk * dkk) /
141 178992 : denom;
142 : }
143 : }
144 :
145 : Real
146 132384 : SolidMechanicsPlasticTensile::dyieldFunction_dintnl(const RankTwoTensor & /*stress*/,
147 : Real intnl) const
148 : {
149 132384 : return -dtensile_strength(intnl);
150 : }
151 :
152 : RankTwoTensor
153 426768 : SolidMechanicsPlasticTensile::flowPotential(const RankTwoTensor & stress, Real intnl) const
154 : {
155 : // This plasticity is associative so
156 426768 : return dyieldFunction_dstress(stress, intnl);
157 : }
158 :
159 : RankFourTensor
160 132384 : SolidMechanicsPlasticTensile::dflowPotential_dstress(const RankTwoTensor & stress,
161 : Real /*intnl*/) const
162 : {
163 132384 : Real mean_stress = stress.trace() / 3.0;
164 132384 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
165 132384 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
166 132384 : if (sin3Lode <= _sin3tt)
167 : {
168 : // the non-edge-smoothed version
169 : std::vector<Real> eigvals;
170 : std::vector<RankTwoTensor> deigvals;
171 : std::vector<RankFourTensor> d2eigvals;
172 90000 : stress.dsymmetricEigenvalues(eigvals, deigvals);
173 90000 : stress.d2symmetricEigenvalues(d2eigvals);
174 :
175 90000 : Real denom = std::sqrt(smooth(stress) + Utility::pow<2>(eigvals[2] - mean_stress));
176 : Real denom3 = Utility::pow<3>(denom);
177 90000 : RankTwoTensor numer_part = deigvals[2] - dmean_stress;
178 : RankTwoTensor numer_full =
179 90000 : 0.5 * dsmooth(stress) * dmean_stress + (eigvals[2] - mean_stress) * numer_part;
180 90000 : Real d2smooth_over_denom = d2smooth(stress) / denom;
181 :
182 90000 : RankFourTensor dr_dstress = (eigvals[2] - mean_stress) * d2eigvals[2] / denom;
183 360000 : for (unsigned i = 0; i < 3; ++i)
184 1080000 : for (unsigned j = 0; j < 3; ++j)
185 3240000 : for (unsigned k = 0; k < 3; ++k)
186 9720000 : for (unsigned l = 0; l < 3; ++l)
187 : {
188 7290000 : dr_dstress(i, j, k, l) +=
189 7290000 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
190 7290000 : dr_dstress(i, j, k, l) += numer_part(i, j) * numer_part(k, l) / denom;
191 7290000 : dr_dstress(i, j, k, l) -= numer_full(i, j) * numer_full(k, l) / denom3;
192 : }
193 90000 : return dr_dstress;
194 90000 : }
195 : else
196 : {
197 : // the edge-smoothed version
198 42384 : RankTwoTensor dsin3Lode = stress.dsin3Lode(_lode_cutoff);
199 42384 : Real kk = _aaa + _bbb * sin3Lode + _ccc * Utility::pow<2>(sin3Lode);
200 42384 : RankTwoTensor dkk = (_bbb + 2.0 * _ccc * sin3Lode) * dsin3Lode;
201 42384 : RankFourTensor d2kk = (_bbb + 2.0 * _ccc * sin3Lode) * stress.d2sin3Lode(_lode_cutoff);
202 169536 : for (unsigned i = 0; i < 3; ++i)
203 508608 : for (unsigned j = 0; j < 3; ++j)
204 1525824 : for (unsigned k = 0; k < 3; ++k)
205 4577472 : for (unsigned l = 0; l < 3; ++l)
206 3433104 : d2kk(i, j, k, l) += 2.0 * _ccc * dsin3Lode(i, j) * dsin3Lode(k, l);
207 :
208 42384 : Real sibar2 = stress.secondInvariant();
209 42384 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
210 42384 : RankFourTensor d2sibar2 = stress.d2secondInvariant();
211 :
212 42384 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
213 : Real denom3 = Utility::pow<3>(denom);
214 42384 : Real d2smooth_over_denom = d2smooth(stress) / denom;
215 : RankTwoTensor numer_full =
216 42384 : 0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * kk * kk + sibar2 * kk * dkk;
217 :
218 42384 : RankFourTensor dr_dstress = (0.5 * d2sibar2 * Utility::pow<2>(kk) + sibar2 * kk * d2kk) / denom;
219 169536 : for (unsigned i = 0; i < 3; ++i)
220 508608 : for (unsigned j = 0; j < 3; ++j)
221 1525824 : for (unsigned k = 0; k < 3; ++k)
222 4577472 : for (unsigned l = 0; l < 3; ++l)
223 : {
224 3433104 : dr_dstress(i, j, k, l) +=
225 3433104 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
226 3433104 : dr_dstress(i, j, k, l) +=
227 3433104 : (dsibar2(i, j) * dkk(k, l) * kk + dkk(i, j) * dsibar2(k, l) * kk +
228 3433104 : sibar2 * dkk(i, j) * dkk(k, l)) /
229 : denom;
230 3433104 : dr_dstress(i, j, k, l) -= numer_full(i, j) * numer_full(k, l) / denom3;
231 : }
232 42384 : return dr_dstress;
233 : }
234 : }
235 :
236 : RankTwoTensor
237 132384 : SolidMechanicsPlasticTensile::dflowPotential_dintnl(const RankTwoTensor & /*stress*/,
238 : Real /*intnl*/) const
239 : {
240 132384 : return RankTwoTensor();
241 : }
242 :
243 : Real
244 348432 : SolidMechanicsPlasticTensile::tensile_strength(const Real internal_param) const
245 : {
246 348432 : return _strength.value(internal_param);
247 : }
248 :
249 : Real
250 132384 : SolidMechanicsPlasticTensile::dtensile_strength(const Real internal_param) const
251 : {
252 132384 : return _strength.derivative(internal_param);
253 : }
254 :
255 : Real
256 1039968 : SolidMechanicsPlasticTensile::smooth(const RankTwoTensor & stress) const
257 : {
258 1039968 : Real smoother2 = _small_smoother2;
259 1039968 : if (_tip_scheme == "cap")
260 : {
261 84480 : Real x = stress.trace() / 3.0 - _cap_start;
262 : Real p = 0;
263 84480 : if (x > 0)
264 72384 : p = x * (1 - std::exp(-_cap_rate * x));
265 84480 : smoother2 += Utility::pow<2>(p);
266 : }
267 1039968 : return smoother2;
268 : }
269 :
270 : Real
271 691536 : SolidMechanicsPlasticTensile::dsmooth(const RankTwoTensor & stress) const
272 : {
273 : Real dsmoother2 = 0;
274 691536 : if (_tip_scheme == "cap")
275 : {
276 58944 : Real x = stress.trace() / 3.0 - _cap_start;
277 : Real p = 0;
278 : Real dp_dx = 0;
279 58944 : if (x > 0)
280 : {
281 50304 : p = x * (1 - std::exp(-_cap_rate * x));
282 50304 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
283 : }
284 58944 : dsmoother2 += 2.0 * p * dp_dx;
285 : }
286 691536 : return dsmoother2;
287 : }
288 :
289 : Real
290 132384 : SolidMechanicsPlasticTensile::d2smooth(const RankTwoTensor & stress) const
291 : {
292 : Real d2smoother2 = 0;
293 132384 : if (_tip_scheme == "cap")
294 : {
295 12768 : Real x = stress.trace() / 3.0 - _cap_start;
296 : Real p = 0;
297 : Real dp_dx = 0;
298 : Real d2p_dx2 = 0;
299 12768 : if (x > 0)
300 : {
301 11040 : p = x * (1 - std::exp(-_cap_rate * x));
302 11040 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
303 11040 : d2p_dx2 = 2.0 * _cap_rate * std::exp(-_cap_rate * x) -
304 11040 : x * Utility::pow<2>(_cap_rate) * std::exp(-_cap_rate * x);
305 : }
306 12768 : d2smoother2 += 2.0 * Utility::pow<2>(dp_dx) + 2.0 * p * d2p_dx2;
307 : }
308 132384 : return d2smoother2;
309 : }
310 :
311 : std::string
312 0 : SolidMechanicsPlasticTensile::modelName() const
313 : {
314 0 : return "Tensile";
315 : }
|
