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