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 "SolidMechanicsPlasticDruckerPrager.h"
11 : #include "RankFourTensor.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", SolidMechanicsPlasticDruckerPrager);
15 : registerMooseObjectRenamed("SolidMechanicsApp",
16 : TensorMechanicsPlasticDruckerPrager,
17 : "01/01/2025 00:00",
18 : SolidMechanicsPlasticDruckerPrager);
19 :
20 : InputParameters
21 468 : SolidMechanicsPlasticDruckerPrager::validParams()
22 : {
23 468 : InputParameters params = SolidMechanicsPlasticModel::validParams();
24 : MooseEnum mc_interpolation_scheme("outer_tip=0 inner_tip=1 lode_zero=2 inner_edge=3 native=4",
25 936 : "lode_zero");
26 936 : params.addParam<MooseEnum>(
27 : "mc_interpolation_scheme",
28 : mc_interpolation_scheme,
29 : "Scheme by which the Drucker-Prager cohesion, friction angle and dilation angle are set from "
30 : "the Mohr-Coulomb parameters mc_cohesion, mc_friction_angle and mc_dilation_angle. Consider "
31 : "the DP and MC yield surfaces on the deviatoric (octahedral) plane. Outer_tip: the DP "
32 : "circle touches the outer tips of the MC hex. Inner_tip: the DP circle touches the inner "
33 : "tips of the MC hex. Lode_zero: the DP circle intersects the MC hex at lode angle=0. "
34 : "Inner_edge: the DP circle is the largest circle that wholly fits inside the MC hex. "
35 : "Native: The DP cohesion, friction angle and dilation angle are set equal to the mc_ "
36 : "parameters entered.");
37 936 : params.addRequiredParam<UserObjectName>(
38 : "mc_cohesion",
39 : "A SolidMechanicsHardening UserObject that defines hardening of the "
40 : "Mohr-Coulomb cohesion. Physically this should not be negative.");
41 936 : params.addRequiredParam<UserObjectName>(
42 : "mc_friction_angle",
43 : "A SolidMechanicsHardening UserObject that defines hardening of the "
44 : "Mohr-Coulomb friction angle (in radians). Physically this should be "
45 : "between 0 and Pi/2.");
46 936 : params.addRequiredParam<UserObjectName>(
47 : "mc_dilation_angle",
48 : "A SolidMechanicsHardening UserObject that defines hardening of the "
49 : "Mohr-Coulomb dilation angle (in radians). Usually the dilation angle "
50 : "is not greater than the friction angle, and it is between 0 and Pi/2.");
51 468 : params.addClassDescription(
52 : "Non-associative Drucker Prager plasticity with no smoothing of the cone tip.");
53 468 : return params;
54 468 : }
55 :
56 234 : SolidMechanicsPlasticDruckerPrager::SolidMechanicsPlasticDruckerPrager(
57 234 : const InputParameters & parameters)
58 : : SolidMechanicsPlasticModel(parameters),
59 234 : _mc_cohesion(getUserObject<SolidMechanicsHardeningModel>("mc_cohesion")),
60 234 : _mc_phi(getUserObject<SolidMechanicsHardeningModel>("mc_friction_angle")),
61 234 : _mc_psi(getUserObject<SolidMechanicsHardeningModel>("mc_dilation_angle")),
62 468 : _mc_interpolation_scheme(getParam<MooseEnum>("mc_interpolation_scheme")),
63 234 : _zero_cohesion_hardening(_mc_cohesion.modelName().compare("Constant") == 0),
64 234 : _zero_phi_hardening(_mc_phi.modelName().compare("Constant") == 0),
65 468 : _zero_psi_hardening(_mc_psi.modelName().compare("Constant") == 0)
66 : {
67 702 : if (_mc_phi.value(0.0) < 0.0 || _mc_psi.value(0.0) < 0.0 ||
68 702 : _mc_phi.value(0.0) > libMesh::pi / 2.0 || _mc_psi.value(0.0) > libMesh::pi / 2.0)
69 0 : mooseError("SolidMechanicsPlasticDruckerPrager: MC friction and dilation angles must lie in "
70 : "[0, Pi/2]");
71 234 : if (_mc_phi.value(0) < _mc_psi.value(0.0))
72 0 : mooseError("SolidMechanicsPlasticDruckerPrager: MC friction angle must not be less than MC "
73 : "dilation angle");
74 234 : if (_mc_cohesion.value(0.0) < 0)
75 0 : mooseError("SolidMechanicsPlasticDruckerPrager: MC cohesion should not be negative");
76 :
77 234 : initializeAandB(0.0, _aaa, _bbb);
78 234 : initializeB(0.0, dilation, _bbb_flow);
79 234 : }
80 :
81 : Real
82 0 : SolidMechanicsPlasticDruckerPrager::yieldFunction(const RankTwoTensor & stress, Real intnl) const
83 : {
84 : Real aaa;
85 : Real bbb;
86 0 : bothAB(intnl, aaa, bbb);
87 0 : return std::sqrt(stress.secondInvariant()) + stress.trace() * bbb - aaa;
88 : }
89 :
90 : RankTwoTensor
91 0 : SolidMechanicsPlasticDruckerPrager::df_dsig(const RankTwoTensor & stress, Real bbb) const
92 : {
93 0 : return 0.5 * stress.dsecondInvariant() / std::sqrt(stress.secondInvariant()) +
94 0 : stress.dtrace() * bbb;
95 : }
96 :
97 : RankTwoTensor
98 0 : SolidMechanicsPlasticDruckerPrager::dyieldFunction_dstress(const RankTwoTensor & stress,
99 : Real intnl) const
100 : {
101 : Real bbb;
102 0 : onlyB(intnl, friction, bbb);
103 0 : return df_dsig(stress, bbb);
104 : }
105 :
106 : Real
107 60576 : SolidMechanicsPlasticDruckerPrager::dyieldFunction_dintnl(const RankTwoTensor & stress,
108 : Real intnl) const
109 : {
110 : Real daaa;
111 : Real dbbb;
112 60576 : dbothAB(intnl, daaa, dbbb);
113 60576 : return stress.trace() * dbbb - daaa;
114 : }
115 :
116 : RankTwoTensor
117 0 : SolidMechanicsPlasticDruckerPrager::flowPotential(const RankTwoTensor & stress, Real intnl) const
118 : {
119 : Real bbb_flow;
120 0 : onlyB(intnl, dilation, bbb_flow);
121 0 : return df_dsig(stress, bbb_flow);
122 : }
123 :
124 : RankFourTensor
125 0 : SolidMechanicsPlasticDruckerPrager::dflowPotential_dstress(const RankTwoTensor & stress,
126 : Real /*intnl*/) const
127 : {
128 0 : RankFourTensor dr_dstress;
129 0 : dr_dstress = 0.5 * stress.d2secondInvariant() / std::sqrt(stress.secondInvariant());
130 0 : dr_dstress += -0.5 * 0.5 * stress.dsecondInvariant().outerProduct(stress.dsecondInvariant()) /
131 0 : std::pow(stress.secondInvariant(), 1.5);
132 0 : return dr_dstress;
133 : }
134 :
135 : RankTwoTensor
136 0 : SolidMechanicsPlasticDruckerPrager::dflowPotential_dintnl(const RankTwoTensor & stress,
137 : Real intnl) const
138 : {
139 : Real dbbb;
140 0 : donlyB(intnl, dilation, dbbb);
141 0 : return stress.dtrace() * dbbb;
142 : }
143 :
144 : std::string
145 0 : SolidMechanicsPlasticDruckerPrager::modelName() const
146 : {
147 0 : return "DruckerPrager";
148 : }
149 :
150 : void
151 1023744 : SolidMechanicsPlasticDruckerPrager::bothAB(Real intnl, Real & aaa, Real & bbb) const
152 : {
153 1023744 : if (_zero_cohesion_hardening && _zero_phi_hardening)
154 : {
155 986368 : aaa = _aaa;
156 986368 : bbb = _bbb;
157 986368 : return;
158 : }
159 37376 : initializeAandB(intnl, aaa, bbb);
160 : }
161 :
162 : void
163 1483744 : SolidMechanicsPlasticDruckerPrager::onlyB(Real intnl, int fd, Real & bbb) const
164 : {
165 1483744 : if (_zero_phi_hardening && (fd == friction))
166 : {
167 59680 : bbb = _bbb;
168 59680 : return;
169 : }
170 1424064 : if (_zero_psi_hardening && (fd == dilation))
171 : {
172 1353536 : bbb = _bbb_flow;
173 1353536 : return;
174 : }
175 70528 : initializeB(intnl, fd, bbb);
176 : }
177 :
178 : void
179 948624 : SolidMechanicsPlasticDruckerPrager::donlyB(Real intnl, int fd, Real & dbbb) const
180 : {
181 948624 : if (_zero_phi_hardening && (fd == friction))
182 : {
183 0 : dbbb = 0;
184 0 : return;
185 : }
186 948624 : if (_zero_psi_hardening && (fd == dilation))
187 : {
188 905872 : dbbb = 0;
189 905872 : return;
190 : }
191 42752 : const Real s = (fd == friction) ? std::sin(_mc_phi.value(intnl)) : std::sin(_mc_psi.value(intnl));
192 85504 : const Real ds = (fd == friction) ? std::cos(_mc_phi.value(intnl)) * _mc_phi.derivative(intnl)
193 42752 : : std::cos(_mc_psi.value(intnl)) * _mc_psi.derivative(intnl);
194 42752 : switch (_mc_interpolation_scheme)
195 : {
196 0 : case 0: // outer_tip
197 0 : dbbb = 2.0 / std::sqrt(3.0) * (ds / (3.0 - s) + s * ds / Utility::pow<2>(3.0 - s));
198 0 : break;
199 0 : case 1: // inner_tip
200 0 : dbbb = 2.0 / std::sqrt(3.0) * (ds / (3.0 + s) - s * ds / Utility::pow<2>(3.0 + s));
201 0 : break;
202 42752 : case 2: // lode_zero
203 42752 : dbbb = ds / 3.0;
204 42752 : break;
205 : case 3: // inner_edge
206 0 : dbbb = ds / std::sqrt(9.0 + 3.0 * Utility::pow<2>(s)) -
207 0 : 3 * s * s * ds / std::pow(9.0 + 3.0 * Utility::pow<2>(s), 1.5);
208 0 : break;
209 0 : case 4: // native
210 : const Real c =
211 0 : (fd == friction) ? std::cos(_mc_phi.value(intnl)) : std::cos(_mc_psi.value(intnl));
212 : const Real dc = (fd == friction)
213 0 : ? -std::sin(_mc_phi.value(intnl)) * _mc_phi.derivative(intnl)
214 0 : : -std::sin(_mc_psi.value(intnl)) * _mc_psi.derivative(intnl);
215 0 : dbbb = ds / c - s * dc / Utility::pow<2>(c);
216 0 : break;
217 : }
218 : }
219 :
220 : void
221 748448 : SolidMechanicsPlasticDruckerPrager::dbothAB(Real intnl, Real & daaa, Real & dbbb) const
222 : {
223 748448 : if (_zero_cohesion_hardening && _zero_phi_hardening)
224 : {
225 721888 : daaa = 0;
226 721888 : dbbb = 0;
227 : return;
228 : }
229 :
230 26560 : const Real C = _mc_cohesion.value(intnl);
231 26560 : const Real dC = _mc_cohesion.derivative(intnl);
232 26560 : const Real cosphi = std::cos(_mc_phi.value(intnl));
233 26560 : const Real dcosphi = -std::sin(_mc_phi.value(intnl)) * _mc_phi.derivative(intnl);
234 26560 : const Real sinphi = std::sin(_mc_phi.value(intnl));
235 26560 : const Real dsinphi = std::cos(_mc_phi.value(intnl)) * _mc_phi.derivative(intnl);
236 26560 : switch (_mc_interpolation_scheme)
237 : {
238 0 : case 0: // outer_tip
239 0 : daaa = 2.0 * std::sqrt(3.0) *
240 0 : (dC * cosphi / (3.0 - sinphi) + C * dcosphi / (3.0 - sinphi) +
241 0 : C * cosphi * dsinphi / Utility::pow<2>(3.0 - sinphi));
242 0 : dbbb = 2.0 / std::sqrt(3.0) *
243 0 : (dsinphi / (3.0 - sinphi) + sinphi * dsinphi / Utility::pow<2>(3.0 - sinphi));
244 0 : break;
245 0 : case 1: // inner_tip
246 0 : daaa = 2.0 * std::sqrt(3.0) *
247 0 : (dC * cosphi / (3.0 + sinphi) + C * dcosphi / (3.0 + sinphi) -
248 0 : C * cosphi * dsinphi / Utility::pow<2>(3.0 + sinphi));
249 0 : dbbb = 2.0 / std::sqrt(3.0) *
250 0 : (dsinphi / (3.0 + sinphi) - sinphi * dsinphi / Utility::pow<2>(3.0 + sinphi));
251 0 : break;
252 26560 : case 2: // lode_zero
253 26560 : daaa = dC * cosphi + C * dcosphi;
254 26560 : dbbb = dsinphi / 3.0;
255 26560 : break;
256 0 : case 3: // inner_edge
257 0 : daaa = 3.0 * dC * cosphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi)) +
258 0 : 3.0 * C * dcosphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi)) -
259 0 : 3.0 * C * cosphi * 3.0 * sinphi * dsinphi /
260 0 : std::pow(9.0 + 3.0 * Utility::pow<2>(sinphi), 1.5);
261 0 : dbbb = dsinphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi)) -
262 0 : 3.0 * sinphi * sinphi * dsinphi / std::pow(9.0 + 3.0 * Utility::pow<2>(sinphi), 1.5);
263 0 : break;
264 0 : case 4: // native
265 0 : daaa = dC;
266 0 : dbbb = dsinphi / cosphi - sinphi * dcosphi / Utility::pow<2>(cosphi);
267 0 : break;
268 : }
269 : }
270 :
271 : void
272 37610 : SolidMechanicsPlasticDruckerPrager::initializeAandB(Real intnl, Real & aaa, Real & bbb) const
273 : {
274 37610 : const Real C = _mc_cohesion.value(intnl);
275 37610 : const Real cosphi = std::cos(_mc_phi.value(intnl));
276 37610 : const Real sinphi = std::sin(_mc_phi.value(intnl));
277 37610 : switch (_mc_interpolation_scheme)
278 : {
279 24 : case 0: // outer_tip
280 24 : aaa = 2.0 * std::sqrt(3.0) * C * cosphi / (3.0 - sinphi);
281 24 : bbb = 2.0 * sinphi / std::sqrt(3.0) / (3.0 - sinphi);
282 24 : break;
283 24 : case 1: // inner_tip
284 24 : aaa = 2.0 * std::sqrt(3.0) * C * cosphi / (3.0 + sinphi);
285 24 : bbb = 2.0 * sinphi / std::sqrt(3.0) / (3.0 + sinphi);
286 24 : break;
287 37514 : case 2: // lode_zero
288 37514 : aaa = C * cosphi;
289 37514 : bbb = sinphi / 3.0;
290 37514 : break;
291 24 : case 3: // inner_edge
292 24 : aaa = 3.0 * C * cosphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi));
293 24 : bbb = sinphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi));
294 24 : break;
295 24 : case 4: // native
296 24 : aaa = C;
297 24 : bbb = sinphi / cosphi;
298 24 : break;
299 : }
300 37610 : }
301 :
302 : void
303 70762 : SolidMechanicsPlasticDruckerPrager::initializeB(Real intnl, int fd, Real & bbb) const
304 : {
305 70762 : const Real s = (fd == friction) ? std::sin(_mc_phi.value(intnl)) : std::sin(_mc_psi.value(intnl));
306 70762 : switch (_mc_interpolation_scheme)
307 : {
308 24 : case 0: // outer_tip
309 24 : bbb = 2.0 * s / std::sqrt(3.0) / (3.0 - s);
310 24 : break;
311 24 : case 1: // inner_tip
312 24 : bbb = 2.0 * s / std::sqrt(3.0) / (3.0 + s);
313 24 : break;
314 70666 : case 2: // lode_zero
315 70666 : bbb = s / 3.0;
316 70666 : break;
317 24 : case 3: // inner_edge
318 24 : bbb = s / std::sqrt(9.0 + 3.0 * Utility::pow<2>(s));
319 24 : break;
320 24 : case 4: // native
321 : const Real c =
322 24 : (fd == friction) ? std::cos(_mc_phi.value(intnl)) : std::cos(_mc_psi.value(intnl));
323 24 : bbb = s / c;
324 24 : break;
325 : }
326 70762 : }
|
