This class uses the stress update material in an anisotropic return mapping. More...

#include <HillElastoPlasticityStressUpdate.h>

Inheritance diagram for HillElastoPlasticityStressUpdateTempl< is_ad >:

Public Member Functions
	HillElastoPlasticityStressUpdateTempl (const InputParameters &parameters)

virtual Real	computeStrainEnergyRateDensity (const GenericMaterialProperty< RankTwoTensor, is_ad > &stress, const GenericMaterialProperty< RankTwoTensor, is_ad > &strain_rate) override

Static Public Member Functions
static InputParameters	validParams ()

Protected Types
enum class	Axis { X , Y , Z }

Protected Member Functions
virtual void	initQpStatefulProperties () override

virtual void	computeStressInitialize (const GenericDenseVector< is_ad > &stress_dev, const GenericDenseVector< is_ad > &stress, const GenericRankFourTensor< is_ad > &elasticity_tensor) override

void	computeElasticityTensorEigenDecomposition ()
	Compute eigendecomposition of element-wise elasticity tensor.

virtual GenericReal< is_ad >	computeResidual (const GenericDenseVector< is_ad > &effective_trial_stress, const GenericDenseVector< is_ad > &stress_new, const GenericReal< is_ad > &scalar) override

virtual GenericReal< is_ad >	computeDerivative (const GenericDenseVector< is_ad > &effective_trial_stress, const GenericDenseVector< is_ad > &stress_new, const GenericReal< is_ad > &scalar) override

virtual Real	computeReferenceResidual (const GenericDenseVector< is_ad > &effective_trial_stress, const GenericDenseVector< is_ad > &stress_new, const GenericReal< is_ad > &residual, const GenericReal< is_ad > &scalar_effective_inelastic_strain) override

virtual void	propagateQpStatefulProperties () override

bool	requiresIsotropicTensor () override
	Does the model require the elasticity tensor to be isotropic? No, this class does anisotropic elasto-plasticity

Real	computeHardeningDerivative ()

GenericReal< is_ad >	computeHardeningValue (const GenericReal< is_ad > &scalar, const GenericReal< is_ad > &omega)

void	computeHillTensorEigenDecomposition (const DenseMatrix< GenericReal< is_ad > > &hill_tensor)
	Compute eigendecomposition of Hill's tensor for anisotropic plasticity.

virtual void	computeStrainFinalize (GenericRankTwoTensor< is_ad > &inelasticStrainIncrement, const GenericRankTwoTensor< is_ad > &stress, const GenericDenseVector< is_ad > &stress_dev, const GenericReal< is_ad > &delta_gamma) override
	Perform any necessary steps to finalize strain increment after return mapping iterations.

virtual void	computeStressFinalize (const GenericRankTwoTensor< is_ad > &inelasticStrainIncrement, const GenericReal< is_ad > &delta_gamma, GenericRankTwoTensor< is_ad > &stress, const GenericDenseVector< is_ad > &stress_dev, const GenericRankTwoTensor< is_ad > &stress_old, const GenericRankFourTensor< is_ad > &elasticity_tensor) override
	Perform any necessary steps to finalize state after return mapping iterations.

GenericReal< is_ad >	computeOmega (const GenericReal< is_ad > &delta_gamma, const GenericDenseVector< is_ad > &stress_trial)

void	computeDeltaDerivatives (const GenericReal< is_ad > &delta_gamma, const GenericDenseVector< is_ad > &stress_trial, const GenericReal< is_ad > &sy_alpha, GenericReal< is_ad > &omega, GenericReal< is_ad > &omega_gamma, GenericReal< is_ad > &sy_gamma)

void	paramError (const std::string &param, Args... args) const

virtual GenericReal< is_ad >	computeStressDerivative (const Real, const Real) override
	Calculate the derivative of the strain increment with respect to the updated stress.

Protected Attributes
	usingTransientInterfaceMembers

GenericReal< is_ad >	_qsigma
	Square of the q function for orthotropy.

GenericDenseVector< is_ad >	_eigenvalues_hill

GenericDenseMatrix< is_ad >	_eigenvectors_hill

GenericDenseMatrix< is_ad >	_anisotropic_elastic_tensor

const std::string	_elasticity_tensor_name
	Name of the elasticity tensor material property.

const GenericMaterialProperty< RankFourTensor, is_ad > &	_elasticity_tensor
	Anisotropic elasticity tensor material property.

const Real	_hardening_constant

const Real	_hardening_exponent

GenericMaterialProperty< Real, is_ad > &	_hardening_variable

const MaterialProperty< Real > &	_hardening_variable_old

GenericMaterialProperty< DenseMatrix< Real >, is_ad > &	_elasticity_eigenvectors

GenericMaterialProperty< DenseVector< Real >, is_ad > &	_elasticity_eigenvalues

GenericMaterialProperty< DenseMatrix< Real >, is_ad > &	_b_eigenvectors

GenericMaterialProperty< DenseVector< Real >, is_ad > &	_b_eigenvalues

GenericMaterialProperty< DenseMatrix< Real >, is_ad > &	_alpha_matrix

GenericMaterialProperty< DenseVector< Real >, is_ad > &	_sigma_tilde

GenericMaterialProperty< DenseVector< Real >, is_ad > &	_sigma_tilde_rotated

GenericReal< is_ad >	_hardening_derivative

GenericReal< is_ad >	_yield_condition

GenericReal< is_ad >	_yield_stress

const MaterialProperty< DenseMatrix< Real > > &	_hill_tensor
	Hill tensor, when global axes do not (somehow) align with those of the material Example: Large rotation due to rigid body and/or large deformation kinematics.

MaterialProperty< RankTwoTensor > &	_rotation_matrix

MaterialProperty< RankTwoTensor > &	_rotation_matrix_transpose

const MaterialProperty< RankTwoTensor > &	_rotation_matrix_old

const MaterialProperty< RankTwoTensor > &	_rotation_matrix_transpose_old

const bool	_local_cylindrical_csys

enum HillElastoPlasticityStressUpdateTempl::Axis	_axis

const Elem *const &	_current_elem

const MooseArray< Point > &	_q_point

unsigned int	_qp

GenericMaterialProperty< RankTwoTensor, is_ad > &	_plasticity_strain
	Plasticity strain tensor material property.

const MaterialProperty< RankTwoTensor > &	_plasticity_strain_old

Detailed Description

template<bool is_ad>
class HillElastoPlasticityStressUpdateTempl< is_ad >

This class uses the stress update material in an anisotropic return mapping.

This class is one of the generalized radial return constitutive models based on Hill's criterion; it assumes and anisotropic elasticity tensor and an anisotropic plastic yield surface. Constitutive models that combine creep and plasticity can be used.

Definition at line 22 of file HillElastoPlasticityStressUpdate.h.

Member Enumeration Documentation

◆ Axis

template<bool is_ad>

enum class HillElastoPlasticityStressUpdateTempl::Axis

strongprotected

Enumerator
X
Y
Z

Definition at line 166 of file HillElastoPlasticityStressUpdate.h.

166{ X, Y, Z } _axis;

HillElastoPlasticityStressUpdateTempl::Axis::X

@ X

HillElastoPlasticityStressUpdateTempl::Axis::Z

@ Z

HillElastoPlasticityStressUpdateTempl::Axis::Y

@ Y

HillElastoPlasticityStressUpdateTempl::_axis

enum HillElastoPlasticityStressUpdateTempl::Axis _axis

Constructor & Destructor Documentation

◆ HillElastoPlasticityStressUpdateTempl()

template<bool is_ad>

HillElastoPlasticityStressUpdateTempl< is_ad >::HillElastoPlasticityStressUpdateTempl ( const InputParameters & parameters )

Definition at line 47 of file HillElastoPlasticityStressUpdate.C.

  : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>(parameters),
    _qsigma(0.0),
    _eigenvalues_hill(6),
    _eigenvectors_hill(6, 6),
    _anisotropic_elastic_tensor(6, 6),
    _elasticity_tensor_name(this->_base_name + "elasticity_tensor"),
    _elasticity_tensor(
        this->template getGenericMaterialProperty<RankFourTensor, is_ad>(_elasticity_tensor_name)),
    _hardening_constant(this->template getParam<Real>("hardening_constant")),
    _hardening_exponent(this->template getParam<Real>("hardening_exponent")),
    _hardening_variable(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
                                                                           "hardening_variable")),
    _hardening_variable_old(
        this->template getMaterialPropertyOld<Real>(this->_base_name + "hardening_variable")),
    _elasticity_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
        this->_base_name + "elasticity_eigenvectors")),
    _elasticity_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
        this->_base_name + "elasticity_eigenvalues")),
    _b_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
        this->_base_name + "b_eigenvectors")),
    _b_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
        this->_base_name + "b_eigenvalues")),
    _alpha_matrix(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(this->_base_name +
                                                                                  "alpha_matrix")),
    _sigma_tilde(this->template declareGenericProperty<DenseVector<Real>, is_ad>(this->_base_name +
                                                                                 "sigma_tilde")),
    _sigma_tilde_rotated(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
        this->_base_name + "sigma_tilde_rotated")),
    _hardening_derivative(0.0),
    _yield_condition(1.0),
    _yield_stress(this->template getParam<Real>("yield_stress")),
    _hill_tensor(this->template getMaterialPropertyByName<DenseMatrix<Real>>(this->_base_name +
                                                                             "hill_tensor")),
    _rotation_matrix(
        this->template declareProperty<RankTwoTensor>(this->_base_name + "rotation_matrix")),
    _rotation_matrix_transpose(this->template declareProperty<RankTwoTensor>(
        this->_base_name + "rotation_matrix_transpose")),
    _rotation_matrix_old(
        this->template getMaterialPropertyOld<RankTwoTensor>(this->_base_name + "rotation_matrix")),
    _rotation_matrix_transpose_old(this->template getMaterialPropertyOld<RankTwoTensor>(
        this->_base_name + "rotation_matrix_transpose")),
    _local_cylindrical_csys(this->template getParam<bool>("local_cylindrical_csys")),
    _axis(this->template getParam<MooseEnum>("axis").template getEnum<Axis>())
{
  if (_local_cylindrical_csys && !parameters.isParamSetByUser("axis"))
  {
    this->paramError(
        "local_cylindrical_csys",
        "\nIf parameter local_cylindrical_csys is provided then parameter axis should be also "
        "provided.");
  }
 
  if (!_local_cylindrical_csys && parameters.isParamSetByUser("axis"))
  {
    this->paramError("axis",
                     "\nIf parameter axis is provided then the parameter local_cylindrical_csys "
                     "should be set to true");
  }
}

Member Function Documentation

◆ computeDeltaDerivatives()

template<bool is_ad>

void HillElastoPlasticityStressUpdateTempl< is_ad >::computeDeltaDerivatives	(	const GenericReal< is_ad > &	delta_gamma,
		const GenericDenseVector< is_ad > &	stress_trial,
		const GenericReal< is_ad > &	sy_alpha,
		GenericReal< is_ad > &	omega,
		GenericReal< is_ad > &	omega_gamma,
		GenericReal< is_ad > &	sy_gamma
	)

protected

Definition at line 465 of file HillElastoPlasticityStressUpdate.C.

{
  omega_gamma = 0.0;
  sy_gamma = 0.0;
 
  GenericDenseVector<is_ad> K_deltaGamma(6);
 
  if (_local_cylindrical_csys)
    omega = computeOmega(delta_gamma, _sigma_tilde_rotated[_qp]);
  else
    omega = computeOmega(delta_gamma, _sigma_tilde[_qp]);
 
  GenericDenseVector<is_ad> K(6);
  for (unsigned int i = 0; i < 6; i++)
    K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
 
  for (unsigned int i = 0; i < 6; i++)
    K_deltaGamma(i) =
        -2.0 * _b_eigenvalues[_qp](i) * K(i) / (1 + _b_eigenvalues[_qp](i) * delta_gamma);
 
  for (unsigned int i = 0; i < 6; i++)
  {
    if (_local_cylindrical_csys)
      omega_gamma += K_deltaGamma(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
    else
      omega_gamma += K_deltaGamma(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
  }
 
  omega_gamma /= 4.0 * omega;
  sy_gamma = 2.0 * sy_alpha * (omega + delta_gamma * omega_gamma);
}

◆ computeDerivative()

template<bool is_ad>

GenericReal< is_ad > HillElastoPlasticityStressUpdateTempl< is_ad >::computeDerivative	(	const GenericDenseVector< is_ad > &	effective_trial_stress,
		const GenericDenseVector< is_ad > &	stress_new,
		const GenericReal< is_ad > &	scalar
	)

overrideprotectedvirtual

Definition at line 435 of file HillElastoPlasticityStressUpdate.C.

{
  using std::pow;
 
  // If in elastic regime, return unit derivative
  if (_yield_condition <= 0.0)
    return 1.0;
 
  GenericReal<is_ad> omega = computeOmega(delta_gamma, stress_sigma);
  _hardening_derivative = computeHardeningDerivative();
 
  GenericReal<is_ad> hardeningVariable = computeHardeningValue(delta_gamma, omega);
  GenericReal<is_ad> sy =
      _hardening_constant * pow(hardeningVariable + 1.0e-30, _hardening_exponent) + _yield_stress;
  GenericReal<is_ad> sy_alpha = _hardening_derivative;
 
  GenericReal<is_ad> omega_gamma;
  GenericReal<is_ad> sy_gamma;
 
  computeDeltaDerivatives(delta_gamma, stress_sigma, sy_alpha, omega, omega_gamma, sy_gamma);
  GenericReal<is_ad> residual_derivative = 1 / omega * (sy_gamma - 1 / omega * omega_gamma * sy);
 
  return residual_derivative;
}

◆ computeElasticityTensorEigenDecomposition()

template<bool is_ad>

void HillElastoPlasticityStressUpdateTempl< is_ad >::computeElasticityTensorEigenDecomposition ( )

protected

Compute eigendecomposition of element-wise elasticity tensor.

Definition at line 218 of file HillElastoPlasticityStressUpdate.C.

{
  // Helper method to compute qp matrices required for the elasto-plasticity algorithm.
  using std::sqrt;
 
  GenericRankFourTensor<is_ad> elasticity_tensor = _elasticity_tensor[_qp];
 
  if (_local_cylindrical_csys)
  {
    elasticity_tensor.rotate(_rotation_matrix_transpose[_qp]);
  }
 
  ElasticityTensorTools::toVoigtNotation<is_ad>(_anisotropic_elastic_tensor, elasticity_tensor);
 
  const unsigned int dimension = _anisotropic_elastic_tensor.n();
 
  AnisotropyMatrixReal A = AnisotropyMatrixReal::Zero();
  for (unsigned int index_i = 0; index_i < dimension; index_i++)
    for (unsigned int index_j = 0; index_j < dimension; index_j++)
      A(index_i, index_j) = MetaPhysicL::raw_value(_anisotropic_elastic_tensor(index_i, index_j));
 
  Eigen::SelfAdjointEigenSolver<AnisotropyMatrixReal> es(A);
 
  auto lambda = es.eigenvalues();
  auto v = es.eigenvectors();
 
  for (unsigned int index_i = 0; index_i < dimension; index_i++)
    _elasticity_eigenvalues[_qp](index_i) = lambda(index_i);
 
  for (unsigned int index_i = 0; index_i < dimension; index_i++)
    for (unsigned int index_j = 0; index_j < dimension; index_j++)
      _elasticity_eigenvectors[_qp](index_i, index_j) = v(index_i, index_j);
 
  // Compute sqrt(Delta_c) * QcT * A * Qc * sqrt(Delta_c)
  GenericDenseMatrix<is_ad> sqrt_Delta(6, 6);
  for (unsigned int i = 0; i < 6; i++)
  {
    sqrt_Delta(i, i) = sqrt(_elasticity_eigenvalues[_qp](i));
  }
 
  GenericDenseMatrix<is_ad> eigenvectors_elasticity_transpose(6, 6);
  _elasticity_eigenvectors[_qp].get_transpose(eigenvectors_elasticity_transpose);
 
  GenericDenseMatrix<is_ad> b_matrix(6, 6);
  b_matrix = _hill_tensor[_qp];
 
  // Right multiply by matrix of eigenvectors transpose
  b_matrix.right_multiply(_elasticity_eigenvectors[_qp]);
  b_matrix.right_multiply(sqrt_Delta);
  b_matrix.left_multiply(eigenvectors_elasticity_transpose);
  b_matrix.left_multiply(sqrt_Delta);
 
  AnisotropyMatrixReal B_eigen = AnisotropyMatrixReal::Zero();
  for (unsigned int index_i = 0; index_i < dimension; index_i++)
    for (unsigned int index_j = 0; index_j < dimension; index_j++)
      B_eigen(index_i, index_j) = MetaPhysicL::raw_value(b_matrix(index_i, index_j));
 
  if (isBlockDiagonal(B_eigen))
  {
    Eigen::SelfAdjointEigenSolver<AnisotropyMatrixRealBlock> es_b_left(B_eigen.block<3, 3>(0, 0));
    Eigen::SelfAdjointEigenSolver<AnisotropyMatrixRealBlock> es_b_right(B_eigen.block<3, 3>(3, 3));
 
    auto lambda_b_left = es_b_left.eigenvalues();
    auto v_b_left = es_b_left.eigenvectors();
 
    auto lambda_b_right = es_b_right.eigenvalues();
    auto v_b_right = es_b_right.eigenvectors();
 
    _b_eigenvalues[_qp](0) = lambda_b_left(0);
    _b_eigenvalues[_qp](1) = lambda_b_left(1);
    _b_eigenvalues[_qp](2) = lambda_b_left(2);
    _b_eigenvalues[_qp](3) = lambda_b_right(0);
    _b_eigenvalues[_qp](4) = lambda_b_right(1);
    _b_eigenvalues[_qp](5) = lambda_b_right(2);
 
    _b_eigenvectors[_qp](0, 0) = v_b_left(0, 0);
    _b_eigenvectors[_qp](0, 1) = v_b_left(0, 1);
    _b_eigenvectors[_qp](0, 2) = v_b_left(0, 2);
    _b_eigenvectors[_qp](1, 0) = v_b_left(1, 0);
    _b_eigenvectors[_qp](1, 1) = v_b_left(1, 1);
    _b_eigenvectors[_qp](1, 2) = v_b_left(1, 2);
    _b_eigenvectors[_qp](2, 0) = v_b_left(2, 0);
    _b_eigenvectors[_qp](2, 1) = v_b_left(2, 1);
    _b_eigenvectors[_qp](2, 2) = v_b_left(2, 2);
    _b_eigenvectors[_qp](3, 3) = v_b_right(0, 0);
    _b_eigenvectors[_qp](3, 4) = v_b_right(0, 1);
    _b_eigenvectors[_qp](3, 5) = v_b_right(0, 2);
    _b_eigenvectors[_qp](4, 3) = v_b_right(1, 0);
    _b_eigenvectors[_qp](4, 4) = v_b_right(1, 1);
    _b_eigenvectors[_qp](4, 5) = v_b_right(1, 2);
    _b_eigenvectors[_qp](5, 3) = v_b_right(2, 0);
    _b_eigenvectors[_qp](5, 4) = v_b_right(2, 1);
    _b_eigenvectors[_qp](5, 5) = v_b_right(2, 2);
  }
  else
  {
    Eigen::SelfAdjointEigenSolver<AnisotropyMatrixReal> es_b(B_eigen);
 
    auto lambda_b = es_b.eigenvalues();
    auto v_b = es_b.eigenvectors();
    for (unsigned int index_i = 0; index_i < dimension; index_i++)
      _b_eigenvalues[_qp](index_i) = lambda_b(index_i);
 
    for (unsigned int index_i = 0; index_i < dimension; index_i++)
      for (unsigned int index_j = 0; index_j < dimension; index_j++)
        _b_eigenvectors[_qp](index_i, index_j) = v_b(index_i, index_j);
  }
 
  _alpha_matrix[_qp] = sqrt_Delta;
  _alpha_matrix[_qp].right_multiply(_b_eigenvectors[_qp]);
  _alpha_matrix[_qp].left_multiply(_elasticity_eigenvectors[_qp]);
}

◆ computeHardeningDerivative()

template<bool is_ad>

Real HillElastoPlasticityStressUpdateTempl< is_ad >::computeHardeningDerivative ( )

protected

Definition at line 513 of file HillElastoPlasticityStressUpdate.C.

{
  using std::pow;
  return _hardening_constant * _hardening_exponent *
         MetaPhysicL::raw_value(pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent - 1));
}

◆ computeHardeningValue()

template<bool is_ad>

GenericReal< is_ad > HillElastoPlasticityStressUpdateTempl< is_ad >::computeHardeningValue	(	const GenericReal< is_ad > &	scalar,
		const GenericReal< is_ad > &	omega
	)

protected

Definition at line 505 of file HillElastoPlasticityStressUpdate.C.

{
  return _hardening_variable_old[_qp] + 2.0 * delta_gamma * omega;
}

◆ computeHillTensorEigenDecomposition()

template<bool is_ad>

void HillElastoPlasticityStressUpdateTempl< is_ad >::computeHillTensorEigenDecomposition ( const DenseMatrix< GenericReal< is_ad > > & hill_tensor )

protected

Compute eigendecomposition of Hill's tensor for anisotropic plasticity.

Parameters

hill_tensor 6x6 matrix representing fourth order Hill's tensor describing anisotropy

◆ computeOmega()

template<bool is_ad>

GenericReal< is_ad > HillElastoPlasticityStressUpdateTempl< is_ad >::computeOmega	(	const GenericReal< is_ad > &	delta_gamma,
		const GenericDenseVector< is_ad > &	stress_trial
	)

protected

Definition at line 333 of file HillElastoPlasticityStressUpdate.C.

{
  GenericDenseVector<is_ad> K(6);
  GenericReal<is_ad> omega = 0.0;
 
  for (unsigned int i = 0; i < 6; i++)
  {
    K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
 
    if (_local_cylindrical_csys)
      omega += K(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
    else
      omega += K(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
  }
  omega *= 0.5;
 
  if (omega == 0.0)
    omega = 1.0e-36;
 
  using std::sqrt;
  return sqrt(omega);
}

◆ computeReferenceResidual()

template<bool is_ad>

Real HillElastoPlasticityStressUpdateTempl< is_ad >::computeReferenceResidual	(	const GenericDenseVector< is_ad > &	effective_trial_stress,
		const GenericDenseVector< is_ad > &	stress_new,
		const GenericReal< is_ad > &	residual,
		const GenericReal< is_ad > &	scalar_effective_inelastic_strain
	)

overrideprotectedvirtual

Definition at line 359 of file HillElastoPlasticityStressUpdate.C.

{
  // Equation is already normalized.
  return 1.0;
}

◆ computeResidual()

template<bool is_ad>

GenericReal< is_ad > HillElastoPlasticityStressUpdateTempl< is_ad >::computeResidual	(	const GenericDenseVector< is_ad > &	effective_trial_stress,
		const GenericDenseVector< is_ad > &	stress_new,
		const GenericReal< is_ad > &	scalar
	)

overrideprotectedvirtual

Definition at line 371 of file HillElastoPlasticityStressUpdate.C.

{
  using std::pow;
 
  // If in elastic regime, just return
  if (_yield_condition <= 0.0)
    return 0.0;
 
  GenericReal<is_ad> omega;
 
  if (_local_cylindrical_csys)
  {
    // Get stress_tilde_rotated
    GenericRankTwoTensor<is_ad> stress;
    RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_sigma, stress);
 
    stress.rotate(_rotation_matrix_transpose[_qp]);
 
    GenericDenseVector<is_ad> stress_sigma_rotated(6);
    RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_sigma_rotated);
 
    GenericDenseVector<is_ad> stress_tilde_rotated(6);
    GenericDenseMatrix<is_ad> alpha_temp_rotated(_alpha_matrix[_qp]);
    alpha_temp_rotated.lu_solve(stress_sigma_rotated, stress_tilde_rotated);
 
    // Material property used in computeStressFinalize
    _sigma_tilde_rotated[_qp] = stress_tilde_rotated;
    omega = computeOmega(delta_gamma, stress_tilde_rotated);
  }
 
  else
  {
    // Get stress_tilde
    GenericDenseVector<is_ad> stress_tilde(6);
    GenericDenseMatrix<is_ad> alpha_temp(_alpha_matrix[_qp]);
    alpha_temp.lu_solve(stress_sigma, stress_tilde);
 
    // Material property used in computeStressFinalize
    _sigma_tilde[_qp] = stress_tilde;
    omega = computeOmega(delta_gamma, stress_tilde);
  }
 
  // Hardening variable is \alpha isotropic hardening for now.
  _hardening_variable[_qp] = computeHardeningValue(delta_gamma, omega);
 
  // A small value of 1.0e-30 is added to the hardening variable to avoid numerical issues
  // related to hardening variable becoming negative early in the iteration leading to non-
  // convergence. Note that pow(x,n) requires x to be positive if n is less than 1.
 
  GenericReal<is_ad> s_y =
      _hardening_constant * pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent) +
      _yield_stress;
 
  GenericReal<is_ad> residual = 0.0;
  residual = s_y / omega - 1.0;
 
  return residual;
}

◆ computeStrainEnergyRateDensity()

template<bool is_ad>

Real HillElastoPlasticityStressUpdateTempl< is_ad >::computeStrainEnergyRateDensity	(	const GenericMaterialProperty< RankTwoTensor, is_ad > &	stress,
		const GenericMaterialProperty< RankTwoTensor, is_ad > &	strain_rate
	)

overridevirtual

Definition at line 616 of file HillElastoPlasticityStressUpdate.C.

{
  mooseError("computeStrainEnergyRateDensity not implemented for anisotropic plasticity.");
  return 0.0;
}

◆ computeStrainFinalize()

template<bool is_ad>

void HillElastoPlasticityStressUpdateTempl< is_ad >::computeStrainFinalize	(	GenericRankTwoTensor< is_ad > &	inelasticStrainIncrement,
		const GenericRankTwoTensor< is_ad > &	stress,
		const GenericDenseVector< is_ad > &	stress_dev,
		const GenericReal< is_ad > &	delta_gamma
	)

overrideprotectedvirtual

Perform any necessary steps to finalize strain increment after return mapping iterations.

Parameters

inelasticStrainIncrement	Inelastic strain increment
stress	Cauchy stresss tensor
stress_dev	Deviatoric partt of the Cauchy stresss tensor
delta_gamma	Generalized radial return's plastic multiplier

Reimplemented from AnisotropicReturnPlasticityStressUpdateBaseTempl< is_ad >.

Definition at line 522 of file HillElastoPlasticityStressUpdate.C.

{
  using std::sqrt;
 
  // e^P = delta_gamma * hill_tensor * stress
  GenericDenseVector<is_ad> inelasticStrainIncrement_vector(6);
  GenericDenseVector<is_ad> hill_stress(6);
  GenericDenseVector<is_ad> stress_vector(6);
 
  if (_local_cylindrical_csys)
  {
    GenericRankTwoTensor<is_ad> stress_rotated(stress);
    stress_rotated.rotate(_rotation_matrix_transpose[_qp]);
 
    RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress_rotated, stress_vector);
  }
  else
  {
    RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_vector);
  }
 
  _hill_tensor[_qp].vector_mult(hill_stress, stress_vector);
  hill_stress.scale(delta_gamma);
  inelasticStrainIncrement_vector = hill_stress;
 
  inelasticStrainIncrement(0, 0) = inelasticStrainIncrement_vector(0);
  inelasticStrainIncrement(1, 1) = inelasticStrainIncrement_vector(1);
  inelasticStrainIncrement(2, 2) = inelasticStrainIncrement_vector(2);
  inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
      inelasticStrainIncrement_vector(3) / 2.0;
  inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
      inelasticStrainIncrement_vector(4) / 2.0;
  inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
      inelasticStrainIncrement_vector(5) / 2.0;
 
  if (_local_cylindrical_csys)
    inelasticStrainIncrement.rotate(_rotation_matrix[_qp]);
 
  // Calculate equivalent plastic strain
  GenericDenseVector<is_ad> Mepsilon(6);
  _hill_tensor[_qp].vector_mult(Mepsilon, inelasticStrainIncrement_vector);
  GenericReal<is_ad> eq_plastic_strain_inc = Mepsilon.dot(inelasticStrainIncrement_vector);
 
  if (eq_plastic_strain_inc > 0.0)
    eq_plastic_strain_inc = sqrt(eq_plastic_strain_inc);
 
  _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp] + eq_plastic_strain_inc;
 
  AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
      inelasticStrainIncrement, stress, stress_dev, delta_gamma);
}

◆ computeStressDerivative()

template<bool is_ad>

virtual GenericReal< is_ad > AnisotropicReturnPlasticityStressUpdateBaseTempl< is_ad >::computeStressDerivative	(	const Real	,
		const Real
	)

inlineoverrideprotectedvirtualinherited

Calculate the derivative of the strain increment with respect to the updated stress.

Parameters

effective_trial_stress	Effective trial stress
scalar	Inelastic strain increment magnitude being solved for

Definition at line 47 of file AnisotropicReturnPlasticityStressUpdateBase.h.

  {
    return 0.0;
  }

◆ computeStressFinalize()

template<bool is_ad>

void HillElastoPlasticityStressUpdateTempl< is_ad >::computeStressFinalize	(	const GenericRankTwoTensor< is_ad > &	inelasticStrainIncrement,
		const GenericReal< is_ad > &	delta_gamma,
		GenericRankTwoTensor< is_ad > &	stress,
		const GenericDenseVector< is_ad > &	stress_dev,
		const GenericRankTwoTensor< is_ad > &	stress_old,
		const GenericRankFourTensor< is_ad > &	elasticity_tensor
	)

overrideprotectedvirtual

Perform any necessary steps to finalize state after return mapping iterations.

Parameters

inelasticStrainIncrement	Inelastic strain increment
delta_gamma	Generalized radial return's plastic multiplier
stress	Cauchy stresss tensor
stress_dev	Deviatoric partt of the Cauchy stresss tensor

Definition at line 580 of file HillElastoPlasticityStressUpdate.C.

{
  // Need to compute this iteration's stress tensor based on the scalar variable
  // For deviatoric
  // sigma(n+1) = {Alpha [I + delta_gamma*Delta_b]^(-1) A^-1}  sigma(trial)
 
  if (_yield_condition <= 0.0)
    return;
  GenericDenseMatrix<is_ad> inv_matrix(6, 6);
  inv_matrix.zero();
 
  for (unsigned int i = 0; i < 6; i++)
    inv_matrix(i, i) = 1 / (1 + delta_gamma * _b_eigenvalues[_qp](i));
 
  _alpha_matrix[_qp].right_multiply(inv_matrix);
 
  GenericDenseVector<is_ad> stress_output(6);
  if (_local_cylindrical_csys)
    _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde_rotated[_qp]);
  else
    _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde[_qp]);
 
  RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_output, stress_new);
 
  if (_local_cylindrical_csys)
    stress_new.rotate(_rotation_matrix[_qp]);
}

◆ computeStressInitialize()

template<bool is_ad>

void HillElastoPlasticityStressUpdateTempl< is_ad >::computeStressInitialize	(	const GenericDenseVector< is_ad > &	stress_dev,
		const GenericDenseVector< is_ad > &	stress,
		const GenericRankFourTensor< is_ad > &	elasticity_tensor
	)

overrideprotectedvirtual

Definition at line 190 of file HillElastoPlasticityStressUpdate.C.

{
  _hardening_variable[_qp] = _hardening_variable_old[_qp];
  _plasticity_strain[_qp] = _plasticity_strain_old[_qp];
  _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
 
  _elasticity_eigenvectors[_qp].resize(6, 6);
  _elasticity_eigenvalues[_qp].resize(6);
 
  _b_eigenvectors[_qp].resize(6, 6);
  _b_eigenvalues[_qp].resize(6);
 
  _alpha_matrix[_qp].resize(6, 6);
  _sigma_tilde[_qp].resize(6);
  _sigma_tilde_rotated[_qp].resize(6);
 
  // Algebra needed for the generalized return mapping of anisotropic elastoplasticity
  computeElasticityTensorEigenDecomposition();
 
  _yield_condition = 1.0; // Some positive value
  _yield_condition = -computeResidual(stress_dev, stress, 0.0);
}

◆ initQpStatefulProperties()

template<bool is_ad>

void HillElastoPlasticityStressUpdateTempl< is_ad >::initQpStatefulProperties ( )

overrideprotectedvirtual

Reimplemented from AnisotropicReturnPlasticityStressUpdateBaseTempl< is_ad >.

Definition at line 111 of file HillElastoPlasticityStressUpdate.C.

{
  RealVectorValue normal_vector(0, 0, 0);
  RealVectorValue axial_vector(0, 0, 0);
 
  if (_local_cylindrical_csys)
  {
    if (_axis == Axis::X)
    {
      normal_vector(0) = 0.0;
      normal_vector(1) = _q_point[_qp](1);
      normal_vector(2) = _q_point[_qp](2);
      axial_vector(0) = 1.0;
    }
    else if (_axis == Axis::Y)
    {
      normal_vector(0) = _q_point[_qp](0);
      normal_vector(1) = 0.0;
      normal_vector(2) = _q_point[_qp](2);
      axial_vector(1) = 1.0;
    }
    else if (_axis == Axis::Z)
    {
      normal_vector(0) = _q_point[_qp](0);
      normal_vector(1) = _q_point[_qp](1);
      normal_vector(2) = 0.0;
      axial_vector(2) = 1.0;
    }
    else
      mooseError("\nInvalid value for axis parameter provided in HillElastoPlasticityStressUpdate");
  }
 
  if (_local_cylindrical_csys)
  {
    Real nv_norm = normal_vector.norm();
    Real av_norm = axial_vector.norm();
 
    if (!(MooseUtils::absoluteFuzzyEqual(nv_norm, 0.0)))
      normal_vector /= nv_norm;
    else
    {
      mooseError("The normal vector cannot be a zero vector in "
                 "HillElastoPlasticityStressUpdate");
    }
 
    if (!(MooseUtils::absoluteFuzzyEqual(av_norm, 0.0)))
      axial_vector /= av_norm;
    else
    {
      mooseError("The axial vector cannot be a zero vector in "
                 "HillElastoPlasticityStressUpdate");
    }
 
    RankTwoScalarTools::setRotationMatrix(
        normal_vector, axial_vector, _rotation_matrix[_qp], false);
    _rotation_matrix_transpose[_qp] = _rotation_matrix[_qp].transpose();
  }
 
  AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::initQpStatefulProperties();
}

◆ paramError()

template<bool is_ad>

void Material::paramError	(	const std::string &	param,
		Args...	args
	)		const

protected

Referenced by HillElastoPlasticityStressUpdateTempl< is_ad >::HillElastoPlasticityStressUpdateTempl().

◆ propagateQpStatefulProperties()

template<bool is_ad>

void HillElastoPlasticityStressUpdateTempl< is_ad >::propagateQpStatefulProperties ( )

overrideprotectedvirtual

Reimplemented from AnisotropicReturnPlasticityStressUpdateBaseTempl< is_ad >.

Definition at line 174 of file HillElastoPlasticityStressUpdate.C.

{
  _hardening_variable[_qp] = _hardening_variable_old[_qp];
  _plasticity_strain[_qp] = _plasticity_strain_old[_qp];
 
  if (_local_cylindrical_csys)
  {
    _rotation_matrix[_qp] = _rotation_matrix_old[_qp];
    _rotation_matrix_transpose[_qp] = _rotation_matrix_transpose_old[_qp];
  }
 
  AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::propagateQpStatefulProperties();
}

◆ requiresIsotropicTensor()

template<bool is_ad>

bool HillElastoPlasticityStressUpdateTempl< is_ad >::requiresIsotropicTensor ( )

inlineoverrideprotected

Does the model require the elasticity tensor to be isotropic? No, this class does anisotropic elasto-plasticity

Definition at line 77 of file HillElastoPlasticityStressUpdate.h.

77{ return false; }

◆ validParams()

template<bool is_ad>

InputParameters HillElastoPlasticityStressUpdateTempl< is_ad >::validParams ( )

static

Definition at line 19 of file HillElastoPlasticityStressUpdate.C.

{
  InputParameters params = AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::validParams();
  params.addClassDescription(
      "This class uses the generalized radial return for anisotropic elasto-plasticity model."
      "This class can be used in conjunction with other creep and plasticity materials for "
      "more complex simulations.");
 
  params.addRequiredParam<Real>("hardening_constant",
                                "Hardening constant (H) for anisotropic plasticity");
  params.addParam<Real>(
      "hardening_exponent", 1.0, "Hardening exponent (n) for anisotropic plasticity");
  params.addRequiredParam<Real>("yield_stress",
                                "Yield stress (constant value) for anisotropic plasticity");
  params.addParam<bool>(
      "local_cylindrical_csys",
      false,
      "Compute inelastic strain increment in local cylindrical coordinate system");
 
  params.addParam<MooseEnum>("axis",
                             MooseEnum("x y z", "y"),
                             "The axis of cylindrical component about which to rotate when "
                             "computing inelastic strain increment in local coordinate system");
 
  return params;
}

Member Data Documentation

◆ _alpha_matrix

template<bool is_ad>

GenericMaterialProperty<DenseMatrix<Real>, is_ad>& HillElastoPlasticityStressUpdateTempl< is_ad >::_alpha_matrix

protected

Definition at line 149 of file HillElastoPlasticityStressUpdate.h.

◆ _anisotropic_elastic_tensor

template<bool is_ad>

GenericDenseMatrix<is_ad> HillElastoPlasticityStressUpdateTempl< is_ad >::_anisotropic_elastic_tensor

protected

Definition at line 130 of file HillElastoPlasticityStressUpdate.h.

◆ _axis

template<bool is_ad>

enum HillElastoPlasticityStressUpdateTempl::Axis HillElastoPlasticityStressUpdateTempl< is_ad >::_axis

protected

◆ _b_eigenvalues

template<bool is_ad>

GenericMaterialProperty<DenseVector<Real>, is_ad>& HillElastoPlasticityStressUpdateTempl< is_ad >::_b_eigenvalues

protected

Definition at line 147 of file HillElastoPlasticityStressUpdate.h.

◆ _b_eigenvectors

template<bool is_ad>

GenericMaterialProperty<DenseMatrix<Real>, is_ad>& HillElastoPlasticityStressUpdateTempl< is_ad >::_b_eigenvectors

protected

Definition at line 146 of file HillElastoPlasticityStressUpdate.h.

◆ _current_elem

template<bool is_ad>

const Elem *const & Material::_current_elem

protected

◆ _eigenvalues_hill

template<bool is_ad>

GenericDenseVector<is_ad> HillElastoPlasticityStressUpdateTempl< is_ad >::_eigenvalues_hill

protected

Definition at line 128 of file HillElastoPlasticityStressUpdate.h.

◆ _eigenvectors_hill

template<bool is_ad>

GenericDenseMatrix<is_ad> HillElastoPlasticityStressUpdateTempl< is_ad >::_eigenvectors_hill

protected

Definition at line 129 of file HillElastoPlasticityStressUpdate.h.

◆ _elasticity_eigenvalues

template<bool is_ad>

GenericMaterialProperty<DenseVector<Real>, is_ad>& HillElastoPlasticityStressUpdateTempl< is_ad >::_elasticity_eigenvalues

protected

Definition at line 144 of file HillElastoPlasticityStressUpdate.h.

◆ _elasticity_eigenvectors

template<bool is_ad>

GenericMaterialProperty<DenseMatrix<Real>, is_ad>& HillElastoPlasticityStressUpdateTempl< is_ad >::_elasticity_eigenvectors

protected

Definition at line 143 of file HillElastoPlasticityStressUpdate.h.

◆ _elasticity_tensor

template<bool is_ad>

const GenericMaterialProperty<RankFourTensor, is_ad>& HillElastoPlasticityStressUpdateTempl< is_ad >::_elasticity_tensor

protected

Anisotropic elasticity tensor material property.

Definition at line 135 of file HillElastoPlasticityStressUpdate.h.

◆ _elasticity_tensor_name

template<bool is_ad>

const std::string HillElastoPlasticityStressUpdateTempl< is_ad >::_elasticity_tensor_name

protected

Name of the elasticity tensor material property.

Definition at line 133 of file HillElastoPlasticityStressUpdate.h.

◆ _hardening_constant

template<bool is_ad>

const Real HillElastoPlasticityStressUpdateTempl< is_ad >::_hardening_constant

protected

Definition at line 137 of file HillElastoPlasticityStressUpdate.h.

◆ _hardening_derivative

template<bool is_ad>

GenericReal<is_ad> HillElastoPlasticityStressUpdateTempl< is_ad >::_hardening_derivative

protected

Definition at line 153 of file HillElastoPlasticityStressUpdate.h.

◆ _hardening_exponent

template<bool is_ad>

const Real HillElastoPlasticityStressUpdateTempl< is_ad >::_hardening_exponent

protected

Definition at line 138 of file HillElastoPlasticityStressUpdate.h.

◆ _hardening_variable

template<bool is_ad>

GenericMaterialProperty<Real, is_ad>& HillElastoPlasticityStressUpdateTempl< is_ad >::_hardening_variable

protected

Definition at line 140 of file HillElastoPlasticityStressUpdate.h.

◆ _hardening_variable_old

template<bool is_ad>

const MaterialProperty<Real>& HillElastoPlasticityStressUpdateTempl< is_ad >::_hardening_variable_old

protected

Definition at line 141 of file HillElastoPlasticityStressUpdate.h.

◆ _hill_tensor

template<bool is_ad>

const MaterialProperty<DenseMatrix<Real> >& HillElastoPlasticityStressUpdateTempl< is_ad >::_hill_tensor

protected

Hill tensor, when global axes do not (somehow) align with those of the material Example: Large rotation due to rigid body and/or large deformation kinematics.

Definition at line 159 of file HillElastoPlasticityStressUpdate.h.