www.mooseframework.org
StressDivergence.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 "StressDivergence.h"
8 
9 // MOOSE includes
10 #include "Assembly.h"
11 #include "Material.h"
12 #include "MooseMesh.h"
13 #include "MooseVariable.h"
14 #include "SymmElasticityTensor.h"
15 #include "SystemBase.h"
16 
17 #include "libmesh/quadrature.h"
18 
19 template <>
20 InputParameters
22 {
23  InputParameters params = validParams<Kernel>();
24  params.addRequiredParam<unsigned int>("component",
25  "An integer corresponding to the direction "
26  "the variable this kernel acts in. (0 for x, "
27  "1 for y, 2 for z)");
28  params.addCoupledVar("disp_x", "The x displacement");
29  params.addCoupledVar("disp_y", "The y displacement");
30  params.addCoupledVar("disp_z", "The z displacement");
31  params.addCoupledVar("temp", "The temperature");
32  params.addParam<Real>("zeta", 0.0, "Stiffness dependent Rayleigh damping coefficient");
33  params.addParam<Real>("alpha", 0.0, "alpha parameter required for HHT time integration");
34  params.addParam<std::string>(
35  "appended_property_name", "", "Name appended to material properties to make them unique");
36  params.addParam<bool>("volumetric_locking_correction",
37  true,
38  "Set to false to turn off volumetric locking correction");
39 
40  params.set<bool>("use_displaced_mesh") = true;
41 
42  return params;
43 }
44 
45 StressDivergence::StressDivergence(const InputParameters & parameters)
46  : Kernel(parameters),
47  _stress_older(getMaterialPropertyOlder<SymmTensor>(
48  "stress" + getParam<std::string>("appended_property_name"))),
49  _stress_old(getMaterialPropertyOld<SymmTensor>(
50  "stress" + getParam<std::string>("appended_property_name"))),
51  _stress(getMaterialProperty<SymmTensor>("stress" +
52  getParam<std::string>("appended_property_name"))),
53  _Jacobian_mult(getMaterialProperty<SymmElasticityTensor>(
54  "Jacobian_mult" + getParam<std::string>("appended_property_name"))),
55  _d_stress_dT(getMaterialProperty<SymmTensor>("d_stress_dT" +
56  getParam<std::string>("appended_property_name"))),
57  _component(getParam<unsigned int>("component")),
58  _xdisp_coupled(isCoupled("disp_x")),
59  _ydisp_coupled(isCoupled("disp_y")),
60  _zdisp_coupled(isCoupled("disp_z")),
61  _temp_coupled(isCoupled("temp")),
62  _xdisp_var(_xdisp_coupled ? coupled("disp_x") : 0),
63  _ydisp_var(_ydisp_coupled ? coupled("disp_y") : 0),
64  _zdisp_var(_zdisp_coupled ? coupled("disp_z") : 0),
65  _temp_var(_temp_coupled ? coupled("temp") : 0),
66  _zeta(getParam<Real>("zeta")),
67  _alpha(getParam<Real>("alpha")),
68  _avg_grad_test(_test.size(), std::vector<Real>(3, 0.0)),
69  _avg_grad_phi(_phi.size(), std::vector<Real>(3, 0.0)),
70  _volumetric_locking_correction(getParam<bool>("volumetric_locking_correction"))
71 {
72 }
73 
74 void
76 {
77  DenseVector<Number> & re = _assembly.residualBlock(_var.number());
78  _local_re.resize(re.size());
79  _local_re.zero();
80 
82  {
83  // calculate volume averaged value of shape function derivative
84  _avg_grad_test.resize(_test.size());
85  for (_i = 0; _i < _test.size(); _i++)
86  {
87  _avg_grad_test[_i].resize(3);
88  _avg_grad_test[_i][_component] = 0.0;
89  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
90  _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
91 
92  _avg_grad_test[_i][_component] /= _current_elem_volume;
93  }
94  }
95 
96  precalculateResidual();
97  for (_i = 0; _i < _test.size(); _i++)
98  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
99  _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
100 
101  re += _local_re;
102 
103  if (_has_save_in)
104  {
105  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
106  for (const auto & var : _save_in)
107  var->sys().solution().add_vector(_local_re, var->dofIndices());
108  }
109 }
110 
111 Real
113 {
114  Real residual = 0.0;
115  if ((_dt > 0) && ((_zeta != 0) || (_alpha != 0)))
116  {
117  residual = _stress[_qp].rowDot(_component, _grad_test[_i][_qp]) *
118  (1 + _alpha + (1 + _alpha) * _zeta / _dt) -
119  (_alpha + (1 + 2 * _alpha) * _zeta / _dt) *
120  _stress_old[_qp].rowDot(_component, _grad_test[_i][_qp]) +
121  (_alpha * _zeta / _dt) * _stress_older[_qp].rowDot(_component, _grad_test[_i][_qp]);
122 
123  // volumetric locking correction
125  residual += (_stress[_qp].trace() * (1 + _alpha + (1 + _alpha) * _zeta / _dt) -
126  (_alpha + (1 + 2 * _alpha) * _zeta / _dt) * _stress_old[_qp].trace() +
127  (_alpha * _zeta / _dt) * _stress_older[_qp].trace()) /
128  3.0 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
129  }
130  else
131  {
132  residual += _stress[_qp].rowDot(_component, _grad_test[_i][_qp]);
133 
134  // volumetric locking correction
136  residual += _stress[_qp].trace() / 3.0 *
137  (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
138  }
139  return residual;
140 }
141 
142 void
144 {
145  DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), _var.number());
146  _local_ke.resize(ke.m(), ke.n());
147  _local_ke.zero();
148 
150  {
151  // calculate volume averaged value of shape function derivative
152  _avg_grad_test.resize(_test.size());
153  for (_i = 0; _i < _test.size(); _i++)
154  {
155  _avg_grad_test[_i].resize(3);
156  _avg_grad_test[_i][_component] = 0.0;
157  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
158  _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
159 
160  _avg_grad_test[_i][_component] /= _current_elem_volume;
161  }
162 
163  _avg_grad_phi.resize(_phi.size());
164  for (_i = 0; _i < _phi.size(); _i++)
165  {
166  _avg_grad_phi[_i].resize(3);
167  for (unsigned int component = 0; component < _mesh.dimension(); component++)
168  {
169  _avg_grad_phi[_i][component] = 0.0;
170  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
171  _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
172 
173  _avg_grad_phi[_i][component] /= _current_elem_volume;
174  }
175  }
176  }
177 
178  for (_i = 0; _i < _test.size(); _i++)
179  for (_j = 0; _j < _phi.size(); _j++)
180  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
181  _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian();
182 
183  ke += _local_ke;
184 
185  if (_has_diag_save_in)
186  {
187  unsigned int rows = ke.m();
188  DenseVector<Number> diag(rows);
189  for (unsigned int i = 0; i < rows; i++)
190  diag(i) = _local_ke(i, i);
191 
192  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
193  for (const auto & var : _diag_save_in)
194  var->sys().solution().add_vector(diag, var->dofIndices());
195  }
196 }
197 
198 Real
200 {
201  Real sum_C3x3 = _Jacobian_mult[_qp].sum_3x3();
202  RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum_3x1();
203 
204  Real jacobian = 0.0;
205  // B^T_i * C * B_j
206  jacobian += _Jacobian_mult[_qp].stiffness(
207  _component, _component, _grad_test[_i][_qp], _grad_phi[_j][_qp]); // B^T_i * C *B_j
208 
210  {
211  // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
212  // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j + Bvol^T_i * C * B_j + B^T_i * C * Bvol_j
213 
214  // Bvol^T_i * C * Bvol_j
215  jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
216  (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 9.0;
217 
218  // B^T_i * C * Bvol_j
219  jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
220  (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 3.0;
221 
222  // Bvol^T_i * C * B_j
223  SymmTensor phi;
224  if (_component == 0)
225  {
226  phi.xx() = _grad_phi[_j][_qp](0);
227  phi.xy() = _grad_phi[_j][_qp](1);
228  phi.xz() = _grad_phi[_j][_qp](2);
229  }
230  else if (_component == 1)
231  {
232  phi.yy() = _grad_phi[_j][_qp](1);
233  phi.xy() = _grad_phi[_j][_qp](0);
234  phi.yz() = _grad_phi[_j][_qp](2);
235  }
236  else if (_component == 2)
237  {
238  phi.zz() = _grad_phi[_j][_qp](2);
239  phi.xz() = _grad_phi[_j][_qp](0);
240  phi.yz() = _grad_phi[_j][_qp](1);
241  }
242 
243  SymmTensor tmp(_Jacobian_mult[_qp] * phi);
244 
245  jacobian += (tmp.xx() + tmp.yy() + tmp.zz()) *
246  (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
247  }
248 
249  if (_dt > 0)
250  return jacobian * (1 + _alpha + _zeta / _dt);
251  else
252  return jacobian;
253 }
254 
255 void
257 {
258  if (jvar == _var.number())
259  computeJacobian();
260  else
261  {
263  {
264  // calculate volume averaged value of shape function derivative
265  _avg_grad_test.resize(_test.size());
266  for (_i = 0; _i < _test.size(); _i++)
267  {
268  _avg_grad_test[_i].resize(3);
269  _avg_grad_test[_i][_component] = 0.0;
270  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
271  _avg_grad_test[_i][_component] +=
272  _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
273 
274  _avg_grad_test[_i][_component] /= _current_elem_volume;
275  }
276 
277  _avg_grad_phi.resize(_phi.size());
278  for (_i = 0; _i < _phi.size(); _i++)
279  {
280  _avg_grad_phi[_i].resize(3);
281  for (unsigned int component = 0; component < _mesh.dimension(); component++)
282  {
283  _avg_grad_phi[_i][component] = 0.0;
284  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
285  _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
286 
287  _avg_grad_phi[_i][component] /= _current_elem_volume;
288  }
289  }
290  }
291 
292  DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
293 
294  for (_i = 0; _i < _test.size(); _i++)
295  for (_j = 0; _j < _phi.size(); _j++)
296  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
297  ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar);
298  }
299 }
300 
301 Real
303 {
304  unsigned int coupled_component = 0;
305 
306  bool active = false;
307 
308  if (_xdisp_coupled && jvar == _xdisp_var)
309  {
310  coupled_component = 0;
311  active = true;
312  }
313  else if (_ydisp_coupled && jvar == _ydisp_var)
314  {
315  coupled_component = 1;
316  active = true;
317  }
318  else if (_zdisp_coupled && jvar == _zdisp_var)
319  {
320  coupled_component = 2;
321  active = true;
322  }
323 
324  if (active)
325  {
326  Real sum_C3x3 = _Jacobian_mult[_qp].sum_3x3();
327  RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum_3x1();
328  Real jacobian = 0.0;
329  jacobian += _Jacobian_mult[_qp].stiffness(
330  _component, coupled_component, _grad_test[_i][_qp], _grad_phi[_j][_qp]); // B^T_i * C *B_j
331 
333  {
334  // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
335  // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j + Bvol^T_i * C * B_j + B^T_i * C *
336  // Bvol_j
337 
338  // Bvol^T_i * C * Bvol_j
339  jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
340  (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
341  9.0;
342 
343  // B^T_i * C * Bvol_j
344  jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
345  (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
346  3.0;
347 
348  // Bvol^T_i * C * B_i
349  SymmTensor phi;
350  if (coupled_component == 0)
351  {
352  phi.xx() = _grad_phi[_j][_qp](0);
353  phi.xy() = _grad_phi[_j][_qp](1);
354  phi.xz() = _grad_phi[_j][_qp](2);
355  }
356  else if (coupled_component == 1)
357  {
358  phi.yy() = _grad_phi[_j][_qp](1);
359  phi.xy() = _grad_phi[_j][_qp](0);
360  phi.yz() = _grad_phi[_j][_qp](2);
361  }
362  else if (coupled_component == 2)
363  {
364  phi.zz() = _grad_phi[_j][_qp](2);
365  phi.xz() = _grad_phi[_j][_qp](0);
366  phi.yz() = _grad_phi[_j][_qp](1);
367  }
368 
369  SymmTensor tmp(_Jacobian_mult[_qp] * phi);
370 
371  jacobian += (tmp.xx() + tmp.yy() + tmp.zz()) *
372  (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
373  }
374 
375  if (_dt > 0)
376  return jacobian * (1 + _alpha + _zeta / _dt);
377  else
378  return jacobian;
379  }
380 
381  if (_temp_coupled && jvar == _temp_var)
382  return _d_stress_dT[_qp].rowDot(_component, _grad_test[_i][_qp]) * _phi[_j][_qp];
383 
384  return 0;
385 }
This class defines a basic set of capabilities any elasticity tensor should have. ...
virtual Real computeQpJacobian()
const bool _temp_coupled
virtual void computeJacobian()
const bool _xdisp_coupled
const unsigned int _zdisp_var
virtual Real computeQpOffDiagJacobian(unsigned int jvar)
Real xx() const
Definition: SymmTensor.h:129
Real component(const SymmTensor &symm_tensor, unsigned int index)
std::vector< std::vector< Real > > _avg_grad_test
const MaterialProperty< SymmTensor > & _d_stress_dT
StressDivergence(const InputParameters &parameters)
const unsigned int _xdisp_var
std::vector< std::vector< Real > > _avg_grad_phi
const bool _ydisp_coupled
virtual void computeResidual()
const MaterialProperty< SymmTensor > & _stress
InputParameters validParams< StressDivergence >()
const unsigned int _ydisp_var
virtual Real computeQpResidual()
bool _volumetric_locking_correction
const MaterialProperty< SymmTensor > & _stress_old
const unsigned int _temp_var
const bool _zdisp_coupled
const MaterialProperty< SymmElasticityTensor > & _Jacobian_mult
const MaterialProperty< SymmTensor > & _stress_older
virtual void computeOffDiagJacobian(unsigned int jvar)
const unsigned int _component