www.mooseframework.org
CappedDruckerPragerCosseratStressUpdate.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 
9 template <>
10 InputParameters
12 {
13  InputParameters params = validParams<CappedDruckerPragerStressUpdate>();
14  params.addClassDescription("Capped Drucker-Prager plasticity stress calculator for the Cosserat "
15  "situation where the host medium (ie, the limit where all Cosserat "
16  "effects are zero) is isotropic. Note that the return-map flow rule "
17  "uses an isotropic elasticity tensor built with the 'host' properties "
18  "defined by the user.");
19  params.addRequiredRangeCheckedParam<Real>("host_youngs_modulus",
20  "host_youngs_modulus>0",
21  "Young's modulus for the isotropic host medium");
22  params.addRequiredRangeCheckedParam<Real>("host_poissons_ratio",
23  "host_poissons_ratio>=0 & host_poissons_ratio<0.5",
24  "Poisson's ratio for the isotropic host medium");
25  return params;
26 }
27 
29  const InputParameters & parameters)
30  : CappedDruckerPragerStressUpdate(parameters),
31  _shear(getParam<Real>("host_youngs_modulus") /
32  (2.0 * (1.0 + getParam<Real>("host_poissons_ratio"))))
33 {
34  const Real young = getParam<Real>("host_youngs_modulus");
35  const Real poisson = getParam<Real>("host_poissons_ratio");
36  const Real lambda = young * poisson / ((1.0 + poisson) * (1.0 - 2.0 * poisson));
37  _Ehost.fillFromInputVector({lambda, _shear}, RankFourTensor::symmetric_isotropic);
38 }
39 
40 void
41 CappedDruckerPragerCosseratStressUpdate::setEppEqq(const RankFourTensor & /*Eijkl*/,
42  Real & Epp,
43  Real & Eqq) const
44 {
45  Epp = _Ehost.sum3x3();
46  Eqq = _shear;
47 }
48 
49 void
51  Real p_ok,
52  Real q_ok,
53  Real /*gaE*/,
54  const std::vector<Real> & /*intnl*/,
55  const yieldAndFlow & /*smoothed_q*/,
56  const RankFourTensor & /*Eijkl*/,
57  RankTwoTensor & stress) const
58 {
59  // symm_stress is the symmetric part of the stress tensor.
60  // symm_stress = (s_ij+s_ji)/2 + de_ij tr(stress) / 3
61  // = q / q_trial * (s_ij^trial+s_ji^trial)/2 + de_ij p / 3
62  // = q / q_trial * (symm_stress_ij^trial - de_ij tr(stress^trial) / 3) + de_ij p / 3
63  const Real p_trial = stress_trial.trace();
64  RankTwoTensor symm_stress = RankTwoTensor(RankTwoTensor::initIdentity) / 3.0 *
65  (p_ok - (_in_q_trial == 0.0 ? 0.0 : p_trial * q_ok / _in_q_trial));
66  if (_in_q_trial > 0)
67  symm_stress += q_ok / _in_q_trial * 0.5 * (stress_trial + stress_trial.transpose());
68  stress = symm_stress + 0.5 * (stress_trial - stress_trial.transpose());
69 }
70 
71 void
73  const RankTwoTensor & /*stress_trial*/,
74  Real /*p_trial*/,
75  Real /*q_trial*/,
76  const RankTwoTensor & stress,
77  Real /*p*/,
78  Real q,
79  Real gaE,
80  const yieldAndFlow & smoothed_q,
81  const RankFourTensor & Eijkl,
82  bool compute_full_tangent_operator,
83  RankFourTensor & cto) const
84 {
85  if (!compute_full_tangent_operator)
86  {
87  cto = Eijkl;
88  return;
89  }
90 
91  RankFourTensor EAijkl;
92  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
93  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
94  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
95  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
96  {
97  cto(i, j, k, l) = 0.5 * (Eijkl(i, j, k, l) + Eijkl(j, i, k, l));
98  EAijkl(i, j, k, l) = 0.5 * (Eijkl(i, j, k, l) - Eijkl(j, i, k, l));
99  }
100 
101  const RankTwoTensor s_over_q =
102  (q == 0.0 ? RankTwoTensor()
103  : (0.5 * (stress + stress.transpose()) -
104  stress.trace() * RankTwoTensor(RankTwoTensor::initIdentity) / 3.0) /
105  q);
106  const RankTwoTensor E_s_over_q = Eijkl.innerProductTranspose(s_over_q); // not symmetric in kl
107  const RankTwoTensor Ekl =
108  RankTwoTensor(RankTwoTensor::initIdentity).initialContraction(Eijkl); // symmetric in kl
109 
110  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
111  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
112  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
113  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
114  {
115  cto(i, j, k, l) -= (i == j) * (1.0 / 3.0) *
116  (Ekl(k, l) * (1.0 - _dp_dpt) + 0.5 * E_s_over_q(k, l) * (-_dp_dqt));
117  cto(i, j, k, l) -=
118  s_over_q(i, j) * (Ekl(k, l) * (-_dq_dpt) + 0.5 * E_s_over_q(k, l) * (1.0 - _dq_dqt));
119  }
120 
121  if (smoothed_q.dg[1] != 0.0)
122  {
123  const RankFourTensor Tijab = _Ehost * (gaE / _Epp) * smoothed_q.dg[1] * d2qdstress2(stress);
124  RankFourTensor inv = RankFourTensor(RankFourTensor::initIdentitySymmetricFour) + Tijab;
125  try
126  {
127  inv = inv.transposeMajor().invSymm();
128  }
129  catch (const MooseException & e)
130  {
131  // Cannot form the inverse, so probably at some degenerate place in stress space.
132  // Just return with the "best estimate" of the cto.
133  mooseWarning("CappedDruckerPragerCosseratStressUpdate: Cannot invert 1+T in consistent "
134  "tangent operator computation at quadpoint ",
135  _qp,
136  " of element ",
137  _current_elem->id());
138  return;
139  }
140  cto = (cto.transposeMajor() * inv).transposeMajor();
141  }
142  cto += EAijkl;
143 }
const Real _shear
Shear modulus for the host medium.
CappedDruckerPragerCosseratStressUpdate(const InputParameters &parameters)
Real _dp_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Real _dq_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
virtual void consistentTangentOperator(const RankTwoTensor &stress_trial, Real p_trial, Real q_trial, const RankTwoTensor &stress, Real p, Real q, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, RankFourTensor &cto) const override
Calculates the consistent tangent operator.
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
InputParameters validParams< CappedDruckerPragerCosseratStressUpdate >()
InputParameters validParams< CappedDruckerPragerStressUpdate >()
RankFourTensor _Ehost
Isotropic elasticity tensor for the host medium.
virtual void setStressAfterReturn(const RankTwoTensor &stress_trial, Real p_ok, Real q_ok, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const override
Sets stress from the admissible parameters.
Real _dq_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Real _Epp
elasticity tensor in p direction
virtual RankFourTensor d2qdstress2(const RankTwoTensor &stress) const override
d2(q)/d(stress)/d(stress) Derived classes must override this
CappedDruckerPragerStressUpdate performs the return-map algorithm and associated stress updates for p...
virtual void setEppEqq(const RankFourTensor &Eijkl, Real &Epp, Real &Eqq) const override
Set Epp and Eqq based on the elasticity tensor Derived classes must override this.
Real _dp_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
static constexpr unsigned _tensor_dimensionality
Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) ...