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