Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "InertialTorque.h"
      11             : 
      12             : #include "MooseVariable.h"
      13             : #include "PermutationTensor.h"
      14             : 
      15             : registerMooseObject("TensorMechanicsApp", InertialTorque);
      16             : 
      17             : InputParameters
      18          42 : InertialTorque::validParams()
      19             : {
      20          42 :   InputParameters params = TimeKernel::validParams();
      21          42 :   params.addClassDescription("Kernel for inertial torque: density * displacement x acceleration");
      22          42 :   params.set<bool>("use_displaced_mesh") = false;
      23          84 :   params.addRequiredRangeCheckedParam<unsigned int>(
      24             :       "component",
      25             :       "component<3",
      26             :       "The component of the Cosserat rotation Variable prescribed to "
      27             :       "this Kernel (0 for x, 1 for y, 2 for z)");
      28          84 :   params.addRequiredCoupledVar("displacements", "The 3 displacement variables");
      29          84 :   params.addRequiredCoupledVar("velocities", "The 3 velocity variables");
      30          84 :   params.addRequiredCoupledVar("accelerations", "The 3 acceleration variables");
      31          84 :   params.addRequiredParam<Real>("beta", "beta parameter for Newmark Time integration");
      32          84 :   params.addRequiredParam<Real>("gamma", "gamma parameter for Newmark Time integration");
      33          84 :   params.addParam<MaterialPropertyName>(
      34             :       "density", "density", "Name of Material Property that provides the density");
      35          84 :   params.addParam<MaterialPropertyName>("eta",
      36          84 :                                         0.0,
      37             :                                         "Name of material property or a constant real "
      38             :                                         "number defining the eta parameter for the "
      39             :                                         "Rayleigh damping.");
      40          84 :   params.addParam<Real>("alpha",
      41          84 :                         0,
      42             :                         "alpha parameter for mass dependent numerical damping induced "
      43             :                         "by HHT time integration scheme");
      44          42 :   return params;
      45           0 : }
      46             : 
      47          21 : InertialTorque::InertialTorque(const InputParameters & parameters)
      48             :   : TimeKernel(parameters),
      49          21 :     _density(getMaterialProperty<Real>("density")),
      50          42 :     _beta(getParam<Real>("beta")),
      51          42 :     _gamma(getParam<Real>("gamma")),
      52          42 :     _eta(getMaterialProperty<Real>("eta")),
      53          42 :     _alpha(getParam<Real>("alpha")),
      54          42 :     _component(getParam<unsigned int>("component")),
      55          21 :     _ndisp(coupledComponents("displacements")),
      56          21 :     _disp_num(coupledIndices("displacements")),
      57          21 :     _disp(coupledValues("displacements")),
      58          21 :     _disp_old(coupledValuesOld("displacements")),
      59          21 :     _vel_old(coupledValuesOld("velocities")),
      60          21 :     _accel_old(coupledValuesOld("accelerations")),
      61          21 :     _accel(_ndisp),
      62          21 :     _vel(_ndisp),
      63          21 :     _daccel(_ndisp),
      64          42 :     _dvel(_ndisp)
      65             : {
      66          63 :   if (_ndisp != 3 || coupledComponents("velocities") != 3 ||
      67          63 :       coupledComponents("accelerations") != 3)
      68           0 :     mooseError(
      69             :         "InertialTorque: This Kernel is only defined for 3-dimensional simulations so 3 "
      70             :         "displacement variables, 3 velocity variables and 3 acceleration variables are needed");
      71          21 : }
      72             : 
      73             : Real
      74       15744 : InertialTorque::computeQpResidual()
      75             : {
      76       62976 :   for (unsigned int j = 0; j < _ndisp; ++j)
      77             :   {
      78             :     // Newmark and damping
      79       47232 :     _accel[j] = 1.0 / _beta *
      80       47232 :                 ((((*_disp[j])[_qp] - (*_disp_old[j])[_qp]) / (_dt * _dt)) -
      81       47232 :                  (*_vel_old[j])[_qp] / _dt - (*_accel_old[j])[_qp] * (0.5 - _beta));
      82       47232 :     _vel[j] = (*_vel_old[j])[_qp] + (_dt * (1 - _gamma)) * (*_accel_old[j])[_qp] +
      83       47232 :               _gamma * _dt * _accel[j];
      84       47232 :     _accel[j] =
      85       47232 :         _accel[j] + _vel[j] * _eta[_qp] * (1 + _alpha) - _alpha * _eta[_qp] * (*_vel_old[j])[_qp];
      86             :   }
      87             : 
      88             :   Real the_sum = 0.0;
      89       62976 :   for (unsigned int j = 0; j < _ndisp; ++j)
      90      188928 :     for (unsigned int k = 0; k < _ndisp; ++k)
      91      141696 :       the_sum += PermutationTensor::eps(_component, j, k) * (*_disp[j])[_qp] * _accel[k];
      92       15744 :   return _test[_i][_qp] * _density[_qp] * the_sum;
      93             : }
      94             : 
      95             : Real
      96       13312 : InertialTorque::computeQpJacobian()
      97             : {
      98       13312 :   return InertialTorque::computeQpOffDiagJacobian(_var.number());
      99             : }
     100             : 
     101             : Real
     102       19456 : InertialTorque::computeQpOffDiagJacobian(unsigned int jvar)
     103             : {
     104       77824 :   for (unsigned int j = 0; j < _ndisp; ++j)
     105             :   {
     106       58368 :     _accel[j] = 1.0 / _beta *
     107       58368 :                 ((((*_disp[j])[_qp] - (*_disp_old[j])[_qp]) / (_dt * _dt)) -
     108       58368 :                  (*_vel_old[j])[_qp] / _dt - (*_accel_old[j])[_qp] * (0.5 - _beta));
     109       58368 :     _daccel[j] = 1.0 / _beta / _dt / _dt;
     110       58368 :     _vel[j] = (*_vel_old[j])[_qp] + (_dt * (1 - _gamma)) * (*_accel_old[j])[_qp] +
     111       58368 :               _gamma * _dt * _accel[j];
     112       58368 :     _dvel[j] = _gamma * _dt * _daccel[j];
     113       58368 :     _accel[j] =
     114       58368 :         _accel[j] + _vel[j] * _eta[_qp] * (1 + _alpha) - _alpha * _eta[_qp] * (*_vel_old[j])[_qp];
     115       58368 :     _daccel[j] = _daccel[j] + _dvel[j] * _eta[_qp] * (1 + _alpha);
     116             :   }
     117             : 
     118       45056 :   for (unsigned v = 0; v < _ndisp; ++v)
     119             :   {
     120       41984 :     if (v == _component)
     121       16384 :       continue; // PermutationTensor will kill this
     122       25600 :     if (jvar == _disp_num[v])
     123             :     {
     124             :       Real the_sum = 0.0;
     125       65536 :       for (unsigned int k = 0; k < _ndisp; ++k)
     126       49152 :         the_sum += PermutationTensor::eps(_component, v, k) * _accel[k];
     127       65536 :       for (unsigned int j = 0; j < _ndisp; ++j)
     128       49152 :         the_sum += PermutationTensor::eps(_component, j, v) * (*_disp[j])[_qp] * _daccel[v];
     129       16384 :       return _test[_i][_qp] * _density[_qp] * the_sum * _phi[_j][_qp];
     130             :     }
     131             :   }
     132             :   return 0.0;
     133             : }