×

In this work, we apply multiscale modeling techniques associated with
Integrated Computational Materials Engineering (ICME) in order to capture the
plastic behavior of pure chromium (Cr). We used Density Functional Theory (DFT)
in the open-source software Quantum Espresso (QE) to generate energy versus
volume (EV), energy versus atomic separation (EA), and generalized stacking
fault energy (GSFE) curves along with the lattice parameter, bulk modulus, and
cohesive energy at equilibrium. The Modified Embedded Atom Method (MEAM)
Parameter Calibration tool along with the DFT data was utilized to create a
MEAM potential for Cr. Both sets of simulations were verified by comparison to
values from literature, and a sensitivity analysis was conducted for the MEAM
potential. These lower scale simulations are upscaled to Dislocation Dynamics
(DD) computations. Information is bridged from atomistics to the microscale by
tracking dislocation mobility and quantifying the stress-strain behavior. A
convergence study is conducted on different sized crystal structures and
dislocation mobility calculations with the Large-scale Atomic/Molecular
Massively Parallel Simulator (LAMMPS) program. The mobility is calculated from
seven simulations at varying applied shear stresses from 250 MPa to 3000 MPa,
or 2500 to 30000 bar, respectively. These calculations were run with three sets
of MEAM parameters to find an upper and lower bound for our data. The
post-processing program OVITO was used to visualize the centrosymmetry
parameter, thermal equilibrium, and track dislocation mobility. Then, we
generated position versus time plots for each simulation to calculate the
dislocation velocity, dislocation mobility, and drag coefficient. Three sets of
stress-strain curves were found from a single Frank-Read source in Multiscale
Dislocation Dynamics Plasticity (MDDP) computations.

Keywords: Chromium (Cr), Modified Embedded Atom Method (MEAM), Density
Functional Theory (DFT), Integrated Computational Materials Engineering (ICME),
Generalized Stacking Fault Energy (GSFE), Molecular Dynamics (MD), Dislocation
Dynamics (DD), Large-scale Atomic/Molecular Massively Parallel Simulator
(LAMMPS), Multiscale Dislocation Dynamics Plasticity (MDDP)

Author(s): Cagle, M. S., Fonville, T. R., Kazandjian, S. L., and Sprow, B. T.

**Step 1:**Downscaling, where we define the desired higher-level “effects,” or information needed.**Step 2:**Modeling and simulation at lower length scales (including calibration and verification)**Step 3:**Upscaling, where we pass information up to meet the requirements set at the higher length scale in step 1. Iteration is required to ensure a proper connection.**Step 4:**Validation at the higher length scale through experiment or simulation.

According to the ICME methodology

Figure 1. Initial Energy (eV) versus lattice parameter (A) curve for BCC chromium.

Figure 2. Initial Energy (eV) versus Volume (V) curve for BCC chromium.

Figure 3. Lattice parameter versus K-point grid for BCC chromium.

Figure 4. Bulk modulus versus K-point grid for BCC chromium.

Figure 5. Minimum cohesive energy versus K-point grid for BCC chromium.

Figure 6. CPU time versus K point grid for BCC chromium.

Figure 7. Lattice parameter versus energy cutoff for BCC chromium.

Figure 8. Bulk modulus versus energy cutoff for BCC chromium.

Figure 9. Minimum cohesive energy versus energy cutoff for BCC chromium.

Figure 10. CPU time versus energy cutoff for BCC chromium.

Figure 11. Lattice parameter versus number of plotted points for BCC chromium with the Birch equation of state.

Figure 12. Bulk modulus verses number of plotted points for BCC chromium with the Birch equation of state.

Figure 13. Minimum cohesive energy versus number of plotted points for BCC chromium with the Birch equation of state.

Figure 14. Lattice parameter versus number of plotted points for BCC chromium with the Murnaghan equation of state.

Figure 15. Bulk modulus verses number of plotted points for BCC chromium with the Murnaghan equation of state.

Figure 16. Minimum cohesive energy versus number of plotted points for BCC chromium with the Murnaghan equation of state.

Figure 17. Converged energy versus lattice parameter (E-A) plot for BCC chromium.

Figure 18. Converged energy versus volume (E-V) plot for BCC chromium.

Figure 19. Converged energy versus lattice parameter (E-A) plot for FCC chromium.

Figure 20. Converged energy versus volume (E-V) plot for FCC chromium.

Figure 21. Converged energy versus lattice parameter (E-A) plot for HCP chromium.

Figure 22. Converged energy versus volume (E-V) plot for HCP chromium.

Table 1: Input Parameters for GSFE Curve Python Script

Name | Slip Plane | Slip Direction |
---|---|---|

Full | (110) | [-111] |

Partial | (110) | [001] |

Longpartial | (110) | [-110] |

Table 2: Input Parameters for GSFE Curve Python Script

Parameter | Value |
---|---|

Lattice Parameter | 2.8497783 Å |

Energy Cutoff | 60 Ry |

Number of K-points | 16 |

Element Weight | 51.9961 atomic units^{[4]} |

Smear | 0.06 |

Vacuum Size | 20.0 |

Number of Stacking Layers | 10 |

Figure 23. Change in stacking fault energy (γ) versus displacement normalized by lattice parameter (2.8497783 Å) for the longpartial and full slip systems.

In Figure 24, we compare data from our full slip system GSFE curve with DFT from Yang et al. [5]. The axes in Figure 24 are the same as in Figure 23 except both curves are scaled such that the displacement spans zero to one. Our DFT calculations show close agreement with Yang et al.

Figure 24. GSFE curves from this study and DFT calculations from Yang et al. ^{[5]} along the full slip system show close agreement. Both curves are scaled such that displacement spans zero to 1.

Figure 25. BCC, FCC, and HCP E-V curves after calibration to DFT data.

Table 3: Comparison of experimental and MEAM single crystal chromium elastic stiffness constants (GPa)^{[5]}.

Elastic Constant | Experimental | MEAM |
---|---|---|

C11 | 391.0 | 391.0 |

C12 | 89.6 | 89.6 |

C44 | 103.2 | 103.2 |

Figure 26. Calibrated GSFE versus shift plot.

Table 4: MEAM model parameter values after calibration to DFT simulation data.

Name | Calibrated Value | Name | Calibrated Value |
---|---|---|---|

alat | 2.9891 | esub | 4.3477 |

asub | 0.6146 | alpha | 5.5799 |

b0 | 5.5535 | b1 | 1.0000 |

b2 | 5.0000 | b3 | 1.0000 |

t0 | 1.0000 | t1 | 1.0000 |

t2 | 0.0000 | t3 | 10.6000 |

Cmax | 2.5418 | Cmin | 0.9839 |

attract | 0.0000 | repuls | 0.0000 |

To study the effect of number of atoms on dislocation mobility, we conducted six LAMMPS simulations with a constant applied shear stress of 1500 MPa (15000 bar). We varied the number of unit cells in the X and Y directions to produce the different volumes. We held the number of unit cells in Z constant at 3 unit cells because the edge dislocation is assumed to be infinitely long in this direction. Therefore, Z dimension can be accurately modeled with only a few unit cells. The number of unit cells in X and Y were selected and varied consistently to prevent edge and size effects. The X direction was the largest dimension because the dislocation translates in the X direction. In Table 5 below, we show the atom position sizes and their resulting dislocation velocity.

Table 5: Atom position convergence study parameters and results

Simulation | (X, Y, Z) units | # of atoms | Velocity (nm/ps) |
---|---|---|---|

1 | 11 x 10 x 3 | 330 | 0.416 |

2 | 31 x 14 x 3 | 1302 | 0.243 |

3 | 41 x 20 x 3 | 2460 | 0.133 |

4 | 55 x 36 x 3 | 5940 | 0.102 |

5 | 73 x 34 x 3 | 7446 | 0.104 |

6 | 97 x 46 x 3 | 13386 | 0.104 |

Previously, we calibrated a MEAM parameter with respect to energy vs. lattice parameter, energy vs. volume, elastic constants, and the Generalized Stacking Fault Energy (GSFE) curve. In this section, we use the calibrated MEAM parameter as a mean value and generate a set of upper and lower parameters by varying the GSFE curve. We found the t3 parameter had the largest influence on the GSFE curve; therefore, we conducted a parametric study on this parameter to find a GSFE curve that has a peak value approximately 10% higher and lower than the mean. We show the parametric study and results from the MEAM Parameter Calibration (MPC) tool in Figure 27 below.

Figure 27. Parametric study on the Generalized Stacking Fault Energy Curve (GSFE) with the t3 parameter ranging from -1 to 1.

In this section we study the effect of applied shear stress on dislocation velocity for the three sets of MEAM parameters to determine the dislocation mobility. We conducted seven simulations for each MEAM parameter set to study the effect of applied shear stress on dislocation velocity. The applied shear stresses ranged from 250 MPa to 3000 MPa, or 2500 to 30000 bar, respectively. From the position vs. time results, we can calculate dislocation velocity as the slope of each line. Then we convert angstroms to nanometer by multiplying each result by 0.1. In Figure 28 below, we show the results for all 21 simulations velocity versus applied shear stress.

Figure 28. Velocity vs. shear stress for mean, upper, and lower MEAM parameters.

Table 6: Mobility and Peierls Stresses for three MEAM parameters.

MEAM Parameter | Mobility (1/Pa*s) | Peierls Stress (MPa) | Drag Coefficient (Pa*s) |
---|---|---|---|

Lower | 726 | 204.34 | 0.00138 |

Mean | 277 | 167.41 | 0.0036 |

Upper | 70.7 | 174.82 | 0.0141 |

Table 7: Chromium input parameters for MDDP.

Parameter | Value |
---|---|

Crystal Size (Burger's vectors) | 10000 X 10000 X 10000 |

X-direction | [1 0 0] |

Y-direction | [0 1 0] |

Z-direction | [0 0 1] |

Strain Rate | 1e5 |

Loading Direction | X-axis |

Density (kg/m3) | 7190 |

Shear Modulus (GPa) | 110 |

Poisson’s Ratio | 0.31 |

Burger’s Vector Length (m) | 2.47E-10 |

Temperature (K) | 300 |

Stacking Fault Energy (J/m2) | 1.338 |

Lower Bound Dislocation Mobility (1/Pa*s) | 726 |

Mean Dislocation Mobility (1/Pa*s) | 277 |

Upper Bound Dislocation Mobility (1/Pa*s) | 70.7 |

Figure 29: Initial state for lower bound.

Figure 30: Dislocation movement for lower bound after 2 frames.

Figure 31: Dislocation movement for lower bound after 4 frames.

Figure 32: Dislocation movement for lower bound after 6 frames.

Figure 33: Initial state for mean.

Figure 34: Dislocation movement for mean after 2 frames.

Figure 35: Dislocation movement for mean after 4 frames.

Figure 36: Dislocation movement for mean after 6 frames.

Figure 37: Initial state for lower bound

Figure 38: Dislocation movement for lower bound after 2 frames

Figure 39: Dislocation movement for lower bound after 4 frames

Figure 40: Dislocation movement for lower bound after 6 frames

Table 8: Chromium input parameters for multiple Frank Read Sources in MDDP.

Parameter | Value |
---|---|

Crystal Size (Burger's vectors) | 10000 X 10000 X 10000 |

X-direction | [1 0 0] |

Y-direction | [0 1 0] |

Z-direction | [0 0 1] |

Strain Rate | 1e5 |

Loading Direction | X-axis |

# Frank Read Source Pairs | 10 |

Density (kg/m3) | 7190 |

Shear Modulus (GPa) | 110 |

Poisson’s Ratio | 0.31 |

Burger’s Vector Length (m) | 2.47E-10 |

Temperature (K) | 300 |

Stacking Fault Energy (J/m2) | 1.338 |

Lower Bound Dislocation Mobility (1/Pa*s) | 726 |

Mean Dislocation Mobility (1/Pa*s) | 277 |

Upper Bound Dislocation Mobility (1/Pa*s) | 70.7 |

Where K is hardening, α is a constant relating to strength (assumed to be 0.35

K0 is the starting value of kappa versus time, and Ks is the asymptotic value in which the function converges. The hardening rate (h0) is the fitted slope in Excel of the decaying exponential graph. Figures 41, 42, and 43 shows the fitted curves to the lower bound, mean, and upper bound data from MDDP. Table 9 lists the resulting values.

Figure 41: Fitted Voce Equation to Lower Bound MDDP Data.

Figure 42: Fitted Voce Equation to Mean MDDP Data

Figure 43: Fitted Voce Equation to Upper Bound MDDP Data

Table 9: K0, KS, and h0 Values for Lower Bound, Mean, and Upper Bound Mobility in MDDP.

Parameter | Value |
---|---|

Lower Bound K0 (MPa) | 52 |

Mean K0 (MPa) | 45 |

Upper Bound K0 (MPa) | 37 |

Lower Bound KS (MPa) | 150 |

Mean KS (MPa) | 600 |

Upper Bound KS (MPa) | 700 |

Lower Bound h0 (GPa) | 27.5 |

Mean h0 (GPa) | 25 |

Upper Bound h0 (GPa) | 10 |

Figure 44: Initial strength (K0) value versus # of FRS in BCC chromium.

Figure 45: Saturation strength (Ks) value versus # of FRS in BCC chromium

Figure 46: Hardening rate (h0) value versus # of FRS in BCC chromium

Figure 47. Stress-strain response in lower bound, mean, and upper bound

Figure 48. Stress-strain response from upper bound, mean, and lower bound dislocation dynamics simulation.

Figure 49.Tensile stress-strain curves of a 500 grain, single-element simulation for the upper, mean, and lower Voce hardening law parameters.

Figure 50.Compressive stress-strain curves of a 500 grain, single-element simulation for the upper, mean, and lower Voce hardening law parameters.

Figure 51.Torsion (simple shear) stress-strain curves of a 500 grain, single-element simulation for the upper, mean, and lower Voce hardening law parameters.

Figure 52.Tensile stress-strain curves of the 500 grain simulations and single grain simulations for the upper, mean, and lower Voce hardening law parameters.

Figure 53.Compressive stress-strain curves of the 500 grain simulations and single grain simulations for the upper, mean, and lower Voce hardening law parameters.

Figure 54.Torsional stress-strain curves of the 500 grain simulations and single grain simulations for the upper, mean, and lower Voce hardening law parameters.

Table 10. ISV dataset input parameters

Dataset Parameter | Value |
---|---|

Units | 5 = MPa |

Strain Rate (1/s) | 1 |

Tinit (initial temperature) | 300 K |

G (Shear Modulus (MPa)) | 110000 |

Bulk (Bulk Modulus (MPa)) | 251000 |

Tmelt (melting temperature (K)) | 2180 |

Table 11. ISV material constants

Material Constant | Lower ISV | Mean ISV | Upper ISV |
---|---|---|---|

Yield Stress (C3) | 158.52 | 148.98 | 115.22 |

Transition (C5) | 1.036 | 1.036 | 1.036 |

Kinematic anisotropic hardening (C9) | 0 | 0 | 127.78 |

Kinematic dynamic recovery (C7) | 0.00217 | 0.00223 | 0.00223 |

Isotropic hardening modulus (C15) | 227.95 | 264.91 | 66.85 |

Isotropic dynamic recovery (C13) | 0.0034 | 0.0034 | 0.5146 |

Figure 55.Mean Compression ISV model compared to the verification single element analysis and original.

Figure 56.Mean Tension ISV model compared to the verification single element analysis and original.

Figure 57.Mean Torsion ISV model compared to the verification single element analysis and original.

where E is energy, ρ

where Z is the number of first neighbors, E

The MEAM embedding function comprises the aforementioned terms and the adjustable parameter A as shown in Equation 8 below:

The background electron density is the combination of the spherically symmetric (ρ

where the terms α, β, γ represent summations over the x, y, and z coordinate directions. Atomic electron density ρ

The equations for the partial electron densities are coordinate invariant and equal to zero for cubic symmetric crystals. We represent the angular contributions from the partial electron densities in simplified form in Equation 10 (with constants t

In order to obtain accurate material properties and avoid imaginary electron densities for Γ<-1, we use the functional form of G (shown in Equation 11) where G(0)=1/2. When G(0)=1, the background electron density is non angular dependent and is equivalent to the EAM expression for ρ

where κ

- ↑Integrated Computational Materials Engineering (ICME) for Metals, 1st ed. John Wiley & Sons, Ltd, 2018. doi: 10.1002/9781119018377.
- ↑M. F. Horstemeyer et al., “Hierarchical Bridging Between Ab Initio and Atomistic Level Computations: Calibrating the Modified Embedded Atom Method (MEAM) Potential (Part A),” JOM, vol. 67, no. 1, pp. 143–147, Dec. 2014, doi: 10.1007/s11837-014-1244-0.
- ↑M. F. Horstemeyer, Integrated Computational Materials Engineering (ICME) for Metals: Using Multiscale Modeling to Invigorate Engineering Design with Science. Wiley, 2012. [Online]. Available: http://books.google.com/books?id=KCRLnYj1YiYC
- ↑J. Meija et al., “Atomic weights of the elements 2013 (IUPAC Technical Report),” Pure and Applied Chemistry, vol. 88, no. 3, pp. 265–291, Mar. 2016, doi: 10.1515/pac-2015-0305.
- ↑K. Yang, L. Lang, H. Deng, F. Gao, and W. Hu, “Modified analytic embedded atom method potential for chromium,” Modelling Simul. Mater. Sci. Eng., vol. 26, no. 6, p. 065001, Jun. 2018, doi: 10.1088/1361-651X/aaca48.
- ↑J. M. Hughes et al., “Hierarchical Bridging Between Ab Initio and Atomistic Level Computations: Sensitivity and Uncertainty Analysis for the Modified Embedded-Atom Method (MEAM) Potential (Part B),” JOM, vol. 67, no. 1, pp. 148–153, Dec. 2014, doi: 10.1007/s11837-014-1205-7.
- ↑B.-J. Lee, W.-S. Ko, H.-K. Kim, and E.-H. Kim, “The modified embedded-atom method interatomic potentials and recent progress in atomistic simulations,” Calphad, vol. 34, no. 4, pp. 510–522, Dec. 2010, doi: 10.1016/j.calphad.2010.10.007.
- ↑M. I. Baskes, “Determination of modified embedded atom method parameters for nickel,” Materials Chemistry and Physics, vol. 50, no. 2, pp. 152–158, Sep. 1997, doi: 10.1016/S0254-0584(97)80252-0.
- ↑D. P. Rao Palaparti, B. K. Choudhary, and T. Jayakumar, “Influence of Temperature on Tensile Stress–Strain and Work Hardening Behaviour of Forged Thick Section 9Cr–1Mo Ferritic Steel,” Trans Indian Inst Met, vol. 65, no. 5, pp. 413–423, Oct. 2012, doi: 10.1007/s12666-012-0146-5.