www.mooseframework.org
StressDivergenceRZTensors.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 
9 #include "Assembly.h"
10 #include "ElasticityTensorTools.h"
11 #include "libmesh/quadrature.h"
12 
13 template <>
14 InputParameters
16 {
17  InputParameters params = validParams<StressDivergenceTensors>();
18  params.addClassDescription(
19  "Calculate stress divergence for an axisymmetric problem in cylinderical coordinates.");
20  params.addRequiredParam<unsigned int>(
21  "component",
22  "An integer corresponding to the direction the variable this kernel acts in. (0 "
23  "for x, 1 for y, 2 for z; note in this kernel disp_x refers to the radial "
24  "displacement and disp_y refers to the axial displacement.)");
25  params.set<bool>("use_displaced_mesh") = true;
26  return params;
27 }
28 
29 StressDivergenceRZTensors::StressDivergenceRZTensors(const InputParameters & parameters)
30  : StressDivergenceTensors(parameters)
31 {
32 }
33 
34 void
36 {
37  if (getBlockCoordSystem() != Moose::COORD_RZ)
38  mooseError("The coordinate system in the Problem block must be set to RZ for axisymmetric "
39  "geometries.");
40 }
41 
42 Real
44 {
45  Real div = 0.0;
46  if (_component == 0)
47  {
48  div = _grad_test[_i][_qp](0) * _stress[_qp](0, 0) +
49  +(_test[_i][_qp] / _q_point[_qp](0)) * _stress[_qp](2, 2) +
50  +_grad_test[_i][_qp](1) * _stress[_qp](0, 1); // stress_{rz}
51 
52  // volumetric locking correction
54  div += (_avg_grad_test[_i][0] - _grad_test[_i][_qp](0) - _test[_i][_qp] / _q_point[_qp](0)) *
55  (_stress[_qp].trace()) / 3.0;
56  }
57  else if (_component == 1)
58  {
59  div = _grad_test[_i][_qp](1) * _stress[_qp](1, 1) +
60  +_grad_test[_i][_qp](0) * _stress[_qp](1, 0); // stress_{zr}
61 
62  // volumetric locking correction
64  div += (_avg_grad_test[_i][1] - _grad_test[_i][_qp](1)) * (_stress[_qp].trace()) / 3.0;
65  }
66  else
67  mooseError("Invalid component for this AxisymmetricRZ problem.");
68 
69  return div;
70 }
71 
72 Real
74 {
76 }
77 
78 Real
80 {
81  for (unsigned int i = 0; i < _ndisp; ++i)
82  {
83  if (jvar == _disp_var[i])
84  return calculateJacobian(_component, i);
85  }
86 
87  if (_temp_coupled && jvar == _temp_var)
88  {
89  Real jac = 0.0;
90  if (_component == 0)
91  {
92  for (unsigned k = 0; k < LIBMESH_DIM; ++k)
93  for (unsigned l = 0; l < LIBMESH_DIM; ++l)
94  jac -= (_grad_test[_i][_qp](0) * _Jacobian_mult[_qp](0, 0, k, l) +
95  _test[_i][_qp] / _q_point[_qp](0) * _Jacobian_mult[_qp](2, 2, k, l) +
96  _grad_test[_i][_qp](1) * _Jacobian_mult[_qp](0, 1, k, l)) *
97  (*_deigenstrain_dT)[_qp](k, l);
98  return jac * _phi[_j][_qp];
99  }
100  else if (_component == 1)
101  {
102  for (unsigned k = 0; k < LIBMESH_DIM; ++k)
103  for (unsigned l = 0; l < LIBMESH_DIM; ++l)
104  jac -= (_grad_test[_i][_qp](1) * _Jacobian_mult[_qp](1, 1, k, l) +
105  _grad_test[_i][_qp](0) * _Jacobian_mult[_qp](1, 0, k, l)) *
106  (*_deigenstrain_dT)[_qp](k, l);
107  return jac * _phi[_j][_qp];
108  }
109  }
110 
111  return 0.0;
112 }
113 
114 Real
115 StressDivergenceRZTensors::calculateJacobian(unsigned int ivar, unsigned int jvar)
116 {
117  // B^T_i * C * B_j
118  RealGradient test, test_z, phi, phi_z;
119  Real first_term = 0.0;
120  if (ivar == 0) // Case grad_test for x, requires contributions from stress_xx, stress_xy, and
121  // stress_zz
122  {
123  test(0) = _grad_test[_i][_qp](0);
124  test(1) = _grad_test[_i][_qp](1);
125  test_z(2) = _test[_i][_qp] / _q_point[_qp](0);
126  }
127  else // Case grad_test for y
128  {
129  test(0) = _grad_test[_i][_qp](0);
130  test(1) = _grad_test[_i][_qp](1);
131  }
132 
133  if (jvar == 0)
134  {
135  phi(0) = _grad_phi[_j][_qp](0);
136  phi(1) = _grad_phi[_j][_qp](1);
137  phi_z(2) = _phi[_j][_qp] / _q_point[_qp](0);
138  }
139  else
140  {
141  phi(0) = _grad_phi[_j][_qp](0);
142  phi(1) = _grad_phi[_j][_qp](1);
143  }
144 
145  if (ivar == 0 &&
146  jvar == 0) // Case when both phi and test are functions of x and z; requires four terms
147  {
148  const Real first_sum = ElasticityTensorTools::elasticJacobian(
149  _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_x and phi_x
150  const Real second_sum = ElasticityTensorTools::elasticJacobian(
151  _Jacobian_mult[_qp], 2, 2, test_z, phi_z); // test_z and phi_z
152  const Real mixed_sum1 = ElasticityTensorTools::elasticJacobian(
153  _Jacobian_mult[_qp], ivar, 2, test, phi_z); // test_x and phi_z
154  const Real mixed_sum2 = ElasticityTensorTools::elasticJacobian(
155  _Jacobian_mult[_qp], 2, jvar, test_z, phi); // test_z and phi_x
156 
157  first_term = first_sum + second_sum + mixed_sum1 + mixed_sum2;
158  }
159  else if (ivar == 0 && jvar == 1)
160  {
161  const Real first_sum = ElasticityTensorTools::elasticJacobian(
162  _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_x and phi_y
163  const Real mixed_sum2 = ElasticityTensorTools::elasticJacobian(
164  _Jacobian_mult[_qp], 2, jvar, test_z, phi); // test_z and phi_y
165 
166  first_term = first_sum + mixed_sum2;
167  }
168  else if (ivar == 1 && jvar == 0)
169  {
170  const Real second_sum = ElasticityTensorTools::elasticJacobian(
171  _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_y and phi_x
172  const Real mixed_sum1 = ElasticityTensorTools::elasticJacobian(
173  _Jacobian_mult[_qp], ivar, 2, test, phi_z); // test_y and phi_z
174 
175  first_term = second_sum + mixed_sum1;
176  }
177  else if (ivar == 1 && jvar == 1)
179  _Jacobian_mult[_qp], ivar, jvar, test, phi); // test_y and phi_y
180  else
181  mooseError("Invalid component in Jacobian Calculation");
182 
183  Real val = 0.0;
184  // volumetric locking correction
185  // K = Bbar^T_i * C * Bbar^T_j where Bbar = B + Bvol
186  // 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
188  {
189  RealGradient new_test(2, 0.0);
190  RealGradient new_phi(2, 0.0);
191 
192  new_test(0) = _grad_test[_i][_qp](0) + _test[_i][_qp] / _q_point[_qp](0);
193  new_test(1) = _grad_test[_i][_qp](1);
194  new_phi(0) = _grad_phi[_j][_qp](0) + _phi[_j][_qp] / _q_point[_qp](0);
195  new_phi(1) = _grad_phi[_j][_qp](1);
196 
197  // Bvol^T_i * C * Bvol_j
198  val += _Jacobian_mult[_qp].sum3x3() * (_avg_grad_test[_i][ivar] - new_test(ivar)) *
199  (_avg_grad_phi[_j][jvar] - new_phi(jvar)) / 3.0;
200 
201  // B^T_i * C * Bvol_j
202  RealGradient sum_3x1 = _Jacobian_mult[_qp].sum3x1();
203  if (ivar == 0 && jvar == 0)
204  val += (sum_3x1(0) * test(0) + sum_3x1(2) * test_z(2)) * (_avg_grad_phi[_j][0] - new_phi(0));
205  else if (ivar == 0 && jvar == 1)
206  val += (sum_3x1(0) * test(0) + sum_3x1(2) * test_z(2)) * (_avg_grad_phi[_j][1] - new_phi(1));
207  else if (ivar == 1 && jvar == 0)
208  val += sum_3x1(1) * test(1) * (_avg_grad_phi[_j][0] - new_phi(0));
209  else
210  val += sum_3x1(1) * test(1) * (_avg_grad_phi[_j][1] - new_phi(1));
211 
212  // Bvol^T_i * C * B_j
213  // val = trace (C * B_j) *(avg_grad_test[_i][ivar] - new_test(ivar))
214  if (jvar == 0)
215  for (unsigned int i = 0; i < 3; ++i)
216  val +=
217  (_Jacobian_mult[_qp](i, i, 0, 0) * phi(0) + _Jacobian_mult[_qp](i, i, 0, 1) * phi(1) +
218  _Jacobian_mult[_qp](i, i, 2, 2) * phi_z(2)) *
219  (_avg_grad_test[_i][ivar] - new_test(ivar));
220  else if (jvar == 1)
221  for (unsigned int i = 0; i < 3; ++i)
222  val +=
223  (_Jacobian_mult[_qp](i, i, 0, 1) * phi(0) + _Jacobian_mult[_qp](i, i, 1, 1) * phi(1)) *
224  (_avg_grad_test[_i][ivar] - new_test(ivar));
225  }
226 
227  return val / 3.0 + first_term;
228 }
229 
230 void
232 {
233  // calculate volume averaged value of shape function derivative
234  _avg_grad_test.resize(_test.size());
235  for (_i = 0; _i < _test.size(); ++_i)
236  {
237  _avg_grad_test[_i].resize(2);
238  _avg_grad_test[_i][_component] = 0.0;
239  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
240  {
241  if (_component == 0)
242  _avg_grad_test[_i][_component] +=
243  (_grad_test[_i][_qp](_component) + _test[_i][_qp] / _q_point[_qp](0)) * _JxW[_qp] *
244  _coord[_qp];
245  else
246  _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
247  }
248  _avg_grad_test[_i][_component] /= _current_elem_volume;
249  }
250 }
251 
252 void
254 {
255  _avg_grad_phi.resize(_phi.size());
256  for (_i = 0; _i < _phi.size(); ++_i)
257  {
258  _avg_grad_phi[_i].resize(2);
259  for (unsigned int component = 0; component < 2; ++component)
260  {
261  _avg_grad_phi[_i][component] = 0.0;
262  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
263  {
264  if (component == 0)
265  _avg_grad_phi[_i][component] +=
266  (_grad_phi[_i][_qp](component) + _phi[_i][_qp] / _q_point[_qp](0)) * _JxW[_qp] *
267  _coord[_qp];
268  else
269  _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
270  }
271  _avg_grad_phi[_i][component] /= _current_elem_volume;
272  }
273  }
274 }
std::vector< std::vector< Real > > _avg_grad_phi
Gradient of phi function averaged over the element. Used in volumetric locking correction calculation...
bool _volumetric_locking_correction
Flag for volumetric locking correction.
unsigned int _ndisp
Coupled displacement variables.
std::vector< std::vector< Real > > _avg_grad_test
Gradient of test function averaged over the element. Used in volumetric locking correction calculatio...
virtual void initialSetup() override
Real component(const SymmTensor &symm_tensor, unsigned int index)
InputParameters validParams< StressDivergenceRZTensors >()
Real elasticJacobian(const RankFourTensor &r4t, unsigned int i, unsigned int k, const RealGradient &grad_test, const RealGradient &grad_phi)
This is used for the standard kernel stress_ij*d(test)/dx_j, when varied wrt u_k Jacobian entry: d(st...
StressDivergenceRZTensors(const InputParameters &parameters)
std::vector< unsigned int > _disp_var
virtual Real computeQpResidual() override
virtual Real computeQpJacobian() override
const MaterialProperty< RankFourTensor > & _Jacobian_mult
Real calculateJacobian(unsigned int ivar, unsigned int jvar)
virtual void computeAverageGradientTest() override
virtual void computeAverageGradientPhi() override
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
const MaterialProperty< RankTwoTensor > & _stress
InputParameters validParams< StressDivergenceTensors >()
StressDivergenceTensors mostly copies from StressDivergence.