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 "ADDynamicStressDivergenceTensors.h"
      11             : #include "ElasticityTensorTools.h"
      12             : 
      13             : registerMooseObject("TensorMechanicsApp", ADDynamicStressDivergenceTensors);
      14             : 
      15             : InputParameters
      16          81 : ADDynamicStressDivergenceTensors::validParams()
      17             : {
      18          81 :   InputParameters params = ADStressDivergenceTensors::validParams();
      19          81 :   params.addClassDescription(
      20             :       "Residual due to stress related Rayleigh damping and HHT time integration terms");
      21         162 :   params.addParam<MaterialPropertyName>("zeta",
      22         162 :                                         0.0,
      23             :                                         "Name of material property or a constant real "
      24             :                                         "number defining the zeta parameter for the "
      25             :                                         "Rayleigh damping.");
      26         162 :   params.addParam<Real>("alpha", 0, "alpha parameter for HHT time integration");
      27         162 :   params.addParam<bool>("static_initialization",
      28         162 :                         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          81 :   return params;
      34           0 : }
      35             : 
      36          27 : ADDynamicStressDivergenceTensors::ADDynamicStressDivergenceTensors(
      37          27 :     const InputParameters & parameters)
      38             :   : ADStressDivergenceTensors(parameters),
      39          27 :     _stress_older(getMaterialPropertyOlder<RankTwoTensor>(_base_name + "stress")),
      40          54 :     _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
      41          54 :     _zeta(getMaterialProperty<Real>("zeta")),
      42          54 :     _alpha(getParam<Real>("alpha")),
      43          81 :     _static_initialization(getParam<bool>("static_initialization"))
      44             : {
      45          27 : }
      46             : 
      47             : ADReal
      48     6481920 : ADDynamicStressDivergenceTensors::computeQpResidual()
      49             : {
      50             :   /**
      51             :    *This kernel needs to be used only if either Rayleigh damping or numerical damping through HHT
      52             :    *time integration scheme needs to be added to the problem through the stiffness dependent damping
      53             :    * parameter _zeta or HHT parameter _alpha, respectively.
      54             :    *
      55             :    * The residual of _zeta*K*[(1+_alpha)vel-_alpha vel_old]+ alpha K [ u - uold] + K u is required
      56             :    * = _zeta*[(1+_alpha)d/dt (Div sigma)-alpha d/dt(Div sigma_old)] +alpha [Div sigma - Div
      57             :    *sigma_old]+ Div sigma
      58             :    * = _zeta*[(1+alpha)(Div sigma - Div sigma_old)/dt - alpha (Div sigma_old - Div sigma_older)/dt]
      59             :    *   + alpha [Div sigma - Div sigma_old] +Div sigma
      60             :    * = [(1+_alpha)*_zeta/dt +_alpha+1]* Div sigma - [(1+2_alpha)*_zeta/dt + _alpha] Div sigma_old +
      61             :    *_alpha*_zeta/dt Div sigma_older
      62             :    */
      63             : 
      64             :   ADReal residual;
      65     6481920 :   if (_static_initialization && _t == _dt)
      66             :   {
      67             :     // If static initialization is true, then in the first step residual is only Ku which is
      68             :     // stress.grad(test).
      69           0 :     residual = _stress[_qp].row(_component) * _grad_test[_i][_qp];
      70             : 
      71           0 :     if (_volumetric_locking_correction)
      72             :       residual +=
      73           0 :           _stress[_qp].trace() / 3.0 * (_avg_grad_test[_i] - _grad_test[_i][_qp](_component));
      74             :   }
      75     6481920 :   else if (_dt > 0)
      76             :   {
      77             :     residual =
      78     6481920 :         _stress[_qp].row(_component) * _grad_test[_i][_qp] *
      79     6481920 :             (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt) -
      80    12963840 :         (_alpha + (1.0 + 2.0 * _alpha) * _zeta[_qp] / _dt) * _stress_old[_qp].row(_component) *
      81     6481920 :             _grad_test[_i][_qp] +
      82    19445760 :         (_alpha * _zeta[_qp] / _dt) * _stress_older[_qp].row(_component) * _grad_test[_i][_qp];
      83             : 
      84     6481920 :     if (_volumetric_locking_correction)
      85           0 :       residual += (_stress[_qp].trace() * (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt) -
      86           0 :                    (_alpha + (1.0 + 2.0 * _alpha) * _zeta[_qp] / _dt) * _stress_old[_qp].trace() +
      87           0 :                    (_alpha * _zeta[_qp] / _dt) * _stress_older[_qp].trace()) /
      88           0 :                   3.0 * (_avg_grad_test[_i] - _grad_test[_i][_qp](_component));
      89             :   }
      90             :   else
      91           0 :     residual = 0.0;
      92             : 
      93     6481920 :   return residual;
      94             : }