Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       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("SolidMechanicsApp", InertialTorque);
      16             : 
      17             : InputParameters
      18          92 : InertialTorque::validParams()
      19             : {
      20          92 :   InputParameters params = TimeKernel::validParams();
      21          92 :   params.addClassDescription("Kernel for inertial torque: density * displacement x acceleration");
      22          92 :   params.set<bool>("use_displaced_mesh") = false;
      23         184 :   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         184 :   params.addRequiredCoupledVar("displacements", "The 3 displacement variables");
      29         184 :   params.addRequiredCoupledVar("velocities", "The 3 velocity variables");
      30         184 :   params.addRequiredCoupledVar("accelerations", "The 3 acceleration variables");
      31         184 :   params.addRequiredParam<Real>("beta", "beta parameter for Newmark Time integration");
      32         184 :   params.addRequiredParam<Real>("gamma", "gamma parameter for Newmark Time integration");
      33         184 :   params.addParam<MaterialPropertyName>(
      34             :       "density", "density", "Name of Material Property that provides the density");
      35         184 :   params.addParam<MaterialPropertyName>("eta",
      36         184 :                                         0.0,
      37             :                                         "Name of material property or a constant real "
      38             :                                         "number defining the eta parameter for the "
      39             :                                         "Rayleigh damping.");
      40         184 :   params.addParam<Real>("alpha",
      41         184 :                         0,
      42             :                         "alpha parameter for mass dependent numerical damping induced "
      43             :                         "by HHT time integration scheme");
      44          92 :   return params;
      45           0 : }
      46             : 
      47          46 : InertialTorque::InertialTorque(const InputParameters & parameters)
      48             :   : TimeKernel(parameters),
      49          46 :     _density(getMaterialProperty<Real>("density")),
      50          92 :     _beta(getParam<Real>("beta")),
      51          92 :     _gamma(getParam<Real>("gamma")),
      52          92 :     _eta(getMaterialProperty<Real>("eta")),
      53          92 :     _alpha(getParam<Real>("alpha")),
      54          92 :     _component(getParam<unsigned int>("component")),
      55          46 :     _ndisp(coupledComponents("displacements")),
      56          46 :     _disp_num(coupledIndices("displacements")),
      57          46 :     _disp(coupledValues("displacements")),
      58          46 :     _disp_old(coupledValuesOld("displacements")),
      59          46 :     _vel_old(coupledValuesOld("velocities")),
      60          46 :     _accel_old(coupledValuesOld("accelerations")),
      61          46 :     _accel(_ndisp),
      62          46 :     _vel(_ndisp),
      63          46 :     _daccel(_ndisp),
      64          92 :     _dvel(_ndisp)
      65             : {
      66         138 :   if (_ndisp != 3 || coupledComponents("velocities") != 3 ||
      67         138 :       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          46 : }
      72             : 
      73             : Real
      74       29184 : InertialTorque::computeQpResidual()
      75             : {
      76      116736 :   for (unsigned int j = 0; j < _ndisp; ++j)
      77             :   {
      78             :     // Newmark and damping
      79       87552 :     _accel[j] = 1.0 / _beta *
      80       87552 :                 ((((*_disp[j])[_qp] - (*_disp_old[j])[_qp]) / (_dt * _dt)) -
      81       87552 :                  (*_vel_old[j])[_qp] / _dt - (*_accel_old[j])[_qp] * (0.5 - _beta));
      82       87552 :     _vel[j] = (*_vel_old[j])[_qp] + (_dt * (1 - _gamma)) * (*_accel_old[j])[_qp] +
      83       87552 :               _gamma * _dt * _accel[j];
      84       87552 :     _accel[j] =
      85       87552 :         _accel[j] + _vel[j] * _eta[_qp] * (1 + _alpha) - _alpha * _eta[_qp] * (*_vel_old[j])[_qp];
      86             :   }
      87             : 
      88             :   Real the_sum = 0.0;
      89      116736 :   for (unsigned int j = 0; j < _ndisp; ++j)
      90      350208 :     for (unsigned int k = 0; k < _ndisp; ++k)
      91      262656 :       the_sum += PermutationTensor::eps(_component, j, k) * (*_disp[j])[_qp] * _accel[k];
      92       29184 :   return _test[_i][_qp] * _density[_qp] * the_sum;
      93             : }
      94             : 
      95             : Real
      96       32256 : InertialTorque::computeQpJacobian()
      97             : {
      98       32256 :   return InertialTorque::computeQpOffDiagJacobian(_var.number());
      99             : }
     100             : 
     101             : Real
     102       44544 : InertialTorque::computeQpOffDiagJacobian(unsigned int jvar)
     103             : {
     104      178176 :   for (unsigned int j = 0; j < _ndisp; ++j)
     105             :   {
     106      133632 :     _accel[j] = 1.0 / _beta *
     107      133632 :                 ((((*_disp[j])[_qp] - (*_disp_old[j])[_qp]) / (_dt * _dt)) -
     108      133632 :                  (*_vel_old[j])[_qp] / _dt - (*_accel_old[j])[_qp] * (0.5 - _beta));
     109      133632 :     _daccel[j] = 1.0 / _beta / _dt / _dt;
     110      133632 :     _vel[j] = (*_vel_old[j])[_qp] + (_dt * (1 - _gamma)) * (*_accel_old[j])[_qp] +
     111      133632 :               _gamma * _dt * _accel[j];
     112      133632 :     _dvel[j] = _gamma * _dt * _daccel[j];
     113      133632 :     _accel[j] =
     114      133632 :         _accel[j] + _vel[j] * _eta[_qp] * (1 + _alpha) - _alpha * _eta[_qp] * (*_vel_old[j])[_qp];
     115      133632 :     _daccel[j] = _daccel[j] + _dvel[j] * _eta[_qp] * (1 + _alpha);
     116             :   }
     117             : 
     118      101376 :   for (unsigned v = 0; v < _ndisp; ++v)
     119             :   {
     120       95232 :     if (v == _component)
     121       38400 :       continue; // PermutationTensor will kill this
     122       56832 :     if (jvar == _disp_num[v])
     123             :     {
     124             :       Real the_sum = 0.0;
     125      153600 :       for (unsigned int k = 0; k < _ndisp; ++k)
     126      115200 :         the_sum += PermutationTensor::eps(_component, v, k) * _accel[k];
     127      153600 :       for (unsigned int j = 0; j < _ndisp; ++j)
     128      115200 :         the_sum += PermutationTensor::eps(_component, j, v) * (*_disp[j])[_qp] * _daccel[v];
     129       38400 :       return _test[_i][_qp] * _density[_qp] * the_sum * _phi[_j][_qp];
     130             :     }
     131             :   }
     132             :   return 0.0;
     133             : }