www.mooseframework.org
CappedWeakPlaneCosseratStressUpdate.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 #include "libmesh/utility.h"
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<CappedWeakPlaneStressUpdate>();
16  params.addClassDescription("Capped weak-plane plasticity Cosserat stress calculator");
17  return params;
18 }
19 
21  const InputParameters & parameters)
22  : CappedWeakPlaneStressUpdate(parameters)
23 {
24 }
25 
26 void
28  Real p_ok,
29  Real q_ok,
30  Real gaE,
31  const std::vector<Real> & /*intnl*/,
32  const yieldAndFlow & smoothed_q,
33  const RankFourTensor & Eijkl,
34  RankTwoTensor & stress) const
35 {
36  stress = stress_trial;
37  stress(2, 2) = p_ok;
38  // stress_xx and stress_yy are sitting at their trial-stress values
39  // so need to bring them back via Poisson's ratio
40  stress(0, 0) -= Eijkl(2, 2, 0, 0) * gaE / _Epp * smoothed_q.dg[0];
41  stress(1, 1) -= Eijkl(2, 2, 1, 1) * gaE / _Epp * smoothed_q.dg[0];
42  if (_in_q_trial == 0.0)
43  stress(0, 2) = stress(1, 2) = 0.0;
44  else
45  {
46  stress(0, 2) = _in_trial02 * q_ok / _in_q_trial;
47  stress(1, 2) = _in_trial12 * q_ok / _in_q_trial;
48  }
49 }
50 
51 void
53  const RankTwoTensor & /*stress_trial*/,
54  Real /*p_trial*/,
55  Real q_trial,
56  const RankTwoTensor & /*stress*/,
57  Real /*p*/,
58  Real q,
59  Real gaE,
60  const yieldAndFlow & smoothed_q,
61  const RankFourTensor & Eijkl,
62  bool compute_full_tangent_operator,
63  RankFourTensor & cto) const
64 {
65  cto = Eijkl;
66  if (!compute_full_tangent_operator)
67  return;
68 
69  const Real Ezzzz = Eijkl(2, 2, 2, 2);
70  const Real Exzxz = Eijkl(0, 2, 0, 2);
71  const Real tanpsi = _tan_psi.value(_intnl[_qp][0]);
72  const Real dintnl0_dq = -1.0 / Exzxz;
73  const Real dintnl0_dqt = 1.0 / Exzxz;
74  const Real dintnl1_dp = -1.0 / Ezzzz;
75  const Real dintnl1_dpt = 1.0 / Ezzzz;
76  const Real dintnl1_dq =
77  tanpsi / Exzxz - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dq / Exzxz;
78  const Real dintnl1_dqt =
79  -tanpsi / Exzxz - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dqt / Exzxz;
80 
81  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
82  {
83  const Real dpt_depii = Eijkl(2, 2, i, i);
84  cto(2, 2, i, i) = _dp_dpt * dpt_depii;
85  const Real poisson_effect =
86  Eijkl(2, 2, 0, 0) / Ezzzz *
87  (_dgaE_dpt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dpt +
88  gaE * smoothed_q.d2g[0][1] * _dq_dpt +
89  gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dq * _dq_dpt) +
90  gaE * smoothed_q.d2g_di[0][1] *
91  (dintnl1_dpt + dintnl1_dp * _dp_dpt + dintnl1_dq * _dq_dpt)) *
92  dpt_depii;
93  cto(0, 0, i, i) -= poisson_effect;
94  cto(1, 1, i, i) -= poisson_effect;
95  if (q_trial > 0.0)
96  {
97  cto(0, 2, i, i) = _in_trial02 / q_trial * _dq_dpt * dpt_depii;
98  cto(1, 2, i, i) = _in_trial12 / q_trial * _dq_dpt * dpt_depii;
99  }
100  }
101 
102  const Real poisson_effect =
103  -Eijkl(2, 2, 0, 0) / Ezzzz *
104  (_dgaE_dqt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dqt +
105  gaE * smoothed_q.d2g[0][1] * _dq_dqt +
106  gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dqt + dintnl0_dq * _dq_dqt) +
107  gaE * smoothed_q.d2g_di[0][1] * (dintnl1_dqt + dintnl1_dp * _dp_dqt + dintnl1_dq * _dq_dqt));
108 
109  const Real dqt_dep02 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(0, 2, 0, 2);
110  cto(2, 2, 0, 2) = _dp_dqt * dqt_dep02;
111  cto(0, 0, 0, 2) = cto(1, 1, 0, 2) = poisson_effect * dqt_dep02;
112  if (q_trial > 0.0)
113  {
114  // for q_trial=0, Jacobian_mult is just given by the elastic case
115  cto(0, 2, 0, 2) = Eijkl(0, 2, 0, 2) * q / q_trial +
116  _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep02;
117  cto(1, 2, 0, 2) = _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep02;
118  }
119 
120  const Real dqt_dep20 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(0, 2, 2, 0);
121  cto(2, 2, 2, 0) = _dp_dqt * dqt_dep20;
122  cto(0, 0, 2, 0) = cto(1, 1, 2, 0) = poisson_effect * dqt_dep20;
123  if (q_trial > 0.0)
124  {
125  // for q_trial=0, Jacobian_mult is just given by the elastic case
126  cto(0, 2, 2, 0) = Eijkl(0, 2, 2, 0) * q / q_trial +
127  _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
128  cto(1, 2, 2, 0) = _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
129  }
130 
131  const Real dqt_dep12 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(1, 2, 1, 2);
132  cto(2, 2, 1, 2) = _dp_dqt * dqt_dep12;
133  cto(0, 0, 1, 2) = cto(1, 1, 1, 2) = poisson_effect * dqt_dep12;
134  if (q_trial > 0.0)
135  {
136  // for q_trial=0, Jacobian_mult is just given by the elastic case
137  cto(0, 2, 1, 2) = _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep12;
138  cto(1, 2, 1, 2) = Eijkl(1, 2, 1, 2) * q / q_trial +
139  _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep12;
140  }
141 
142  const Real dqt_dep21 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(1, 2, 2, 1);
143  cto(2, 2, 2, 1) = _dp_dqt * dqt_dep21;
144  cto(0, 0, 2, 1) = cto(1, 1, 2, 1) = poisson_effect * dqt_dep21;
145  if (q_trial > 0.0)
146  {
147  // for q_trial=0, Jacobian_mult is just given by the elastic case
148  cto(0, 2, 2, 1) = _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
149  cto(1, 2, 2, 1) = Eijkl(1, 2, 2, 1) * q / q_trial +
150  _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
151  }
152 }
153 
154 RankTwoTensor
155 CappedWeakPlaneCosseratStressUpdate::dqdstress(const RankTwoTensor & stress) const
156 {
157  RankTwoTensor deriv = RankTwoTensor();
158  const Real q = std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
159  if (q > 0.0)
160  {
161  deriv(0, 2) = stress(0, 2) / q;
162  deriv(1, 2) = stress(1, 2) / q;
163  }
164  else
165  {
166  // derivative is not defined here. For now i'll set:
167  deriv(0, 2) = 1.0;
168  deriv(1, 2) = 1.0;
169  }
170  return deriv;
171 }
172 
173 RankFourTensor
174 CappedWeakPlaneCosseratStressUpdate::d2qdstress2(const RankTwoTensor & stress) const
175 {
176  RankFourTensor d2 = RankFourTensor();
177 
178  const Real q = std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
179  if (q == 0.0)
180  return d2;
181 
182  const Real dq02 = stress(0, 2) / q;
183  const Real dq12 = stress(1, 2) / q;
184 
185  d2(0, 2, 0, 2) = 1.0 / q - dq02 * dq02 / q;
186  d2(0, 2, 1, 2) = -dq02 * dq12 / q;
187  d2(1, 2, 0, 2) = -dq12 * dq02 / q;
188  d2(1, 2, 1, 2) = 1.0 / q - dq12 * dq12 / q;
189 
190  return d2;
191 }
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.
InputParameters validParams< CappedWeakPlaneStressUpdate >()
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 Real derivative(Real intnl) const
InputParameters validParams< CappedWeakPlaneCosseratStressUpdate >()
virtual RankFourTensor d2qdstress2(const RankTwoTensor &stress) const override
d2(q)/d(stress)/d(stress) Derived classes must override this
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
Real _dgaE_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Real _dgaE_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
Real _in_trial02
trial value of stress(0, 2)
Real _dq_dqt
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.
const TensorMechanicsHardeningModel & _tan_psi
Hardening model for tan(psi)
Real _Epp
elasticity tensor in p direction
virtual Real value(Real intnl) const
CappedWeakPlaneCosseratStressUpdate(const InputParameters &parameters)
CappedWeakPlaneStressUpdate performs the return-map algorithm and associated stress updates for plast...
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) ...
virtual RankTwoTensor dqdstress(const RankTwoTensor &stress) const override
d(q)/d(stress) Derived classes must override this
Real _in_trial12
trial value of stress(1, 2)