www.mooseframework.org
InertialTorque.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 "InertialTorque.h"
8 
9 #include "MooseVariable.h"
10 #include "PermutationTensor.h"
11 
12 template <>
13 InputParameters
15 {
16  InputParameters params = validParams<TimeKernel>();
17  params.addClassDescription("Kernel for interial torque: density * displacement x acceleration");
18  params.set<bool>("use_displaced_mesh") = false;
19  params.addRequiredRangeCheckedParam<unsigned int>(
20  "component",
21  "component<3",
22  "The component of the Cosserat rotation Variable prescribed to "
23  "this Kernel (0 for x, 1 for y, 2 for z)");
24  params.addRequiredCoupledVar("displacements", "The 3 displacement variables");
25  params.addRequiredCoupledVar("velocities", "The 3 velocity variables");
26  params.addRequiredCoupledVar("accelerations", "The 3 acceleration variables");
27  params.addRequiredParam<Real>("beta", "beta parameter for Newmark Time integration");
28  params.addRequiredParam<Real>("gamma", "gamma parameter for Newmark Time integration");
29  params.addParam<MaterialPropertyName>(
30  "density", "density", "Name of Material Property that provides the density");
31  params.addParam<MaterialPropertyName>("eta",
32  0.0,
33  "Name of material property or a constant real "
34  "number defining the eta parameter for the "
35  "Rayleigh damping.");
36  params.addParam<Real>("alpha",
37  0,
38  "alpha parameter for mass dependent numerical damping induced "
39  "by HHT time integration scheme");
40  return params;
41 }
42 
43 InertialTorque::InertialTorque(const InputParameters & parameters)
44  : TimeKernel(parameters),
45  _density(getMaterialProperty<Real>("density")),
46  _beta(getParam<Real>("beta")),
47  _gamma(getParam<Real>("gamma")),
48  _eta(getMaterialProperty<Real>("eta")),
49  _alpha(getParam<Real>("alpha")),
50  _component(getParam<unsigned int>("component")),
51  _ndisp(coupledComponents("displacements")),
52  _disp_num(_ndisp),
53  _disp(_ndisp),
54  _disp_old(_ndisp),
55  _vel_old(_ndisp),
56  _accel_old(_ndisp),
57  _accel(_ndisp),
58  _vel(_ndisp),
59  _daccel(_ndisp),
60  _dvel(_ndisp)
61 {
62  if (_ndisp != 3 || coupledComponents("velocities") != 3 ||
63  coupledComponents("accelerations") != 3)
64  mooseError(
65  "InertialTorque: This Kernel is only defined for 3-dimensional simulations so 3 "
66  "displacement variables, 3 velocity variables and 3 acceleration variables are needed");
67  for (unsigned i = 0; i < _ndisp; ++i)
68  {
69  _disp_num[i] = coupled("displacements", i);
70  _disp[i] = &coupledValue("displacements", i);
71  _disp_old[i] = &coupledValueOld("displacements", i);
72  _vel_old[i] = &coupledValueOld("velocities", i);
73  _accel_old[i] = &coupledValueOld("accelerations", i);
74  }
75 }
76 
77 Real
79 {
80  for (unsigned int j = 0; j < _ndisp; ++j)
81  {
82  // Newmark and damping
83  _accel[j] = 1.0 / _beta * ((((*_disp[j])[_qp] - (*_disp_old[j])[_qp]) / (_dt * _dt)) -
84  (*_vel_old[j])[_qp] / _dt - (*_accel_old[j])[_qp] * (0.5 - _beta));
85  _vel[j] = (*_vel_old[j])[_qp] + (_dt * (1 - _gamma)) * (*_accel_old[j])[_qp] +
86  _gamma * _dt * _accel[j];
87  _accel[j] =
88  _accel[j] + _vel[j] * _eta[_qp] * (1 + _alpha) - _alpha * _eta[_qp] * (*_vel_old[j])[_qp];
89  }
90 
91  Real the_sum = 0.0;
92  for (unsigned int j = 0; j < _ndisp; ++j)
93  for (unsigned int k = 0; k < _ndisp; ++k)
94  the_sum += PermutationTensor::eps(_component, j, k) * (*_disp[j])[_qp] * _accel[k];
95  return _test[_i][_qp] * _density[_qp] * the_sum;
96 }
97 
98 Real
100 {
101  return InertialTorque::computeQpOffDiagJacobian(_var.number());
102 }
103 
104 Real
106 {
107  for (unsigned int j = 0; j < _ndisp; ++j)
108  {
109  _accel[j] = 1.0 / _beta * ((((*_disp[j])[_qp] - (*_disp_old[j])[_qp]) / (_dt * _dt)) -
110  (*_vel_old[j])[_qp] / _dt - (*_accel_old[j])[_qp] * (0.5 - _beta));
111  _daccel[j] = 1.0 / _beta / _dt / _dt;
112  _vel[j] = (*_vel_old[j])[_qp] + (_dt * (1 - _gamma)) * (*_accel_old[j])[_qp] +
113  _gamma * _dt * _accel[j];
114  _dvel[j] = _gamma * _dt * _daccel[j];
115  _accel[j] =
116  _accel[j] + _vel[j] * _eta[_qp] * (1 + _alpha) - _alpha * _eta[_qp] * (*_vel_old[j])[_qp];
117  _daccel[j] = _daccel[j] + _dvel[j] * _eta[_qp] * (1 + _alpha);
118  }
119 
120  for (unsigned v = 0; v < _ndisp; ++v)
121  {
122  if (v == _component)
123  continue; // PermutationTensor will kill this
124  if (jvar == _disp_num[v])
125  {
126  Real the_sum = 0.0;
127  for (unsigned int k = 0; k < _ndisp; ++k)
128  the_sum += PermutationTensor::eps(_component, v, k) * _accel[k];
129  for (unsigned int j = 0; j < _ndisp; ++j)
130  the_sum += PermutationTensor::eps(_component, j, v) * (*_disp[j])[_qp] * _daccel[v];
131  return _test[_i][_qp] * _density[_qp] * the_sum * _phi[_j][_qp];
132  }
133  }
134  return 0.0;
135 }
std::vector< Real > _dvel
Derivative of velocity with respect to displacement.
const MaterialProperty< Real > & _eta
Rayleigh-damping eta parameter.
std::vector< const VariableValue * > _disp_old
Old value of displacements.
InertialTorque(const InputParameters &parameters)
const unsigned _ndisp
Number of displacement variables. This must be 3.
const Real _beta
Newmark beta parameter.
std::vector< Real > _accel
Acceleration (instantiating this vector avoids re-creating a new vector every residual calculation) ...
const Real _alpha
HHT alpha parameter.
virtual Real computeQpJacobian() override
std::vector< const VariableValue * > _accel_old
Old value of accelerations.
virtual Real computeQpResidual() override
const unsigned _component
Component of the cross-product desired.
std::vector< Real > _daccel
Derivative of acceleration with respect to displacement.
const Real _gamma
Newmark gamma parameter.
const MaterialProperty< Real > & _density
density
std::vector< const VariableValue * > _disp
Displacements.
std::vector< unsigned > _disp_num
MOOSE internal variable numbers corresponding to the displacments.
std::vector< const VariableValue * > _vel_old
Old value of velocities.
std::vector< Real > _vel
Velocity (instantiating this vector avoids re-creating a new vector every residual calculation) ...
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
InputParameters validParams< InertialTorque >()