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 "SolidMechanicsPlasticMohrCoulomb.h"
11 : #include "RankFourTensor.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", SolidMechanicsPlasticMohrCoulomb);
15 : registerMooseObjectRenamed("SolidMechanicsApp",
16 : TensorMechanicsPlasticMohrCoulomb,
17 : "01/01/2025 00:00",
18 : SolidMechanicsPlasticMohrCoulomb);
19 :
20 : InputParameters
21 276 : SolidMechanicsPlasticMohrCoulomb::validParams()
22 : {
23 276 : InputParameters params = SolidMechanicsPlasticModel::validParams();
24 552 : params.addRequiredParam<UserObjectName>(
25 : "cohesion",
26 : "A SolidMechanicsHardening UserObject that defines hardening of the cohesion. "
27 : "Physically the cohesion should not be negative.");
28 552 : params.addRequiredParam<UserObjectName>(
29 : "friction_angle",
30 : "A SolidMechanicsHardening UserObject that defines hardening of the "
31 : "friction angle (in radians). Physically the friction angle should be "
32 : "between 0 and 90deg.");
33 552 : params.addRequiredParam<UserObjectName>(
34 : "dilation_angle",
35 : "A SolidMechanicsHardening UserObject that defines hardening of the "
36 : "dilation angle (in radians). Usually the dilation angle is not greater "
37 : "than the friction angle, and it is between 0 and 90deg.");
38 828 : params.addRangeCheckedParam<Real>(
39 : "mc_edge_smoother",
40 552 : 25.0,
41 : "mc_edge_smoother>=0 & mc_edge_smoother<=30",
42 : "Smoothing parameter: the edges of the cone are smoothed by the given amount.");
43 552 : MooseEnum tip_scheme("hyperbolic cap", "hyperbolic");
44 552 : params.addParam<MooseEnum>(
45 : "tip_scheme", tip_scheme, "Scheme by which the pyramid's tip will be smoothed.");
46 552 : params.addRequiredRangeCheckedParam<Real>("mc_tip_smoother",
47 : "mc_tip_smoother>=0",
48 : "Smoothing parameter: the cone vertex at mean = "
49 : "cohesion*cot(friction_angle), will be smoothed by "
50 : "the given amount. Typical value is 0.1*cohesion");
51 552 : params.addParam<Real>(
52 : "cap_start",
53 552 : 0.0,
54 : "For the 'cap' tip_scheme, smoothing is performed in the stress_mean > cap_start region");
55 828 : params.addRangeCheckedParam<Real>("cap_rate",
56 552 : 0.0,
57 : "cap_rate>=0",
58 : "For the 'cap' tip_scheme, this controls how quickly the cap "
59 : "degenerates to a hemisphere: small values mean a slow "
60 : "degeneration to a hemisphere (and zero means the 'cap' will "
61 : "be totally inactive). Typical value is 1/tensile_strength");
62 552 : params.addParam<Real>("mc_lode_cutoff",
63 : "If the second invariant of stress is less than this "
64 : "amount, the Lode angle is assumed to be zero. This is "
65 : "to gaurd against precision-loss problems, and this "
66 : "parameter should be set small. Default = "
67 : "0.00001*((yield_Function_tolerance)^2)");
68 276 : params.addClassDescription("Non-associative Mohr-Coulomb plasticity with hardening/softening");
69 :
70 276 : return params;
71 276 : }
72 :
73 142 : SolidMechanicsPlasticMohrCoulomb::SolidMechanicsPlasticMohrCoulomb(
74 142 : const InputParameters & parameters)
75 : : SolidMechanicsPlasticModel(parameters),
76 142 : _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
77 142 : _phi(getUserObject<SolidMechanicsHardeningModel>("friction_angle")),
78 142 : _psi(getUserObject<SolidMechanicsHardeningModel>("dilation_angle")),
79 284 : _tip_scheme(getParam<MooseEnum>("tip_scheme")),
80 284 : _small_smoother2(Utility::pow<2>(getParam<Real>("mc_tip_smoother"))),
81 284 : _cap_start(getParam<Real>("cap_start")),
82 284 : _cap_rate(getParam<Real>("cap_rate")),
83 284 : _tt(getParam<Real>("mc_edge_smoother") * libMesh::pi / 180.0),
84 142 : _costt(std::cos(_tt)),
85 142 : _sintt(std::sin(_tt)),
86 142 : _cos3tt(std::cos(3 * _tt)),
87 142 : _sin3tt(std::sin(3 * _tt)),
88 142 : _cos6tt(std::cos(6 * _tt)),
89 142 : _sin6tt(std::sin(6 * _tt)),
90 146 : _lode_cutoff(parameters.isParamValid("mc_lode_cutoff") ? getParam<Real>("mc_lode_cutoff")
91 142 : : 1.0E-5 * Utility::pow<2>(_f_tol))
92 :
93 : {
94 142 : if (_lode_cutoff < 0)
95 2 : mooseError("mc_lode_cutoff must not be negative");
96 :
97 : // With arbitary UserObjects, it is impossible to check everything, and
98 : // I think this is the best I can do
99 140 : if (phi(0) < 0 || psi(0) < 0 || phi(0) > libMesh::pi / 2.0 || psi(0) > libMesh::pi / 2.0)
100 0 : mooseError("Mohr-Coulomb friction and dilation angles must lie in [0, Pi/2]");
101 140 : if (phi(0) < psi(0))
102 2 : mooseError("Mohr-Coulomb friction angle must not be less than Mohr-Coulomb dilation angle");
103 138 : if (cohesion(0) < 0)
104 0 : mooseError("Mohr-Coulomb cohesion must not be negative");
105 :
106 : // check Abbo et al's convexity constraint (Eqn c.18 in their paper)
107 : // With an arbitrary UserObject, it is impossible to check for all angles
108 : // I think the following is the best we can do
109 138 : Real sin_angle = std::sin(std::max(phi(0), psi(0)));
110 138 : sin_angle = std::max(sin_angle, std::sin(std::max(phi(1E6), psi(1E6))));
111 138 : Real rhs = std::sqrt(3) * (35 * std::sin(_tt) + 14 * std::sin(5 * _tt) - 5 * std::sin(7 * _tt)) /
112 138 : 16 / Utility::pow<5>(std::cos(_tt)) / (11 - 10 * std::cos(2 * _tt));
113 138 : if (rhs <= sin_angle)
114 4 : mooseError("Mohr-Coulomb edge smoothing angle is too small and a non-convex yield surface will "
115 : "result. Please choose a larger value");
116 134 : }
117 :
118 : Real
119 1294800 : SolidMechanicsPlasticMohrCoulomb::yieldFunction(const RankTwoTensor & stress, Real intnl) const
120 : {
121 1294800 : Real mean_stress = stress.trace() / 3.0;
122 1294800 : Real sinphi = std::sin(phi(intnl));
123 1294800 : Real cosphi = std::cos(phi(intnl));
124 1294800 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
125 1294800 : if (std::abs(sin3Lode) <= _sin3tt)
126 : {
127 : // the non-edge-smoothed version
128 : std::vector<Real> eigvals;
129 344624 : stress.symmetricEigenvalues(eigvals);
130 689248 : return mean_stress * sinphi +
131 344624 : std::sqrt(smooth(stress) +
132 344624 : 0.25 * Utility::pow<2>(eigvals[2] - eigvals[0] +
133 344624 : (eigvals[2] + eigvals[0] - 2 * mean_stress) * sinphi)) -
134 344624 : cohesion(intnl) * cosphi;
135 344624 : }
136 : else
137 : {
138 : // the edge-smoothed version
139 : Real aaa, bbb, ccc;
140 950176 : abbo(sin3Lode, sinphi, aaa, bbb, ccc);
141 950176 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
142 950176 : Real sibar2 = stress.secondInvariant();
143 950176 : return mean_stress * sinphi + std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk)) -
144 950176 : cohesion(intnl) * cosphi;
145 : }
146 : }
147 :
148 : RankTwoTensor
149 2140968 : SolidMechanicsPlasticMohrCoulomb::df_dsig(const RankTwoTensor & stress, const Real sin_angle) const
150 : {
151 2140968 : Real mean_stress = stress.trace() / 3.0;
152 2140968 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
153 2140968 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
154 2140968 : if (std::abs(sin3Lode) <= _sin3tt)
155 : {
156 : // the non-edge-smoothed version
157 : std::vector<Real> eigvals;
158 : std::vector<RankTwoTensor> deigvals;
159 512104 : stress.dsymmetricEigenvalues(eigvals, deigvals);
160 512104 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2.0 * mean_stress) * sin_angle;
161 : RankTwoTensor dtmp =
162 512104 : deigvals[2] - deigvals[0] + (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * sin_angle;
163 512104 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
164 512104 : return dmean_stress * sin_angle +
165 512104 : (0.5 * dsmooth(stress) * dmean_stress + 0.25 * tmp * dtmp) / denom;
166 512104 : }
167 : else
168 : {
169 : // the edge-smoothed version
170 : Real aaa, bbb, ccc;
171 1628864 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
172 1628864 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
173 1628864 : RankTwoTensor dkk = (bbb + 2 * ccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff);
174 1628864 : Real sibar2 = stress.secondInvariant();
175 1628864 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
176 1628864 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
177 1628864 : return dmean_stress * sin_angle + (0.5 * dsmooth(stress) * dmean_stress +
178 1628864 : 0.5 * dsibar2 * Utility::pow<2>(kk) + sibar2 * kk * dkk) /
179 1628864 : denom;
180 : }
181 : }
182 :
183 : RankTwoTensor
184 504888 : SolidMechanicsPlasticMohrCoulomb::dyieldFunction_dstress(const RankTwoTensor & stress,
185 : Real intnl) const
186 : {
187 504888 : Real sinphi = std::sin(phi(intnl));
188 504888 : return df_dsig(stress, sinphi);
189 : }
190 :
191 : Real
192 504888 : SolidMechanicsPlasticMohrCoulomb::dyieldFunction_dintnl(const RankTwoTensor & stress,
193 : Real intnl) const
194 : {
195 504888 : Real sin_angle = std::sin(phi(intnl));
196 504888 : Real cos_angle = std::cos(phi(intnl));
197 504888 : Real dsin_angle = cos_angle * dphi(intnl);
198 504888 : Real dcos_angle = -sin_angle * dphi(intnl);
199 :
200 504888 : Real mean_stress = stress.trace() / 3.0;
201 504888 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
202 504888 : if (std::abs(sin3Lode) <= _sin3tt)
203 : {
204 : // the non-edge-smoothed version
205 : std::vector<Real> eigvals;
206 118928 : stress.symmetricEigenvalues(eigvals);
207 118928 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2 * mean_stress) * sin_angle;
208 118928 : Real dtmp = (eigvals[2] + eigvals[0] - 2 * mean_stress) * dsin_angle;
209 118928 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
210 118928 : return mean_stress * dsin_angle + 0.25 * tmp * dtmp / denom - dcohesion(intnl) * cos_angle -
211 118928 : cohesion(intnl) * dcos_angle;
212 118928 : }
213 : else
214 : {
215 : // the edge-smoothed version
216 : Real aaa, bbb, ccc;
217 385960 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
218 : Real daaa, dbbb, dccc;
219 385960 : dabbo(sin3Lode, sin_angle, daaa, dbbb, dccc);
220 385960 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
221 385960 : Real dkk = (daaa + dbbb * sin3Lode + dccc * Utility::pow<2>(sin3Lode)) * dsin_angle;
222 385960 : Real sibar2 = stress.secondInvariant();
223 385960 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
224 385960 : return mean_stress * dsin_angle + sibar2 * kk * dkk / denom - dcohesion(intnl) * cos_angle -
225 385960 : cohesion(intnl) * dcos_angle;
226 : }
227 : }
228 :
229 : RankTwoTensor
230 1636080 : SolidMechanicsPlasticMohrCoulomb::flowPotential(const RankTwoTensor & stress, Real intnl) const
231 : {
232 1636080 : Real sinpsi = std::sin(psi(intnl));
233 1636080 : return df_dsig(stress, sinpsi);
234 : }
235 :
236 : RankFourTensor
237 504888 : SolidMechanicsPlasticMohrCoulomb::dflowPotential_dstress(const RankTwoTensor & stress,
238 : Real intnl) const
239 : {
240 504888 : RankFourTensor dr_dstress;
241 504888 : Real sin_angle = std::sin(psi(intnl));
242 504888 : Real mean_stress = stress.trace() / 3.0;
243 504888 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
244 504888 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
245 504888 : if (std::abs(sin3Lode) <= _sin3tt)
246 : {
247 : // the non-edge-smoothed version
248 : std::vector<Real> eigvals;
249 : std::vector<RankTwoTensor> deigvals;
250 : std::vector<RankFourTensor> d2eigvals;
251 118928 : stress.dsymmetricEigenvalues(eigvals, deigvals);
252 118928 : stress.d2symmetricEigenvalues(d2eigvals);
253 :
254 118928 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2.0 * mean_stress) * sin_angle;
255 : RankTwoTensor dtmp =
256 118928 : deigvals[2] - deigvals[0] + (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * sin_angle;
257 118928 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
258 : Real denom3 = Utility::pow<3>(denom);
259 118928 : Real d2smooth_over_denom = d2smooth(stress) / denom;
260 118928 : RankTwoTensor numer = dsmooth(stress) * dmean_stress + 0.5 * tmp * dtmp;
261 :
262 118928 : dr_dstress = 0.25 * tmp *
263 237856 : (d2eigvals[2] - d2eigvals[0] + (d2eigvals[2] + d2eigvals[0]) * sin_angle) / denom;
264 :
265 475712 : for (unsigned i = 0; i < 3; ++i)
266 1427136 : for (unsigned j = 0; j < 3; ++j)
267 4281408 : for (unsigned k = 0; k < 3; ++k)
268 12844224 : for (unsigned l = 0; l < 3; ++l)
269 : {
270 9633168 : dr_dstress(i, j, k, l) +=
271 9633168 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
272 9633168 : dr_dstress(i, j, k, l) += 0.25 * dtmp(i, j) * dtmp(k, l) / denom;
273 9633168 : dr_dstress(i, j, k, l) -= 0.25 * numer(i, j) * numer(k, l) / denom3;
274 : }
275 118928 : }
276 : else
277 : {
278 : // the edge-smoothed version
279 : Real aaa, bbb, ccc;
280 385960 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
281 385960 : RankTwoTensor dsin3Lode = stress.dsin3Lode(_lode_cutoff);
282 385960 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
283 385960 : RankTwoTensor dkk = (bbb + 2 * ccc * sin3Lode) * dsin3Lode;
284 385960 : RankFourTensor d2kk = (bbb + 2 * ccc * sin3Lode) * stress.d2sin3Lode(_lode_cutoff);
285 1543840 : for (unsigned i = 0; i < 3; ++i)
286 4631520 : for (unsigned j = 0; j < 3; ++j)
287 13894560 : for (unsigned k = 0; k < 3; ++k)
288 41683680 : for (unsigned l = 0; l < 3; ++l)
289 31262760 : d2kk(i, j, k, l) += 2 * ccc * dsin3Lode(i, j) * dsin3Lode(k, l);
290 :
291 385960 : Real sibar2 = stress.secondInvariant();
292 385960 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
293 385960 : RankFourTensor d2sibar2 = stress.d2secondInvariant();
294 :
295 385960 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
296 : Real denom3 = Utility::pow<3>(denom);
297 385960 : Real d2smooth_over_denom = d2smooth(stress) / denom;
298 : RankTwoTensor numer_full =
299 385960 : 0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * kk * kk + sibar2 * kk * dkk;
300 :
301 385960 : dr_dstress = (0.5 * d2sibar2 * Utility::pow<2>(kk) + sibar2 * kk * d2kk) / denom;
302 1543840 : for (unsigned i = 0; i < 3; ++i)
303 4631520 : for (unsigned j = 0; j < 3; ++j)
304 13894560 : for (unsigned k = 0; k < 3; ++k)
305 41683680 : for (unsigned l = 0; l < 3; ++l)
306 : {
307 31262760 : dr_dstress(i, j, k, l) +=
308 31262760 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
309 31262760 : dr_dstress(i, j, k, l) +=
310 31262760 : (dsibar2(i, j) * dkk(k, l) * kk + dkk(i, j) * dsibar2(k, l) * kk +
311 31262760 : sibar2 * dkk(i, j) * dkk(k, l)) /
312 : denom;
313 31262760 : dr_dstress(i, j, k, l) -= numer_full(i, j) * numer_full(k, l) / denom3;
314 : }
315 : }
316 504888 : return dr_dstress;
317 : }
318 :
319 : RankTwoTensor
320 504888 : SolidMechanicsPlasticMohrCoulomb::dflowPotential_dintnl(const RankTwoTensor & stress,
321 : Real intnl) const
322 : {
323 504888 : Real sin_angle = std::sin(psi(intnl));
324 504888 : Real dsin_angle = std::cos(psi(intnl)) * dpsi(intnl);
325 :
326 504888 : Real mean_stress = stress.trace() / 3.0;
327 504888 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
328 504888 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
329 :
330 504888 : if (std::abs(sin3Lode) <= _sin3tt)
331 : {
332 : // the non-edge-smoothed version
333 : std::vector<Real> eigvals;
334 : std::vector<RankTwoTensor> deigvals;
335 118928 : stress.dsymmetricEigenvalues(eigvals, deigvals);
336 118928 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2.0 * mean_stress) * sin_angle;
337 118928 : Real dtmp_dintnl = (eigvals[2] + eigvals[0] - 2 * mean_stress) * dsin_angle;
338 : RankTwoTensor dtmp_dstress =
339 118928 : deigvals[2] - deigvals[0] + (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * sin_angle;
340 : RankTwoTensor d2tmp_dstress_dintnl =
341 118928 : (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * dsin_angle;
342 118928 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
343 118928 : return dmean_stress * dsin_angle + 0.25 * dtmp_dintnl * dtmp_dstress / denom +
344 118928 : 0.25 * tmp * d2tmp_dstress_dintnl / denom -
345 118928 : 0.5 * (dsmooth(stress) * dmean_stress + 0.5 * tmp * dtmp_dstress) * 0.25 * tmp *
346 118928 : dtmp_dintnl / Utility::pow<3>(denom);
347 118928 : }
348 : else
349 : {
350 : // the edge-smoothed version
351 : Real aaa, bbb, ccc;
352 385960 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
353 385960 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
354 :
355 : Real daaa, dbbb, dccc;
356 385960 : dabbo(sin3Lode, sin_angle, daaa, dbbb, dccc);
357 385960 : Real dkk_dintnl = (daaa + dbbb * sin3Lode + dccc * Utility::pow<2>(sin3Lode)) * dsin_angle;
358 385960 : RankTwoTensor dkk_dstress = (bbb + 2 * ccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff);
359 : RankTwoTensor d2kk_dstress_dintnl =
360 385960 : (dbbb + 2 * dccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff) * dsin_angle;
361 :
362 385960 : Real sibar2 = stress.secondInvariant();
363 385960 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
364 385960 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
365 :
366 385960 : return dmean_stress * dsin_angle +
367 385960 : (dsibar2 * kk * dkk_dintnl + sibar2 * dkk_dintnl * dkk_dstress +
368 385960 : sibar2 * kk * d2kk_dstress_dintnl) /
369 : denom -
370 385960 : (0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * Utility::pow<2>(kk) +
371 771920 : sibar2 * kk * dkk_dstress) *
372 385960 : sibar2 * kk * dkk_dintnl / Utility::pow<3>(denom);
373 : }
374 : }
375 :
376 : Real
377 1799826 : SolidMechanicsPlasticMohrCoulomb::cohesion(const Real internal_param) const
378 : {
379 1799826 : return _cohesion.value(internal_param);
380 : }
381 :
382 : Real
383 504888 : SolidMechanicsPlasticMohrCoulomb::dcohesion(const Real internal_param) const
384 : {
385 504888 : return _cohesion.derivative(internal_param);
386 : }
387 :
388 : Real
389 4104960 : SolidMechanicsPlasticMohrCoulomb::phi(const Real internal_param) const
390 : {
391 4104960 : return _phi.value(internal_param);
392 : }
393 :
394 : Real
395 1009776 : SolidMechanicsPlasticMohrCoulomb::dphi(const Real internal_param) const
396 : {
397 1009776 : return _phi.derivative(internal_param);
398 : }
399 :
400 : Real
401 3151440 : SolidMechanicsPlasticMohrCoulomb::psi(const Real internal_param) const
402 : {
403 3151440 : return _psi.value(internal_param);
404 : }
405 :
406 : Real
407 504888 : SolidMechanicsPlasticMohrCoulomb::dpsi(const Real internal_param) const
408 : {
409 504888 : return _psi.derivative(internal_param);
410 : }
411 :
412 : void
413 3736920 : SolidMechanicsPlasticMohrCoulomb::abbo(
414 : const Real sin3lode, const Real sin_angle, Real & aaa, Real & bbb, Real & ccc) const
415 : {
416 3736920 : Real tmp1 = (sin3lode >= 0 ? _costt - sin_angle * _sintt / std::sqrt(3.0)
417 1469856 : : _costt + sin_angle * _sintt / std::sqrt(3.0));
418 3736920 : Real tmp2 = (sin3lode >= 0 ? _sintt + sin_angle * _costt / std::sqrt(3.0)
419 1469856 : : -_sintt + sin_angle * _costt / std::sqrt(3.0));
420 :
421 3736920 : ccc = -_cos3tt * tmp1;
422 3736920 : ccc += (sin3lode >= 0 ? -3 * _sin3tt * tmp2 : 3 * _sin3tt * tmp2);
423 3736920 : ccc /= 18 * Utility::pow<3>(_cos3tt);
424 :
425 3736920 : bbb = (sin3lode >= 0 ? _sin6tt * tmp1 : -_sin6tt * tmp1);
426 3736920 : bbb -= 6 * _cos6tt * tmp2;
427 3736920 : bbb /= 18 * Utility::pow<3>(_cos3tt);
428 :
429 3736920 : aaa = (sin3lode >= 0 ? -sin_angle * _sintt / std::sqrt(3.0) - bbb * _sin3tt
430 1469856 : : sin_angle * _sintt / std::sqrt(3.0) + bbb * _sin3tt);
431 3736920 : aaa += -ccc * Utility::pow<2>(_sin3tt) + _costt;
432 3736920 : }
433 :
434 : void
435 771920 : SolidMechanicsPlasticMohrCoulomb::dabbo(
436 : const Real sin3lode, const Real /*sin_angle*/, Real & daaa, Real & dbbb, Real & dccc) const
437 : {
438 771920 : Real dtmp1 = (sin3lode >= 0 ? -_sintt / std::sqrt(3.0) : _sintt / std::sqrt(3.0));
439 771920 : Real dtmp2 = _costt / std::sqrt(3.0);
440 :
441 771920 : dccc = -_cos3tt * dtmp1;
442 771920 : dccc += (sin3lode >= 0 ? -3 * _sin3tt * dtmp2 : 3 * _sin3tt * dtmp2);
443 771920 : dccc /= 18 * Utility::pow<3>(_cos3tt);
444 :
445 771920 : dbbb = (sin3lode >= 0 ? _sin6tt * dtmp1 : -_sin6tt * dtmp1);
446 771920 : dbbb -= 6 * _cos6tt * dtmp2;
447 771920 : dbbb /= 18 * Utility::pow<3>(_cos3tt);
448 :
449 771920 : daaa = (sin3lode >= 0 ? -_sintt / std::sqrt(3.0) - dbbb * _sin3tt
450 298656 : : _sintt / std::sqrt(3.0) + dbbb * _sin3tt);
451 771920 : daaa += -dccc * Utility::pow<2>(_sin3tt);
452 771920 : }
453 :
454 : Real
455 4950432 : SolidMechanicsPlasticMohrCoulomb::smooth(const RankTwoTensor & stress) const
456 : {
457 4950432 : Real smoother2 = _small_smoother2;
458 4950432 : if (_tip_scheme == "cap")
459 : {
460 611616 : Real x = stress.trace() / 3.0 - _cap_start;
461 : Real p = 0;
462 611616 : if (x > 0)
463 548640 : p = x * (1 - std::exp(-_cap_rate * x));
464 611616 : smoother2 += Utility::pow<2>(p);
465 : }
466 4950432 : return smoother2;
467 : }
468 :
469 : Real
470 3150744 : SolidMechanicsPlasticMohrCoulomb::dsmooth(const RankTwoTensor & stress) const
471 : {
472 : Real dsmoother2 = 0;
473 3150744 : if (_tip_scheme == "cap")
474 : {
475 392184 : Real x = stress.trace() / 3.0 - _cap_start;
476 : Real p = 0;
477 : Real dp_dx = 0;
478 392184 : if (x > 0)
479 : {
480 351120 : p = x * (1 - std::exp(-_cap_rate * x));
481 351120 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
482 : }
483 392184 : dsmoother2 += 2 * p * dp_dx;
484 : }
485 3150744 : return dsmoother2;
486 : }
487 :
488 : Real
489 504888 : SolidMechanicsPlasticMohrCoulomb::d2smooth(const RankTwoTensor & stress) const
490 : {
491 : Real d2smoother2 = 0;
492 504888 : if (_tip_scheme == "cap")
493 : {
494 62736 : Real x = stress.trace() / 3.0 - _cap_start;
495 : Real p = 0;
496 : Real dp_dx = 0;
497 : Real d2p_dx2 = 0;
498 62736 : if (x > 0)
499 : {
500 55392 : p = x * (1 - std::exp(-_cap_rate * x));
501 55392 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
502 55392 : d2p_dx2 = 2 * _cap_rate * std::exp(-_cap_rate * x) -
503 55392 : x * Utility::pow<2>(_cap_rate) * std::exp(-_cap_rate * x);
504 : }
505 62736 : d2smoother2 += 2 * Utility::pow<2>(dp_dx) + 2 * p * d2p_dx2;
506 : }
507 504888 : return d2smoother2;
508 : }
509 :
510 : std::string
511 0 : SolidMechanicsPlasticMohrCoulomb::modelName() const
512 : {
513 0 : return "MohrCoulomb";
514 : }
|
