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