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