www.mooseframework.org
AbaqusCreepMaterial.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 "AbaqusCreepMaterial.h"
8 
9 #include "SymmTensor.h"
10 #include "Factory.h"
11 
12 #include <dlfcn.h>
13 #define QUOTE(macro) stringifyName(macro)
14 
15 template <>
16 InputParameters
18 {
19  InputParameters params = validParams<SolidModel>();
20  params.addRequiredParam<FileName>("plugin",
21  "The path to the compiled dynamic library for the "
22  "plugin you want to use (without -opt.plugin or "
23  "-dbg.plugin)");
24  params.addRequiredParam<Real>("youngs_modulus", "Young's Modulus");
25  params.addRequiredParam<Real>("poissons_ratio", "Poissons Ratio");
26  params.addRequiredParam<Real>("num_state_vars",
27  "The number of state variables this CREEP routine will use");
28  params.addRequiredParam<unsigned int>(
29  "integration_flag", "The creep integration method: Explicit = 0 and Implicit = 1");
30  params.addRequiredParam<unsigned int>(
31  "solve_definition", "Creep/Swell Explicit/Implicit Integration Definition to use: 1 - 5");
32  params.addParam<unsigned int>("routine_flag",
33  0,
34  "The flag determining when the routine is "
35  "called: Start of increment = 0 and End of "
36  "Increment = 1");
37  return params;
38 }
39 
40 AbaqusCreepMaterial::AbaqusCreepMaterial(const InputParameters & parameters)
41  : SolidModel(parameters),
42  _plugin(getParam<FileName>("plugin")),
43  _youngs_modulus(getParam<Real>("youngs_modulus")),
44  _poissons_ratio(getParam<Real>("poissons_ratio")),
45  _num_state_vars(getParam<Real>("num_state_vars")),
46  _integration_flag(getParam<unsigned int>("integration_flag")),
47  _solve_definition(getParam<unsigned int>("solve_definition")),
48  _routine_flag(getParam<unsigned int>("routine_flag")),
49  _state_var(declareProperty<std::vector<Real>>("state_var")),
50  _state_var_old(getMaterialPropertyOld<std::vector<Real>>("state_var")),
51  _trial_stress(declareProperty<SymmTensor>("trial_stress")),
52  _trial_stress_old(getMaterialPropertyOld<SymmTensor>("trial_stress")),
53  _dev_trial_stress(declareProperty<SymmTensor>("dev_trial_stress")),
54  _dev_trial_stress_old(getMaterialPropertyOld<SymmTensor>("dev_trial_stress")),
55  _ets(declareProperty<Real>("effective_trial_stress")),
56  _ets_old(getMaterialPropertyOld<Real>("effective_trial_stress")),
57  _stress(declareProperty<SymmTensor>("stress")),
58  _stress_old(getMaterialPropertyOld<SymmTensor>("stress")),
59  _creep_inc(declareProperty<Real>("creep_inc")),
60  _creep_inc_old(getMaterialPropertyOld<Real>("creep_inc")),
61  _total_creep(declareProperty<Real>("total_creep")),
62  _total_creep_old(getMaterialPropertyOld<Real>("total_creep")),
63  _swell_inc(declareProperty<Real>("swell_inc")),
64  _swell_inc_old(getMaterialPropertyOld<Real>("swell_inc")),
65  _total_swell(declareProperty<Real>("total_swell")),
66  _total_swell_old(getMaterialPropertyOld<Real>("total_swell"))
67 {
68 #if defined(METHOD)
69  _plugin += std::string("-") + QUOTE(METHOD) + ".plugin";
70 #endif
71 
72  _STATEV = new Real[_num_state_vars];
73 
74  // Set subroutine values
78 
79  // Calculate constants needed for elasticity tensor
80  _ebulk3 = _youngs_modulus / (1. - 2. * _poissons_ratio);
82  _eg = _eg2 / 2.;
83  _eg3 = 3. * _eg;
84  _elam = (_ebulk3 - _eg2) / 3.;
85 
89 
90  // Open the library
91  _handle = dlopen(_plugin.c_str(), RTLD_LAZY);
92 
93  if (!_handle)
94  {
95  std::ostringstream error;
96  error << "Cannot open library: " << dlerror() << '\n';
97  mooseError(error.str());
98  }
99 
100  // Reset errors
101  dlerror();
102 
103  // Snag the function pointer from the library
104  {
105  void * pointer = dlsym(_handle, "creep_");
106  _creep = *reinterpret_cast<creep_t *>(&pointer);
107  }
108 
109  // Catch errors
110  const char * dlsym_error = dlerror();
111  if (dlsym_error)
112  {
113  dlclose(_handle);
114  std::ostringstream error;
115  error << "Cannot load symbol 'creep_': " << dlsym_error << '\n';
116  mooseError(error.str());
117  }
118 }
119 
121 {
122  delete _STATEV;
123 
124  dlclose(_handle);
125 }
126 
127 void
129 {
130  for (unsigned qp(0); qp < n_points; ++qp)
131  {
132  // Initialize state variable vector
133  _state_var[qp].resize(_num_state_vars);
134  for (unsigned int i = 0; i < _num_state_vars; i++)
135  {
136  _state_var[qp][i] = 0.0;
137  }
138  }
139 }
140 
141 // void AbaqusCreepMaterial::modifyStrain(const unsigned int qp,
142 // const Real /*scale_factor*/,
143 // SymmTensor & strain_increment,
144 // SymmTensor & /*dstrain_increment_dT*/)
145 void
147 {
148  // Recover "old" state variables
149  for (unsigned int i = 0; i < _num_state_vars; i++)
150  _STATEV[i] = _state_var_old[_qp][i];
151 
152  // Initialize DECRA and DESWA arrays
153  for (unsigned int i = 0; i < 5; i++)
154  {
155  _DECRA[i] = 0.0;
156  _DESWA[i] = 0.0;
157  }
158 
159  // Calculate stress array components
163 
167 
171 
175 
176  // Calculate trial stress and deviatoric trial stress
177  SymmTensor trial_stress(_stress_component[0],
182  _stress_component[5]);
183  _trial_stress[_qp] = trial_stress;
184  _trial_stress[_qp] += _stress_old[_qp];
185  _dev_trial_stress[_qp] = _trial_stress[_qp];
186  _dev_trial_stress[_qp].addDiag(-_trial_stress[_qp].trace() / 3.0);
187 
188  // Calculate effective trial stress (_ets)
189  Real dts_squared = _dev_trial_stress[_qp].doubleContraction(_dev_trial_stress[_qp]);
190 
191  if (dts_squared >= 0.)
192  _ets[_qp] = std::sqrt(1.5 * dts_squared);
193  else
194  mooseError("Attempted to take square root of a negative number!\n");
195 
196  // Calculate gradient of dev stress potential (grad_dts)
197  // grad_dts = d(qtild)/d(dev_trial_stress)
198  SymmTensor delta_dts = _dev_trial_stress[_qp] - _dev_trial_stress_old[_qp];
199 
200  Real delta_ets = _ets[_qp] - _ets_old[_qp];
201 
202  Real grad_dts_potential[6];
203  for (unsigned int i = 0; i < 6; i++)
204  {
205  if (delta_dts.component(i) == 0.0)
206  grad_dts_potential[i] = 0.0;
207  else
208  grad_dts_potential[i] = delta_ets / delta_dts.component(i);
209  }
210 
211  SymmTensor grad_dts(grad_dts_potential[0],
212  grad_dts_potential[1],
213  grad_dts_potential[2],
214  grad_dts_potential[3],
215  grad_dts_potential[4],
216  grad_dts_potential[5]);
217 
218  // Pass variables in for information
219  _KSTEP = _t_step; // Step number
220  _TIME[0] = _dt; // Value of step time at the end of the increment - Check
221  _TIME[1] = _t; // Value of total time at the end of the increment - Check
222  _DTIME = _dt; // Time increment
223  _EC[0] = _total_creep_old[_qp]; // Metal and Drucker-Prager creep at the start of the increment
224  _EC[1] = _total_creep[_qp]; // Metal and Drucker-Prager creep at the end of the increment
225  _ESW[0] = _total_swell_old[_qp]; // Metal swell at the start of the increment
226  _ESW[1] = _total_swell[_qp]; // Metal swell at the end of the increment
227  _QTILD = _ets[_qp]; // Von mises equivalent stress
228  _P = -_trial_stress[_qp].trace(); // Equivalent pressure stress
229 
230  // Connection to extern statement
231  _creep(_DECRA,
232  _DESWA,
233  _STATEV,
234  &_SERD,
235  _EC,
236  _ESW,
237  &_P,
238  &_QTILD,
239  &_TEMP,
240  &_DTEMP,
241  _PREDEF,
242  _DPRED,
243  _TIME,
244  &_DTIME,
245  &_CMNAME,
246  &_LEXIMP,
247  &_LEND,
248  _COORDS,
249  &_NSTATV,
250  &_NOEL,
251  &_NPT,
252  &_LAYER,
253  &_KSPT,
254  &_KSTEP,
255  &_KINC);
256 
257  // Update state variables
258  for (unsigned int i = 0; i < _num_state_vars; i++)
259  _state_var[_qp][i] = _STATEV[i];
260 
261  // Solve for Incremental creep (creep_inc_used) and/or Incremental Swell (swell_inc_used) based on
262  // definition
263  // NOTE: Below are equations for metal creep and/or time-dependent volumetric swelling materials
264  // only
265  // Drucker-Prager, Capped Drucker-Prager, and Gasket materials have not been included
266  _creep_inc[_qp] = _DECRA[0];
267  _total_creep[_qp] = _creep_inc[_qp];
268  _total_creep[_qp] += _total_creep_old[_qp];
269 
270  _swell_inc[_qp] = _DESWA[0];
271  _total_swell[_qp] = _swell_inc[_qp];
272  _total_swell[_qp] += _total_swell_old[_qp];
273 
274  Real p = -_trial_stress[_qp].trace();
275  Real pold = -_trial_stress_old[_qp].trace();
276 
277  Real creep_inc_used = 0.0;
278  Real swell_inc_used = 0.0;
279 
280  if (_integration_flag == 0 && _solve_definition == 1)
281  {
282  creep_inc_used = _DECRA[0];
283  swell_inc_used = _DESWA[0];
284  }
285  else if (_integration_flag == 1 && _solve_definition == 2)
286  {
287  creep_inc_used =
288  (_DECRA[1] * (_total_creep[_qp] - _total_creep_old[_qp])) + _creep_inc_old[_qp];
289  swell_inc_used =
290  (_DESWA[1] * (_total_creep[_qp] - _total_creep_old[_qp])) + _swell_inc_old[_qp];
291  }
292  else if (_integration_flag == 1 && _solve_definition == 3)
293  {
294  creep_inc_used =
295  (_DECRA[2] * (_total_swell[_qp] - _total_swell_old[_qp])) + _creep_inc_old[_qp];
296  swell_inc_used =
297  (_DESWA[2] * (_total_swell[_qp] - _total_swell_old[_qp])) + _swell_inc_old[_qp];
298  }
299  else if (_integration_flag == 1 && _solve_definition == 4)
300  {
301  creep_inc_used = (_DECRA[3] * (p - pold)) + _creep_inc_old[_qp];
302  swell_inc_used = (_DESWA[3] * (p - pold)) + _swell_inc_old[_qp];
303  }
304  else if (_integration_flag == 1 && _solve_definition == 5)
305  {
306  creep_inc_used = (_DECRA[4] * (_ets[_qp] - _ets_old[_qp])) + _creep_inc_old[_qp];
307  swell_inc_used = (_DESWA[4] * (_ets[_qp] - _ets_old[_qp])) + _swell_inc_old[_qp];
308  }
309 
310  // Calculate Incremental Creep Strain (total_effects)
311  // Incremental creep strain = ((1/3)*(swell_inc_used)*R) + (creep_inc_used*grad_dts)
312  // R = The matrix with the anisotropic swelling ratios in the diagonal if anisotropic swelling is
313  // defined; Otherwise R = Identity
314  SymmTensor R(1., 1., 1., 0., 0., 0.);
315  SymmTensor total_effects = (R * (swell_inc_used / 3.)) + (grad_dts * (creep_inc_used));
316 
317  // Modify strain increment
318  _strain_increment += total_effects;
319 
323 
327 
331 
335 
336  // Update Stress
337  SymmTensor stressnew(_stress_component[0],
342  _stress_component[5]);
343  _stress[_qp] = stressnew;
344  _stress[_qp] += _stress_old[_qp];
345 }
MaterialProperty< SymmTensor > & _trial_stress
MaterialProperty< Real > & _ets
MaterialProperty< std::vector< Real > > & _state_var
MaterialProperty< SymmTensor > & _stress
const MaterialProperty< SymmTensor > & _stress_old
virtual void initStatefulProperties(unsigned n_points)
SolidModel is the base class for all this module&#39;s solid mechanics material models.
Definition: SolidModel.h:31
const MaterialProperty< Real > & _creep_inc_old
MaterialProperty< Real > & _swell_inc
const MaterialProperty< Real > & _total_swell_old
AbaqusCreepMaterial(const InputParameters &parameters)
const MaterialProperty< Real > & _ets_old
InputParameters validParams< AbaqusCreepMaterial >()
void(* creep_t)(Real DECRA[], Real DESWA[], Real STATEV[], Real *SERD, Real EC[], Real ESW[], Real *P, Real *QTILD, Real *TEMP, Real *DTEMP, Real PREDEF[], Real DPRED[], Real TIME[], Real *DTIME, Real *CMNAME, Real *LEXIMP, Real *LEND, Real COORDS[], Real *NSTATV, int *NOEL, int *NPT, int *LAYER, int *KSPT, int *KSTEP, int *KINC)
MaterialProperty< Real > & _total_creep
const MaterialProperty< Real > & _swell_inc_old
const MaterialProperty< SymmTensor > & _trial_stress_old
MaterialProperty< Real > & _total_swell
InputParameters validParams< SolidModel >()
Definition: SolidModel.C:27
MaterialProperty< Real > & _creep_inc
unsigned int _num_state_vars
SymmTensor _strain_increment
Definition: SolidModel.h:143
const MaterialProperty< std::vector< Real > > & _state_var_old
const MaterialProperty< SymmTensor > & _dev_trial_stress_old
unsigned int _integration_flag
unsigned int _solve_definition
void computeStress()
Compute the stress (sigma += deltaSigma)
const MaterialProperty< Real > & _total_creep_old
Real component(unsigned int i) const
Definition: SymmTensor.h:97
MaterialProperty< SymmTensor > & _dev_trial_stress