StressDivergence.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 "StressDivergence.h"
 
// MOOSE includes
#include "Assembly.h"
#include "Material.h"
#include "MooseMesh.h"
#include "MooseVariable.h"
#include "SymmElasticityTensor.h"
#include "SystemBase.h"
 
#include "libmesh/quadrature.h"
 
registerMooseObjectDeprecated("SolidMechanicsApp", StressDivergence, "07/30/2020 24:00");
 
template <>
InputParameters
validParams<StressDivergence>()
{
  InputParameters params = validParams<Kernel>();
  params.addRequiredParam<unsigned int>("component",
                                        "An integer corresponding to the direction "
                                        "the variable this kernel acts in. (0 for x, "
                                        "1 for y, 2 for z)");
  params.addCoupledVar("disp_x", "The x displacement");
  params.addCoupledVar("disp_y", "The y displacement");
  params.addCoupledVar("disp_z", "The z displacement");
  params.addCoupledVar("temp", "The temperature");
  params.addParam<Real>("zeta", 0.0, "Stiffness dependent Rayleigh damping coefficient");
  params.addParam<Real>("alpha", 0.0, "alpha parameter required for HHT time integration");
  params.addParam<std::string>(
      "appended_property_name", "", "Name appended to material properties to make them unique");
  params.addParam<bool>("volumetric_locking_correction",
                        true,
                        "Set to false to turn off volumetric locking correction");
 
  params.set<bool>("use_displaced_mesh") = true;
 
  return params;
}
 
StressDivergence::StressDivergence(const InputParameters & parameters)
  : Kernel(parameters),
    _stress_older(getMaterialPropertyOlder<SymmTensor>(
        "stress" + getParam<std::string>("appended_property_name"))),
    _stress_old(getMaterialPropertyOld<SymmTensor>(
        "stress" + getParam<std::string>("appended_property_name"))),
    _stress(getMaterialProperty<SymmTensor>("stress" +
                                            getParam<std::string>("appended_property_name"))),
    _Jacobian_mult(getMaterialProperty<SymmElasticityTensor>(
        "Jacobian_mult" + getParam<std::string>("appended_property_name"))),
    _d_stress_dT(getMaterialProperty<SymmTensor>("d_stress_dT" +
                                                 getParam<std::string>("appended_property_name"))),
    _component(getParam<unsigned int>("component")),
    _xdisp_coupled(isCoupled("disp_x")),
    _ydisp_coupled(isCoupled("disp_y")),
    _zdisp_coupled(isCoupled("disp_z")),
    _temp_coupled(isCoupled("temp")),
    _xdisp_var(_xdisp_coupled ? coupled("disp_x") : 0),
    _ydisp_var(_ydisp_coupled ? coupled("disp_y") : 0),
    _zdisp_var(_zdisp_coupled ? coupled("disp_z") : 0),
    _temp_var(_temp_coupled ? coupled("temp") : 0),
    _zeta(getParam<Real>("zeta")),
    _alpha(getParam<Real>("alpha")),
    _avg_grad_test(_test.size(), std::vector<Real>(3, 0.0)),
    _avg_grad_phi(_phi.size(), std::vector<Real>(3, 0.0)),
    _volumetric_locking_correction(getParam<bool>("volumetric_locking_correction"))
{
  mooseDeprecated(name(), ": StressDivergence is deprecated.\n\
The solid_mechanics module will be removed from MOOSE on July 31, 2020.\n\
Please update your input files to utilize the tensor_mechanics equivalents of\n\
models based on solid_mechanics. A detailed migration guide that was developed\n\
for BISON, but which is generally applicable to any MOOSE model is available at:\n\
https://mooseframework.org/bison/tutorials/mechanics_conversion/overview.html");
}
 
void
StressDivergence::computeResidual()
{
  prepareVectorTag(_assembly, _var.number());
 
  if (_volumetric_locking_correction)
  {
    // calculate volume averaged value of shape function derivative
    _avg_grad_test.resize(_test.size());
    for (_i = 0; _i < _test.size(); _i++)
    {
      _avg_grad_test[_i].resize(3);
      _avg_grad_test[_i][_component] = 0.0;
      for (_qp = 0; _qp < _qrule->n_points(); _qp++)
        _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
 
      _avg_grad_test[_i][_component] /= _current_elem_volume;
    }
  }
 
  precalculateResidual();
  for (_i = 0; _i < _test.size(); _i++)
    for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
 
  accumulateTaggedLocalResidual();
 
  if (_has_save_in)
  {
    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
    for (const auto & var : _save_in)
      var->sys().solution().add_vector(_local_re, var->dofIndices());
  }
}
 
Real
StressDivergence::computeQpResidual()
{
  Real residual = 0.0;
  if ((_dt > 0) && ((_zeta != 0) || (_alpha != 0)))
  {
    residual = _stress[_qp].rowDot(_component, _grad_test[_i][_qp]) *
                   (1 + _alpha + (1 + _alpha) * _zeta / _dt) -
               (_alpha + (1 + 2 * _alpha) * _zeta / _dt) *
                   _stress_old[_qp].rowDot(_component, _grad_test[_i][_qp]) +
               (_alpha * _zeta / _dt) * _stress_older[_qp].rowDot(_component, _grad_test[_i][_qp]);
 
    // volumetric locking correction
    if (_volumetric_locking_correction)
      residual += (_stress[_qp].trace() * (1 + _alpha + (1 + _alpha) * _zeta / _dt) -
                   (_alpha + (1 + 2 * _alpha) * _zeta / _dt) * _stress_old[_qp].trace() +
                   (_alpha * _zeta / _dt) * _stress_older[_qp].trace()) /
                  3.0 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
  }
  else
  {
    residual += _stress[_qp].rowDot(_component, _grad_test[_i][_qp]);
 
    // volumetric locking correction
    if (_volumetric_locking_correction)
      residual += _stress[_qp].trace() / 3.0 *
                  (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component));
  }
  return residual;
}
 
void
StressDivergence::computeJacobian()
{
  prepareMatrixTag(_assembly, _var.number(), _var.number());
 
  if (_volumetric_locking_correction)
  {
    // calculate volume averaged value of shape function derivative
    _avg_grad_test.resize(_test.size());
    for (_i = 0; _i < _test.size(); _i++)
    {
      _avg_grad_test[_i].resize(3);
      _avg_grad_test[_i][_component] = 0.0;
      for (_qp = 0; _qp < _qrule->n_points(); _qp++)
        _avg_grad_test[_i][_component] += _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
 
      _avg_grad_test[_i][_component] /= _current_elem_volume;
    }
 
    _avg_grad_phi.resize(_phi.size());
    for (_i = 0; _i < _phi.size(); _i++)
    {
      _avg_grad_phi[_i].resize(3);
      for (unsigned int component = 0; component < _mesh.dimension(); component++)
      {
        _avg_grad_phi[_i][component] = 0.0;
        for (_qp = 0; _qp < _qrule->n_points(); _qp++)
          _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
 
        _avg_grad_phi[_i][component] /= _current_elem_volume;
      }
    }
  }
 
  for (_i = 0; _i < _test.size(); _i++)
    for (_j = 0; _j < _phi.size(); _j++)
      for (_qp = 0; _qp < _qrule->n_points(); _qp++)
        _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian();
 
  accumulateTaggedLocalMatrix();
 
  if (_has_diag_save_in)
  {
    unsigned int rows = _local_ke.m();
    DenseVector<Number> diag(rows);
    for (unsigned int i = 0; i < rows; i++)
      diag(i) = _local_ke(i, i);
 
    Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
    for (const auto & var : _diag_save_in)
      var->sys().solution().add_vector(diag, var->dofIndices());
  }
}
 
Real
StressDivergence::computeQpJacobian()
{
  Real sum_C3x3 = _Jacobian_mult[_qp].sum_3x3();
  RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum_3x1();
 
  Real jacobian = 0.0;
  // B^T_i * C * B_j
  jacobian += _Jacobian_mult[_qp].stiffness(
      _component, _component, _grad_test[_i][_qp], _grad_phi[_j][_qp]); // B^T_i * C *B_j
 
  if (_volumetric_locking_correction)
  {
    // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
    // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j +  Bvol^T_i * C * B_j + B^T_i * C * Bvol_j
 
    // Bvol^T_i * C * Bvol_j
    jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
                (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 9.0;
 
    // B^T_i * C * Bvol_j
    jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
                (_avg_grad_phi[_j][_component] - _grad_phi[_j][_qp](_component)) / 3.0;
 
    // Bvol^T_i * C * B_j
    SymmTensor phi;
    if (_component == 0)
    {
      phi.xx() = _grad_phi[_j][_qp](0);
      phi.xy() = _grad_phi[_j][_qp](1);
      phi.xz() = _grad_phi[_j][_qp](2);
    }
    else if (_component == 1)
    {
      phi.yy() = _grad_phi[_j][_qp](1);
      phi.xy() = _grad_phi[_j][_qp](0);
      phi.yz() = _grad_phi[_j][_qp](2);
    }
    else if (_component == 2)
    {
      phi.zz() = _grad_phi[_j][_qp](2);
      phi.xz() = _grad_phi[_j][_qp](0);
      phi.yz() = _grad_phi[_j][_qp](1);
    }
 
    SymmTensor tmp(_Jacobian_mult[_qp] * phi);
 
    jacobian += (tmp.xx() + tmp.yy() + tmp.zz()) *
                (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
  }
 
  if (_dt > 0)
    return jacobian * (1 + _alpha + _zeta / _dt);
  else
    return jacobian;
}
 
void
StressDivergence::computeOffDiagJacobian(MooseVariableFEBase & jvar)
{
  size_t jvar_num = jvar.number();
  if (jvar_num == _var.number())
    computeJacobian();
  else
  {
    // This (undisplaced) jvar could potentially yield the wrong phi size if this object is acting
    // on the displaced mesh
    auto phi_size = _sys.getVariable(_tid, jvar.number()).dofIndices().size();
 
    if (_volumetric_locking_correction)
    {
      // calculate volume averaged value of shape function derivative
      _avg_grad_test.resize(_test.size());
      for (_i = 0; _i < _test.size(); _i++)
      {
        _avg_grad_test[_i].resize(3);
        _avg_grad_test[_i][_component] = 0.0;
        for (_qp = 0; _qp < _qrule->n_points(); _qp++)
          _avg_grad_test[_i][_component] +=
              _grad_test[_i][_qp](_component) * _JxW[_qp] * _coord[_qp];
 
        _avg_grad_test[_i][_component] /= _current_elem_volume;
      }
 
      _avg_grad_phi.resize(phi_size);
      for (_i = 0; _i < phi_size; _i++)
      {
        _avg_grad_phi[_i].resize(3);
        for (unsigned int component = 0; component < _mesh.dimension(); component++)
        {
          _avg_grad_phi[_i][component] = 0.0;
          for (_qp = 0; _qp < _qrule->n_points(); _qp++)
            _avg_grad_phi[_i][component] += _grad_phi[_i][_qp](component) * _JxW[_qp] * _coord[_qp];
 
          _avg_grad_phi[_i][component] /= _current_elem_volume;
        }
      }
    }
 
    prepareMatrixTag(_assembly, _var.number(), jvar_num);
 
    for (_i = 0; _i < _test.size(); _i++)
      for (_j = 0; _j < phi_size; _j++)
        for (_qp = 0; _qp < _qrule->n_points(); _qp++)
          _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar_num);
 
    accumulateTaggedLocalMatrix();
  }
}
 
Real
StressDivergence::computeQpOffDiagJacobian(unsigned int jvar)
{
  unsigned int coupled_component = 0;
 
  bool active = false;
 
  if (_xdisp_coupled && jvar == _xdisp_var)
  {
    coupled_component = 0;
    active = true;
  }
  else if (_ydisp_coupled && jvar == _ydisp_var)
  {
    coupled_component = 1;
    active = true;
  }
  else if (_zdisp_coupled && jvar == _zdisp_var)
  {
    coupled_component = 2;
    active = true;
  }
 
  if (active)
  {
    Real sum_C3x3 = _Jacobian_mult[_qp].sum_3x3();
    RealGradient sum_C3x1 = _Jacobian_mult[_qp].sum_3x1();
    Real jacobian = 0.0;
    jacobian += _Jacobian_mult[_qp].stiffness(
        _component, coupled_component, _grad_test[_i][_qp], _grad_phi[_j][_qp]); // B^T_i * C *B_j
 
    if (_volumetric_locking_correction)
    {
      // jacobian = Bbar^T_i * C * Bbar_j where Bbar = B + Bvol
      // jacobian = B^T_i * C * B_j + Bvol^T_i * C * Bvol_j +  Bvol^T_i * C * B_j + B^T_i * C *
      // Bvol_j
 
      // Bvol^T_i * C * Bvol_j
      jacobian += sum_C3x3 * (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) *
                  (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
                  9.0;
 
      // B^T_i * C * Bvol_j
      jacobian += sum_C3x1(_component) * _grad_test[_i][_qp](_component) *
                  (_avg_grad_phi[_j][coupled_component] - _grad_phi[_j][_qp](coupled_component)) /
                  3.0;
 
      // Bvol^T_i * C * B_i
      SymmTensor phi;
      if (coupled_component == 0)
      {
        phi.xx() = _grad_phi[_j][_qp](0);
        phi.xy() = _grad_phi[_j][_qp](1);
        phi.xz() = _grad_phi[_j][_qp](2);
      }
      else if (coupled_component == 1)
      {
        phi.yy() = _grad_phi[_j][_qp](1);
        phi.xy() = _grad_phi[_j][_qp](0);
        phi.yz() = _grad_phi[_j][_qp](2);
      }
      else if (coupled_component == 2)
      {
        phi.zz() = _grad_phi[_j][_qp](2);
        phi.xz() = _grad_phi[_j][_qp](0);
        phi.yz() = _grad_phi[_j][_qp](1);
      }
 
      SymmTensor tmp(_Jacobian_mult[_qp] * phi);
 
      jacobian += (tmp.xx() + tmp.yy() + tmp.zz()) *
                  (_avg_grad_test[_i][_component] - _grad_test[_i][_qp](_component)) / 3.0;
    }
 
    if (_dt > 0)
      return jacobian * (1 + _alpha + _zeta / _dt);
    else
      return jacobian;
  }
 
  if (_temp_coupled && jvar == _temp_var)
    return _d_stress_dT[_qp].rowDot(_component, _grad_test[_i][_qp]) * _phi[_j][_qp];
 
  return 0;
}