CP with threshold

[amassaf]

Hello,

I want to implement a slip rate for the crystal plasticity with a backstress and a threshold of this form:

I have created a new class based on the existing CrystalPlasticityKalidindiUpdate class that changes the calculateSlipRate and state variable evolution accordingly, the application compiles all is well when I try to run it on one single crystal, but when I try to run it with a finer mesh or with a polycrystal I keep getting convergence error with this printed message :

Time Step 1, time = 10, dt = 10

Performing automatic scaling calculation


ComputeMultipleCrystalPlasticityStress: Constitutive failure
To recover, the solution will fail and then be re-attempted with a reduced time step.

 Solve Did NOT Converge!
Aborting as solve did not converge

I do not know what exactly leads to this constitutive failure as there is no other error message even with print_state_variable_convergence_error_messages = true. If anyone has any idea about where the convergence issue might come from, I would appreciate the help, thank you :)

[GiudGiud]

Hello

The inner solves of the computeMultipleCrystalPlasticity Stress are failing
For example, it could fail if this new calculateSlipRate in your model is not returning true.

    _convergence_failed = !_models[i]->calculateSlipRate();

this is an example implementation

bool
CrystalPlasticityHCPDislocationSlipBeyerleinUpdate::calculateSlipRate()
{
  for (const auto i : make_range(_number_slip_systems))
  {
    Real driving_force = std::abs(_tau[_qp][i] / _slip_resistance[_qp][i]);
    if (driving_force < _zero_tol)
      _slip_increment[_qp][i] = 0.0;
    else
    {
      _slip_increment[_qp][i] =
          _reference_strain_rate * std::pow(driving_force, (1.0 / _rate_sensitivity_exponent));
      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;
}

does your return true? Or is your implementation detecting an error and returning false?

[GiudGiud]

@dewenyushu @sapitts what are your thoughts on this?
Any technique to get this model to avid the constitutive failure?

[sapitts]

@amassaf would you post your entire input file for the run that uses the CrystalPlasticityKalidindiUpdate model here please? Or share a github branch? I'll have an easier time giving relevant suggestions if I have the context from the input file

[amassaf]

Hi, here is the input file :

[Mesh]
  [msh]
    type = GeneratedMeshGenerator
    dim = 3
    xmax = 1
    ymax = 1
    nx = 1
    ny = 1
    nz = 1
  []
  [top_node]
    type = ExtraNodesetGenerator
    coord = '0 1 0'
    input = msh
    new_boundary = top_node
  []
  [bottom_node]
    type = ExtraNodesetGenerator
    coord = '0 0 0'
    input = top_node
    new_boundary = bottom_node
  []
  [bottom_node2]
    type = ExtraNodesetGenerator
    coord = '1 0 1'
    input = bottom_node
    new_boundary = bottom_node2
  []

[]  
  
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]  

[AuxVariables]
  [./resid_x]
  [../]
  [./resid_y]
  [../] 
  [./resid_z]
  [../]

  [./pk2]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./fp_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./rotout]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./e_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./gss]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./slip_increment]
   order = CONSTANT
   family = MONOMIAL
  [../]
[]

[Modules/TensorMechanics/Master/all]
  strain = FINITE
  add_variables = true
  generate_output = stress_yy
  save_in = 'resid_x resid_y resid_z'
[]

[AuxKernels]
[./fp_yy]
    type = RankTwoAux
    variable = fp_yy
    rank_two_tensor = plastic_deformation_gradient
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  [../]
  [./pk2]
   type = RankTwoAux
   variable = pk2
   rank_two_tensor = second_piola_kirchhoff_stress
   index_j = 1
   index_i = 1
   execute_on = timestep_end
  [../]
  [./e_yy]
    type = RankTwoAux
    variable = e_yy
    rank_two_tensor = total_lagrangian_strain
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  [../]
  [./gss]
    type = MaterialStdVectorAux
    variable = gss
    property = slip_resistance
    index = 0
    execute_on = timestep_end
  [../]
  [./slip_inc]
   type = MaterialStdVectorAux
   variable = slip_increment
   property = slip_increment
   index = 0
   execute_on = timestep_end
  [../]
[]

[BCs]
  [./symmy]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0
  [../]
  [./symmx]
    type = DirichletBC
    variable = disp_x
    boundary = bottom_node
    value = 0
  [../]
  [./symmz]
    type = DirichletBC
    variable = disp_z
    boundary = bottom_node
    value = 0
  [../]
[./symmz2]
    type = DirichletBC
    variable = disp_z
    boundary = bottom_node2
    value = 0
  [../]
  [./tdisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = '-1e-5*t'
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeElasticityTensorCP
    C_ijkl = '48.7e3  11.9e3  11.9e3  48.7e3 11.9e3  48.7e3   13.1e3 13.1e3 13.1e3'
    fill_method = symmetric9
  [../]

  [./stress]
    type = ComputeMultipleCrystalPlasticityStress
    crystal_plasticity_models = 'trial_xtalpl'
    tan_mod_type = exact
    maximum_substep_iteration = 5
    rtol = 1e-7 
    maxiter_state_variable = 10 
    #use_line_search = true
    #print_state_variable_convergence_error_messages = true
  [../]
  [./trial_xtalpl]
    type = CrystalPlasticityKalidindiUpdate
    number_slip_systems = 12 
    slip_sys_file_name = 2input.txt
    ao = 2e-3
    gss_a = 2
    gss_initial = 9
    h = 80
    r = 1.4
    xm = 2
    t_sat = 80.
    slip_increment_tolerance = 0.1 
    stol = 0.1 
    resistance_tol = 0.1
    #print_state_variable_convergence_error_messages = true

  [../]
[]

[Postprocessors]
    [./top_x]
    type = NodalExtremeValue
    variable = disp_x
    boundary = top
    value_type = max
  [../]
[./top_y]
    type = NodalExtremeValue
    variable = disp_y
    boundary = top
    value_type = max
  [../]
  [./resid_x]
    type = NodalSum
    variable = resid_x
    boundary = bottom
  [../]
  [./resid_y]
    type = NodalSum
    variable = resid_y
    boundary = bottom
  [../]

  [./stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  [../]
  [./pk2]
   type = ElementAverageValue
   variable = pk2
  [../]
  [./fp_yy]
    type = ElementAverageValue
    variable = fp_yy
  [../]
  [./e_yy]
    type = ElementAverageValue
    variable = e_yy
  [../]
  [./gss]
    type = ElementAverageValue
    variable = gss
  [../]
  [./slip_increment]
   type = ElementAverageValue
   variable = slip_increment
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
   type = Transient
  solve_type = 'newton'
 petsc_options_iname = ' -pc_hypre_type -pc_factor_mat_solver_package  -pc_type  '
 petsc_options_value = ' boomeramg mumps hypre  '
reuse_preconditioner_max_linear_its = 100
scheme = 'bdf2'
  automatic_scaling = true
  dt = 10.0
  nl_abs_tol = 1e-6
  nl_rel_tol = 1e-4
  l_tol = 1e-3
  nl_max_its = 25
  end_time = 5000

[]

[Outputs]
  exodus = true
  [./csv]
    type = CSV
  [../]
[]

[bmurgas-LANL]

Hi @amassaf,
I think the time step you are using are too large: 10 seconds.
Can you try with a smaller time step?
Also, what kind of test are you trying to replicate?

[./tdisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = '-1e-5*t'
  [../]

You can try a compression strain rate by doing:

[./BC_Loading]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = '-0.00001*y*t'
  [../]

Hope this helps.
Brayan

[amassaf]

why '-0.00001yt' ? isn't y=1 ? or is it another y? even with a dt=1 I get convergence issues.