Abstract

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.

Introduction

Horstemeyer ^[1] proposed the ICME paradigm as a means to capture the multi-scale phenomenon of a material by bridging information between length scales from the electron scale to the structural (continuum) scale. Multiscale modeling with the ICME paradigm offers the accuracy required to predict failure, reduce design costs, and reduce time-to-market with the simulation-based design paradigm. The fundamental steps required to achieve multiscale modeling with the ICME paradigm are as follows:

• 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.

In this work, we select pure Cr as our material of interest where we wish to model the plastic behavior at the structural scale. First, we downscaled (Step 1) the end-goal of modeling plastic behavior to the electronic scale, where we identified EV and GSFE curves as the quantities necessary to bridge information up (upscaling) to the nanoscale. We calculated the EV and GSFE curves with ab initio DFT simulations that we ran in QE. At the nanoscale, we calibrated the MEAM potential, derived from the Embedded Atom Method (EAM), using the EV and GSFE curves. The MEAM potential accounts for the angular forces between atoms and provides electronic structure and dislocation mobility information that we will need to conduct dislocation dynamics simulations at the microscale. MEAM potentials are applicable with many kinds of crystal lattices including the primary crystal structure of chromium that is body centered cubic (BCC). We used the MEAM Parameter Calibration (MPC) program with MATLAB interface created by Horstemeyer et al. ^[2] in combination with the molecular dynamics (MD) code known as LAMMPS to calibrate our MEAM potential. For both the DFT and MEAM simulations, we verify the results with values from the literature. In the next section, we give an overview of the theoretical models used in this work.

According to the ICME methodology ^[3], capturing plastic deformation of Cr at the macroscale requires the linking of atomistic and DD simulations at lower length scales. A convergence study is performed to determine the appropriate atomic size for our DD simulations using the calibrated MEAM potential and an applied shear stress of 1500 MPa (15000 bar). Then, we vary the MEAM potential to create an upper and lower bound. We used the three MEAM potentials to run Molecular Dynamics (MD) simulations within the LAMMPS program. Seven simulations were performed per MEAM potential with an applied shear stress ranging from 250 to 3000 MPa (2500 to 30000 bar). We visualized the simulation results with the open-source visualization tool “Ovito.” In Ovito, we can study the centrosymmetry parameter, crystal structures, thermal equilibrium, and dislocation behavior. The dislocation behavior was captured and used to create position vs. time plots that we use to calculate dislocation velocity and dislocation mobility. The dislocation mobility was then passed on to run MDDP calculations. We can use the MDDP calculations with a single Frank-Read source to describe the stress-strain behavior. Three stress-strain curves are generated from the MDDP calculations utilizing the three MEAM parameters used to calculate an upper and lower bound for the dislocation mobility.

Materials and Methods

Density Functional Theory (DFT)

Energy vs. Volume (E-V) Curve

We plotted initial energy versus lattice parameter (E-A) and energy versus volume (E-V) curves from points calculated with Quantum Espresso (QE). The range of values was 80 – 120% of the assumed lattice constant of chromium (2.880). We set the initial K-point grid to 3 and assumed a constant energy cutoff value of 30 Rydberg. We show the initial E-A and E-V curves below in Figures 1 and 2, respectively.

KPOINT Convergence

Next, we performed a convergence study on a single atom of BCC chromium to determine which K-point grid offered a good tradeoff between lattice parameter, bulk modulus, minimum cohesive energy, and CPU time. We varied the K-points from 1 to 10 and plotted the results against derived parameters. We held energy cutoff constant at 30 Rydberg through each K-point convergence study. We show the K-point versus lattice parameter, bulk modulus, minimum cohesive energy, and CPU time plots below in Figures 3-6. We determine a K-point grid of 16 constituted convergence of lattice parameter (± 0.00086 A), bulk modulus (± 1.963958 Kbar), and minimum cohesive energy (± 0.000753 eV) without requiring more CPU time. Using our converged K-point grid, we conducted a second convergence study on a single atom of BCC chromium to determine which energy cutoff value offered a good tradeoff between lattice parameter, bulk modulus, and CPU time. We varied the energy cutoff values from 20 to 100 Rydberg and then plotted the results against the derived parameters. We show the energy cutoff versus lattice parameter, bulk modulus, minimum cohesive energy, and CPU time plots below in Figures 7-10.

Equation of State (EOS)

We determined an energy cutoff of 60 Rydberg constituted convergence of lattice parameter (± 0.000144 A), bulk modulus (± 1.623263 Kbar), and minimum cohesive energy (± 0.000896 eV) without requiring more CPU time. Next, we compared two equations of state (EOS), namely, Birch and Murnaghan. We use an EOS to derive the lattice parameter, bulk modulus, and minimum cohesive energy. We show the results for lattice parameter, bulk modulus, and minimum cohesive energy using the Birch EOS below in Figures 11-13. We do not find good convergence for our lattice parameter, bulk modulus, or minimum cohesive energy using the Birch EOS. We show the results for lattice parameter, bulk modulus, and minimum cohesive energy using the Murnaghan EOS below in Figures 14-16. Again, we do not find good convergence for our lattice parameter, bulk modulus, or minimum cohesive energy using the Murnaghan EOS. However, we select the Birch EOS for because the results for the three parameters were closer to experimental values found in the literature.

Final Converged E-A and E-V Plots

Finally, using our converged K-point grid (16 points), and converged energy cutoff (60 Rydberg) values, we created a set of energy versus lattice parameter (E-A) and energy versus volume (E-V) curves for the primary chromium crystal structure (BCC), as well as the secondary (FCC) and tertiary (HCP) crystal structures. We show the final E-A and E-V curves for BCC, FCC, and HCP chromium below in Figures 17-22.

Generalized Stacking Fault Energy (GSFE) Curves

We ran DFT simulations in QE to generate the GSFE curves for pure chromium using a Python script available here. This script can generate GSFE curves along the three BCC slip systems presented in Table 1.

Table 2 shows the inputs to the code including the converged lattice parameter, energy cutoff, and number of k-points from the previous energy-volume curve study. We did not conduct a convergence study for the GSFE curves assuming that the converged values from previous DFT simulations produced reliable results. These calculations also employed the same pseudopotential as the previous DFT calculations.

GSFE Curves

We generated GSFE curves for pure chromium along three BCC slip systems: (110)[-111] (full), (110)[001] (partial), and (110) [-110] (longpartial). We computed sixteen GSFE data points for all three slip systems. Data for the longpartial slip system only spanned half of the slip distance; therefore, we assumed this curve was symmetric based on the symmetry of the BCC structure. We show the GSFE curves (change in energy versus upper atom displacement) in Figure 23 below. The Y-axis represents change in energy (mJ/m²), and the X-axis represents displacement normalized by lattice parameter (2.8497783 Å).

Based on the GSFE curves, the full slip system presents the lowest energy barrier to slip. Therefore, slip is most likely to occur along the full slip system instead of either the partial or longpartial 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.^[5].

Modified Embedded Atom Method Potential Calibration

MEAM potential calibration has been a manual process since the creation of the MEAM model. Horstemeyer et al.^[2] presented a simple MEAM potential calibration methodology with the MPC software tool. In this work, we follow a similar approach where we calibrated BCC, FCC, and HCP energy versus lattice spacing using energy versus atomic volume (E-V) data produced in QE. After loading the E-V curve results from DFT, we calibrate the MEAM model parameter values of equilibrium lattice constant (alat), energy per unit volume (esub), and the constant related to the bulk modulus (alpha) for LAMMPS to match the minimum BCC energy as shown in Figure 25 below.

Sensitivity Analysis

We calculate values far away from equilibrium using the attract and repuls parameters proposed by Hughes et al.^[6] that represent the cubic attraction and repulsion terms in the UEOS. We use single crystal elastic moduli experimental values from Yang et al.^[5] (Table 3) to calibrate the MEAM potential elastic constants by varying the screening parameter, Cmin, scaling factor for embedding energy, asub, and exponential decay factor for atomic density, b0. These variables were shown to have the largest effect on the E-V curve calibration based on parametric sensitivity studies.

Furthermore, calibration to the GSFE curve from DFT is achieved by adjusting the t3 parameter (weighting factor for atomic density). In Figure 26 below, we compare the GSFE curves obtained from DFT and the MPC tool and in Table 4, we summarize the MEAM potential constants after calibration to the DFT results.

Molecular Dynamics

Atomic Position Generation

We generated the atom positions files for dislocation mobility simulations via a Fortran routine available on Atomistic Dislocation Generation. To setup the code, we selected the body centered cubic (BCC) crystal structure to match the crystal structure of Cr. Then, we selected to study an edge dislocation for a Periodic Array of Dislocations (PAD). The atom positions application does not allow for custom lattice parameters with BCC crystal structures. For this reason, we selected the pre-programmed iron, which has the closest lattice parameter (2.85 Å) to chromium. However, the first step of each LAMMPS simulation equilibrates the system using the user-provided MEAM potential, which moves the atoms to the lowest energy state. This allows for a slight tolerance in the initial lattice parameter; simulations can be visually inspected after equilibrium to check that the initial lattice parameter does not adversely affect the resulting structure as was the case in our simulations.

Dislocation Mobility

Dislocation Velocity vs. Applied Shear Stress

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.

We can see the change in dislocation velocity is below 5% after simulation number 4, or 5940 atoms, therefore we consider the result to be converged. For all future calculations, we select simulation 5 as our baseline atom position file having the dimensions of 73 x 34 x 3 and a total of 7446 atoms.

MEAM Parameter Bounds

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.

We can see a good variation in the GSFE curve by ranging the t3 parameter from -1 to 1, where the original t3 parameter was approximately -0.435. From this parametric study, we could further refine the parameters and we selected t3=-1.3 for our new lower GSFE curve, and t3=0.7 for our new upper GSFE curve. The three MEAM parameters were then combined with the atom positions file to run the DD simulations in LAMMPS.

Dislocation Velocity vs. Applied Shear Stress

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.

We can calculate the dislocation mobility from the linear region of the velocity vs. applied shear stress curve where the slope is equivalent to the burgers vector divided by the drag coefficient. We know that dislocation mobility is equivalent to the inverse of the drag coefficient ^[3], therefore, we must only divide the slope of the linear region by the burgers vector to obtain the dislocation mobility. From these results shown in Figure 28, we take the linear region to lie between 0 and 1000 MPa. We can then find the slope of this line and then divide the slope by the burgers vector (2.47e-10) to calculate the mobility. We can also find where this line intersects the X-Axis, which is the Peierls stress. In Table 6 below, we show the dislocation mobility, Peierls Stress, and drag coefficient (inverse of mobility) for each of our MEAM parameter sets.

Dislocation Dynamics

MDDP with Single Frank-Read Source

We generated a single Frank Read Source (FRS) pair within BCC chromium using Multiscale Dislocation Dynamics Plasticity (MDDP). This software illustrates dislocation movement throughout a crystal. The input files included values shown in Table 7. Dislocation mobility was calculated from previous LAMMPS simulations for a lower, mean, and upper MEAM parameter set.

We illustrated dislocation loops from the Frank Read Source pairs using TecPlot 360. Several intervals are depicted in Figures 29 – 40. We ran MDDP until the FRS had completed one loop.

MDDP with multiple Frank-Read Sources

We generated a multiple Frank Read Source (FRS) pairs within BCC chromium using Multiscale Dislocation Dynamics Plasticity (MDDP). This software illustrates dislocation movement throughout a crystal. The input files included values shown in Table 8. Dislocation mobility was calculated from previous LAMMPS simulations for a lower, mean, and upper MEAM parameter set.

Determination of Voce Hardening Law

We calculated the Voce hardening law parameters by using Equation 1^[5].



Where K is hardening, α is a constant relating to strength (assumed to be 0.35^[4]), μ is the shear modulus (110 GPa for chromium), b is the Burger’s vector (2.47E-10 m for chromium), and ρ is the dislocation density. We determined dislocation densities from MDDP and plugged those values into Equation 1 to calculate kappa (κ). We plotted K values against time, and determined initial strength (κ₀), saturation strength (κ_S), and initial hardening rate (h₀) from the decaying exponential function (Equation 2^[5]).



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.

To investigate the effect of different numbers of Frank Read Source (FRS) pairs on the Voce hardening relationship, we ran several simulations with MDDP. Figures 44-46 depict the change in K0, Ks, and h0 parameters with respect to varying FRS pairs. K0 seems to increase steadily while Ks remains the same. The h0 value seems to converge after 8 FRS pairs. We attempted a simulation with zero FRS; however, MDDP was not able to run. With no dislocation source initially, a single crystal has theoretical infinite strength and produces no stress-strain response.

Crystal Plasticity

Single Crystal Stress Strain Behavior

The Voce hardening law material constants are upscaled to the mesoscale via MDDP computations using lower scale mobility information. We simulate hardening of Cr within a single Crystal Plasticity Finite Element Method (CPFEM) ABAQUS model. By understanding dislocation interaction at the microscale, we can determine hardening and plastic deformation at higher length scales. The mechanical response for each MEAM parameter set was nearly identical. Figure 47 depicts the stress-strain data provided from MDDP. Each MDDP calculation ran for 30 time steps. The plotted data shows mostly linear elastic response, but dislocations are present within the crystal.

Figure 48 shows the stress-strain response from ten Frank Read Source pairs within BCC chromium from three different mobility inputs. The upper bound resulted in the largest yield point, which is expected from the lower dislocation mobility. The mean and lower bound simulations yield much sooner due to higher dislocation mobility. The yield point is much higher than chromium’s yield point at macroscale due to size effects at the lower length scales.

A user-material subroutine is implemented in ABAQUS to perform a single C3D8R brick element simulation outputting plots of the normal stress (S33) and logarithmic strain (LE33) components along the z-axis for uniaxial tension and compression. For the case of simple shear, the Von Mises stress is plotted against the logarithmic strain (LE13) component acting in the z coordinate direction on the plane whose outward normal is parallel to the x-axis. The dislocation drag coefficient and hardening parameters for Cr are specified within ABAQUS input files along with other quantities such as Euler angles dictating the crystal orientation. Computations are performed in uniaxial tension, uniaxial compression, and simple shear for lower bound, upper bound, and mean DD simulation results. In the results section we show the stress-strain behaviors for single-crystal chromium evaluated at various levels of plastic strain.

Polycrystal Stress Strain Behavior

In this section we used Abaqus to generate true stress-true strain curves from three different hardening laws for tension, compression, and torsion (simple shear) at a strain rate of roughly 1/s. Theses curves were generated for both a single crystal, single element simulations and a polycrystalline (500 grains), single element simulations. We used an elastic modulus of 210000 GPa and Poisson’s ratio of 0.30 for all the polycrystalline simulations and the single-crystal torsion simulations. Figure 49, Figure 50, and Figure 51 present the tension, compression, and torsion stress-strain curves, respectively. Tensile and compressive stress-strain curves are presented as axial stresses and strains, but torsional stress-strain curves are presented as Von Mises equivalent stresses and strains. Furthermore, the shear strain in the torsion (simple shear) simulations was assumed to be roughly equivalent to the Von Mises equivalent strain.

For every stress state, the lower set of Voce hardening parameters produced the highest yield stress and vice versa. Both the tension and compression simulations were carried out to relatively high strains; however, the physical behavior of Chromium is not necessarily expected to be accurately captured at these high strains. We compared the polycrystalline simulations to the single-grain simulations in Figure 52, Figure 53, and Figure 54.

The single crystal stress-strain curves yielded at lower stresses in tension and compression and higher stresses in torsion. This behavior may be due to the orientation of the single-grain; more single-grain simulations should be run with various Euler angles to understand this factor. The Hall Petch effect may also play a role in the polycrystalline simulations increased strength. In the next section we calibrate an Internal State Variable (ISV) model to the polycrystalline results to model the macroscopic stress-strain behavior of pure Chromium.

Internal State Variable Model

Calibration

In this section we use the stress-strain behavior from the polycrystalline single element CPFEM to calibrate an ISV model with the DMGfit tool and then verify our results with another set of single element analysis in Abaqus. A tutorial for using the MSU DMGfit tool to calibrate an ISV model to experimental data can be found here. With the tool, we can find useful material constants for the ISV model that are used in a finite element calculation as a User Material subroutine (UMAT) in Abaqus. In this work, we fit three ISV models, one for the lower, mean, and upper sets of data we have been collecting. Ideally, the ISV model would capture temperature effects, strain rate dependencies, damage, hardening, and recovery for a range of loading conditions. However, in this work, we are limited to a single strain rate, a constant temperature, and disregard cyclic hardening/recovery and damage. Upon opening the DMGfit tool, we load in a set of tension, compression, and torsion data points and then list out our material dataset parameters. The only parameters we adjusted in the dataset are listed below in Table 10. The period was calculated automatically from the assigned strain rate.

Finally, we imported the stress strain data from the polycrystal single element analysis for calibration. This included data for compression, tension, and torsion up to 30% strain. As we mentioned earlier, in this work we only calibrate the ISV model for the yield stress, hardening, recovery, and loading conditions; therefore, the only parameters we optimized were the Yield Stress (C3), Transition (C5), Isotropic Hardening Modulus (C15), and Isotropic Dynamic Recovery (C13). We list the calibrated parameters for these variables in Table 11 below.

Verification

To verify the ISV model calibration, we updated the input decks from the original polycrystalline CPFEM single element analysis with the new set of dependent variables and material constants. Additionally, to avoid hourglassing errors in the analysis, we changed the element type from C3D8R to a fully integrated C3D8 element. We submitted nine simulations to verify the three ISV models for tension, compression, and torsional loading. In Figures 55, 56, and 57, we show the overlaid stress-strain results for the original single element analysis, the ISV model fit, and the single element analysis with the ISV user material subroutine.

The results were similar for the lower and upper datasets and are not reported here. From Figure 57, we see we need to work on the ISV fit for the torsion including material constants Ca and Cb. However, we are satisfied with the results and believe we could move forward to validate the model by experimentally testing a piece of pure Chromium for comparison to ISV model.

Theoretical Models

MEAM

We use atomistic analysis methods, suggested by Lee et al.^[7], to describe the energy of a system related to interatomic forces and the relative position of atoms. According to the EAM developed by Baskes ^[8], the total energy of a system is composed of the embedding function Fi and pair interaction φ terms shown in Equation 3 below:



where E is energy, ρ_i is background electron density, φ is atomic pair interaction, and R_ij is distance. We describe the atomic pair interactions (φ) at a distance R in Equation 4 below:



where Z is the number of first neighbors, E^u (R) is the background electronic density, and ρ^-0 (R) is energy per atom for the reference structure. We use the universal equation of state (UEOS) or first principles computations to obtain the atomic energy (E^u (R)). We calculate the UEOS parameters of cohesive energy (E_c), nearest neighbor distance (r_e), atomic volume (Ω), and bulk modulus (B) at equilibrium in the reference structure in Equations 5, 6, and 7, respectively.





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 (ρ⁽⁰⁾) and angular (ρ⁽¹⁾, ρ⁽²⁾, ρ⁽³⁾) partial electron density terms:






where the terms α, β, γ represent summations over the x, y, and z coordinate directions. Atomic electron density ρ^α(h) decreases exponentially with increasing distance rⁱ between an atom i and the site of interest as illustrated in Equation 10 that contains constant decay lengths, β^h.



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^h) and the equation for background electron density in Equation 11.




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 ρ^[8].

Voce Law

We define the Voce law for materials exhibiting isotropic work-hardening or stress saturation at elevated stresses or strains in Equation 2 below. This equation represents how the flow stress evolves with plastic work.



where κ_S is the saturation strength. The initial strength and hardening rate are given by &kapps;₀ and h₀, respectively. Moreover, the constant plastic shear strain rate, C, is computed with DD^{[3], [9]}.

Conclusions

In this work, we present DFT and MEAM simulations for pure chromium. We started by generating EA and EV curves and computing the bulk modulus, lattice parameter, and minimum cohesive energy. These values converged at 8 K point grid and energy cutoff of 60 Rydberg. A bulk modulus of 251 GPa, lattice parameter of 2.845 Å, and minimum cohesive energy of -4.10 eV was found for BCC chromium. DFT simulation data is in good agreement with existing literature. First principles computations were linked to atomistics through the calibration of MEAM model constants to DFT results and elastic moduli found in literature. Information obtained using atomistic methods is passed on to dislocation dynamics simulations at the microscale. To visualize dislocation movement within chromium, we produced a structure in LAMMPS with several atoms and a dislocation present. Once we had a converged system without too much image effects, we sheared the top and bottom rows of atoms and received dislocation position and time information. From those results, we calculated the dislocation velocity and calculated the mobility from three MEAM parameter sets. The stress-strain response and dislocation loop movement from MDDP gave us crucial information to feed a crystal plasticity model. After analyzing dislocation formation and movement within a crystal given a mechanical stimulus, we create our crystal plasticity model of Cr through hardening law development.

References

1. ↑Integrated Computational Materials Engineering (ICME) for Metals, 1st ed. John Wiley & Sons, Ltd, 2018. doi: 10.1002/9781119018377.
2. ↑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.
3. ↑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
4. ↑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.
5. ↑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.
6. ↑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.
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.
8. ↑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.
9. ↑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.