www.mooseframework.org
StressDivergenceTensors.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 
10 // MOOSE includes
11 #include "ElasticityTensorTools.h"
12 #include "Material.h"
13 #include "MooseMesh.h"
14 #include "MooseVariable.h"
15 #include "SystemBase.h"
16 
17 #include "libmesh/quadrature.h"
18 
19 template <>
20 InputParameters
22 {
23  InputParameters params = validParams<ALEKernel>();
24  params.addClassDescription("Stress divergence kernel for the Cartesian coordinate system");
25  params.addRequiredParam<unsigned int>("component",
26  "An integer corresponding to the direction "
27  "the variable this kernel acts in. (0 for x, "
28  "1 for y, 2 for z)");
29  params.addRequiredCoupledVar("displacements",
30  "The string of displacements suitable for the problem statement");
31  params.addCoupledVar("temperature",
32  "The name of the temperature variable used in the "
33  "ComputeThermalExpansionEigenstrain. (Not required for "
34  "simulations without temperature coupling.)");
35  params.addParam<std::string>(
36  "thermal_eigenstrain_name",
37  "thermal_eigenstrain",
38  "The eigenstrain_name used in the ComputeThermalExpansionEigenstrain.");
39  params.addParam<std::string>("base_name", "Material property base name");
40  params.set<bool>("use_displaced_mesh") = false;
41  params.addParam<bool>(
42  "use_finite_deform_jacobian", false, "Jacobian for corotational finite strain");
43  params.addParam<bool>("volumetric_locking_correction",
44  false,
45  "Set to false to turn off volumetric locking correction");
46  return params;
47 }
48 
49 StressDivergenceTensors::StressDivergenceTensors(const InputParameters & parameters)
50  : ALEKernel(parameters),
51  _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
52  _use_finite_deform_jacobian(getParam<bool>("use_finite_deform_jacobian")),
53  _stress(getMaterialPropertyByName<RankTwoTensor>(_base_name + "stress")),
54  _Jacobian_mult(getMaterialPropertyByName<RankFourTensor>(_base_name + "Jacobian_mult")),
55  _component(getParam<unsigned int>("component")),
56  _ndisp(coupledComponents("displacements")),
57  _disp_var(_ndisp),
58  _temp_coupled(isCoupled("temperature")),
59  _temp_var(_temp_coupled ? coupled("temperature") : 0),
60  _deigenstrain_dT(_temp_coupled
61  ? &getMaterialPropertyDerivative<RankTwoTensor>(
62  getParam<std::string>("thermal_eigenstrain_name"),
63  getVar("temperature", 0)->name())
64  : nullptr),
65  _avg_grad_test(_test.size(), std::vector<Real>(3, 0.0)),
66  _avg_grad_phi(_phi.size(), std::vector<Real>(3, 0.0)),
67  _volumetric_locking_correction(getParam<bool>("volumetric_locking_correction"))
68 {
69  for (unsigned int i = 0; i < _ndisp; ++i)
70  _disp_var[i] = coupled("displacements", i);
71 
72  // Checking for consistency between mesh size and length of the provided displacements vector
73  if (_ndisp != _mesh.dimension())
74  mooseError("The number of displacement variables supplied must match the mesh dimension.");
75 
77  {
79  &getMaterialProperty<RankTwoTensor>(_base_name + "deformation_gradient");
81  &getMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient");
82  _rotation_increment = &getMaterialProperty<RankTwoTensor>(_base_name + "rotation_increment");
83  }
84 
85  // Error if volumetic locking correction is turned on for 1D problems
86  if (_ndisp == 1 && _volumetric_locking_correction)
87  mooseError("Volumetric locking correction should be set to false for 1-D problems.");
88 }
89 
90 void
92 {
93  if (getBlockCoordSystem() != Moose::COORD_XYZ)
94  mooseError(
95  "The coordinate system in the Problem block must be set to XYZ for cartesian geometries.");
96 }
97 
98 void
100 {
101  DenseVector<Number> & re = _assembly.residualBlock(_var.number());
102  _local_re.resize(re.size());
103  _local_re.zero();
104 
107 
108  precalculateResidual();
109  for (_i = 0; _i < _test.size(); ++_i)
110  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
111  _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
112 
113  re += _local_re;
114 
115  if (_has_save_in)
116  {
117  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
118  for (const auto & var : _save_in)
119  var->sys().solution().add_vector(_local_re, var->dofIndices());
120  }
121 }
122 
123 Real
125 {
126  Real residual = _stress[_qp].row(_component) * _grad_test[_i][_qp];
127  // volumetric locking correction
129  residual += _stress[_qp].trace() / 3.0 *
130  (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
131 
132  return residual;
133 }
134 
135 void
137 {
139  {
140  _finite_deform_Jacobian_mult.resize(_qrule->n_points());
141 
142  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
144 
146  }
147  else
148  {
150  {
153  }
154  Kernel::computeJacobian();
155  }
156 }
157 
158 void
160 {
162  {
163  _finite_deform_Jacobian_mult.resize(_qrule->n_points());
164 
165  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
167 
169  }
170  else
171  {
173  {
176  }
177  Kernel::computeOffDiagJacobian(jvar);
178  }
179 }
180 
181 Real
183 {
186  _component,
187  _component,
188  _grad_test[_i][_qp],
189  _grad_phi_undisplaced[_j][_qp]);
190 
191  Real sum_C3x3 = _Jacobian_mult[_qp].sum3x3();
192  RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum3x1();
193 
194  Real jacobian = 0.0;
195  // B^T_i * C * B_j
197  _Jacobian_mult[_qp], _component, _component, _grad_test[_i][_qp], _grad_phi[_j][_qp]);
198 
200  {
201  // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
202  // 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
203 
204  // Bvol^T_i * C * Bvol_j
205  jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
206  (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 9.0;
207 
208  // B^T_i * C * Bvol_j
209  jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
210  (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 3.0;
211 
212  // Bvol^T_i * C * B_j
213  RankTwoTensor phi;
214  if (_component == 0)
215  {
216  phi(0, 0) = _grad_phi[_j][_qp](0);
217  phi(0, 1) = phi(1, 0) = _grad_phi[_j][_qp](1);
218  phi(0, 2) = phi(2, 0) = _grad_phi[_j][_qp](2);
219  }
220  else if (_component == 1)
221  {
222  phi(1, 1) = _grad_phi[_j][_qp](1);
223  phi(0, 1) = phi(1, 0) = _grad_phi[_j][_qp](0);
224  phi(1, 2) = phi(2, 1) = _grad_phi[_j][_qp](2);
225  }
226  else if (_component == 2)
227  {
228  phi(2, 2) = _grad_phi[_j][_qp](2);
229  phi(0, 2) = phi(2, 0) = _grad_phi[_j][_qp](0);
230  phi(1, 2) = phi(2, 1) = _grad_phi[_j][_qp](1);
231  }
232 
233  jacobian += (_Jacobian_mult[_qp] * phi).trace() *
234  (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
235  }
236  return jacobian;
237 }
238 
239 Real
241 {
242  // off-diagonal Jacobian with respect to a coupled displacement component
243  for (unsigned int coupled_component = 0; coupled_component < _ndisp; ++coupled_component)
244  if (jvar == _disp_var[coupled_component])
245  {
248  _component,
249  coupled_component,
250  _grad_test[_i][_qp],
251  _grad_phi_undisplaced[_j][_qp]);
252 
253  const Real sum_C3x3 = _Jacobian_mult[_qp].sum3x3();
254  const RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum3x1();
255  Real jacobian = 0.0;
256 
257  // B^T_i * C * B_j
259  _component,
260  coupled_component,
261  _grad_test[_i][_qp],
262  _grad_phi[_j][_qp]);
263 
265  {
266  // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
267  // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j + Bvol^T_i * C * B_j + B^T_i * C *
268  // Bvol_j
269 
270  // Bvol^T_i * C * Bvol_j
271  jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
272  (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
273  9.0;
274 
275  // B^T_i * C * Bvol_j
276  jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
277  (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
278  3.0;
279 
280  // Bvol^T_i * C * B_i
281  RankTwoTensor phi;
282  for (unsigned int i = 0; i < 3; ++i)
283  phi(coupled_component, i) = _grad_phi[_j][_qp](i);
284 
285  jacobian += (_Jacobian_mult[_qp] * phi).trace() *
286  (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
287  }
288 
289  return jacobian;
290  }
291 
292  // off-diagonal Jacobian with respect to a coupled temperature variable
293  if (_temp_coupled && jvar == _temp_var)
294  return -((_Jacobian_mult[_qp] * (*_deigenstrain_dT)[_qp]) *
295  _grad_test[_i][_qp])(_component)*_phi[_j][_qp];
296 
297  return 0.0;
298 }
299 
300 void
302 {
303  const RankTwoTensor I(RankTwoTensor::initIdentity);
304  const RankFourTensor II_ijkl = I.mixedProductIkJl(I);
305 
306  // Bring back to unrotated config
307  const RankTwoTensor unrotated_stress =
308  (*_rotation_increment)[_qp].transpose() * _stress[_qp] * (*_rotation_increment)[_qp];
309 
310  // Incremental deformation gradient Fhat
311  const RankTwoTensor Fhat =
312  (*_deformation_gradient)[_qp] * (*_deformation_gradient_old)[_qp].inverse();
313  const RankTwoTensor Fhatinv = Fhat.inverse();
314 
315  const RankTwoTensor rot_times_stress = (*_rotation_increment)[_qp] * unrotated_stress;
316  const RankFourTensor dstress_drot =
317  I.mixedProductIkJl(rot_times_stress) + I.mixedProductJkIl(rot_times_stress);
318  const RankFourTensor rot_rank_four =
319  (*_rotation_increment)[_qp].mixedProductIkJl((*_rotation_increment)[_qp]);
320  const RankFourTensor drot_dUhatinv = Fhat.mixedProductIkJl(I);
321 
322  const RankTwoTensor A = I - Fhatinv;
323 
324  // Ctilde = Chat^-1 - I
325  const RankTwoTensor Ctilde = A * A.transpose() - A - A.transpose();
326  const RankFourTensor dCtilde_dFhatinv =
327  -I.mixedProductIkJl(A) - I.mixedProductJkIl(A) + II_ijkl + I.mixedProductJkIl(I);
328 
329  // Second order approximation of Uhat - consistent with strain increment definition
330  // const RankTwoTensor Uhat = I - 0.5 * Ctilde - 3.0/8.0 * Ctilde * Ctilde;
331 
332  RankFourTensor dUhatinv_dCtilde =
333  0.5 * II_ijkl - 1.0 / 8.0 * (I.mixedProductIkJl(Ctilde) + Ctilde.mixedProductIkJl(I));
334  RankFourTensor drot_dFhatinv = drot_dUhatinv * dUhatinv_dCtilde * dCtilde_dFhatinv;
335 
336  drot_dFhatinv -= Fhat.mixedProductIkJl((*_rotation_increment)[_qp].transpose());
337  _finite_deform_Jacobian_mult[_qp] = dstress_drot * drot_dFhatinv;
338 
339  const RankFourTensor dstrain_increment_dCtilde =
340  -0.5 * II_ijkl + 0.25 * (I.mixedProductIkJl(Ctilde) + Ctilde.mixedProductIkJl(I));
342  rot_rank_four * _Jacobian_mult[_qp] * dstrain_increment_dCtilde * dCtilde_dFhatinv;
343  _finite_deform_Jacobian_mult[_qp] += Fhat.mixedProductJkIl(_stress[_qp]);
344 
345  const RankFourTensor dFhat_dFhatinv = -Fhat.mixedProductIkJl(Fhat.transpose());
346  const RankTwoTensor dJ_dFhatinv = dFhat_dFhatinv.innerProductTranspose(Fhat.ddet());
347 
348  // Component from Jacobian derivative
349  _finite_deform_Jacobian_mult[_qp] += _stress[_qp].outerProduct(dJ_dFhatinv);
350 
351  // Derivative of Fhatinv w.r.t. undisplaced coordinates
352  const RankTwoTensor Finv = (*_deformation_gradient)[_qp].inverse();
353  const RankFourTensor dFhatinv_dGradu = -Fhatinv.mixedProductIkJl(Finv.transpose());
354  _finite_deform_Jacobian_mult[_qp] = _finite_deform_Jacobian_mult[_qp] * dFhatinv_dGradu;
355 }
356 
357 void
359 {
360  // Calculate volume averaged value of shape function derivative
361  _avg_grad_test.resize(_test.size());
362  for (_i = 0; _i < _test.size(); ++_i)
363  {
364  _avg_grad_test[_i].resize(3);
365  _avg_grad_test[_i][_component] = 0.0;
366  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
367  _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
368 
369  _avg_grad_test[_i][_component] /= _current_elem_volume;
370  }
371 }
372 
373 void
375 {
376  // Calculate volume average derivatives for phi
377  _avg_grad_phi.resize(_phi.size());
378  for (_i = 0; _i < _phi.size(); ++_i)
379  {
380  _avg_grad_phi[_i].resize(3);
381  for (unsigned int component = 0; component < _mesh.dimension(); ++component)
382  {
383  _avg_grad_phi[_i][component] = 0.0;
384  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
385  _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
386 
387  _avg_grad_phi[_i][component] /= _current_elem_volume;
388  }
389  }
390 }
std::vector< std::vector< Real > > _avg_grad_phi
Gradient of phi function averaged over the element. Used in volumetric locking correction calculation...
virtual void computeOffDiagJacobian(unsigned int jvar)
Definition: ALEKernel.C:40
bool _volumetric_locking_correction
Flag for volumetric locking correction.
unsigned int _ndisp
Coupled displacement variables.
virtual void computeJacobian()
Definition: ALEKernel.C:33
std::vector< std::vector< Real > > _avg_grad_test
Gradient of test function averaged over the element. Used in volumetric locking correction calculatio...
virtual void computeOffDiagJacobian(unsigned int jvar) override
Real component(const SymmTensor &symm_tensor, unsigned int index)
InputParameters validParams< StressDivergenceTensors >()
const MaterialProperty< RankTwoTensor > * _deformation_gradient_old
virtual void computeJacobian() override
virtual Real computeQpJacobian() override
const MaterialProperty< RankTwoTensor > * _rotation_increment
InputParameters validParams< ALEKernel >()
Definition: ALEKernel.C:15
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...
std::vector< unsigned int > _disp_var
std::vector< RankFourTensor > _finite_deform_Jacobian_mult
StressDivergenceTensors(const InputParameters &parameters)
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
virtual void initialSetup() override
const MaterialProperty< RankTwoTensor > * _deformation_gradient
const MaterialProperty< RankFourTensor > & _Jacobian_mult
virtual void computeFiniteDeformJacobian()
const VariablePhiGradient & _grad_phi_undisplaced
Shape and test functions on the undisplaced mesh.
Definition: ALEKernel.h:35
virtual void computeResidual() override
virtual Real computeQpResidual() override
const MaterialProperty< RankTwoTensor > & _stress