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 "SolidMechanicsPlasticWeakPlaneShear.h"
11 : #include "RankFourTensor.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", SolidMechanicsPlasticWeakPlaneShear);
15 : registerMooseObjectRenamed("SolidMechanicsApp",
16 : TensorMechanicsPlasticWeakPlaneShear,
17 : "01/01/2025 00:00",
18 : SolidMechanicsPlasticWeakPlaneShear);
19 :
20 : InputParameters
21 146 : SolidMechanicsPlasticWeakPlaneShear::validParams()
22 : {
23 146 : InputParameters params = SolidMechanicsPlasticModel::validParams();
24 292 : params.addRequiredParam<UserObjectName>(
25 : "cohesion",
26 : "A SolidMechanicsHardening UserObject that defines hardening of the cohesion. "
27 : "Physically the cohesion should not be negative.");
28 292 : params.addRequiredParam<UserObjectName>("tan_friction_angle",
29 : "A SolidMechanicsHardening UserObject that defines "
30 : "hardening of tan(friction angle). Physically the "
31 : "friction angle should be between 0 and 90deg.");
32 292 : params.addRequiredParam<UserObjectName>(
33 : "tan_dilation_angle",
34 : "A SolidMechanicsHardening UserObject that defines hardening of the "
35 : "tan(dilation angle). Usually the dilation angle is not greater than "
36 : "the friction angle, and it is between 0 and 90deg.");
37 292 : MooseEnum tip_scheme("hyperbolic cap", "hyperbolic");
38 292 : params.addParam<MooseEnum>(
39 : "tip_scheme", tip_scheme, "Scheme by which the cone's tip will be smoothed.");
40 292 : params.addRequiredRangeCheckedParam<Real>(
41 : "smoother",
42 : "smoother>=0",
43 : "For the 'hyperbolic' tip_scheme, the cone vertex at shear-stress "
44 : "= 0 will be smoothed by the given amount. For the 'cap' "
45 : "tip_scheme, additional smoothing will occur. Typical value is "
46 : "0.1*cohesion");
47 292 : params.addParam<Real>(
48 : "cap_start",
49 292 : 0.0,
50 : "For the 'cap' tip_scheme, smoothing is performed in the stress_zz > cap_start region");
51 438 : params.addRangeCheckedParam<Real>("cap_rate",
52 292 : 0.0,
53 : "cap_rate>=0",
54 : "For the 'cap' tip_scheme, this controls how quickly the cap "
55 : "degenerates to a hemisphere: small values mean a slow "
56 : "degeneration to a hemisphere (and zero means the 'cap' will "
57 : "be totally inactive). Typical value is 1/cohesion");
58 146 : params.addClassDescription("Non-associative finite-strain weak-plane shear perfect plasticity. "
59 : "Here cohesion = 1, tan(phi) = 1 = tan(psi)");
60 :
61 146 : return params;
62 146 : }
63 :
64 74 : SolidMechanicsPlasticWeakPlaneShear::SolidMechanicsPlasticWeakPlaneShear(
65 74 : const InputParameters & parameters)
66 : : SolidMechanicsPlasticModel(parameters),
67 74 : _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
68 74 : _tan_phi(getUserObject<SolidMechanicsHardeningModel>("tan_friction_angle")),
69 74 : _tan_psi(getUserObject<SolidMechanicsHardeningModel>("tan_dilation_angle")),
70 148 : _tip_scheme(getParam<MooseEnum>("tip_scheme")),
71 148 : _small_smoother2(Utility::pow<2>(getParam<Real>("smoother"))),
72 148 : _cap_start(getParam<Real>("cap_start")),
73 222 : _cap_rate(getParam<Real>("cap_rate"))
74 : {
75 : // With arbitary UserObjects, it is impossible to check everything, and
76 : // I think this is the best I can do
77 74 : if (tan_phi(0) < 0 || tan_psi(0) < 0)
78 0 : mooseError("Weak-Plane-Shear friction and dilation angles must lie in [0, Pi/2]");
79 74 : if (tan_phi(0) < tan_psi(0))
80 2 : mooseError(
81 : "Weak-Plane-Shear friction angle must not be less than Weak-Plane-Shear dilation angle");
82 72 : if (cohesion(0) < 0)
83 0 : mooseError("Weak-Plane-Shear cohesion must not be negative");
84 72 : }
85 :
86 : Real
87 271440 : SolidMechanicsPlasticWeakPlaneShear::yieldFunction(const RankTwoTensor & stress, Real intnl) const
88 : {
89 : // note that i explicitly symmeterise in preparation for Cosserat
90 271440 : return std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
91 271440 : Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0) + smooth(stress)) +
92 271440 : stress(2, 2) * tan_phi(intnl) - cohesion(intnl);
93 : }
94 :
95 : RankTwoTensor
96 363520 : SolidMechanicsPlasticWeakPlaneShear::df_dsig(const RankTwoTensor & stress,
97 : Real _tan_phi_or_psi) const
98 : {
99 363520 : RankTwoTensor deriv; // the constructor zeroes this
100 :
101 363520 : Real tau = std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
102 363520 : Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0) + smooth(stress));
103 : // note that i explicitly symmeterise in preparation for Cosserat
104 363520 : if (tau == 0.0)
105 : {
106 : // the derivative is not defined here, but i want to set it nonzero
107 : // because otherwise the return direction might be too crazy
108 0 : deriv(0, 2) = deriv(2, 0) = 0.5;
109 0 : deriv(1, 2) = deriv(2, 1) = 0.5;
110 : }
111 : else
112 : {
113 363520 : deriv(0, 2) = deriv(2, 0) = 0.25 * (stress(0, 2) + stress(2, 0)) / tau;
114 363520 : deriv(1, 2) = deriv(2, 1) = 0.25 * (stress(1, 2) + stress(2, 1)) / tau;
115 363520 : deriv(2, 2) = 0.5 * dsmooth(stress) / tau;
116 : }
117 363520 : deriv(2, 2) += _tan_phi_or_psi;
118 363520 : return deriv;
119 : }
120 :
121 : RankTwoTensor
122 98528 : SolidMechanicsPlasticWeakPlaneShear::dyieldFunction_dstress(const RankTwoTensor & stress,
123 : Real intnl) const
124 : {
125 98528 : return df_dsig(stress, tan_phi(intnl));
126 : }
127 :
128 : Real
129 98528 : SolidMechanicsPlasticWeakPlaneShear::dyieldFunction_dintnl(const RankTwoTensor & stress,
130 : Real intnl) const
131 : {
132 98528 : return stress(2, 2) * dtan_phi(intnl) - dcohesion(intnl);
133 : }
134 :
135 : RankTwoTensor
136 264992 : SolidMechanicsPlasticWeakPlaneShear::flowPotential(const RankTwoTensor & stress, Real intnl) const
137 : {
138 264992 : return df_dsig(stress, tan_psi(intnl));
139 : }
140 :
141 : RankFourTensor
142 98528 : SolidMechanicsPlasticWeakPlaneShear::dflowPotential_dstress(const RankTwoTensor & stress,
143 : Real /*intnl*/) const
144 : {
145 98528 : RankFourTensor dr_dstress;
146 98528 : Real tau = std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
147 98528 : Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0) + smooth(stress));
148 98528 : if (tau == 0.0)
149 : return dr_dstress;
150 :
151 : // note that i explicitly symmeterise
152 98528 : RankTwoTensor dtau;
153 98528 : dtau(0, 2) = dtau(2, 0) = 0.25 * (stress(0, 2) + stress(2, 0)) / tau;
154 98528 : dtau(1, 2) = dtau(2, 1) = 0.25 * (stress(1, 2) + stress(2, 1)) / tau;
155 98528 : dtau(2, 2) = 0.5 * dsmooth(stress) / tau;
156 :
157 394112 : for (unsigned i = 0; i < 3; ++i)
158 1182336 : for (unsigned j = 0; j < 3; ++j)
159 3547008 : for (unsigned k = 0; k < 3; ++k)
160 10641024 : for (unsigned l = 0; l < 3; ++l)
161 7980768 : dr_dstress(i, j, k, l) = -dtau(i, j) * dtau(k, l) / tau;
162 :
163 : // note that i explicitly symmeterise
164 98528 : dr_dstress(0, 2, 0, 2) += 0.25 / tau;
165 98528 : dr_dstress(0, 2, 2, 0) += 0.25 / tau;
166 98528 : dr_dstress(2, 0, 0, 2) += 0.25 / tau;
167 98528 : dr_dstress(2, 0, 2, 0) += 0.25 / tau;
168 98528 : dr_dstress(1, 2, 1, 2) += 0.25 / tau;
169 98528 : dr_dstress(1, 2, 2, 1) += 0.25 / tau;
170 98528 : dr_dstress(2, 1, 1, 2) += 0.25 / tau;
171 98528 : dr_dstress(2, 1, 2, 1) += 0.25 / tau;
172 98528 : dr_dstress(2, 2, 2, 2) += 0.5 * d2smooth(stress) / tau;
173 :
174 : return dr_dstress;
175 : }
176 :
177 : RankTwoTensor
178 98528 : SolidMechanicsPlasticWeakPlaneShear::dflowPotential_dintnl(const RankTwoTensor & /*stress*/,
179 : Real intnl) const
180 : {
181 98528 : RankTwoTensor dr_dintnl;
182 98528 : dr_dintnl(2, 2) = dtan_psi(intnl);
183 98528 : return dr_dintnl;
184 : }
185 :
186 : Real
187 287320 : SolidMechanicsPlasticWeakPlaneShear::cohesion(const Real internal_param) const
188 : {
189 287320 : return _cohesion.value(internal_param);
190 : }
191 :
192 : Real
193 98528 : SolidMechanicsPlasticWeakPlaneShear::dcohesion(const Real internal_param) const
194 : {
195 98528 : return _cohesion.derivative(internal_param);
196 : }
197 :
198 : Real
199 398932 : SolidMechanicsPlasticWeakPlaneShear::tan_phi(const Real internal_param) const
200 : {
201 398932 : return _tan_phi.value(internal_param);
202 : }
203 :
204 : Real
205 98528 : SolidMechanicsPlasticWeakPlaneShear::dtan_phi(const Real internal_param) const
206 : {
207 98528 : return _tan_phi.derivative(internal_param);
208 : }
209 :
210 : Real
211 293956 : SolidMechanicsPlasticWeakPlaneShear::tan_psi(const Real internal_param) const
212 : {
213 293956 : return _tan_psi.value(internal_param);
214 : }
215 :
216 : Real
217 98528 : SolidMechanicsPlasticWeakPlaneShear::dtan_psi(const Real internal_param) const
218 : {
219 98528 : return _tan_psi.derivative(internal_param);
220 : }
221 :
222 : Real
223 733488 : SolidMechanicsPlasticWeakPlaneShear::smooth(const RankTwoTensor & stress) const
224 : {
225 733488 : Real smoother2 = _small_smoother2;
226 733488 : if (_tip_scheme == "cap")
227 : {
228 39168 : Real x = stress(2, 2) - _cap_start;
229 : Real p = 0;
230 39168 : if (x > 0)
231 39168 : p = x * (1 - std::exp(-_cap_rate * x));
232 39168 : smoother2 += Utility::pow<2>(p);
233 : }
234 733488 : return smoother2;
235 : }
236 :
237 : Real
238 462048 : SolidMechanicsPlasticWeakPlaneShear::dsmooth(const RankTwoTensor & stress) const
239 : {
240 : Real dsmoother2 = 0;
241 462048 : if (_tip_scheme == "cap")
242 : {
243 26880 : Real x = stress(2, 2) - _cap_start;
244 : Real p = 0;
245 : Real dp_dx = 0;
246 26880 : if (x > 0)
247 : {
248 26880 : const Real exp_cap_rate_x = std::exp(-_cap_rate * x);
249 26880 : p = x * (1 - exp_cap_rate_x);
250 26880 : dp_dx = (1 - exp_cap_rate_x) + x * _cap_rate * exp_cap_rate_x;
251 : }
252 26880 : dsmoother2 += 2 * p * dp_dx;
253 : }
254 462048 : return dsmoother2;
255 : }
256 :
257 : Real
258 98528 : SolidMechanicsPlasticWeakPlaneShear::d2smooth(const RankTwoTensor & stress) const
259 : {
260 : Real d2smoother2 = 0;
261 98528 : if (_tip_scheme == "cap")
262 : {
263 6144 : Real x = stress(2, 2) - _cap_start;
264 : Real p = 0;
265 : Real dp_dx = 0;
266 : Real d2p_dx2 = 0;
267 6144 : if (x > 0)
268 : {
269 6144 : const Real exp_cap_rate_x = std::exp(-_cap_rate * x);
270 6144 : p = x * (1.0 - exp_cap_rate_x);
271 6144 : dp_dx = (1.0 - exp_cap_rate_x) + x * _cap_rate * exp_cap_rate_x;
272 6144 : d2p_dx2 = 2.0 * _cap_rate * exp_cap_rate_x - x * Utility::pow<2>(_cap_rate) * exp_cap_rate_x;
273 : }
274 6144 : d2smoother2 += 2.0 * Utility::pow<2>(dp_dx) + 2.0 * p * d2p_dx2;
275 : }
276 98528 : return d2smoother2;
277 : }
278 :
279 : void
280 28848 : SolidMechanicsPlasticWeakPlaneShear::activeConstraints(const std::vector<Real> & f,
281 : const RankTwoTensor & stress,
282 : Real intnl,
283 : const RankFourTensor & Eijkl,
284 : std::vector<bool> & act,
285 : RankTwoTensor & returned_stress) const
286 : {
287 : act.assign(1, false);
288 28848 : returned_stress = stress;
289 :
290 28848 : if (f[0] <= _f_tol)
291 472 : return;
292 :
293 : // in the following i will derive returned_stress for the case smoother=0
294 :
295 28816 : Real tanpsi = tan_psi(intnl);
296 28816 : Real tanphi = tan_phi(intnl);
297 :
298 : // norm is the normal to the yield surface
299 : // with f having psi (dilation angle) instead of phi:
300 : // norm(0) = df/dsig(2,0) = df/dsig(0,2)
301 : // norm(1) = df/dsig(2,1) = df/dsig(1,2)
302 : // norm(2) = df/dsig(2,2)
303 28816 : std::vector<Real> norm(3, 0.0);
304 28816 : const Real tau = std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
305 28816 : Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0));
306 28816 : if (tau > 0.0)
307 : {
308 28376 : norm[0] = 0.25 * (stress(0, 2) + stress(2, 0)) / tau;
309 28376 : norm[1] = 0.25 * (stress(1, 2) + stress(2, 1)) / tau;
310 : }
311 : else
312 : {
313 440 : returned_stress(2, 2) = cohesion(intnl) / tanphi;
314 : act[0] = true;
315 : return;
316 : }
317 28376 : norm[2] = tanpsi;
318 :
319 : // to get the flow directions, we have to multiply norm by Eijkl.
320 : // I assume that E(0,2,0,2) = E(1,2,1,2), and E(2,2,0,2) = 0 = E(0,2,1,2), etc
321 : // with the usual symmetry. This makes finding the returned_stress
322 : // much easier.
323 : // returned_stress = stress - alpha*n
324 : // where alpha is chosen so that f = 0
325 28376 : Real alpha = f[0] / (Eijkl(0, 2, 0, 2) + Eijkl(2, 2, 2, 2) * tanpsi * tanphi);
326 :
327 28376 : if (1 - alpha * Eijkl(0, 2, 0, 2) / tau >= 0)
328 : {
329 : // returning to the "surface" of the cone
330 13008 : returned_stress(2, 2) = stress(2, 2) - alpha * Eijkl(2, 2, 2, 2) * norm[2];
331 13008 : returned_stress(0, 2) = returned_stress(2, 0) =
332 13008 : stress(0, 2) - alpha * 2.0 * Eijkl(0, 2, 0, 2) * norm[0];
333 13008 : returned_stress(1, 2) = returned_stress(2, 1) =
334 13008 : stress(1, 2) - alpha * 2.0 * Eijkl(1, 2, 1, 2) * norm[1];
335 : }
336 : else
337 : {
338 : // returning to the "tip" of the cone
339 15368 : returned_stress(2, 2) = cohesion(intnl) / tanphi;
340 15368 : returned_stress(0, 2) = returned_stress(2, 0) = returned_stress(1, 2) = returned_stress(2, 1) =
341 : 0.0;
342 : }
343 28376 : returned_stress(0, 0) =
344 28376 : stress(0, 0) - Eijkl(0, 0, 2, 2) * (stress(2, 2) - returned_stress(2, 2)) / Eijkl(2, 2, 2, 2);
345 28376 : returned_stress(1, 1) =
346 28376 : stress(1, 1) - Eijkl(1, 1, 2, 2) * (stress(2, 2) - returned_stress(2, 2)) / Eijkl(2, 2, 2, 2);
347 :
348 : act[0] = true;
349 : }
350 :
351 : std::string
352 12 : SolidMechanicsPlasticWeakPlaneShear::modelName() const
353 : {
354 12 : return "WeakPlaneShear";
355 : }
|
