www.mooseframework.org
TensileStressUpdate.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 /****************************************************************/
7 #include "TensileStressUpdate.h"
8 #include "libmesh/utility.h"
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<MultiParameterPlasticityStressUpdate>();
15  params.addRequiredParam<UserObjectName>(
16  "tensile_strength",
17  "A TensorMechanicsHardening UserObject that defines hardening of the tensile strength");
18  params.addParam<bool>("perfect_guess",
19  true,
20  "Provide a guess to the Newton-Raphson proceedure "
21  "that is the result from perfect plasticity. With "
22  "severe hardening/softening this may be "
23  "suboptimal.");
24  params.addClassDescription(
25  "Associative, smoothed, tensile (Rankine) plasticity with hardening/softening");
26  return params;
27 }
28 
29 TensileStressUpdate::TensileStressUpdate(const InputParameters & parameters)
30  : MultiParameterPlasticityStressUpdate(parameters, 3, 3, 1),
31  _strength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
32  _perfect_guess(getParam<bool>("perfect_guess")),
33  _eigvecs(RankTwoTensor())
34 {
35 }
36 
37 void
38 TensileStressUpdate::computeStressParams(const RankTwoTensor & stress,
39  std::vector<Real> & stress_params) const
40 {
41  // stress_params[0] = smallest eigenvalue, stress_params[2] = largest eigenvalue
42  stress.symmetricEigenvalues(stress_params);
43 }
44 
45 std::vector<RankTwoTensor>
46 TensileStressUpdate::dstress_param_dstress(const RankTwoTensor & stress) const
47 {
48  std::vector<Real> sp;
49  std::vector<RankTwoTensor> dsp;
50  stress.dsymmetricEigenvalues(sp, dsp);
51  return dsp;
52 }
53 
54 std::vector<RankFourTensor>
55 TensileStressUpdate::d2stress_param_dstress(const RankTwoTensor & stress) const
56 {
57  std::vector<RankFourTensor> d2;
58  stress.d2symmetricEigenvalues(d2);
59  return d2;
60 }
61 
62 void
63 TensileStressUpdate::preReturnMapV(const std::vector<Real> & /*trial_stress_params*/,
64  const RankTwoTensor & stress_trial,
65  const std::vector<Real> & /*intnl_old*/,
66  const std::vector<Real> & /*yf*/,
67  const RankFourTensor & /*Eijkl*/)
68 {
69  std::vector<Real> eigvals;
70  stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs);
71 }
72 
73 void
74 TensileStressUpdate::setStressAfterReturnV(const RankTwoTensor & /*stress_trial*/,
75  const std::vector<Real> & stress_params,
76  Real /*gaE*/,
77  const std::vector<Real> & /*intnl*/,
78  const yieldAndFlow & /*smoothed_q*/,
79  const RankFourTensor & /*Eijkl*/,
80  RankTwoTensor & stress) const
81 {
82  // form the diagonal stress
83  stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
84  // rotate to the original frame
85  stress = _eigvecs * stress * (_eigvecs.transpose());
86 }
87 
88 void
89 TensileStressUpdate::yieldFunctionValuesV(const std::vector<Real> & stress_params,
90  const std::vector<Real> & intnl,
91  std::vector<Real> & yf) const
92 {
93  const Real ts = tensile_strength(intnl[0]);
94  // The smoothing strategy means that the last yield function
95  // gives the smallest value when returning to a line/point where,
96  // without smoothing, the yield functions are equal.
97  // Therefore, i use the smallest eigenvalue in this yield function
98  yf[0] = stress_params[2] - ts; // use largest eigenvalue
99  yf[1] = stress_params[1] - ts;
100  yf[2] = stress_params[0] - ts; // use smallest eigenvalue
101 }
102 
103 void
104 TensileStressUpdate::computeAllQV(const std::vector<Real> & stress_params,
105  const std::vector<Real> & intnl,
106  std::vector<yieldAndFlow> & all_q) const
107 {
108  const Real ts = tensile_strength(intnl[0]);
109  const Real dts = dtensile_strength(intnl[0]);
110 
111  // yield functions. See comment in yieldFunctionValuesV
112  all_q[0].f = stress_params[2] - ts;
113  all_q[1].f = stress_params[1] - ts;
114  all_q[2].f = stress_params[0] - ts;
115 
116  // d(yield function)/d(stress_params)
117  for (unsigned yf = 0; yf < _num_yf; ++yf)
118  for (unsigned a = 0; a < _num_sp; ++a)
119  all_q[yf].df[a] = 0.0;
120  all_q[0].df[2] = 1.0;
121  all_q[1].df[1] = 1.0;
122  all_q[2].df[0] = 1.0;
123 
124  // d(yield function)/d(intnl)
125  for (unsigned yf = 0; yf < _num_yf; ++yf)
126  all_q[yf].df_di[0] = -dts;
127 
128  // the flow potential is just the yield function
129  // d(flow potential)/d(stress_params)
130  for (unsigned yf = 0; yf < _num_yf; ++yf)
131  for (unsigned a = 0; a < _num_sp; ++a)
132  all_q[yf].dg[a] = all_q[yf].df[a];
133 
134  // d(flow potential)/d(stress_params)/d(intnl)
135  for (unsigned yf = 0; yf < _num_yf; ++yf)
136  for (unsigned a = 0; a < _num_sp; ++a)
137  all_q[yf].d2g_di[a][0] = 0.0;
138 
139  // d(flow potential)/d(stress_params)/d(stress_params)
140  for (unsigned yf = 0; yf < _num_yf; ++yf)
141  for (unsigned a = 0; a < _num_sp; ++a)
142  for (unsigned b = 0; b < _num_sp; ++b)
143  all_q[yf].d2g[a][b] = 0.0;
144 }
145 
146 void
147 TensileStressUpdate::setEffectiveElasticity(const RankFourTensor & Eijkl)
148 {
149  // Eijkl is required to be isotropic, so we can use the
150  // frame where stress is diagonal
151  for (unsigned a = 0; a < _num_sp; ++a)
152  for (unsigned b = 0; b < _num_sp; ++b)
153  _Eij[a][b] = Eijkl(a, a, b, b);
154  _En = _Eij[2][2];
155  const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]);
156  for (unsigned a = 0; a < _num_sp; ++a)
157  {
158  _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom;
159  for (unsigned b = 0; b < a; ++b)
160  _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom;
161  }
162 }
163 
164 void
165 TensileStressUpdate::initializeVarsV(const std::vector<Real> & trial_stress_params,
166  const std::vector<Real> & intnl_old,
167  std::vector<Real> & stress_params,
168  Real & gaE,
169  std::vector<Real> & intnl) const
170 {
171  if (!_perfect_guess)
172  {
173  for (unsigned i = 0; i < _num_sp; ++i)
174  stress_params[i] = trial_stress_params[i];
175  gaE = 0.0;
176  }
177  else
178  {
179  const Real ts = tensile_strength(intnl_old[0]);
180  stress_params[2] = ts; // largest eigenvalue
181  stress_params[1] = std::min(stress_params[1], ts);
182  stress_params[0] = std::min(stress_params[0], ts);
183  gaE = trial_stress_params[2] - stress_params[2];
184  }
185  setIntnlValuesV(trial_stress_params, stress_params, intnl_old, intnl);
186 }
187 
188 void
189 TensileStressUpdate::setIntnlValuesV(const std::vector<Real> & trial_stress_params,
190  const std::vector<Real> & current_stress_params,
191  const std::vector<Real> & intnl_old,
192  std::vector<Real> & intnl) const
193 {
194  intnl[0] = intnl_old[0] + (trial_stress_params[2] - current_stress_params[2]) / _Eij[2][2];
195 }
196 
197 void
198 TensileStressUpdate::setIntnlDerivativesV(const std::vector<Real> & /*trial_stress_params*/,
199  const std::vector<Real> & /*current_stress_params*/,
200  const std::vector<Real> & /*intnl*/,
201  std::vector<std::vector<Real>> & dintnl) const
202 {
203  dintnl[0][0] = 0.0;
204  dintnl[0][1] = 0.0;
205  dintnl[0][2] = -1.0 / _Eij[2][2];
206 }
207 
208 Real
209 TensileStressUpdate::tensile_strength(const Real internal_param) const
210 {
211  return _strength.value(internal_param);
212 }
213 
214 Real
215 TensileStressUpdate::dtensile_strength(const Real internal_param) const
216 {
217  return _strength.derivative(internal_param);
218 }
219 
220 void
221 TensileStressUpdate::consistentTangentOperatorV(const RankTwoTensor & stress_trial,
222  const std::vector<Real> & trial_stress_params,
223  const RankTwoTensor & /*stress*/,
224  const std::vector<Real> & stress_params,
225  Real /*gaE*/,
226  const yieldAndFlow & /*smoothed_q*/,
227  const RankFourTensor & elasticity_tensor,
228  bool compute_full_tangent_operator,
229  const std::vector<std::vector<Real>> & dvar_dtrial,
230  RankFourTensor & cto)
231 {
232  cto = elasticity_tensor;
233  if (!compute_full_tangent_operator)
234  return;
235 
236  // dvar_dtrial has been computed already, so
237  // d(stress)/d(trial_stress) = d(eigvecs * stress_params * eigvecs.transpose())/d(trial_stress)
238  // eigvecs is a rotation matrix, rot(i, j) = e_j(i) = i^th component of j^th eigenvector
239  // d(rot_ij)/d(stress_kl) = d(e_j(i))/d(stress_kl)
240  // = sum_a 0.5 * e_a(i) * (e_a(k)e_j(l) + e_a(l)e_j(k)) / (la_j - la_a)
241  // = sum_a 0.5 * rot(i,a) * (rot(k,a)rot(l,j) + rot(l,a)*rot(k,j)) / (la_j - la_a)
242  RankFourTensor drot_dstress;
243  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
244  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
245  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
246  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
247  for (unsigned a = 0; a < _num_sp; ++a)
248  {
249  if (trial_stress_params[a] == trial_stress_params[j])
250  continue;
251  drot_dstress(i, j, k, l) += 0.5 * _eigvecs(i, a) * (_eigvecs(k, a) * _eigvecs(l, j) +
252  _eigvecs(l, a) * _eigvecs(k, j)) /
253  (trial_stress_params[j] - trial_stress_params[a]);
254  }
255 
256  const RankTwoTensor eT = _eigvecs.transpose();
257 
258  RankFourTensor dstress_dtrial;
259  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
260  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
261  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
262  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
263  for (unsigned a = 0; a < _num_sp; ++a)
264  dstress_dtrial(i, j, k, l) +=
265  drot_dstress(i, a, k, l) * stress_params[a] * eT(a, j) +
266  _eigvecs(i, a) * stress_params[a] * drot_dstress(j, a, k, l);
267 
268  const std::vector<RankTwoTensor> dsp_trial = dstress_param_dstress(stress_trial);
269  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
270  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
271  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
272  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
273  for (unsigned a = 0; a < _num_sp; ++a)
274  for (unsigned b = 0; b < _num_sp; ++b)
275  dstress_dtrial(i, j, k, l) +=
276  _eigvecs(i, a) * dvar_dtrial[a][b] * dsp_trial[b](k, l) * eT(a, j);
277 
278  cto = dstress_dtrial * elasticity_tensor;
279 }
TensileStressUpdate(const InputParameters &parameters)
virtual void preReturnMapV(const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress_trial, const std::vector< Real > &intnl_old, const std::vector< Real > &yf, const RankFourTensor &Eijkl) override
Derived classes may employ this function to record stuff or do other computations prior to the return...
virtual void consistentTangentOperatorV(const RankTwoTensor &stress_trial, const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, const std::vector< std::vector< Real >> &dvar_dtrial, RankFourTensor &cto) override
Calculates the consistent tangent operator.
InputParameters validParams< TensileStressUpdate >()
InputParameters validParams< MultiParameterPlasticityStressUpdate >()
std::vector< RankFourTensor > d2stress_param_dstress(const RankTwoTensor &stress) const override
d2(stress_param[i])/d(stress)/d(stress) at given stress
virtual Real derivative(Real intnl) const
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
std::vector< std::vector< Real > > _Cij
_Cij[i, j] * _Eij[j, k] = 1 iff j == k
void setEffectiveElasticity(const RankFourTensor &Eijkl) override
Sets _Eij and _En and _Cij.
RankTwoTensor _eigvecs
Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to s...
void setIntnlValuesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const override
Sets the internal parameters based on the trial values of stress_params, their current values...
const TensorMechanicsHardeningModel & _strength
Hardening model for tensile strength.
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
const unsigned _num_yf
Number of yield functions.
void yieldFunctionValuesV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< Real > &yf) const override
Computes the values of the yield functions, given stress_params and intnl parameters.
std::vector< RankTwoTensor > dstress_param_dstress(const RankTwoTensor &stress) const override
d(stress_param[i])/d(stress) at given stress
void setIntnlDerivativesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl, std::vector< std::vector< Real >> &dintnl) const override
Sets the derivatives of internal parameters, based on the trial values of stress_params, their current values, and the current values of the internal parameters.
virtual Real value(Real intnl) const
void computeStressParams(const RankTwoTensor &stress, std::vector< Real > &stress_params) const override
Computes stress_params, given stress.
void initializeVarsV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &stress_params, Real &gaE, std::vector< Real > &intnl) const override
Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm.
virtual Real dtensile_strength(const Real internal_param) const
d(tensile strength)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual void setStressAfterReturnV(const RankTwoTensor &stress_trial, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const override
Sets stress from the admissible parameters.
virtual Real tensile_strength(const Real internal_param) const
tensile strength as a function of residual value, rate, and internal_param
const bool _perfect_guess
Whether to provide an estimate of the returned stress, based on perfect plasticity.
MultiParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates ...
void computeAllQV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< yieldAndFlow > &all_q) const override
Completely fills all_q with correct values.
static constexpr unsigned _tensor_dimensionality
Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) ...
const unsigned _num_sp
Number of stress parameters.