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