slipsystem.C file

[Emmanuel-324]

I am trying to locate the sourcecode that contains the code for different types of crystal structures used in the crystal plasticity model.
Below is a snippet of the relevant code:

[trial_xtalpl_phase0]
    type = CrystalPlasticityKalidindiUpdate
    number_slip_systems = 12
    slip_sys_file_name = input_slip_sys.txt
    crystal_lattice_type = BCC
    resistance_tol = 0.01
    r = 1.4             
    h = 6000            
    t_sat = 598.5        
    gss_a = 1.5         
    ao = 0.001           
    xm = 0.017             
    gss_initial = 400
    base_name = phase0
  []

Below is also the CrystalPlasticityKalidindiUpdate.C

//* 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 "CrystalPlasticityKalidindiUpdate.h"
#include "libmesh/int_range.h"

registerMooseObject("SolidMechanicsApp", CrystalPlasticityKalidindiUpdate);

InputParameters
CrystalPlasticityKalidindiUpdate::validParams()
{
  InputParameters params = CrystalPlasticityStressUpdateBase::validParams();
  params.addClassDescription("Kalidindi version of homogeneous crystal plasticity.");
  params.addParam<Real>("r", 1.0, "Latent hardening coefficient");
  params.addParam<Real>("h", 541.5, "hardening constants");
  params.addParam<Real>("t_sat", 109.8, "saturated slip system strength");
  params.addParam<Real>("gss_a", 2.5, "coefficient for hardening");
  params.addParam<Real>("ao", 0.001, "slip rate coefficient");
  params.addParam<Real>("xm", 0.1, "exponent for slip rate");
  params.addParam<Real>("gss_initial", 60.8, "initial lattice friction strength of the material");

  params.addParam<MaterialPropertyName>(
      "total_twin_volume_fraction",
      "Total twin volume fraction, if twinning is considered in the simulation");

  return params;
}

CrystalPlasticityKalidindiUpdate::CrystalPlasticityKalidindiUpdate(
    const InputParameters & parameters)
  : CrystalPlasticityStressUpdateBase(parameters),
    // Constitutive values
    _r(getParam<Real>("r")),
    _h(getParam<Real>("h")),
    _tau_sat(getParam<Real>("t_sat")),
    _gss_a(getParam<Real>("gss_a")),
    _ao(getParam<Real>("ao")),
    _xm(getParam<Real>("xm")),
    _gss_initial(getParam<Real>("gss_initial")),

    // resize vectors used in the consititutive slip hardening
    _hb(_number_slip_systems, 0.0),
    _slip_resistance_increment(_number_slip_systems, 0.0),

    // resize local caching vectors used for substepping
    _previous_substep_slip_resistance(_number_slip_systems, 0.0),
    _slip_resistance_before_update(_number_slip_systems, 0.0),

    // Twinning contributions, if used
    _include_twinning_in_Lp(parameters.isParamValid("total_twin_volume_fraction")),
    _twin_volume_fraction_total(_include_twinning_in_Lp
                                    ? &getMaterialPropertyOld<Real>("total_twin_volume_fraction")
                                    : nullptr)
{
}

void
CrystalPlasticityKalidindiUpdate::initQpStatefulProperties()
{
  CrystalPlasticityStressUpdateBase::initQpStatefulProperties();
  for (const auto i : make_range(_number_slip_systems))
  {
    _slip_resistance[_qp][i] = _gss_initial;
    _slip_increment[_qp][i] = 0.0;
  }
}

void
CrystalPlasticityKalidindiUpdate::setInitialConstitutiveVariableValues()
{
  // Would also set old dislocation densities here if included in this model
  _slip_resistance[_qp] = _slip_resistance_old[_qp];
  _previous_substep_slip_resistance = _slip_resistance_old[_qp];
}

void
CrystalPlasticityKalidindiUpdate::setSubstepConstitutiveVariableValues()
{
  // Would also set substepped dislocation densities here if included in this model
  _slip_resistance[_qp] = _previous_substep_slip_resistance;
}

bool
CrystalPlasticityKalidindiUpdate::calculateSlipRate()
{
  for (const auto i : make_range(_number_slip_systems))
  {
    _slip_increment[_qp][i] =
        _ao * std::pow(std::abs(_tau[_qp][i] / _slip_resistance[_qp][i]), 1.0 / _xm);
    if (_tau[_qp][i] < 0.0)
      _slip_increment[_qp][i] *= -1.0;

    if (std::abs(_slip_increment[_qp][i]) * _substep_dt > _slip_incr_tol)
    {
      if (_print_convergence_message)
        mooseWarning("Maximum allowable slip increment exceeded ",
                     std::abs(_slip_increment[_qp][i]) * _substep_dt);

      return false;
    }
  }
  return true;
}

void
CrystalPlasticityKalidindiUpdate::calculateEquivalentSlipIncrement(
    RankTwoTensor & equivalent_slip_increment)
{
  if (_include_twinning_in_Lp)
  {
    for (const auto i : make_range(_number_slip_systems))
      equivalent_slip_increment += (1.0 - (*_twin_volume_fraction_total)[_qp]) *
                                   _flow_direction[_qp][i] * _slip_increment[_qp][i] * _substep_dt;
  }
  else // if no twinning volume fraction material property supplied, use base class
    CrystalPlasticityStressUpdateBase::calculateEquivalentSlipIncrement(equivalent_slip_increment);
}

void
CrystalPlasticityKalidindiUpdate::calculateConstitutiveSlipDerivative(
    std::vector<Real> & dslip_dtau)
{
  for (const auto i : make_range(_number_slip_systems))
  {
    if (MooseUtils::absoluteFuzzyEqual(_tau[_qp][i], 0.0))
      dslip_dtau[i] = 0.0;
    else
      dslip_dtau[i] = _ao / _xm *
                      std::pow(std::abs(_tau[_qp][i] / _slip_resistance[_qp][i]), 1.0 / _xm - 1.0) /
                      _slip_resistance[_qp][i];
  }
}

bool
CrystalPlasticityKalidindiUpdate::areConstitutiveStateVariablesConverged()
{
  return isConstitutiveStateVariableConverged(_slip_resistance[_qp],
                                              _slip_resistance_before_update,
                                              _previous_substep_slip_resistance,
                                              _resistance_tol);
}

void
CrystalPlasticityKalidindiUpdate::updateSubstepConstitutiveVariableValues()
{
  // Would also set substepped dislocation densities here if included in this model
  _previous_substep_slip_resistance = _slip_resistance[_qp];
}

void
CrystalPlasticityKalidindiUpdate::cacheStateVariablesBeforeUpdate()
{
  _slip_resistance_before_update = _slip_resistance[_qp];
}

void
CrystalPlasticityKalidindiUpdate::calculateStateVariableEvolutionRateComponent()
{
  for (const auto i : make_range(_number_slip_systems))
  {
    // Clear out increment from the previous iteration
    _slip_resistance_increment[i] = 0.0;

    _hb[i] = _h * std::pow(std::abs(1.0 - _slip_resistance[_qp][i] / _tau_sat), _gss_a);
    const Real hsign = 1.0 - _slip_resistance[_qp][i] / _tau_sat;
    if (hsign < 0.0)
      _hb[i] *= -1.0;
  }

  for (const auto i : make_range(_number_slip_systems))
  {
    for (const auto j : make_range(_number_slip_systems))
    {
      unsigned int iplane, jplane;
      iplane = i / 3;
      jplane = j / 3;

      if (iplane == jplane) // self vs. latent hardening
        _slip_resistance_increment[i] +=
            std::abs(_slip_increment[_qp][j]) * _hb[j]; // q_{ab} = 1.0 for self hardening
      else
        _slip_resistance_increment[i] +=
            std::abs(_slip_increment[_qp][j]) * _r * _hb[j]; // latent hardenign
    }
  }
}

bool
CrystalPlasticityKalidindiUpdate::updateStateVariables()
{
  // Now perform the check to see if the slip system should be updated
  for (const auto i : make_range(_number_slip_systems))
  {
    _slip_resistance_increment[i] *= _substep_dt;
    if (_previous_substep_slip_resistance[i] < _zero_tol && _slip_resistance_increment[i] < 0.0)
      _slip_resistance[_qp][i] = _previous_substep_slip_resistance[i];
    else
      _slip_resistance[_qp][i] =
          _previous_substep_slip_resistance[i] + _slip_resistance_increment[i];

    if (_slip_resistance[_qp][i] < 0.0)
      return false;
  }
  return true;
}

@dschwen @sapitts Can you help me locate the code/ source for the different crystal structures(FCC, BCC, HCP) used in the model ?above.

[sapitts]

yes, of course. That's set by a parameter in the base class, see CrystalPlasticityStressUpdateBase.C#L34 Search through the code in that class to see how and where the lattice type is used.

Please note too that the Kalidindi model is hard coded for an FCC structure, as described in the documentation, Eq 3: https://mooseframework.inl.gov/source/materials/crystal_plasticity/CrystalPlasticityKalidindiUpdate.html

[Emmanuel-324]

@sapitts I tried to add a new crystal structure(specifically the DO22 structure, which is a tetragonal structure). It compiles successfully, but I just wanted to share with you if there is room for improvement or if there is something that I am missing.

Below are relevant snippets of the code:
In the CrystalPlasticityStressUpdateBase.C

void
CrystalPlasticityStressUpdateBase_abs::getSlipSystems()
{
  bool orthonormal_error = false;

  // read in the slip system data from auxiliary text file
  MooseUtils::DelimitedFileReader _reader(_slip_sys_file_name);
  _reader.setFormatFlag(MooseUtils::DelimitedFileReader::FormatFlag::ROWS);
  _reader.read();

  // check the size of the input
  if (_reader.getData().size() != _number_slip_systems)
    paramError(
        "number_slip_systems",
        "The number of rows in the slip system file should match the number of slip system.");
  for (const auto i : make_range(_number_slip_systems))
  {
    // initialize to zero
    _slip_direction[i].zero();
    _slip_plane_normal[i].zero();
  }
  if (_crystal_lattice_type == CrystalLatticeType::DO22)
    transformDO22SlipSystems(_reader);
  else if (_crystal_lattice_type == CrystalLatticeType::HCP)
    transformHexagonalMillerBravaisSlipSystems(_reader);
  else if (_crystal_lattice_type == CrystalLatticeType::BCC ||
           _crystal_lattice_type == CrystalLatticeType::FCC)
  {
    for (const auto i : make_range(_number_slip_systems))
    {
      // directly grab the raw data and scale it by the unit cell dimension
      for (const auto j : index_range(_reader.getData(i)))
      {
        if (j < LIBMESH_DIM)
          _slip_plane_normal[i](j) = _reader.getData(i)[j] / _unit_cell_dimension[j];
        else
          _slip_direction[i](j - LIBMESH_DIM) =
              _reader.getData(i)[j] * _unit_cell_dimension[j - LIBMESH_DIM];
      }
    }
  }
  for (const auto i : make_range(_number_slip_systems))
  {
    // normalize
    _slip_plane_normal[i] /= _slip_plane_normal[i].norm();
    _slip_direction[i] /= _slip_direction[i].norm();

    if (_crystal_lattice_type != CrystalLatticeType::HCP)
    {
      const auto magnitude = _slip_plane_normal[i] * _slip_direction[i];
      if (std::abs(magnitude) > libMesh::TOLERANCE)
      {
        orthonormal_error = true;
        break;
      }
    }
  }
  if (orthonormal_error)
    mooseError("CrystalPlasticityStressUpdateBase_abs Error: The slip system file contains a slip "
               "direction and plane normal pair that are not orthonormal in the Cartesian "
               "coordinate system.");
}

void CrystalPlasticityStressUpdateBase_abs::transformDO22SlipSystems(
    const MooseUtils::DelimitedFileReader & reader)
{
    // Example logic for transforming DO22 slip systems 
    const unsigned int expected_columns = 6; 
    RealVectorValue temporary_slip_direction, temporary_slip_plane;

    for (const auto i : make_range(_number_slip_systems))
    {
        // Read the slip plane and direction data from the file
        for (const auto j : index_range(reader.getData(i)))
        {
            if (j < 3)
                temporary_slip_plane(j) = reader.getData(i)[j];
            else
                temporary_slip_direction(j - 3) = reader.getData(i)[j];
        }
        // Normalize the slip direction and slip plane normal

        _slip_plane_normal[i] = temporary_slip_plane / temporary_slip_plane.norm();
        _slip_direction[i] = temporary_slip_direction / temporary_slip_direction.norm();
    }

}

[GiudGiud]

how often are getSlipSystems() and transformDO22SlipSystems() called? once at initialization right?

[Emmanuel-324]

@GiudGiud Yes.

Also, @sapitts @dschwen do you have any suggestions about how to get different slip system strengths for certain slip systems? For example, out of the 12 slip systems for FCC, I specify every 4 with different slip system strengths in the input script.

So basically each of the 4 groupings has a different slip system in the same input script.

[Emmanuel-324]

@sapitts @dschwen A quick update on the different slip modes I tried to implement in a similar way as the HCP. Below are snippets of it.

InputParameters
CrystalPlasticityKalidindiUpdate_slip::validParams()
{
  InputParameters params = CrystalPlasticityStressUpdateBase_abs::validParams();
  params.addParam<unsigned int>(
    "slip_system_modes",
    1,
    "Number of different types of slip systems, each with its own friction strength.");
  params.addParam<std::vector<unsigned int>>(
    "number_slip_systems_per_mode",
    std::vector<unsigned int>(),
    "The number of slip systems per slip mode.");
  params.addParam<std::vector<Real>>(
    "lattice_friction_per_mode",
    std::vector<Real>(),
    "Lattice friction strength for each slip mode in MPa.");

  params.addClassDescription("Kalidindi version of homogeneous crystal plasticity.");
  params.addParam<Real>("r", 1.0, "Latent hardening coefficient");
  params.addParam<Real>("h", 541.5, "hardening constants");
  params.addParam<Real>("t_sat", 109.8, "saturated slip system strength");
  params.addParam<Real>("gss_a", 2.5, "coefficient for hardening");
  params.addParam<Real>("ao", 0.001, "slip rate coefficient");
  params.addParam<Real>("xm", 0.1, "exponent for slip rate");
  params.addParam<Real>("gss_initial", 60.8, "initial lattice friction strength of the material");

  params.addParam<MaterialPropertyName>(
      "total_twin_volume_fraction",
      "Total twin volume fraction, if twinning is considered in the simulation");

  return params;
}

CrystalPlasticityKalidindiUpdate_slip::CrystalPlasticityKalidindiUpdate_slip(
    const InputParameters & parameters)
  : CrystalPlasticityStressUpdateBase_abs(parameters),
    // Constitutive values
    _slip_system_modes(getParam<unsigned int>("slip_system_modes")),
    _number_slip_systems_per_mode(getParam<std::vector<unsigned int>>("number_slip_systems_per_mode")),
    _lattice_friction(getParam<std::vector<Real>>("lattice_friction_per_mode")),
    _r(getParam<Real>("r")),
    _h(getParam<Real>("h")),
    _tau_sat(getParam<Real>("t_sat")),
    _gss_a(getParam<Real>("gss_a")),
    _ao(getParam<Real>("ao")),
    _xm(getParam<Real>("xm")),
    _gss_initial(getParam<Real>("gss_initial")),

    // resize vectors used in the consititutive slip hardening
    _hb(_number_slip_systems, 0.0),
    _slip_resistance_increment(_number_slip_systems, 0.0),

    // resize local caching vectors used for substepping
    _previous_substep_slip_resistance(_number_slip_systems, 0.0),
    _slip_resistance_before_update(_number_slip_systems, 0.0),

    // Twinning contributions, if used
    _include_twinning_in_Lp(parameters.isParamValid("total_twin_volume_fraction")),
    _twin_volume_fraction_total(_include_twinning_in_Lp
                                    ? &getMaterialPropertyOld<Real>("total_twin_volume_fraction")
                                    : nullptr)
{
    // check that the number of slip systems is equal to the sum of the types of slip system
    if (_number_slip_systems_per_mode.size() != _slip_system_modes)
    paramError("number_slip_systems_per_mode",
               "The size the number of slip systems per mode is not equal to the number of slip "
               "system types.");                                
                                    
    if (_lattice_friction.size() != _slip_system_modes)
    paramError("lattice_friction_per_mode",
               "Please ensure that the size of lattice_friction_per_mode equals the value supplied "
               "for slip_system_modes");

  unsigned int sum = 0;
  for (const auto i : make_range(_slip_system_modes))
    sum += _number_slip_systems_per_mode[i];
  if (sum != _number_slip_systems)
    paramError("slip_system_modes",
               "The number of slip systems and the sum of the slip systems in each of the slip "
               "system modes are not equal");

}

void CrystalPlasticityKalidindiUpdate_slip::initQpStatefulProperties()
{
  CrystalPlasticityStressUpdateBase_abs::initQpStatefulProperties();
  
  // Set initial resistance from lattice friction, which is type dependent
  unsigned int slip_mode = 0;
  unsigned int counter_adjustment = 0;

  for (const auto i : make_range(_number_slip_systems))
  {
    if ((i - counter_adjustment) < _number_slip_systems_per_mode[slip_mode])
      _slip_resistance[_qp][i] = _lattice_friction[slip_mode];
    else
    {
      counter_adjustment += _number_slip_systems_per_mode[slip_mode];
      ++slip_mode;
      _slip_resistance[_qp][i] = _lattice_friction[slip_mode];
    }

    _slip_increment[_qp][i] = 0.0;
  }

}

void CrystalPlasticityKalidindiUpdate_slip::calculateStateVariableEvolutionRateComponent()
{
  for (const auto i : make_range(_number_slip_systems))
  {
    _slip_resistance_increment[i] = 0.0;

    unsigned int slip_mode = 0;
    unsigned int counter_adjustment = 0;
    for (const auto j : make_range(_number_slip_systems))
    {
      if ((j - counter_adjustment) < _number_slip_systems_per_mode[slip_mode])
      {
        _hb[j] = _h * std::pow(std::abs(1.0 - _slip_resistance[_qp][j] / _tau_sat), _gss_a);
      }
      else
      {
        counter_adjustment += _number_slip_systems_per_mode[slip_mode];
        ++slip_mode;
        _hb[j] = _h * std::pow(std::abs(1.0 - _slip_resistance[_qp][j] / _tau_sat), _gss_a);
      }

      Real hsign = 1.0 - _slip_resistance[_qp][j] / _tau_sat;
      if (hsign < 0.0)
        _hb[j] *= -1.0;

      unsigned int iplane = i / 3;
      unsigned int jplane = j / 3;

      if (iplane == jplane) 
        _slip_resistance_increment[i] += std::abs(_slip_increment[_qp][j]) * _hb[j];
      else 
        _slip_resistance_increment[i] += std::abs(_slip_increment[_qp][j]) * _r * _hb[j];
    }
  }
}