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