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 170 : SolidMechanicsPlasticWeakPlaneShear::validParams()
22 : {
23 170 : InputParameters params = SolidMechanicsPlasticModel::validParams();
24 340 : params.addRequiredParam<UserObjectName>(
25 : "cohesion",
26 : "A SolidMechanicsHardening UserObject that defines hardening of the cohesion. "
27 : "Physically the cohesion should not be negative.");
28 340 : 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 340 : 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 340 : MooseEnum tip_scheme("hyperbolic cap", "hyperbolic");
38 340 : params.addParam<MooseEnum>(
39 : "tip_scheme", tip_scheme, "Scheme by which the cone's tip will be smoothed.");
40 340 : 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 340 : params.addParam<Real>(
48 : "cap_start",
49 340 : 0.0,
50 : "For the 'cap' tip_scheme, smoothing is performed in the stress_zz > cap_start region");
51 510 : params.addRangeCheckedParam<Real>("cap_rate",
52 340 : 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 170 : params.addClassDescription("Non-associative finite-strain weak-plane shear perfect plasticity. "
59 : "Here cohesion = 1, tan(phi) = 1 = tan(psi)");
60 :
61 170 : return params;
62 170 : }
63 :
64 86 : SolidMechanicsPlasticWeakPlaneShear::SolidMechanicsPlasticWeakPlaneShear(
65 86 : const InputParameters & parameters)
66 : : SolidMechanicsPlasticModel(parameters),
67 86 : _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
68 86 : _tan_phi(getUserObject<SolidMechanicsHardeningModel>("tan_friction_angle")),
69 86 : _tan_psi(getUserObject<SolidMechanicsHardeningModel>("tan_dilation_angle")),
70 172 : _tip_scheme(getParam<MooseEnum>("tip_scheme")),
71 172 : _small_smoother2(Utility::pow<2>(getParam<Real>("smoother"))),
72 172 : _cap_start(getParam<Real>("cap_start")),
73 258 : _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 86 : 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 86 : 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 84 : if (cohesion(0) < 0)
83 0 : mooseError("Weak-Plane-Shear cohesion must not be negative");
84 84 : }
85 :
86 : Real
87 394848 : SolidMechanicsPlasticWeakPlaneShear::yieldFunction(const RankTwoTensor & stress, Real intnl) const
88 : {
89 : // note that i explicitly symmeterise in preparation for Cosserat
90 394848 : return std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
91 394848 : Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0) + smooth(stress)) +
92 394848 : stress(2, 2) * tan_phi(intnl) - cohesion(intnl);
93 : }
94 :
95 : RankTwoTensor
96 524816 : SolidMechanicsPlasticWeakPlaneShear::df_dsig(const RankTwoTensor & stress,
97 : Real _tan_phi_or_psi) const
98 : {
99 524816 : RankTwoTensor deriv; // the constructor zeroes this
100 :
101 524816 : Real tau = std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
102 524816 : Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0) + smooth(stress));
103 : // note that i explicitly symmeterise in preparation for Cosserat
104 524816 : 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 524816 : deriv(0, 2) = deriv(2, 0) = 0.25 * (stress(0, 2) + stress(2, 0)) / tau;
114 524816 : deriv(1, 2) = deriv(2, 1) = 0.25 * (stress(1, 2) + stress(2, 1)) / tau;
115 524816 : deriv(2, 2) = 0.5 * dsmooth(stress) / tau;
116 : }
117 524816 : deriv(2, 2) += _tan_phi_or_psi;
118 524816 : return deriv;
119 : }
120 :
121 : RankTwoTensor
122 141744 : SolidMechanicsPlasticWeakPlaneShear::dyieldFunction_dstress(const RankTwoTensor & stress,
123 : Real intnl) const
124 : {
125 141744 : return df_dsig(stress, tan_phi(intnl));
126 : }
127 :
128 : Real
129 141744 : SolidMechanicsPlasticWeakPlaneShear::dyieldFunction_dintnl(const RankTwoTensor & stress,
130 : Real intnl) const
131 : {
132 141744 : return stress(2, 2) * dtan_phi(intnl) - dcohesion(intnl);
133 : }
134 :
135 : RankTwoTensor
136 383072 : SolidMechanicsPlasticWeakPlaneShear::flowPotential(const RankTwoTensor & stress, Real intnl) const
137 : {
138 383072 : return df_dsig(stress, tan_psi(intnl));
139 : }
140 :
141 : RankFourTensor
142 141744 : SolidMechanicsPlasticWeakPlaneShear::dflowPotential_dstress(const RankTwoTensor & stress,
143 : Real /*intnl*/) const
144 : {
145 141744 : RankFourTensor dr_dstress;
146 141744 : Real tau = std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
147 141744 : Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0) + smooth(stress));
148 141744 : if (tau == 0.0)
149 : return dr_dstress;
150 :
151 : // note that i explicitly symmeterise
152 141744 : RankTwoTensor dtau;
153 141744 : dtau(0, 2) = dtau(2, 0) = 0.25 * (stress(0, 2) + stress(2, 0)) / tau;
154 141744 : dtau(1, 2) = dtau(2, 1) = 0.25 * (stress(1, 2) + stress(2, 1)) / tau;
155 141744 : dtau(2, 2) = 0.5 * dsmooth(stress) / tau;
156 :
157 566976 : for (unsigned i = 0; i < 3; ++i)
158 1700928 : for (unsigned j = 0; j < 3; ++j)
159 5102784 : for (unsigned k = 0; k < 3; ++k)
160 15308352 : for (unsigned l = 0; l < 3; ++l)
161 11481264 : dr_dstress(i, j, k, l) = -dtau(i, j) * dtau(k, l) / tau;
162 :
163 : // note that i explicitly symmeterise
164 141744 : dr_dstress(0, 2, 0, 2) += 0.25 / tau;
165 141744 : dr_dstress(0, 2, 2, 0) += 0.25 / tau;
166 141744 : dr_dstress(2, 0, 0, 2) += 0.25 / tau;
167 141744 : dr_dstress(2, 0, 2, 0) += 0.25 / tau;
168 141744 : dr_dstress(1, 2, 1, 2) += 0.25 / tau;
169 141744 : dr_dstress(1, 2, 2, 1) += 0.25 / tau;
170 141744 : dr_dstress(2, 1, 1, 2) += 0.25 / tau;
171 141744 : dr_dstress(2, 1, 2, 1) += 0.25 / tau;
172 141744 : dr_dstress(2, 2, 2, 2) += 0.5 * d2smooth(stress) / tau;
173 :
174 : return dr_dstress;
175 : }
176 :
177 : RankTwoTensor
178 141744 : SolidMechanicsPlasticWeakPlaneShear::dflowPotential_dintnl(const RankTwoTensor & /*stress*/,
179 : Real intnl) const
180 : {
181 141744 : RankTwoTensor dr_dintnl;
182 141744 : dr_dintnl(2, 2) = dtan_psi(intnl);
183 141744 : return dr_dintnl;
184 : }
185 :
186 : Real
187 417484 : SolidMechanicsPlasticWeakPlaneShear::cohesion(const Real internal_param) const
188 : {
189 417484 : return _cohesion.value(internal_param);
190 : }
191 :
192 : Real
193 141744 : SolidMechanicsPlasticWeakPlaneShear::dcohesion(const Real internal_param) const
194 : {
195 141744 : return _cohesion.derivative(internal_param);
196 : }
197 :
198 : Real
199 578108 : SolidMechanicsPlasticWeakPlaneShear::tan_phi(const Real internal_param) const
200 : {
201 578108 : return _tan_phi.value(internal_param);
202 : }
203 :
204 : Real
205 141744 : SolidMechanicsPlasticWeakPlaneShear::dtan_phi(const Real internal_param) const
206 : {
207 141744 : return _tan_phi.derivative(internal_param);
208 : }
209 :
210 : Real
211 424588 : SolidMechanicsPlasticWeakPlaneShear::tan_psi(const Real internal_param) const
212 : {
213 424588 : return _tan_psi.value(internal_param);
214 : }
215 :
216 : Real
217 141744 : SolidMechanicsPlasticWeakPlaneShear::dtan_psi(const Real internal_param) const
218 : {
219 141744 : return _tan_psi.derivative(internal_param);
220 : }
221 :
222 : Real
223 1061408 : SolidMechanicsPlasticWeakPlaneShear::smooth(const RankTwoTensor & stress) const
224 : {
225 1061408 : Real smoother2 = _small_smoother2;
226 1061408 : if (_tip_scheme == "cap")
227 : {
228 49416 : Real x = stress(2, 2) - _cap_start;
229 : Real p = 0;
230 49416 : if (x > 0)
231 49416 : p = x * (1 - std::exp(-_cap_rate * x));
232 49416 : smoother2 += Utility::pow<2>(p);
233 : }
234 1061408 : return smoother2;
235 : }
236 :
237 : Real
238 666560 : SolidMechanicsPlasticWeakPlaneShear::dsmooth(const RankTwoTensor & stress) const
239 : {
240 : Real dsmoother2 = 0;
241 666560 : if (_tip_scheme == "cap")
242 : {
243 33912 : Real x = stress(2, 2) - _cap_start;
244 : Real p = 0;
245 : Real dp_dx = 0;
246 33912 : if (x > 0)
247 : {
248 33912 : const Real exp_cap_rate_x = std::exp(-_cap_rate * x);
249 33912 : p = x * (1 - exp_cap_rate_x);
250 33912 : dp_dx = (1 - exp_cap_rate_x) + x * _cap_rate * exp_cap_rate_x;
251 : }
252 33912 : dsmoother2 += 2 * p * dp_dx;
253 : }
254 666560 : return dsmoother2;
255 : }
256 :
257 : Real
258 141744 : SolidMechanicsPlasticWeakPlaneShear::d2smooth(const RankTwoTensor & stress) const
259 : {
260 : Real d2smoother2 = 0;
261 141744 : if (_tip_scheme == "cap")
262 : {
263 7752 : Real x = stress(2, 2) - _cap_start;
264 : Real p = 0;
265 : Real dp_dx = 0;
266 : Real d2p_dx2 = 0;
267 7752 : if (x > 0)
268 : {
269 7752 : const Real exp_cap_rate_x = std::exp(-_cap_rate * x);
270 7752 : p = x * (1.0 - exp_cap_rate_x);
271 7752 : dp_dx = (1.0 - exp_cap_rate_x) + x * _cap_rate * exp_cap_rate_x;
272 7752 : d2p_dx2 = 2.0 * _cap_rate * exp_cap_rate_x - x * Utility::pow<2>(_cap_rate) * exp_cap_rate_x;
273 : }
274 7752 : d2smoother2 += 2.0 * Utility::pow<2>(dp_dx) + 2.0 * p * d2p_dx2;
275 : }
276 141744 : return d2smoother2;
277 : }
278 :
279 : void
280 41392 : 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 41392 : returned_stress = stress;
289 :
290 41392 : if (f[0] <= _f_tol)
291 598 : return;
292 :
293 : // in the following i will derive returned_stress for the case smoother=0
294 :
295 41344 : Real tanpsi = tan_psi(intnl);
296 41344 : 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 41344 : std::vector<Real> norm(3, 0.0);
304 41344 : const Real tau = std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
305 41344 : Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0));
306 41344 : if (tau > 0.0)
307 : {
308 40794 : norm[0] = 0.25 * (stress(0, 2) + stress(2, 0)) / tau;
309 40794 : norm[1] = 0.25 * (stress(1, 2) + stress(2, 1)) / tau;
310 : }
311 : else
312 : {
313 550 : returned_stress(2, 2) = cohesion(intnl) / tanphi;
314 : act[0] = true;
315 : return;
316 : }
317 40794 : 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 40794 : Real alpha = f[0] / (Eijkl(0, 2, 0, 2) + Eijkl(2, 2, 2, 2) * tanpsi * tanphi);
326 :
327 40794 : if (1 - alpha * Eijkl(0, 2, 0, 2) / tau >= 0)
328 : {
329 : // returning to the "surface" of the cone
330 18792 : returned_stress(2, 2) = stress(2, 2) - alpha * Eijkl(2, 2, 2, 2) * norm[2];
331 18792 : returned_stress(0, 2) = returned_stress(2, 0) =
332 18792 : stress(0, 2) - alpha * 2.0 * Eijkl(0, 2, 0, 2) * norm[0];
333 18792 : returned_stress(1, 2) = returned_stress(2, 1) =
334 18792 : 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 22002 : returned_stress(2, 2) = cohesion(intnl) / tanphi;
340 22002 : returned_stress(0, 2) = returned_stress(2, 0) = returned_stress(1, 2) = returned_stress(2, 1) =
341 : 0.0;
342 : }
343 40794 : returned_stress(0, 0) =
344 40794 : stress(0, 0) - Eijkl(0, 0, 2, 2) * (stress(2, 2) - returned_stress(2, 2)) / Eijkl(2, 2, 2, 2);
345 40794 : returned_stress(1, 1) =
346 40794 : 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 41344 : }
350 :
351 : std::string
352 15 : SolidMechanicsPlasticWeakPlaneShear::modelName() const
353 : {
354 15 : return "WeakPlaneShear";
355 : }
|
