Convergence issues with nonlocal model for Al: Crashes despite fine mesh and small time steps

[libooo2023-beep]

Dear DAMASK developers and community,

I am trying to perform a CP-FEM simulation for Aluminum using the nonlocal plasticity model. I have referred to the parameters from C. Kords' thesis (RWTH Aachen, 2013) and the example configuration provided in the documentation.

• The Problem:

1. At Realistic $\rho_0$ ($10^{6} \sim 10^{14} m^{-2}$):The simulation diverges immediately at the very first increment ("cell failed"), even with very fine mesh and small time steps.
2. At Artificially High $\rho_0$ ($1.25 \times 10^{15} m^{-2}$):To stabilize the run, I tried increasing $\rho_0$ to $1.25 \times 10^{15} m^{-2}$.Result on Voronoi & EBSD RVE: The simulation runs partially but crashes around 90% completion (e.g., step 450/500 for a 50s run).Physical Inaccuracy: I compared the stress-strain curve from this high- $\rho_0$ simulation (EBSD based) with my experimental data. The simulation shows excessive hardening (the slope in the plastic region is almost identical to the elastic modulus), confirming that using such high density is not a valid physical solution, but I cannot lower it without causing immediate divergence.

• What I have tried (Debugging attempts): To rule out numerical resolution issues, I have already tried the following, but none solved the problem:

1. Refined Mesh: I significantly increased the grid resolution for the Voronoi tessellation (more voxels per grain), but the divergence
   persists.
2. Physical Inaccuracy: I compared the stress-strain curve from this high-$\rho_0$ simulation (EBSD based) with my experimental data. The simulation shows excessive hardening (the slope in the plastic region is almost identical to the elastic modulus), confirming that using such high density is not a valid physical solution, but I cannot lower it without causing immediate divergence.
3. Extreme Time Stepping: I reduced the time step drasticallly. For a 15s simulation ($10^{-3} s^{-1}$ strain rate), I set 150,000
   increments. Even with this setup, the simulation failed around increment 1000.
4. Flux Check: Setting plastic: flux: false allows the simulation to run, confirming the instability comes from the flux/transport term.

• Simulation Details:
  Material: Aluminum (FCC), parameters based on Kords (2013), downlod from the damask website.

• My Questions:
  Since mesh refinement and small time steps didn't help, are there specific stabilization parameters for the nonlocal model to handle the initial plastic transient? Or is there a way to overcome my current predicament?

• Attachments: To help reproduce the issue, I have attached the following files:
  geometry.vti ( Voronoi versions)
  material.yaml (The configuration causing issues)
  load.yaml (Simple uniaxial tension)
  Curve_Comparison.png (Plot showing the excessive hardening of the simulation vs. Experiment)

geometry.vti:

import damask 
import numpy as np 
size = np.ones(3)*1e-5
cells = [32,32,32]
N_grains = 30
seeds = damask.seeds.from_random(size,N_grains,cells) 
grid = damask.GeomGrid.from_Voronoi_tessellation(cells,size,seeds) 
grid.save('Gemogrid') 
grid

material.yaml:

import damask 
import numpy as np 
config_material = damask.ConfigMaterial() 
config_material['homogenization']['homo'] = {'N_constituents':1, 'mechanical':{'type':'pass'}} 
config_material['phase']['A'] = damask.ConfigMaterial.load('model_al.yaml') 
rng = np.random.default_rng(40280828)
O_A = damask.Rotation.from_random(30,rng_seed=rng)
config_material = config_material.material_add(homogenization='homo',phase='A',O=O_A) 
config_material.save('material.yaml')

model_al.yaml:

lattice: cF
mechanical:
  output: [F, P, F_e, F_p, L_p, O]
  elastic: {type: Hooke, C_11: 106.75e+9, C_12: 60.41e+9, C_44: 28.34e+9}
  plastic:
    type: nonlocal
    references:
      - C. Kords,
        On the role of dislocation transport in the constitutive description of crystal plasticity,
        RWTH Aachen 2013,
        http://publications.rwth-aachen.de/record/229993/files/4862.pdf
    output: [rho_u_ed_pos, rho_b_ed_pos, rho_u_ed_neg, rho_b_ed_neg, rho_u_sc_pos, rho_b_sc_pos, rho_u_sc_neg, rho_b_sc_neg, rho_d_ed, rho_d_sc]
    N_sl: [12]
    b_sl: [2.86e-10]
    V_at: 0.017e-27 # Omega
    d_ed: [1.6e-9]
    d_sc: [10.e-9]
    i_sl: [60]      # k_2 (lambda_0 in Tab. 7.1)
    f_ed_mult: 0.1  # k_1

    rho_u_ed_neg_0: [1.25e9] # 6e10 / (12*4)
    rho_u_ed_pos_0: [1.25e9] # 6e10 / (12*4)
    rho_u_sc_neg_0: [1.25e9] # 6e10 / (12*4)
    rho_u_sc_pos_0: [1.25e9] # 6e10 / (12*4)
    rho_d_ed_0: [0]
    rho_d_sc_0: [0]

    D_0: 7.e-29
    Q_cl: 0.0 # no temperature dependency
    Q_sol: 2.00272e-19 # 1.25 eV
    c_sol: 1.5e-6       # correct unit?
    f_sol: 2.0          # d_obst in multiples of b
    tau_Peierls_ed: [.1e6]
    tau_Peierls_sc: [.1e6]
    w: 10              # w_k in multiple of b
    p_sl: 1
    q_sl: 1
    nu_a: 50.e+9
    B: 1.e-2
    f_ed: 1.0  # k_3
    h_sl-sl: [0, 0, 0.625, 0.07, 0.137, 0.137, 0.122] # Tab. 3.4
    chi_GB: 0.0      # full blocking at GB
    chi_surface: 1.0 # no blocking at surface
    f_F: 0.0         # no line tension correction
    sigma_rho_u: 0   # no random distribution
    short_range_stress_correction: false
    rho_significant: 1.e6

load.yaml:

solver: 
  # mechanical: spectral_polarization
  # mechanical: FEM
  mechanical: spectral_basic
loadstep: 
  - boundary_conditions: 
      mechanical: 
       dot_F: [[1.0e-3, 0, 0], 
               [0, x, 0], 
               [0, 0, x]] 
       dot_P: [[x, x, x], 
               [x, 0, x], 
               [x, x, 0]] 
    discretization: 
      t: 15 
      N: 150
    f_out: 1

Gemogrid_material_load.zip