SolidMechanicsPlasticJ2.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* 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 "SolidMechanicsPlasticJ2.h"
#include "RankFourTensor.h"

registerMooseObject("SolidMechanicsApp", SolidMechanicsPlasticJ2);
registerMooseObjectRenamed("SolidMechanicsApp",
                           TensorMechanicsPlasticJ2,
                           "01/01/2025 00:00",
                           SolidMechanicsPlasticJ2);

InputParameters
SolidMechanicsPlasticJ2::validParams()
{
  InputParameters params = SolidMechanicsPlasticModel::validParams();
  params.addRequiredParam<UserObjectName>(
      "yield_strength",
      "A SolidMechanicsHardening UserObject that defines hardening of the yield strength");
  params.addRangeCheckedParam<unsigned>(
      "max_iterations", 10, "max_iterations>0", "Maximum iterations for custom J2 return map");
  params.addParam<bool>("use_custom_returnMap",
                        true,
                        "Whether to use the custom returnMap "
                        "algorithm.  Set to true if you are using "
                        "isotropic elasticity.");
  params.addParam<bool>("use_custom_cto",
                        true,
                        "Whether to use the custom consistent tangent "
                        "operator computations.  Set to true if you are "
                        "using isotropic elasticity.");
  params.addClassDescription("J2 plasticity, associative, with hardening");

  return params;
}

SolidMechanicsPlasticJ2::SolidMechanicsPlasticJ2(const InputParameters & parameters)
  : SolidMechanicsPlasticModel(parameters),
    _strength(getUserObject<SolidMechanicsHardeningModel>("yield_strength")),
    _max_iters(getParam<unsigned>("max_iterations")),
    _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
    _use_custom_cto(getParam<bool>("use_custom_cto"))
{
}

Real
SolidMechanicsPlasticJ2::yieldFunction(const RankTwoTensor & stress, Real intnl) const
{
  return std::sqrt(3.0 * stress.secondInvariant()) - yieldStrength(intnl);
}

RankTwoTensor
SolidMechanicsPlasticJ2::dyieldFunction_dstress(const RankTwoTensor & stress, Real /*intnl*/) const
{
  Real sII = stress.secondInvariant();
  if (sII == 0.0)
    return RankTwoTensor();
  else
    return 0.5 * std::sqrt(3.0 / sII) * stress.dsecondInvariant();
}

Real
SolidMechanicsPlasticJ2::dyieldFunction_dintnl(const RankTwoTensor & /*stress*/, Real intnl) const
{
  return -dyieldStrength(intnl);
}

RankTwoTensor
SolidMechanicsPlasticJ2::flowPotential(const RankTwoTensor & stress, Real intnl) const
{
  return dyieldFunction_dstress(stress, intnl);
}

RankFourTensor
SolidMechanicsPlasticJ2::dflowPotential_dstress(const RankTwoTensor & stress, Real /*intnl*/) const
{
  Real sII = stress.secondInvariant();
  if (sII == 0)
    return RankFourTensor();

  RankFourTensor dfp = 0.5 * std::sqrt(3.0 / sII) * stress.d2secondInvariant();
  Real pre = -0.25 * std::sqrt(3.0) * std::pow(sII, -1.5);
  RankTwoTensor dII = stress.dsecondInvariant();
  for (unsigned i = 0; i < 3; ++i)
    for (unsigned j = 0; j < 3; ++j)
      for (unsigned k = 0; k < 3; ++k)
        for (unsigned l = 0; l < 3; ++l)
          dfp(i, j, k, l) += pre * dII(i, j) * dII(k, l);
  return dfp;
}

RankTwoTensor
SolidMechanicsPlasticJ2::dflowPotential_dintnl(const RankTwoTensor & /*stress*/,
                                               Real /*intnl*/) const
{
  return RankTwoTensor();
}

Real
SolidMechanicsPlasticJ2::yieldStrength(Real intnl) const
{
  return _strength.value(intnl);
}

Real
SolidMechanicsPlasticJ2::dyieldStrength(Real intnl) const
{
  return _strength.derivative(intnl);
}

std::string
SolidMechanicsPlasticJ2::modelName() const
{
  return "J2";
}

bool
SolidMechanicsPlasticJ2::returnMap(const RankTwoTensor & trial_stress,
                                   Real intnl_old,
                                   const RankFourTensor & E_ijkl,
                                   Real ep_plastic_tolerance,
                                   RankTwoTensor & returned_stress,
                                   Real & returned_intnl,
                                   std::vector<Real> & dpm,
                                   RankTwoTensor & delta_dp,
                                   std::vector<Real> & yf,
                                   bool & trial_stress_inadmissible) const
{
  if (!(_use_custom_returnMap))
    return SolidMechanicsPlasticModel::returnMap(trial_stress,
                                                 intnl_old,
                                                 E_ijkl,
                                                 ep_plastic_tolerance,
                                                 returned_stress,
                                                 returned_intnl,
                                                 dpm,
                                                 delta_dp,
                                                 yf,
                                                 trial_stress_inadmissible);

  yf.resize(1);

  Real yf_orig = yieldFunction(trial_stress, intnl_old);

  yf[0] = yf_orig;

  if (yf_orig < _f_tol)
  {
    // the trial_stress is admissible
    trial_stress_inadmissible = false;
    return true;
  }

  trial_stress_inadmissible = true;
  Real mu = E_ijkl(0, 1, 0, 1);

  // Perform a Newton-Raphson to find dpm when
  // residual = 3*mu*dpm - trial_equivalent_stress + yieldStrength(intnl_old + dpm) = 0
  Real trial_equivalent_stress = yf_orig + yieldStrength(intnl_old);
  Real residual;
  Real jac;
  dpm[0] = 0;
  unsigned int iter = 0;
  do
  {
    residual = 3.0 * mu * dpm[0] - trial_equivalent_stress + yieldStrength(intnl_old + dpm[0]);
    jac = 3.0 * mu + dyieldStrength(intnl_old + dpm[0]);
    dpm[0] += -residual / jac;
    if (iter > _max_iters) // not converging
      return false;
    iter++;
  } while (residual * residual > _f_tol * _f_tol);

  // set the returned values
  yf[0] = 0;
  returned_intnl = intnl_old + dpm[0];
  RankTwoTensor nn = 1.5 * trial_stress.deviatoric() /
                     trial_equivalent_stress; // = dyieldFunction_dstress(trial_stress, intnl_old) =
                                              // the normal to the yield surface, at the trial
                                              // stress
  returned_stress = 2.0 / 3.0 * nn * yieldStrength(returned_intnl);
  returned_stress.addIa(1.0 / 3.0 * trial_stress.trace());
  delta_dp = nn * dpm[0];

  return true;
}

RankFourTensor
SolidMechanicsPlasticJ2::consistentTangentOperator(const RankTwoTensor & trial_stress,
                                                   Real intnl_old,
                                                   const RankTwoTensor & stress,
                                                   Real intnl,
                                                   const RankFourTensor & E_ijkl,
                                                   const std::vector<Real> & cumulative_pm) const
{
  if (!_use_custom_cto)
    return SolidMechanicsPlasticModel::consistentTangentOperator(
        trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);

  Real mu = E_ijkl(0, 1, 0, 1);

  Real h = 3 * mu + dyieldStrength(intnl);
  RankTwoTensor sij = stress.deviatoric();
  Real sII = stress.secondInvariant();
  Real equivalent_stress = std::sqrt(3.0 * sII);
  Real zeta = cumulative_pm[0] / (1.0 + 3.0 * mu * cumulative_pm[0] / equivalent_stress);

  return E_ijkl - 3.0 * mu * mu / sII / h * sij.outerProduct(sij) -
         4.0 * mu * mu * zeta * dflowPotential_dstress(stress, intnl);
}

bool
SolidMechanicsPlasticJ2::useCustomReturnMap() const
{
  return _use_custom_returnMap;
}

bool
SolidMechanicsPlasticJ2::useCustomCTO() const
{
  return _use_custom_cto;
}