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 244 : SolidMechanicsPlasticTensile::validParams()
22 : {
23 244 : InputParameters params = SolidMechanicsPlasticModel::validParams();
24 488 : params.addRequiredParam<UserObjectName>(
25 : "tensile_strength",
26 : "A SolidMechanicsHardening UserObject that defines hardening of the tensile strength");
27 732 : params.addRangeCheckedParam<Real>(
28 : "tensile_edge_smoother",
29 488 : 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 488 : MooseEnum tip_scheme("hyperbolic cap", "hyperbolic");
33 488 : params.addParam<MooseEnum>(
34 : "tip_scheme", tip_scheme, "Scheme by which the pyramid's tip will be smoothed.");
35 488 : 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 488 : params.addParam<Real>(
42 : "cap_start",
43 488 : 0.0,
44 : "For the 'cap' tip_scheme, smoothing is performed in the stress_mean > cap_start region");
45 732 : params.addRangeCheckedParam<Real>("cap_rate",
46 488 : 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 488 : 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 244 : params.addClassDescription(
59 : "Associative tensile plasticity with hardening/softening, and tensile_strength = 1");
60 :
61 244 : return params;
62 244 : }
63 :
64 122 : SolidMechanicsPlasticTensile::SolidMechanicsPlasticTensile(const InputParameters & parameters)
65 : : SolidMechanicsPlasticModel(parameters),
66 122 : _strength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
67 244 : _tip_scheme(getParam<MooseEnum>("tip_scheme")),
68 244 : _small_smoother2(Utility::pow<2>(getParam<Real>("tensile_tip_smoother"))),
69 244 : _cap_start(getParam<Real>("cap_start")),
70 244 : _cap_rate(getParam<Real>("cap_rate")),
71 244 : _tt(getParam<Real>("tensile_edge_smoother") * libMesh::pi / 180.0),
72 122 : _sin3tt(std::sin(3.0 * _tt)),
73 122 : _lode_cutoff(parameters.isParamValid("tensile_lode_cutoff")
74 122 : ? getParam<Real>("tensile_lode_cutoff")
75 122 : : 1.0e-5 * Utility::pow<2>(_f_tol))
76 :
77 : {
78 122 : if (_lode_cutoff < 0)
79 0 : mooseError("tensile_lode_cutoff must not be negative");
80 122 : _ccc = (-std::cos(3.0 * _tt) * (std::cos(_tt) - std::sin(_tt) / std::sqrt(3.0)) -
81 122 : 3.0 * _sin3tt * (std::sin(_tt) + std::cos(_tt) / std::sqrt(3.0))) /
82 122 : (18.0 * Utility::pow<3>(std::cos(3.0 * _tt)));
83 122 : _bbb = (std::sin(6.0 * _tt) * (std::cos(_tt) - std::sin(_tt) / std::sqrt(3.0)) -
84 122 : 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 122 : _aaa = -std::sin(_tt) / std::sqrt(3.0) - _bbb * _sin3tt - _ccc * Utility::pow<2>(_sin3tt) +
87 : std::cos(_tt);
88 122 : }
89 :
90 : Real
91 243312 : SolidMechanicsPlasticTensile::yieldFunction(const RankTwoTensor & stress, Real intnl) const
92 : {
93 243312 : Real mean_stress = stress.trace() / 3.0;
94 243312 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
95 243312 : if (sin3Lode <= _sin3tt)
96 : {
97 : // the non-edge-smoothed version
98 : std::vector<Real> eigvals;
99 166368 : stress.symmetricEigenvalues(eigvals);
100 166368 : return mean_stress + std::sqrt(smooth(stress) + Utility::pow<2>(eigvals[2] - mean_stress)) -
101 166368 : tensile_strength(intnl);
102 : }
103 : else
104 : {
105 : // the edge-smoothed version
106 76944 : Real kk = _aaa + _bbb * sin3Lode + _ccc * Utility::pow<2>(sin3Lode);
107 76944 : Real sibar2 = stress.secondInvariant();
108 76944 : return mean_stress + std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk)) -
109 76944 : tensile_strength(intnl);
110 : }
111 : }
112 :
113 : RankTwoTensor
114 392016 : SolidMechanicsPlasticTensile::dyieldFunction_dstress(const RankTwoTensor & stress,
115 : Real /*intnl*/) const
116 : {
117 392016 : Real mean_stress = stress.trace() / 3.0;
118 392016 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
119 392016 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
120 392016 : if (sin3Lode <= _sin3tt)
121 : {
122 : // the non-edge-smoothed version
123 : std::vector<Real> eigvals;
124 : std::vector<RankTwoTensor> deigvals;
125 267440 : stress.dsymmetricEigenvalues(eigvals, deigvals);
126 267440 : Real denom = std::sqrt(smooth(stress) + Utility::pow<2>(eigvals[2] - mean_stress));
127 267440 : return dmean_stress + (0.5 * dsmooth(stress) * dmean_stress +
128 267440 : (eigvals[2] - mean_stress) * (deigvals[2] - dmean_stress)) /
129 267440 : denom;
130 : }
131 : else
132 : {
133 : // the edge-smoothed version
134 124576 : Real kk = _aaa + _bbb * sin3Lode + _ccc * Utility::pow<2>(sin3Lode);
135 124576 : RankTwoTensor dkk = (_bbb + 2.0 * _ccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff);
136 124576 : Real sibar2 = stress.secondInvariant();
137 124576 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
138 124576 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
139 124576 : return dmean_stress + (0.5 * dsmooth(stress) * dmean_stress +
140 124576 : 0.5 * dsibar2 * Utility::pow<2>(kk) + sibar2 * kk * dkk) /
141 124576 : denom;
142 : }
143 : }
144 :
145 : Real
146 93552 : SolidMechanicsPlasticTensile::dyieldFunction_dintnl(const RankTwoTensor & /*stress*/,
147 : Real intnl) const
148 : {
149 93552 : return -dtensile_strength(intnl);
150 : }
151 :
152 : RankTwoTensor
153 298464 : SolidMechanicsPlasticTensile::flowPotential(const RankTwoTensor & stress, Real intnl) const
154 : {
155 : // This plasticity is associative so
156 298464 : return dyieldFunction_dstress(stress, intnl);
157 : }
158 :
159 : RankFourTensor
160 93552 : SolidMechanicsPlasticTensile::dflowPotential_dstress(const RankTwoTensor & stress,
161 : Real /*intnl*/) const
162 : {
163 93552 : Real mean_stress = stress.trace() / 3.0;
164 93552 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
165 93552 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
166 93552 : 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 63856 : stress.dsymmetricEigenvalues(eigvals, deigvals);
173 63856 : stress.d2symmetricEigenvalues(d2eigvals);
174 :
175 63856 : Real denom = std::sqrt(smooth(stress) + Utility::pow<2>(eigvals[2] - mean_stress));
176 : Real denom3 = Utility::pow<3>(denom);
177 63856 : RankTwoTensor numer_part = deigvals[2] - dmean_stress;
178 : RankTwoTensor numer_full =
179 63856 : 0.5 * dsmooth(stress) * dmean_stress + (eigvals[2] - mean_stress) * numer_part;
180 63856 : Real d2smooth_over_denom = d2smooth(stress) / denom;
181 :
182 63856 : RankFourTensor dr_dstress = (eigvals[2] - mean_stress) * d2eigvals[2] / denom;
183 255424 : for (unsigned i = 0; i < 3; ++i)
184 766272 : for (unsigned j = 0; j < 3; ++j)
185 2298816 : for (unsigned k = 0; k < 3; ++k)
186 6896448 : for (unsigned l = 0; l < 3; ++l)
187 : {
188 5172336 : dr_dstress(i, j, k, l) +=
189 5172336 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
190 5172336 : dr_dstress(i, j, k, l) += numer_part(i, j) * numer_part(k, l) / denom;
191 5172336 : dr_dstress(i, j, k, l) -= numer_full(i, j) * numer_full(k, l) / denom3;
192 : }
193 63856 : return dr_dstress;
194 : }
195 : else
196 : {
197 : // the edge-smoothed version
198 29696 : RankTwoTensor dsin3Lode = stress.dsin3Lode(_lode_cutoff);
199 29696 : Real kk = _aaa + _bbb * sin3Lode + _ccc * Utility::pow<2>(sin3Lode);
200 29696 : RankTwoTensor dkk = (_bbb + 2.0 * _ccc * sin3Lode) * dsin3Lode;
201 29696 : RankFourTensor d2kk = (_bbb + 2.0 * _ccc * sin3Lode) * stress.d2sin3Lode(_lode_cutoff);
202 118784 : for (unsigned i = 0; i < 3; ++i)
203 356352 : for (unsigned j = 0; j < 3; ++j)
204 1069056 : for (unsigned k = 0; k < 3; ++k)
205 3207168 : for (unsigned l = 0; l < 3; ++l)
206 2405376 : d2kk(i, j, k, l) += 2.0 * _ccc * dsin3Lode(i, j) * dsin3Lode(k, l);
207 :
208 29696 : Real sibar2 = stress.secondInvariant();
209 29696 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
210 29696 : RankFourTensor d2sibar2 = stress.d2secondInvariant();
211 :
212 29696 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
213 : Real denom3 = Utility::pow<3>(denom);
214 29696 : Real d2smooth_over_denom = d2smooth(stress) / denom;
215 : RankTwoTensor numer_full =
216 29696 : 0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * kk * kk + sibar2 * kk * dkk;
217 :
218 29696 : RankFourTensor dr_dstress = (0.5 * d2sibar2 * Utility::pow<2>(kk) + sibar2 * kk * d2kk) / denom;
219 118784 : for (unsigned i = 0; i < 3; ++i)
220 356352 : for (unsigned j = 0; j < 3; ++j)
221 1069056 : for (unsigned k = 0; k < 3; ++k)
222 3207168 : for (unsigned l = 0; l < 3; ++l)
223 : {
224 2405376 : dr_dstress(i, j, k, l) +=
225 2405376 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
226 2405376 : dr_dstress(i, j, k, l) +=
227 2405376 : (dsibar2(i, j) * dkk(k, l) * kk + dkk(i, j) * dsibar2(k, l) * kk +
228 2405376 : sibar2 * dkk(i, j) * dkk(k, l)) /
229 : denom;
230 2405376 : dr_dstress(i, j, k, l) -= numer_full(i, j) * numer_full(k, l) / denom3;
231 : }
232 29696 : return dr_dstress;
233 : }
234 : }
235 :
236 : RankTwoTensor
237 93552 : SolidMechanicsPlasticTensile::dflowPotential_dintnl(const RankTwoTensor & /*stress*/,
238 : Real /*intnl*/) const
239 : {
240 93552 : return RankTwoTensor();
241 : }
242 :
243 : Real
244 243312 : SolidMechanicsPlasticTensile::tensile_strength(const Real internal_param) const
245 : {
246 243312 : return _strength.value(internal_param);
247 : }
248 :
249 : Real
250 93552 : SolidMechanicsPlasticTensile::dtensile_strength(const Real internal_param) const
251 : {
252 93552 : return _strength.derivative(internal_param);
253 : }
254 :
255 : Real
256 728880 : SolidMechanicsPlasticTensile::smooth(const RankTwoTensor & stress) const
257 : {
258 728880 : Real smoother2 = _small_smoother2;
259 728880 : if (_tip_scheme == "cap")
260 : {
261 66048 : Real x = stress.trace() / 3.0 - _cap_start;
262 : Real p = 0;
263 66048 : if (x > 0)
264 56640 : p = x * (1 - std::exp(-_cap_rate * x));
265 66048 : smoother2 += Utility::pow<2>(p);
266 : }
267 728880 : return smoother2;
268 : }
269 :
270 : Real
271 485568 : SolidMechanicsPlasticTensile::dsmooth(const RankTwoTensor & stress) const
272 : {
273 : Real dsmoother2 = 0;
274 485568 : if (_tip_scheme == "cap")
275 : {
276 46080 : Real x = stress.trace() / 3.0 - _cap_start;
277 : Real p = 0;
278 : Real dp_dx = 0;
279 46080 : if (x > 0)
280 : {
281 39360 : p = x * (1 - std::exp(-_cap_rate * x));
282 39360 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
283 : }
284 46080 : dsmoother2 += 2.0 * p * dp_dx;
285 : }
286 485568 : return dsmoother2;
287 : }
288 :
289 : Real
290 93552 : SolidMechanicsPlasticTensile::d2smooth(const RankTwoTensor & stress) const
291 : {
292 : Real d2smoother2 = 0;
293 93552 : if (_tip_scheme == "cap")
294 : {
295 9984 : Real x = stress.trace() / 3.0 - _cap_start;
296 : Real p = 0;
297 : Real dp_dx = 0;
298 : Real d2p_dx2 = 0;
299 9984 : if (x > 0)
300 : {
301 8640 : p = x * (1 - std::exp(-_cap_rate * x));
302 8640 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
303 8640 : d2p_dx2 = 2.0 * _cap_rate * std::exp(-_cap_rate * x) -
304 8640 : x * Utility::pow<2>(_cap_rate) * std::exp(-_cap_rate * x);
305 : }
306 9984 : d2smoother2 += 2.0 * Utility::pow<2>(dp_dx) + 2.0 * p * d2p_dx2;
307 : }
308 93552 : return d2smoother2;
309 : }
310 :
311 : std::string
312 0 : SolidMechanicsPlasticTensile::modelName() const
313 : {
314 0 : return "Tensile";
315 : }
|
