www.mooseframework.org
PLC_LSH.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 "PLC_LSH.h"
8 
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<SolidModel>();
16 
17  // Power-law creep material parameters
18  params.addRequiredParam<Real>("coefficient", "Leading coefficent in power-law equation");
19  params.addRequiredParam<Real>("n_exponent", "Exponent on effective stress in power-law equation");
20  params.addParam<Real>("m_exponent", 0.0, "Exponent on time in power-law equation");
21  params.addRequiredParam<Real>("activation_energy", "Activation energy");
22  params.addParam<Real>("gas_constant", 8.3143, "Universal gas constant");
23 
24  // Linear strain hardening parameters
25  params.addRequiredParam<Real>("yield_stress",
26  "The point at which plastic strain begins accumulating");
27  params.addRequiredParam<Real>("hardening_constant", "Hardening slope");
28 
29  // Sub-Newton Iteration control parameters
30  params.addParam<unsigned int>("max_its", 30, "Maximum number of sub-newton iterations");
31  params.addParam<bool>(
32  "output_iteration_info", false, "Set true to output sub-newton iteration information");
33  params.addParam<Real>(
34  "relative_tolerance", 1e-5, "Relative convergence tolerance for sub-newtion iteration");
35  params.addParam<Real>(
36  "absolute_tolerance", 1e-20, "Absolute convergence tolerance for sub-newtion iteration");
37  params.addParam<PostprocessorName>(
38  "output", "", "The reporting postprocessor to use for the max_iterations value.");
39 
40  // Control of combined plasticity-creep iterarion
41  params.addParam<Real>("absolute_stress_tolerance",
42  1e-5,
43  "Convergence tolerance for combined plasticity-creep stress iteration");
44 
45  return params;
46 }
47 
48 PLC_LSH::PLC_LSH(const InputParameters & parameters)
49  : SolidModel(parameters),
50  _coefficient(parameters.get<Real>("coefficient")),
51  _n_exponent(parameters.get<Real>("n_exponent")),
52  _m_exponent(parameters.get<Real>("m_exponent")),
53  _activation_energy(parameters.get<Real>("activation_energy")),
54  _gas_constant(parameters.get<Real>("gas_constant")),
55 
56  _yield_stress(parameters.get<Real>("yield_stress")),
57  _hardening_constant(parameters.get<Real>("hardening_constant")),
58 
59  _max_its(parameters.get<unsigned int>("max_its")),
60  _output_iteration_info(getParam<bool>("output_iteration_info")),
61  _relative_tolerance(parameters.get<Real>("relative_tolerance")),
62  _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
63 
64  _absolute_stress_tolerance(parameters.get<Real>("absolute_stress_tolerance")),
65 
66  _creep_strain(declareProperty<SymmTensor>("creep_strain")),
67  _creep_strain_old(getMaterialPropertyOld<SymmTensor>("creep_strain")),
68 
69  _plastic_strain(declareProperty<SymmTensor>("plastic_strain")),
70  _plastic_strain_old(getMaterialPropertyOld<SymmTensor>("plastic_strain")),
71 
72  _hardening_variable(declareProperty<Real>("hardening_variable")),
73  _hardening_variable_old(getMaterialPropertyOld<Real>("hardening_variable")),
74 
75  _output(getParam<PostprocessorName>("output") != "" ? &getPostprocessorValue("output") : NULL)
76 
77 {
78  if (_yield_stress <= 0)
79  {
80  mooseError("Yield stress must be greater than zero");
81  }
82 }
83 
84 void
86 {
87  _hardening_variable[_qp] = 0;
89 }
90 
91 void
93 {
94  // Given the stretching, compute the stress increment and add it to the old stress. Also update
95  // the creep strain
96  // stress = stressOld + stressIncrement
97  // creep_strain = creep_strainOld + creep_strainIncrement
98 
99  if (_t_step == 0 && !_app.isRestarting())
100  return;
101 
102  if (_output_iteration_info == true)
103  {
104  _console << std::endl
105  << "iteration output for combined creep-plasticity solve:"
106  << " time=" << _t << " temperature=" << _temperature[_qp] << " int_pt=" << _qp
107  << std::endl;
108  }
109 
110  // compute trial stress
112  stress_new += _stress_old;
113 
114  SymmTensor creep_strain_increment;
115  SymmTensor plastic_strain_increment;
116  SymmTensor elastic_strain_increment;
117  SymmTensor stress_new_last(stress_new);
118  Real delS(_absolute_stress_tolerance + 1);
119  Real first_delS(delS);
120  unsigned int counter(0);
121 
122  while (counter < _max_its && delS > _absolute_stress_tolerance &&
123  (delS / first_delS) > _relative_tolerance)
124  {
125  elastic_strain_increment = _strain_increment;
126  elastic_strain_increment -= plastic_strain_increment;
127  stress_new = *elasticityTensor() * elastic_strain_increment;
128  stress_new += _stress_old;
129 
130  elastic_strain_increment = _strain_increment;
131  computeCreep(elastic_strain_increment, creep_strain_increment, stress_new);
132 
133  // now use stress_new to calculate a new effective_trial_stress and determine if
134  // yield has occured and if so, calculate the corresponding plastic strain
135 
136  elastic_strain_increment -= creep_strain_increment;
137 
138  computeLSH(elastic_strain_increment, plastic_strain_increment, stress_new);
139 
140  elastic_strain_increment -= plastic_strain_increment;
141 
142  // now check convergence
143  SymmTensor deltaS(stress_new_last - stress_new);
144  delS = std::sqrt(deltaS.doubleContraction(deltaS));
145  if (counter == 0)
146  {
147  first_delS = delS;
148  }
149  stress_new_last = stress_new;
150 
151  if (_output_iteration_info == true)
152  {
153  _console << "stress_it=" << counter << " rel_delS=" << delS / first_delS
154  << " rel_tol=" << _relative_tolerance << " abs_delS=" << delS
155  << " abs_tol=" << _absolute_stress_tolerance << std::endl;
156  }
157 
158  ++counter;
159  }
160 
161  if (counter == _max_its && delS > _absolute_stress_tolerance &&
162  (delS / first_delS) > _relative_tolerance)
163  {
164  mooseError("Max stress iteration hit during plasticity-creep solve!");
165  }
166 
167  _strain_increment = elastic_strain_increment;
168  _stress[_qp] = stress_new;
169 }
170 
171 void
172 PLC_LSH::computeCreep(const SymmTensor & strain_increment,
173  SymmTensor & creep_strain_increment,
174  SymmTensor & stress_new)
175 {
176  // compute deviatoric trial stress
177  SymmTensor dev_trial_stress(stress_new);
178  dev_trial_stress.addDiag(-dev_trial_stress.trace() / 3.0);
179 
180  // compute effective trial stress
181  Real dts_squared = dev_trial_stress.doubleContraction(dev_trial_stress);
182  Real effective_trial_stress = std::sqrt(1.5 * dts_squared);
183 
184  // Use Newton sub-iteration to determine effective creep strain increment
185 
186  Real exponential(1);
187  if (_has_temp)
188  {
189  exponential = std::exp(-_activation_energy / (_gas_constant * _temperature[_qp]));
190  }
191  Real expTime = std::pow(_t, _m_exponent);
192 
193  Real del_p = 0;
194  unsigned int it = 0;
195  Real creep_residual = 10;
196  Real norm_creep_residual = 10;
197  Real first_norm_creep_residual = 10;
198 
199  while (it < _max_its && norm_creep_residual > _absolute_tolerance &&
200  (norm_creep_residual / first_norm_creep_residual) > _relative_tolerance)
201  {
202 
203  Real phi = _coefficient *
204  std::pow(effective_trial_stress - 3 * _shear_modulus * del_p, _n_exponent) *
205  exponential * expTime;
206  Real dphi_ddelp =
208  std::pow(effective_trial_stress - 3 * _shear_modulus * del_p, _n_exponent - 1) *
209  exponential * expTime;
210 
211  creep_residual = phi - del_p / _dt;
212  norm_creep_residual = std::abs(creep_residual);
213  if (it == 0)
214  {
215  first_norm_creep_residual = norm_creep_residual;
216  }
217 
218  del_p += creep_residual / (1 / _dt - dphi_ddelp);
219 
220  if (_output_iteration_info == true)
221  {
222  _console << "crp_it=" << it << " trl_strs=" << effective_trial_stress << " phi=" << phi
223  << " dphi=" << dphi_ddelp << " del_p=" << del_p
224  << " rel_res=" << norm_creep_residual / first_norm_creep_residual
225  << " rel_tol=" << _relative_tolerance << " abs_res=" << norm_creep_residual
226  << " abs_tol=" << _absolute_tolerance << std::endl;
227  }
228 
229  ++it;
230  }
231 
232  if (it == _max_its && norm_creep_residual > _absolute_tolerance &&
233  (norm_creep_residual / first_norm_creep_residual) > _relative_tolerance)
234  {
235  mooseError("Max sub-newton iteration hit during creep solve!");
236  }
237 
238  // compute creep and elastic strain increments (avoid potential divide by zero - how should this
239  // be done)?
240  if (effective_trial_stress < 0.01)
241  {
242  effective_trial_stress = 0.01;
243  }
244 
245  creep_strain_increment = dev_trial_stress;
246  creep_strain_increment *= (1.5 * del_p / effective_trial_stress);
247 
248  SymmTensor elastic_strain_increment(strain_increment);
249  elastic_strain_increment -= creep_strain_increment;
250 
251  // compute stress increment
252  stress_new = *elasticityTensor() * elastic_strain_increment;
253 
254  // update stress and creep strain
255  stress_new += _stress_old;
256 
257  _creep_strain[_qp] = creep_strain_increment;
258  _creep_strain[_qp] += _creep_strain_old[_qp];
259 }
260 
261 void
262 PLC_LSH::computeLSH(const SymmTensor & strain_increment,
263  SymmTensor & plastic_strain_increment,
264  SymmTensor & stress_new)
265 {
266 
267  // compute deviatoric trial stress
268  SymmTensor dev_trial_stress(stress_new);
269  dev_trial_stress.addDiag(-stress_new.trace() / 3);
270 
271  // effective trial stress
272  Real dts_squared = dev_trial_stress.doubleContraction(dev_trial_stress);
273  Real effective_trial_stress = std::sqrt(1.5 * dts_squared);
274 
275  // determine if yield condition is satisfied
276  Real yield_condition = effective_trial_stress - _hardening_variable_old[_qp] - _yield_stress;
277 
280 
281  if (yield_condition >
282  0) // then use newton iteration to determine effective plastic strain increment
283  {
284  unsigned int it = 0;
285  Real plastic_residual = 0;
286  Real norm_plas_residual = 10;
287  Real first_norm_plas_residual = 10;
288  Real scalar_plastic_strain_increment = 0;
289 
290  while (it < _max_its && norm_plas_residual > _absolute_tolerance &&
291  (norm_plas_residual / first_norm_plas_residual) > _relative_tolerance)
292  {
293  plastic_residual = effective_trial_stress -
294  (3. * _shear_modulus * scalar_plastic_strain_increment) -
295  _hardening_variable[_qp] - _yield_stress;
296  norm_plas_residual = std::abs(plastic_residual);
297  if (it == 0)
298  {
299  first_norm_plas_residual = norm_plas_residual;
300  }
301 
302  scalar_plastic_strain_increment +=
303  plastic_residual / (3. * _shear_modulus + _hardening_constant);
304 
305  _hardening_variable[_qp] =
306  _hardening_variable_old[_qp] + (_hardening_constant * scalar_plastic_strain_increment);
307 
308  if (_output_iteration_info == true)
309  {
310  _console << "pls_it=" << it << " trl_strs=" << effective_trial_stress
311  << " del_p=" << scalar_plastic_strain_increment
312  << " harden=" << _hardening_variable[_qp]
313  << " rel_res=" << norm_plas_residual / first_norm_plas_residual
314  << " rel_tol=" << _relative_tolerance << " abs_res=" << norm_plas_residual
315  << " abs_tol=" << _absolute_tolerance << std::endl;
316  }
317 
318  ++it;
319  }
320 
321  if (it == _max_its && norm_plas_residual > _absolute_tolerance &&
322  (norm_plas_residual / first_norm_plas_residual) > _relative_tolerance)
323  {
324  mooseError("Max sub-newton iteration hit during plasticity increment solve!");
325  }
326 
327  if (effective_trial_stress < 0.01)
328  {
329  effective_trial_stress = 0.01;
330  }
331  plastic_strain_increment = dev_trial_stress;
332  plastic_strain_increment *= (1.5 * scalar_plastic_strain_increment / effective_trial_stress);
333 
334  SymmTensor elastic_strain_increment(strain_increment);
335  elastic_strain_increment -= plastic_strain_increment;
336 
337  // compute stress increment
338  stress_new = *elasticityTensor() * elastic_strain_increment;
339 
340  // update stress and plastic strain
341  stress_new += _stress_old;
342  _plastic_strain[_qp] += plastic_strain_increment;
343 
344  } // end of if statement
345 }
PLC_LSH(const InputParameters &parameters)
Definition: PLC_LSH.C:48
const MaterialProperty< SymmTensor > & _plastic_strain_old
Definition: PLC_LSH.h:52
SymmElasticityTensor * elasticityTensor() const
Definition: SolidModel.h:217
virtual void computeStress()
Compute the stress (sigma += deltaSigma)
Definition: PLC_LSH.C:92
bool _output_iteration_info
Definition: PLC_LSH.h:42
Real _n_exponent
Definition: PLC_LSH.h:33
virtual void initQpStatefulProperties()
Definition: SolidModel.C:692
Real _relative_tolerance
Definition: PLC_LSH.h:43
SolidModel is the base class for all this module&#39;s solid mechanics material models.
Definition: SolidModel.h:31
Real _absolute_stress_tolerance
Definition: PLC_LSH.h:46
virtual void initQpStatefulProperties()
Definition: PLC_LSH.C:85
void computeLSH(const SymmTensor &strain_increment, SymmTensor &plastic_strain_increment, SymmTensor &stress_new)
Definition: PLC_LSH.C:262
Real _gas_constant
Definition: PLC_LSH.h:36
InputParameters validParams< SolidModel >()
Definition: SolidModel.C:27
const bool _has_temp
Definition: SolidModel.h:91
const MaterialProperty< SymmTensor > & _creep_strain_old
Definition: PLC_LSH.h:49
const MaterialProperty< Real > & _hardening_variable_old
Definition: PLC_LSH.h:55
MaterialProperty< SymmTensor > & _creep_strain
Definition: PLC_LSH.h:48
Real _coefficient
Definition: PLC_LSH.h:32
Real trace() const
Definition: SymmTensor.h:95
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
SymmTensor _strain_increment
Definition: SolidModel.h:143
SymmTensor _stress_old
Definition: SolidModel.h:112
Real _yield_stress
Definition: PLC_LSH.h:38
const VariableValue & _temperature
Definition: SolidModel.h:92
unsigned int _max_its
Definition: PLC_LSH.h:41
InputParameters validParams< PLC_LSH >()
Definition: PLC_LSH.C:13
void computeCreep(const SymmTensor &strain_increment, SymmTensor &creep_strain_increment, SymmTensor &stress_new)
Definition: PLC_LSH.C:172
Real _hardening_constant
Definition: PLC_LSH.h:39
MaterialProperty< SymmTensor > & _stress
Definition: SolidModel.h:106
Real _m_exponent
Definition: PLC_LSH.h:34
void addDiag(Real value)
Definition: SymmTensor.h:279
Real doubleContraction(const SymmTensor &rhs) const
Definition: SymmTensor.h:257
MaterialProperty< Real > & _hardening_variable
Definition: PLC_LSH.h:54
Real _activation_energy
Definition: PLC_LSH.h:35
Real _absolute_tolerance
Definition: PLC_LSH.h:44
Real _shear_modulus
Definition: SolidModel.h:71
static unsigned int counter
MaterialProperty< SymmTensor > & _plastic_strain
Definition: PLC_LSH.h:51