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 "DynamicStressDivergenceTensors.h"
      11             : #include "ElasticityTensorTools.h"
      12             : 
      13             : registerMooseObject("SolidMechanicsApp", DynamicStressDivergenceTensors);
      14             : 
      15             : InputParameters
      16        1356 : DynamicStressDivergenceTensors::validParams()
      17             : {
      18        1356 :   InputParameters params = StressDivergenceTensors::validParams();
      19        1356 :   params.addClassDescription(
      20             :       "Residual due to stress related Rayleigh damping and HHT time integration terms ");
      21        2712 :   params.addParam<MaterialPropertyName>("zeta",
      22        2712 :                                         0.0,
      23             :                                         "Name of material property or a constant real "
      24             :                                         "number defining the zeta parameter for the "
      25             :                                         "Rayleigh damping.");
      26        2712 :   params.addParam<Real>("alpha", 0, "alpha parameter for HHT time integration");
      27        2712 :   params.addParam<bool>("static_initialization",
      28        2712 :                         false,
      29             :                         "Set to true to get the system to "
      30             :                         "equilibrium under gravity by running a "
      31             :                         "quasi-static analysis (by solving Ku = F) "
      32             :                         "in the first time step");
      33        1356 :   return params;
      34           0 : }
      35             : 
      36         678 : DynamicStressDivergenceTensors::DynamicStressDivergenceTensors(const InputParameters & parameters)
      37             :   : StressDivergenceTensors(parameters),
      38         678 :     _stress_older(getMaterialPropertyOlderByName<RankTwoTensor>(_base_name + "stress")),
      39        1356 :     _stress_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "stress")),
      40        1356 :     _zeta(getMaterialProperty<Real>("zeta")),
      41        1356 :     _alpha(getParam<Real>("alpha")),
      42        2034 :     _static_initialization(getParam<bool>("static_initialization"))
      43             : {
      44         678 : }
      45             : 
      46             : Real
      47    24292736 : DynamicStressDivergenceTensors::computeQpResidual()
      48             : {
      49             :   /**
      50             :    *This kernel needs to be used only if either Rayleigh damping or numerical damping through HHT
      51             :    *time integration scheme needs to be added to the problem through the stiffness dependent damping
      52             :    * parameter _zeta or HHT parameter _alpha, respectively.
      53             :    *
      54             :    * The residual of _zeta*K*[(1+_alpha)vel-_alpha vel_old]+ alpha K [ u - uold] + K u is required
      55             :    * = _zeta*[(1+_alpha)d/dt (Div sigma)-alpha d/dt(Div sigma_old)] +alpha [Div sigma - Div
      56             :    *sigma_old]+ Div sigma
      57             :    * = _zeta*[(1+alpha)(Div sigma - Div sigma_old)/dt - alpha (Div sigma_old - Div sigma_older)/dt]
      58             :    *   + alpha [Div sigma - Div sigma_old] +Div sigma
      59             :    * = [(1+_alpha)*_zeta/dt +_alpha+1]* Div sigma - [(1+2_alpha)*_zeta/dt + _alpha] Div sigma_old +
      60             :    *_alpha*_zeta/dt Div sigma_older
      61             :    */
      62             : 
      63             :   Real residual = 0.0;
      64    24292736 :   if (_static_initialization && _t == _dt)
      65             :   {
      66             :     // If static initialization is true, then in the first step residual is only Ku which is
      67             :     // stress.grad(test).
      68        3840 :     residual += _stress[_qp].row(_component) * _grad_test[_i][_qp];
      69             : 
      70        3840 :     if (_volumetric_locking_correction)
      71           0 :       residual += _stress[_qp].trace() / 3.0 *
      72           0 :                   (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
      73             :   }
      74    24288896 :   else if (_dt > 0)
      75             :   {
      76    24288896 :     residual +=
      77    24288896 :         _stress[_qp].row(_component) * _grad_test[_i][_qp] *
      78    24288896 :             (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt) -
      79    48577792 :         (_alpha + (1.0 + 2.0 * _alpha) * _zeta[_qp] / _dt) * _stress_old[_qp].row(_component) *
      80    24288896 :             _grad_test[_i][_qp] +
      81    48577792 :         (_alpha * _zeta[_qp] / _dt) * _stress_older[_qp].row(_component) * _grad_test[_i][_qp];
      82             : 
      83    24288896 :     if (_volumetric_locking_correction)
      84           0 :       residual += (_stress[_qp].trace() * (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt) -
      85           0 :                    (_alpha + (1.0 + 2.0 * _alpha) * _zeta[_qp] / _dt) * _stress_old[_qp].trace() +
      86           0 :                    (_alpha * _zeta[_qp] / _dt) * _stress_older[_qp].trace()) /
      87           0 :                   3.0 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
      88             :   }
      89             : 
      90    24292736 :   return residual;
      91             : }
      92             : 
      93             : Real
      94    49298224 : DynamicStressDivergenceTensors::computeQpJacobian()
      95             : {
      96    49298224 :   if (_static_initialization && _t == _dt)
      97        6144 :     return StressDivergenceTensors::computeQpJacobian();
      98    49292080 :   else if (_dt > 0)
      99    49292080 :     return StressDivergenceTensors::computeQpJacobian() *
     100    49292080 :            (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt);
     101             :   else
     102             :     return 0.0;
     103             : }
     104             : 
     105             : Real
     106    59379712 : DynamicStressDivergenceTensors::computeQpOffDiagJacobian(unsigned int jvar)
     107             : {
     108             :   bool active = true;
     109             : 
     110             :   for (unsigned int i = 0; i < _ndisp; ++i)
     111             :     if (jvar == _disp_var[i])
     112             :       active = true;
     113             : 
     114             :   if (active)
     115             :   {
     116    59379712 :     if (_static_initialization && _t == _dt)
     117           0 :       return StressDivergenceTensors::computeQpOffDiagJacobian(jvar);
     118    59379712 :     else if (_dt > 0)
     119    59379712 :       return StressDivergenceTensors::computeQpOffDiagJacobian(jvar) *
     120    59379712 :              (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt);
     121             :     else
     122             :       return 0.0;
     123             :   }
     124             : 
     125             :   return 0;
     126             : }