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