DynamicStressDivergenceTensors.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "DynamicStressDivergenceTensors.h"
#include "ElasticityTensorTools.h"

registerMooseObject("SolidMechanicsApp", DynamicStressDivergenceTensors);

InputParameters
DynamicStressDivergenceTensors::validParams()
{
  InputParameters params = StressDivergenceTensors::validParams();
  params.addClassDescription(
      "Residual due to stress related Rayleigh damping and HHT time integration terms ");
  params.addParam<MaterialPropertyName>("zeta",
                                        0.0,
                                        "Name of material property or a constant real "
                                        "number defining the zeta parameter for the "
                                        "Rayleigh damping.");
  params.addParam<Real>("alpha", 0, "alpha parameter for HHT time integration");
  params.addParam<bool>("static_initialization",
                        false,
                        "Set to true to get the system to "
                        "equilibrium under gravity by running a "
                        "quasi-static analysis (by solving Ku = F) "
                        "in the first time step");
  return params;
}

DynamicStressDivergenceTensors::DynamicStressDivergenceTensors(const InputParameters & parameters)
  : StressDivergenceTensors(parameters),
    _stress_older(getMaterialPropertyOlderByName<RankTwoTensor>(_base_name + "stress")),
    _stress_old(getMaterialPropertyOldByName<RankTwoTensor>(_base_name + "stress")),
    _zeta(getMaterialProperty<Real>("zeta")),
    _alpha(getParam<Real>("alpha")),
    _static_initialization(getParam<bool>("static_initialization"))
{
}

Real
DynamicStressDivergenceTensors::computeQpResidual()
{
  Real residual = 0.0;
  if (_static_initialization && _t == _dt)
  {
    // If static initialization is true, then in the first step residual is only Ku which is
    // stress.grad(test).
    residual += _stress[_qp].row(_component) * _grad_test[_i][_qp];

    if (_volumetric_locking_correction)
      residual += _stress[_qp].trace() / 3.0 *
                  (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
  }
  else if (_dt > 0)
  {
    residual +=
        _stress[_qp].row(_component) * _grad_test[_i][_qp] *
            (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt) -
        (_alpha + (1.0 + 2.0 * _alpha) * _zeta[_qp] / _dt) * _stress_old[_qp].row(_component) *
            _grad_test[_i][_qp] +
        (_alpha * _zeta[_qp] / _dt) * _stress_older[_qp].row(_component) * _grad_test[_i][_qp];

    if (_volumetric_locking_correction)
      residual += (_stress[_qp].trace() * (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt) -
                   (_alpha + (1.0 + 2.0 * _alpha) * _zeta[_qp] / _dt) * _stress_old[_qp].trace() +
                   (_alpha * _zeta[_qp] / _dt) * _stress_older[_qp].trace()) /
                  3.0 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
  }

  return residual;
}

Real
DynamicStressDivergenceTensors::computeQpJacobian()
{
  if (_static_initialization && _t == _dt)
    return StressDivergenceTensors::computeQpJacobian();
  else if (_dt > 0)
    return StressDivergenceTensors::computeQpJacobian() *
           (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt);
  else
    return 0.0;
}

Real
DynamicStressDivergenceTensors::computeQpOffDiagJacobian(unsigned int jvar)
{
  bool active = true;

  for (unsigned int i = 0; i < _ndisp; ++i)
    if (jvar == _disp_var[i])
      active = true;

  if (active)
  {
    if (_static_initialization && _t == _dt)
      return StressDivergenceTensors::computeQpOffDiagJacobian(jvar);
    else if (_dt > 0)
      return StressDivergenceTensors::computeQpOffDiagJacobian(jvar) *
             (1.0 + _alpha + (1.0 + _alpha) * _zeta[_qp] / _dt);
    else
      return 0.0;
  }

  return 0;
}