www.mooseframework.org
TensorMechanicsPlasticDruckerPrager.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 #include "libmesh/utility.h"
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<TensorMechanicsPlasticModel>();
15  MooseEnum mc_interpolation_scheme("outer_tip=0 inner_tip=1 lode_zero=2 inner_edge=3 native=4",
16  "lode_zero");
17  params.addParam<MooseEnum>(
18  "mc_interpolation_scheme",
19  mc_interpolation_scheme,
20  "Scheme by which the Drucker-Prager cohesion, friction angle and dilation angle are set from "
21  "the Mohr-Coulomb parameters mc_cohesion, mc_friction_angle and mc_dilation_angle. Consider "
22  "the DP and MC yield surfaces on the devatoric (octahedral) plane. Outer_tip: the DP circle "
23  "touches the outer tips of the MC hex. Inner_tip: the DP circle touches the inner tips of "
24  "the MC hex. Lode_zero: the DP circle intersects the MC hex at lode angle=0. Inner_edge: "
25  "the DP circle is the largest circle that wholey fits inside the MC hex. Native: The DP "
26  "cohesion, friction angle and dilation angle are set equal to the mc_ parameters entered.");
27  params.addRequiredParam<UserObjectName>(
28  "mc_cohesion",
29  "A TensorMechanicsHardening UserObject that defines hardening of the "
30  "Mohr-Coulomb cohesion. Physically this should not be negative.");
31  params.addRequiredParam<UserObjectName>(
32  "mc_friction_angle",
33  "A TensorMechanicsHardening UserObject that defines hardening of the "
34  "Mohr-Coulomb friction angle (in radians). Physically this should be "
35  "between 0 and Pi/2.");
36  params.addRequiredParam<UserObjectName>(
37  "mc_dilation_angle",
38  "A TensorMechanicsHardening UserObject that defines hardening of the "
39  "Mohr-Coulomb dilation angle (in radians). Usually the dilation angle "
40  "is not greater than the friction angle, and it is between 0 and Pi/2.");
41  params.addClassDescription(
42  "Non-associative Drucker Prager plasticity with no smoothing of the cone tip.");
43  return params;
44 }
45 
47  const InputParameters & parameters)
48  : TensorMechanicsPlasticModel(parameters),
49  _mc_cohesion(getUserObject<TensorMechanicsHardeningModel>("mc_cohesion")),
50  _mc_phi(getUserObject<TensorMechanicsHardeningModel>("mc_friction_angle")),
51  _mc_psi(getUserObject<TensorMechanicsHardeningModel>("mc_dilation_angle")),
52  _mc_interpolation_scheme(getParam<MooseEnum>("mc_interpolation_scheme")),
53  _zero_cohesion_hardening(_mc_cohesion.modelName().compare("Constant") == 0),
54  _zero_phi_hardening(_mc_phi.modelName().compare("Constant") == 0),
55  _zero_psi_hardening(_mc_psi.modelName().compare("Constant") == 0)
56 {
57  if (_mc_phi.value(0.0) < 0.0 || _mc_psi.value(0.0) < 0.0 ||
58  _mc_phi.value(0.0) > libMesh::pi / 2.0 || _mc_psi.value(0.0) > libMesh::pi / 2.0)
59  mooseError("TensorMechanicsPlasticDruckerPrager: MC friction and dilation angles must lie in "
60  "[0, Pi/2]");
61  if (_mc_phi.value(0) < _mc_psi.value(0.0))
62  mooseError("TensorMechanicsPlasticDruckerPrager: MC friction angle must not be less than MC "
63  "dilation angle");
64  if (_mc_cohesion.value(0.0) < 0)
65  mooseError("TensorMechanicsPlasticDruckerPrager: MC cohesion should not be negative");
66 
67  initializeAandB(0.0, _aaa, _bbb);
69 }
70 
71 Real
72 TensorMechanicsPlasticDruckerPrager::yieldFunction(const RankTwoTensor & stress, Real intnl) const
73 {
74  Real aaa;
75  Real bbb;
76  bothAB(intnl, aaa, bbb);
77  return std::sqrt(stress.secondInvariant()) + stress.trace() * bbb - aaa;
78 }
79 
80 RankTwoTensor
81 TensorMechanicsPlasticDruckerPrager::df_dsig(const RankTwoTensor & stress, Real bbb) const
82 {
83  return 0.5 * stress.dsecondInvariant() / std::sqrt(stress.secondInvariant()) +
84  stress.dtrace() * bbb;
85 }
86 
87 RankTwoTensor
89  Real intnl) const
90 {
91  Real bbb;
92  onlyB(intnl, friction, bbb);
93  return df_dsig(stress, bbb);
94 }
95 
96 Real
98  Real intnl) const
99 {
100  Real daaa;
101  Real dbbb;
102  dbothAB(intnl, daaa, dbbb);
103  return stress.trace() * dbbb - daaa;
104 }
105 
106 RankTwoTensor
107 TensorMechanicsPlasticDruckerPrager::flowPotential(const RankTwoTensor & stress, Real intnl) const
108 {
109  Real bbb_flow;
110  onlyB(intnl, dilation, bbb_flow);
111  return df_dsig(stress, bbb_flow);
112 }
113 
114 RankFourTensor
116  Real /*intnl*/) const
117 {
118  RankFourTensor dr_dstress;
119  dr_dstress = 0.5 * stress.d2secondInvariant() / std::sqrt(stress.secondInvariant());
120  dr_dstress += -0.5 * 0.5 * stress.dsecondInvariant().outerProduct(stress.dsecondInvariant()) /
121  std::pow(stress.secondInvariant(), 1.5);
122  return dr_dstress;
123 }
124 
125 RankTwoTensor
127  Real intnl) const
128 {
129  Real dbbb;
130  donlyB(intnl, dilation, dbbb);
131  return stress.dtrace() * dbbb;
132 }
133 
134 std::string
136 {
137  return "DruckerPrager";
138 }
139 
140 void
141 TensorMechanicsPlasticDruckerPrager::bothAB(Real intnl, Real & aaa, Real & bbb) const
142 {
144  {
145  aaa = _aaa;
146  bbb = _bbb;
147  return;
148  }
149  initializeAandB(intnl, aaa, bbb);
150 }
151 
152 void
153 TensorMechanicsPlasticDruckerPrager::onlyB(Real intnl, int fd, Real & bbb) const
154 {
155  if (_zero_phi_hardening && (fd == friction))
156  {
157  bbb = _bbb;
158  return;
159  }
160  if (_zero_psi_hardening && (fd == dilation))
161  {
162  bbb = _bbb_flow;
163  return;
164  }
165  initializeB(intnl, fd, bbb);
166 }
167 
168 void
169 TensorMechanicsPlasticDruckerPrager::donlyB(Real intnl, int fd, Real & dbbb) const
170 {
171  if (_zero_phi_hardening && (fd == friction))
172  {
173  dbbb = 0;
174  return;
175  }
176  if (_zero_psi_hardening && (fd == dilation))
177  {
178  dbbb = 0;
179  return;
180  }
181  const Real s = (fd == friction) ? std::sin(_mc_phi.value(intnl)) : std::sin(_mc_psi.value(intnl));
182  const Real ds = (fd == friction) ? std::cos(_mc_phi.value(intnl)) * _mc_phi.derivative(intnl)
183  : std::cos(_mc_psi.value(intnl)) * _mc_psi.derivative(intnl);
184  switch (_mc_interpolation_scheme)
185  {
186  case 0: // outer_tip
187  dbbb = 2.0 / std::sqrt(3.0) * (ds / (3.0 - s) + s * ds / Utility::pow<2>(3.0 - s));
188  break;
189  case 1: // inner_tip
190  dbbb = 2.0 / std::sqrt(3.0) * (ds / (3.0 + s) - s * ds / Utility::pow<2>(3.0 + s));
191  break;
192  case 2: // lode_zero
193  dbbb = ds / 3.0;
194  break;
195  case 3: // inner_edge
196  dbbb = ds / std::sqrt(9.0 + 3.0 * Utility::pow<2>(s)) -
197  3 * s * s * ds / std::pow(9.0 + 3.0 * Utility::pow<2>(s), 1.5);
198  break;
199  case 4: // native
200  const Real c =
201  (fd == friction) ? std::cos(_mc_phi.value(intnl)) : std::cos(_mc_psi.value(intnl));
202  const Real dc = (fd == friction)
203  ? -std::sin(_mc_phi.value(intnl)) * _mc_phi.derivative(intnl)
204  : -std::sin(_mc_psi.value(intnl)) * _mc_psi.derivative(intnl);
205  dbbb = ds / c - s * dc / Utility::pow<2>(c);
206  break;
207  }
208 }
209 
210 void
211 TensorMechanicsPlasticDruckerPrager::dbothAB(Real intnl, Real & daaa, Real & dbbb) const
212 {
214  {
215  daaa = 0;
216  dbbb = 0;
217  return;
218  }
219 
220  const Real C = _mc_cohesion.value(intnl);
221  const Real dC = _mc_cohesion.derivative(intnl);
222  const Real cosphi = std::cos(_mc_phi.value(intnl));
223  const Real dcosphi = -std::sin(_mc_phi.value(intnl)) * _mc_phi.derivative(intnl);
224  const Real sinphi = std::sin(_mc_phi.value(intnl));
225  const Real dsinphi = std::cos(_mc_phi.value(intnl)) * _mc_phi.derivative(intnl);
226  switch (_mc_interpolation_scheme)
227  {
228  case 0: // outer_tip
229  daaa = 2.0 * std::sqrt(3.0) * (dC * cosphi / (3.0 - sinphi) + C * dcosphi / (3.0 - sinphi) +
230  C * cosphi * dsinphi / Utility::pow<2>(3.0 - sinphi));
231  dbbb = 2.0 / std::sqrt(3.0) *
232  (dsinphi / (3.0 - sinphi) + sinphi * dsinphi / Utility::pow<2>(3.0 - sinphi));
233  break;
234  case 1: // inner_tip
235  daaa = 2.0 * std::sqrt(3.0) * (dC * cosphi / (3.0 + sinphi) + C * dcosphi / (3.0 + sinphi) -
236  C * cosphi * dsinphi / Utility::pow<2>(3.0 + sinphi));
237  dbbb = 2.0 / std::sqrt(3.0) *
238  (dsinphi / (3.0 + sinphi) - sinphi * dsinphi / Utility::pow<2>(3.0 + sinphi));
239  break;
240  case 2: // lode_zero
241  daaa = dC * cosphi + C * dcosphi;
242  dbbb = dsinphi / 3.0;
243  break;
244  case 3: // inner_edge
245  daaa = 3.0 * dC * cosphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi)) +
246  3.0 * C * dcosphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi)) -
247  3.0 * C * cosphi * 3.0 * sinphi * dsinphi /
248  std::pow(9.0 + 3.0 * Utility::pow<2>(sinphi), 1.5);
249  dbbb = dsinphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi)) -
250  3.0 * sinphi * sinphi * dsinphi / std::pow(9.0 + 3.0 * Utility::pow<2>(sinphi), 1.5);
251  break;
252  case 4: // native
253  daaa = dC;
254  dbbb = dsinphi / cosphi - sinphi * dcosphi / Utility::pow<2>(cosphi);
255  break;
256  }
257 }
258 
259 void
260 TensorMechanicsPlasticDruckerPrager::initializeAandB(Real intnl, Real & aaa, Real & bbb) const
261 {
262  const Real C = _mc_cohesion.value(intnl);
263  const Real cosphi = std::cos(_mc_phi.value(intnl));
264  const Real sinphi = std::sin(_mc_phi.value(intnl));
265  switch (_mc_interpolation_scheme)
266  {
267  case 0: // outer_tip
268  aaa = 2.0 * std::sqrt(3.0) * C * cosphi / (3.0 - sinphi);
269  bbb = 2.0 * sinphi / std::sqrt(3.0) / (3.0 - sinphi);
270  break;
271  case 1: // inner_tip
272  aaa = 2.0 * std::sqrt(3.0) * C * cosphi / (3.0 + sinphi);
273  bbb = 2.0 * sinphi / std::sqrt(3.0) / (3.0 + sinphi);
274  break;
275  case 2: // lode_zero
276  aaa = C * cosphi;
277  bbb = sinphi / 3.0;
278  break;
279  case 3: // inner_edge
280  aaa = 3.0 * C * cosphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi));
281  bbb = sinphi / std::sqrt(9.0 + 3.0 * Utility::pow<2>(sinphi));
282  break;
283  case 4: // native
284  aaa = C;
285  bbb = sinphi / cosphi;
286  break;
287  }
288 }
289 
290 void
291 TensorMechanicsPlasticDruckerPrager::initializeB(Real intnl, int fd, Real & bbb) const
292 {
293  const Real s = (fd == friction) ? std::sin(_mc_phi.value(intnl)) : std::sin(_mc_psi.value(intnl));
294  switch (_mc_interpolation_scheme)
295  {
296  case 0: // outer_tip
297  bbb = 2.0 * s / std::sqrt(3.0) / (3.0 - s);
298  break;
299  case 1: // inner_tip
300  bbb = 2.0 * s / std::sqrt(3.0) / (3.0 + s);
301  break;
302  case 2: // lode_zero
303  bbb = s / 3.0;
304  break;
305  case 3: // inner_edge
306  bbb = s / std::sqrt(9.0 + 3.0 * Utility::pow<2>(s));
307  break;
308  case 4: // native
309  const Real c =
310  (fd == friction) ? std::cos(_mc_phi.value(intnl)) : std::cos(_mc_psi.value(intnl));
311  bbb = s / c;
312  break;
313  }
314 }
virtual Real dyieldFunction_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to the internal parameter.
virtual Real yieldFunction(const RankTwoTensor &stress, Real intnl) const override
The following functions are what you should override when building single-plasticity models...
const bool _zero_psi_hardening
True if there is no hardening of dilation angle.
virtual RankTwoTensor df_dsig(const RankTwoTensor &stress, Real bbb) const
Function that&#39;s used in dyieldFunction_dstress and flowPotential.
virtual Real derivative(Real intnl) const
TensorMechanicsPlasticDruckerPrager(const InputParameters &parameters)
InputParameters validParams< TensorMechanicsPlasticDruckerPrager >()
virtual std::string modelName() const override
const TensorMechanicsHardeningModel & _mc_phi
Hardening model for tan(phi)
void bothAB(Real intnl, Real &aaa, Real &bbb) const
Calculates aaa and bbb as a function of the internal parameter intnl.
virtual RankTwoTensor dyieldFunction_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to stress.
const TensorMechanicsHardeningModel & _mc_psi
Hardening model for tan(psi)
const bool _zero_phi_hardening
True if there is no hardening of friction angle.
void donlyB(Real intnl, int fd, Real &dbbb) const
Calculate d(bbb)/d(intnl) or d(bbb_flow)/d(intnl)
InputParameters validParams< TensorMechanicsPlasticModel >()
const bool _zero_cohesion_hardening
True if there is no hardening of cohesion.
void dbothAB(Real intnl, Real &daaa, Real &dbbb) const
Calculates d(aaa)/d(intnl) and d(bbb)/d(intnl) as a function of the internal parameter intnl...
virtual Real value(Real intnl) const
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
void initializeB(Real intnl, int fd, Real &bbb) const
Returns the Drucker-Prager parameters A nice reference on the different relationships between Drucker...
void onlyB(Real intnl, int fd, Real &bbb) const
Calculate bbb or bbb_flow.
const TensorMechanicsHardeningModel & _mc_cohesion
Hardening model for cohesion.
virtual RankTwoTensor flowPotential(const RankTwoTensor &stress, Real intnl) const override
The flow potential.
virtual RankFourTensor dflowPotential_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to stress.
Plastic Model base class The virtual functions written below must be over-ridden in derived classes t...
const MooseEnum _mc_interpolation_scheme
The parameters aaa and bbb are chosen to closely match the Mohr-Coulomb yield surface.
void initializeAandB(Real intnl, Real &aaa, Real &bbb) const
Returns the Drucker-Prager parameters A nice reference on the different relationships between Drucker...
virtual RankTwoTensor dflowPotential_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to the internal parameter.