DynamicStressDivergenceTensors derives from StressDivergenceTensors and adds stress related Rayleigh and HHT time integration terms. More...

#include <DynamicStressDivergenceTensors.h>

Inheritance diagram for DynamicStressDivergenceTensors:

[legend]

Public Types
typedef std::vector< int >	JvarMap

Public Member Functions
	DynamicStressDivergenceTensors (const InputParameters &parameters)

virtual void	computeJacobian () override

virtual void	computeOffDiagJacobian (unsigned int jvar) override

unsigned int	mapJvarToCvar (unsigned int jvar)

int	mapJvarToCvar (unsigned int jvar, const JvarMap &jvar_map)

bool	mapJvarToCvar (unsigned int jvar, unsigned int &cvar)

const JvarMap &	getJvarMap ()

const JvarMap &	getParameterJvarMap (std::string parameter_name)

Static Public Member Functions
static InputParameters	validParams ()

Protected Member Functions
virtual Real	computeQpResidual ()

virtual Real	computeQpJacobian ()

virtual Real	computeQpOffDiagJacobian (unsigned int jvar)

virtual void	initialSetup () override

virtual void	computeResidual () override

virtual void	computeFiniteDeformJacobian ()

virtual void	computeAverageGradientTest ()

virtual void	computeAverageGradientPhi ()

Protected Attributes
const MaterialProperty< Real > &	_zeta

const Real	_alpha

const bool	_static_initialization

const std::string	_base_name
	Base name of the material system that this kernel applies to.

bool	_use_finite_deform_jacobian

const MaterialProperty< RankTwoTensor > &	_stress
	The stress tensor that the divergence operator operates on.

const MaterialProperty< RankFourTensor > &	_Jacobian_mult

std::vector< RankFourTensor >	_finite_deform_Jacobian_mult

const MaterialProperty< RankTwoTensor > *	_deformation_gradient

const MaterialProperty< RankTwoTensor > *	_deformation_gradient_old

const MaterialProperty< RankTwoTensor > *	_rotation_increment

const unsigned int	_component
	An integer corresponding to the direction this kernel acts in.

unsigned int	_ndisp
	Coupled displacement variables.

std::vector< unsigned int >	_disp_var
	Displacement variables IDs.

std::vector< std::vector< const MaterialProperty< RankTwoTensor > * > >	_deigenstrain_dargs
	eigen strain derivatives wrt coupled variables

const bool	_out_of_plane_strain_coupled

const VariableValue *const	_out_of_plane_strain

const unsigned int	_out_of_plane_strain_var

const unsigned int	_out_of_plane_direction

const bool	_use_displaced_mesh
	Whether this object is acting on the displaced mesh.

std::vector< std::vector< Real > >	_avg_grad_test
	Gradient of test function averaged over the element. Used in volumetric locking correction calculation.

std::vector< std::vector< Real > >	_avg_grad_phi
	Gradient of phi function averaged over the element. Used in volumetric locking correction calculation.

bool	_volumetric_locking_correction
	Flag for volumetric locking correction.

const unsigned int	_n_args


const MaterialProperty< RankTwoTensor > &	_stress_older
	The old and older states of the stress tensor that the divergence operator operates on.

const MaterialProperty< RankTwoTensor > &	_stress_old

Private Attributes
const std::size_t	_jvar_max_size

JvarMap	_jvar_map

std::map< std::string, JvarMap >	_jvar_local_map

Detailed Description

DynamicStressDivergenceTensors derives from StressDivergenceTensors and adds stress related Rayleigh and HHT time integration terms.

Definition at line 18 of file DynamicStressDivergenceTensors.h.

Constructor & Destructor Documentation

◆ DynamicStressDivergenceTensors()

DynamicStressDivergenceTensors::DynamicStressDivergenceTensors ( const InputParameters & parameters )

Definition at line 36 of file DynamicStressDivergenceTensors.C.

  : 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"))
{
}

Member Function Documentation

◆ computeAverageGradientPhi()

void StressDivergenceTensors::computeAverageGradientPhi ( )

protectedvirtualinherited

Reimplemented in StressDivergenceRZTensors.

Definition at line 443 of file StressDivergenceTensors.C.

{
  // Calculate volume average derivatives for phi
  _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;
    }
  }
}

Referenced by StressDivergenceTensors::computeJacobian(), and StressDivergenceTensors::computeOffDiagJacobian().

◆ computeAverageGradientTest()

void StressDivergenceTensors::computeAverageGradientTest ( )

protectedvirtualinherited

Reimplemented in StressDivergenceRZTensors.

Definition at line 427 of file StressDivergenceTensors.C.

{
  // 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;
  }
}

Referenced by StressDivergenceTensors::computeJacobian(), StressDivergenceTensors::computeOffDiagJacobian(), and StressDivergenceTensors::computeResidual().

◆ computeFiniteDeformJacobian()

void StressDivergenceTensors::computeFiniteDeformJacobian ( )

protectedvirtualinherited

Definition at line 369 of file StressDivergenceTensors.C.

{
  usingTensorIndices(i_, j_, k_, l_);
  const auto I = RankTwoTensor::Identity();
  const RankFourTensor I2 = I.times<i_, k_, j_, l_>(I);
 
  // Bring back to unrotated config
  const RankTwoTensor unrotated_stress =
      (*_rotation_increment)[_qp].transpose() * _stress[_qp] * (*_rotation_increment)[_qp];
 
  // Incremental deformation gradient Fhat
  const RankTwoTensor Fhat =
      (*_deformation_gradient)[_qp] * (*_deformation_gradient_old)[_qp].inverse();
  const RankTwoTensor Fhatinv = Fhat.inverse();
 
  const RankTwoTensor rot_times_stress = (*_rotation_increment)[_qp] * unrotated_stress;
  const RankFourTensor dstress_drot =
      I.times<i_, k_, j_, l_>(rot_times_stress) + I.times<j_, k_, i_, l_>(rot_times_stress);
  const RankFourTensor rot_rank_four =
      (*_rotation_increment)[_qp].times<i_, k_, j_, l_>((*_rotation_increment)[_qp]);
  const RankFourTensor drot_dUhatinv = Fhat.times<i_, k_, j_, l_>(I);
 
  const RankTwoTensor A = I - Fhatinv;
 
  // Ctilde = Chat^-1 - I
  const RankTwoTensor Ctilde = A * A.transpose() - A - A.transpose();
  const RankFourTensor dCtilde_dFhatinv =
      -I.times<i_, k_, j_, l_>(A) - I.times<j_, k_, i_, l_>(A) + I2 + I.times<j_, k_, i_, l_>(I);
 
  // Second order approximation of Uhat - consistent with strain increment definition
  // const RankTwoTensor Uhat = I - 0.5 * Ctilde - 3.0/8.0 * Ctilde * Ctilde;
 
  RankFourTensor dUhatinv_dCtilde =
      0.5 * I2 - 1.0 / 8.0 * (I.times<i_, k_, j_, l_>(Ctilde) + Ctilde.times<i_, k_, j_, l_>(I));
  RankFourTensor drot_dFhatinv = drot_dUhatinv * dUhatinv_dCtilde * dCtilde_dFhatinv;
 
  drot_dFhatinv -= Fhat.times<i_, k_, j_, l_>((*_rotation_increment)[_qp].transpose());
  _finite_deform_Jacobian_mult[_qp] = dstress_drot * drot_dFhatinv;
 
  const RankFourTensor dstrain_increment_dCtilde =
      -0.5 * I2 + 0.25 * (I.times<i_, k_, j_, l_>(Ctilde) + Ctilde.times<i_, k_, j_, l_>(I));
  _finite_deform_Jacobian_mult[_qp] +=
      rot_rank_four * _Jacobian_mult[_qp] * dstrain_increment_dCtilde * dCtilde_dFhatinv;
  _finite_deform_Jacobian_mult[_qp] += Fhat.times<j_, k_, i_, l_>(_stress[_qp]);
 
  const RankFourTensor dFhat_dFhatinv = -Fhat.times<i_, k_, j_, l_>(Fhat.transpose());
  const RankTwoTensor dJ_dFhatinv = dFhat_dFhatinv.innerProductTranspose(Fhat.ddet());
 
  // Component from Jacobian derivative
  _finite_deform_Jacobian_mult[_qp] += _stress[_qp].times<i_, j_, k_, l_>(dJ_dFhatinv);
 
  // Derivative of Fhatinv w.r.t. undisplaced coordinates
  const RankTwoTensor Finv = (*_deformation_gradient)[_qp].inverse();
  const RankFourTensor dFhatinv_dGradu = -Fhatinv.times<i_, k_, j_, l_>(Finv.transpose());
  _finite_deform_Jacobian_mult[_qp] = _finite_deform_Jacobian_mult[_qp] * dFhatinv_dGradu;
}

Referenced by StressDivergenceTensors::computeJacobian(), and StressDivergenceTensors::computeOffDiagJacobian().

◆ computeJacobian()

void StressDivergenceTensors::computeJacobian ( )

overridevirtualinherited

Definition at line 178 of file StressDivergenceTensors.C.

{
  if (_volumetric_locking_correction)
  {
    computeAverageGradientTest();
    computeAverageGradientPhi();
  }
 
  if (_use_finite_deform_jacobian)
  {
    _finite_deform_Jacobian_mult.resize(_qrule->n_points());
 
    for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      computeFiniteDeformJacobian();
 
    ALEKernel::computeJacobian();
  }
  else
    Kernel::computeJacobian();
}

◆ computeOffDiagJacobian()

void StressDivergenceTensors::computeOffDiagJacobian ( unsigned int jvar )

overridevirtualinherited

Reimplemented from JvarMapKernelInterface< ALEKernel >.

Definition at line 200 of file StressDivergenceTensors.C.

{
  if (_volumetric_locking_correction)
  {
    computeAverageGradientPhi();
    computeAverageGradientTest();
  }
 
  if (_use_finite_deform_jacobian)
  {
    _finite_deform_Jacobian_mult.resize(_qrule->n_points());
 
    for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      computeFiniteDeformJacobian();
 
    ALEKernel::computeOffDiagJacobian(jvar);
  }
  else
    Kernel::computeOffDiagJacobian(jvar);
}

◆ computeQpJacobian()

Real DynamicStressDivergenceTensors::computeQpJacobian ( )

protectedvirtual

Reimplemented from StressDivergenceTensors.

Definition at line 94 of file DynamicStressDivergenceTensors.C.

{
  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;
}

◆ computeQpOffDiagJacobian()

Real DynamicStressDivergenceTensors::computeQpOffDiagJacobian ( unsigned int jvar )

protectedvirtual

Reimplemented from StressDivergenceTensors.

Definition at line 106 of file DynamicStressDivergenceTensors.C.

{
  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;
}

◆ computeQpResidual()

Real DynamicStressDivergenceTensors::computeQpResidual ( )

protectedvirtual

This kernel needs to be used only if either Rayleigh damping or numerical damping through HHT time integration scheme needs to be added to the problem through the stiffness dependent damping parameter _zeta or HHT parameter _alpha, respectively.

The residual of _zeta*K*[(1+_alpha)vel-_alpha vel_old]+ alpha K [ u - uold] + K u is required = _zeta*[(1+_alpha)d/dt (Div sigma)-alpha d/dt(Div sigma_old)] +alpha [Div sigma - Div sigma_old]+ Div sigma = _zeta*[(1+alpha)(Div sigma - Div sigma_old)/dt - alpha (Div sigma_old - Div sigma_older)/dt]

alpha [Div sigma - Div sigma_old] +Div sigma = [(1+_alpha)*_zeta/dt +_alpha+1]* Div sigma - [(1+2_alpha)*_zeta/dt + _alpha] Div sigma_old + _alpha*_zeta/dt Div sigma_older

Reimplemented from StressDivergenceTensors.

Definition at line 47 of file DynamicStressDivergenceTensors.C.

{
  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;
}

◆ computeResidual()

void StressDivergenceTensors::computeResidual ( )

overrideprotectedvirtualinherited

Definition at line 137 of file StressDivergenceTensors.C.

{
  prepareVectorTag(_assembly, _var.number());
 
  if (_volumetric_locking_correction)
    computeAverageGradientTest();
 
  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());
  }
}

◆ initialSetup()

void StressDivergenceTensors::initialSetup ( )

overrideprotectedvirtualinherited

Reimplemented in StressDivergenceRSphericalTensors, and StressDivergenceRZTensors.

Definition at line 124 of file StressDivergenceTensors.C.

{
  // check if any of the eigenstrains provide derivatives wrt variables that are not coupled
  for (auto eigenstrain_name : getParam<std::vector<MaterialPropertyName>>("eigenstrain_names"))
    validateNonlinearCoupling<RankTwoTensor>(eigenstrain_name);
 
  // make sure the coordinate system is cartesioan
  if (getBlockCoordSystem() != Moose::COORD_XYZ)
    mooseError("The coordinate system in the Problem block must be set to XYZ for cartesian "
               "geometries.");
}

◆ validParams()

InputParameters DynamicStressDivergenceTensors::validParams ( )

static

Definition at line 16 of file DynamicStressDivergenceTensors.C.

{
  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;
}