# Multi-Scale Modeling of Pure Vanadium

## Contents |

# Multi-Scale Modeling of Pure Vanadium

## Abstract

Authors: Caleb O. Yenusah ^{a,b}, Abhijith Madabhushi ^{a}, William M. Furr ^{a,b}, Benmbarek Mouhsine ^{a}

a. Department of Mechanical Engineering, Mississippi State University, Mississippi State, MS, USA

b. Center for Advanced Vehicular Systems, Mississippi State University, Mississippi State, MS, USA

## Introduction

Integrated Computational Materials Engineering (ICME) sets a framework for the modeling of materials from the lowest length and time scales to the continuum scale. Bridges between the different scales are established in relation to the nature of the problem, or a continuum scale phenomenon of interest. In this case study, the plasticity of Vanadium is the continuum scale phenomenon of interest. To achieve this goal, Density Functional Theory (DFT) calculations were performed in Quantum Espresso^{[1]} to obtain the energy verse volume curves for BCC, FCC, and HCP Vanadium, along with the coherence energy, elastic constants, generalized stacking fault energy (GSFE), and lattice constants for the stable BCC V phase. To upscale the information generated from the electronic scale to the atomic scale, A second nearest-neighbor (2NN) Modified Embedded Atom Method (MEAM) interatomic potential was generated from these results ^{[2]}. This MEAM potential is used in the atomic scale in molecular dynamics (MD) calculations. The MD calculations can then be used to examine atomic scale phenomena like the movement of dislocations as a function of applied stress, temperature, and other boundary conditions. From the MD simulations, the atomic positions in the lattice under the influence of the dislocation are obtained, from which important variables like the dislocation velocities can be calculated; which in turn can be used to obtain the dislocation mobility. The dislocation mobility will be upscaled to the Mesoscale via Dislocation Dynamics (DD) calculations to study the stress-strain response of the material, especially the hardening behavior.

This exercise is of utmost importance for capturing the critical stress-strain behavior post yielding known as plastic deformation. During the deformation of a material in the plastic regime, dislocation density increases, resulting in the interaction of dislocation, this interaction dictates the hardening behavior of a material. In this study, the Voce hardening rule parameters for a single crystal are determined from DD calculations, which are then used in a Finite Element Code (FEM) called CPFEM^{[3]} to simulate hardening in a single crystal, followed by polycrystalline simulations for three different loading conditions (Tension, Compression, and Torsion) for a realistic quantification of the stress-strain response of the material to different stress states. This stress-strain behavior of the polycrystalline solid is then used to obtain a parameterized internal state variable model (ISV) for the material. The ISV model serves as a bridge to the continuum scale where large systems of engineering significance can be simulated.
During the upscaling process, care is taken to capture the sensitivity of the computed values of the current scale to the values from the previous scale. To accomplish this, a sensitivity analysis is performed by creating three different test cases, starting from the MEAM potential creation. The best fit case is the resultant MEAM potential previously mentioned using the MPC calibration tool^{[4]}. Information from an upper and lower bound potential is also upscaled to the MD simulations to observe the sensitivity. Similar care is taken at the mesoscale to capture the sensitivity of the values from the atomic scale. This is done to capture the uncertainties in the analyses like the uncertainty in the generalized stacking fault energy which is propagated to the MD scale calculations of dislocation mobility, which is upscaled to the mesoscale DD calculations of the hardening behavior, the CPFEM simulation of stress/strain response of single and polycrystalline, and finally, the parameterization of the ISV model.

While the aim of this research endeavor is to obtain a working knowledge of the ICME design paradigm, future work will involve the improvement of the results produced in this work to accurately capture the continuum scale plasticity behavior of pure BCC Vanadium. The following sections describe the methodology involved in this work and the corresponding results obtained for the best fit, upper bound and lower bound cases.

## Methodology

### Electronic Scale Simulations

#### Ab Initio calculations

Ab initio calculations in this work were based on density functional theory (DFT) as implemented in the Quantum Espresso code. Projector augmented-wave (PAW) potential^{[5]} and the Perdew-Burke-Ernzerhof generalized gradient approximation^{[6]} were used. Convergence study was performed in two stages. First to determine the appropriate k-point mesh point needed, the plane wave energy cutoff was fixed at 653.07 eV and k-point mesh were varied. After the converged k-point mesh with respect to the energy cutoff is obtained, a second study was performed using the converged k-point mesh and different energy cutoff ensuring that all results are converged with respect to both parameters. K-point convergence study was also performed for lattice constant and bulk modulus. The results of this study can be found in the Appendix. Energy versus atomic spacing for bcc, fcc, and hcp structures was calculated using the converged energy cutoff and K-point mesh. This is needed for calibration of the MEAM parameters. The lattice constant, bulk modulus, and cohesive energy for the ground state bcc, V was calculated using Murnaghan equation of state (EOS).
The generalized stacking fault energy (GSFE) values (E_{sf}) were calculated using Eq. 1 for calibration of the MEAM parameters for a specific slip system. This ensures that plastic deformation is correctly reproduced by the MEAM potential developed. Due to the fact that vanadium belongs to bcc metal, the (110) plane is the most close-packed plane and the {110}<111> is the most preferable slip system family and should have the smallest GSFE curve. The GSFE curves for (110)[-111] , (110)[001] , and (110)[-110] were plotted. Since the former slip system had the smallest GSFE curve, the MEAM parameters were calibrated to the (110)[-111] slip plane and direction.

where E_{tot} is the total energy of the structure with a stacking fault, N is the number of atoms in the system, ε is the total energy per atom in the bulk, and A is the unit cell area that is perpendicular to the stacking fault^{[7]}.

#### MEAM parameter calibration

MEAM parameter calibration was performed using the MEAM Parameter Calibration software (MPC) developed by Mississippi State University^{[4]}. MEAM parameters were calibrated to the energy verse atomic spacing and GSFE data calculated from DFT and the elastic moduli of bcc V found in the literature^{[8]}. A comparison was made between the DFT calculated energy verse atomic spacing and GSFE curve and results from Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) utilizing the developed MEAM potential.

### Molecular Dynamics Simulations

The MD simulations were performed using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS). A “pad” type structure with BCC Vanadium (V) atoms was created for the purpose of the simulations as shown in the figure below. It was periodic in the X and Z, and non-periodic/shrink-wrapped in the Y direction. An Edge dislocation was introduced in the middle of the cell and constant shear stress was applied on the bottom and top face to propagate the dislocation in the X direction. A constant number of atoms, volume, and energy conditions were prescribed. The temperature of the system was fixed at 300K and equilibrated for 10000ps. The dislocation mobility for the edge dislocation as a function of applied stress and velocity was calculated from the MD simulation results. A convergence study was performed for applied shear stress for the different total number of atoms. The results from the MD simulations serve as the bridging parameters to the mesoscale DD calculations to capture the stress-strain behavior of the material.

### MDDP Simulations

In the next step, Multiscale Dislocation Dynamics Plasticity (MDDP)^{[9]} code was used to create a single Frank-Read source (SFRS) to perform the DD calculations. The code used two input files and an executable (BCCdata.exe) to generate the dislocation input file (DDinput). This input file was modified so that the Frank-Read pair created earlier is modified to a single Frank-Read source. Further modifications to the input file to account for the nearest neighbor atom information was corrected. The results of the MDDP calculations were later extracted and visualized using Tecplot. The process was repeated for a single crystal with multiple Frank-Read sources (MFRS). The material properties used in MDDP calculations are detailed in Table 1. Annealed initial state was assumed for the material with an initial dislocation density in the order of 10^{12}. Preliminary studies were carried out to understand the effect of the number of Frank-Read sources on the initial dislocation density to determine the number of Frank-Read sources needed for a dislocation density in the order of 10^{12}. The MFRS were randomly generated on different planes. For BCC vanadium, all glide planes are of the {110} family. A strain rate of 1.0E5 was used for all simulations. The material properties used in MDDP calculations are detailed in Table 1.

All the values in the above table were obtained from the results of the calculations from the previous scales. The Voce Hardening law parameters were obtained from the results of the MDDP simulations which are upscaled to the mesoscale via the CPFEM simulations.

### CPFEM Simulations

Crystal plasticity finite element method (CPFEM) simulations were run for single crystals and 500 crystals in tension, compression, and shear. The Voce hardening parameters obtained from the MDDP results were used in these simulations, as well as the parameters in Table 2. C11, C12, and C44 are the single crystal elastic constants obtained from previous MPC calibrations, B is the dislocation drag coefficient previously obtained from Lammps simulations.

These values were used in an Abaqus user material subroutine (UMAT) and applied to a single element in tension, compression, and shear. The slip mode examined was <-111> {110} for bcc Vanadium. To evaluate the effects of varying the number of crystals, the best fit case obtained from MDDP simulations was simulated in tension and compared.

### Macroscale ISV Model Calibration

The results of the polycrystalline CPFEM simulations were used to generate the ISV model using the MSU ISV DMGfit 55p v1p1 tool^{[10]}. This was done for the purpose of upscaling the results of the mesoscale onto the macroscale simulations. DMGfit was designed to be an interactive calibration tool to model damage and plasticity using the MSU ISV model. The model is a FORTRAN subroutine implemented in ABAQUS as a UMAT (user material subroutine). The purpose of the calibration tool is to determine ‘material constants’ that will be used as inputs in the finite element simulations. The DMG model itself is specified by 55 constants, which include the dataset and the fitted parameters as shown in the figure below. The upscaled results like bulk and shear modulus and other material properties serve as the input to the calibration tool along with the stress-strain data from CPFEM. The tool uses these inputs and subjects them to an optimization protocol called fminsearch. The relevant Best Fit, Upper Bound and the Lower Bound results of the previous scales like the tension, compression and torsion simulations from the polycrystalline CPFEM simulations are saved as data set files and are subjected to this process and three ISV models are obtained.
The fitted ISV constants are added to an input deck to be used alongside the DMG UMAT to run a single element calculation in ABAQUS for the three loading conditions. The results obtained from the finite element calculation serve as verification of the plasticity-damage calibration.

## Results

### Electronic Scale Simulations

The Shear Modulus, the lattice constant, and the cohesive energy are obtained after the convergence study using the Universal Equations of State (UEOS). Out of the four EOS choices given, BIRCH1 and Keane models were not able to generate results. It seems like the code was unable to converge to a result. From the results generated by BIRCH1 and MURNAGHAN as shown in Table 3, the former yielded a lower lattice constant value but a higher bulk modulus and cohesive energy values. The value generated by MURNAGHAN was closer to Vanadium’s lattice constant as per literature^{[11]} Therefore MURNAGHAN was used for the rest of the analysis performed in this study.

### MEAM Parameter Calibraion

As a prelude to this study, the MEAM potential previously obtained was fine-tuned to accurately capture the behavior of the generalized stacking fault energy (GSFE). The coherence energy, bulk modulus, and lattice constants of the stable BCC V phase were well established in the literature from experiments and first principle calculations. However, the characteristics of the GSFE can have only been reported through first principle calculations ^{[12]}^{[13]}. Therefore, any errors in it are quantified through a sensitivity analysis, since the GSFE is paramount to the plastic behavior of materials. This was performed using three different MEAM potentials, a good fit to the DFT calculations, another that underestimated the GSFE curve by 26% ^{[2]}, and a third one that overestimated the GSFE curve by 22% (Figure 3). These MEAM potentials will hereon be referred to as Best fit, lower bound, and upper bound. It should be noted that the MEAM potential for V from Lee et al [1] underestimated the GSFE curve and this was use used as the lower bound. All the MEAM potential values are shown in Table 4. The values modified in the best fit potential in comparison with the original Lee et al. are A, β(0), β(2), β(3), t(1), t(2), t(3), the attraction and the repulsion parameters. ^{[2]} MEAM potential to better fit the GSFE curve was therefore calculated from DFT. The energy versus atomic spacing for the three MEAM potentials is reported in Figure A1 in the Appendix along with predicted material properties which were compared with experimental results (Table A1).

### LAMMPS Simulations

For the created structure type described in the earlier section, LAMMPS simulations were run to evaluate the dislocation velocity. The LAMMPS output files were then subjected to post-processing using the OVITO visualization software ^{[14]}. The OVITO simulations give the position and time values for the edge dislocation during the simulation. The dislocation velocity is obtained by plotting the position value against time and taking its slope. To select appropriate dimensions for the Pad structure, a convergence study was run to see the behavior of dislocation velocity for applied shear stress of 1500 bar as a function of the number of atoms in the simulation cell. This study gave an optimal number of atoms to run the rest of the simulations. The results of the convergence study and the position vs time plot for the optimal number of atoms (15708) are shown in Figure A3 in the appendix. The value was selected with a view to balancing computational time and the accuracy of the results. A 1-5% range for convergence was considered during the study and the search for convergence was concluded after the last value for the number of atoms shown in Table 3 since a convergence of .457% was attained.

A graph showing the variation in dislocation velocity with respect to the change in the applied shear stress was plotted (Figure 4) and the characteristics were compared with a similar graph from the ICME for Metals textbook (page no: 388, figure (9.7 a))^{[15]}. The figure in the textbook describes the similar effects for an edge location within FCC Aluminum. The velocity scale between the two curves is almost the same however, the velocity of the dislocation in BCC V is lower than FCC Al for the same values of applied stress. Thus, the slope of the FCC Al curve is higher than BCC V. Both curves predict a linear behavior for low applied shear values; however, the linear region in BCC V is between 0-150 MPa which higher than the FCC Al for which the linear region is between 0 – 100 MPa as shown in Figure 4.

#### Calculation of Dislocation Mobility

The dislocation mobility of a metal is a very important property that affects its strength and ductility^{[16]}. The drag-coefficient (B) is a parameter used to quantify the resistance to dislocation motion. Both quantities were evaluated from the velocity vs applied shear stress graph. The drag-coefficient (B) was found by taking the inverse of the slope of the curve in Figure 4 and multiplying it with the Burgers vector of the edge dislocation. Dislocation mobility, in turn, is the inverse of B. The obtained B and dislocation mobility values are shown in Table 5.

#### Sensitivity Analysis

After the evaluation of the dislocation mobility and its drag coefficient (B), a sensitivity analysis was conducted on the previously calibrated MEAM potential by comparing it with two other MEAM potentials. The lower bound potential is generated using the MEAM potential for Vanadium evaluated by Lee et al.^{[2]} It underestimated the maximum GSFE by 26 %. The upper bound was created using the existing Best fit potential and offsetting its maximum GSFE by 22% more than the best fit. The purpose of the sensitivity analysis is to also study the effect of the GSFE on the Mobility of the Edge dislocation. The LAMMPS simulations were performed for the other two MEAM potentials (upper and lower) using 150 MPa applied shear stress value. Figure 5 shows the position vs time curve for all the three cases and their corresponding velocity vs time curves. The slope for the lower bound potential in the position vs time curve was found to be greater than the best fit and the upper bound curve respectively, indicating higher dislocation velocities for the former.

### MDDP Simulations

Figure 6(a-c) shows the dislocation motion under constant applied stress for the best fit MEAM. Figure 6a shows the initial SFRS, the propagation of the dislocation as the strain increases is shown in these figures.

Multiple Frank-Read sources were modeled next to get the dislocation evolution and hardening behavior. Dislocation density evolution for different number of FRS is shown in Figure 7a. The initial dislocation density was 6.30E12, 1.30E13, 1.88E13, and 2.55E13 for 10, 20, 30, and 40 FRS respectively. The 10 FRS was shown for the subsequent calculation for the assumed anneal initial state of the material having dislocation density of the order of 1E12. From Figure 1, forest hardening region was assumed at 7E-8s for 10FRS, while the others have not yet reached a saturation point. The stress/strain response of the single crystal shown in Figure 7b. The stress/strain curve is identical for all the simulations with a yield point between 30MPa to 35MPa. The yield strain was about 3E-2m/m for all initial number of FRS.

#### Voce Hardening

The Voce hardening equation can be written as:

where 𝜅_{𝑠},𝜅_{0},𝑎𝑛𝑑 ℎ_{0} are material properties obtained by correlating the equation with the hardening evolution predicted by DD. Because DD calculations start initially using a random distribution of source segments, the dislocation density calculated at the beginning of the simulation can be ignored and only the forest hardening region, which captures dislocation interactions would be used in the determination of the Voce hardening relationship. The dislocation hardening can be written using the classical Taylor relation:

Here, ρ_{f} is the forest dislocation density, b is the magnitude of the Burgers vector, μ shear modulus, and α is a constant representing an average of the junction strength over all existing configurations. α=0.35 in this work**Cite error: Invalid <ref> tag;
refs with no name must have content**[3]. Using Eq. 2, hardening (κ) as a function of time for various initial FRS is calculated and the results are shown in table XXX along with the fitted Voce equation and material parameters κ_s,κ_0, and h_0. A forest hardening region beginning at about 7E-8s can be observed for 10 FRS. For the others, a forest hardening region could not be identified, and extended simulation times are needed to identify the hardening region. But from a rough fit, κ_0 increases with the number of FRS

### CPFEM Simulations

### Macroscale ISV Model Calibration

## Conclusion

## References

- ↑ P. Giannozzi et al., “QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials,” J. Phys. Condens. Matter, vol. 21, no. 39, p. 395502, Sep. 2009.
- ↑
^{2.0}^{2.1}^{2.2}^{2.3}Lee, Byeong-Joo, et al, "Second nearest-neighbor modified embedded atom method potentials for bcc transition metals," Physical Review B, vol. 64, no. 18, 2001. - ↑ F. Roters, "Advanced material models for the crystal plasticity finite element method: development of a general CPFEM framework.No. RWTH-CONV-144,865th," in Material Science and Engineering, AACHEN, Germany, 2011
- ↑
^{4.0}^{4.1}Carino, Ricolindo L., and Mark F. Horstemeyer. "Case Studies in Using MATLAB to Build Model Calibration Tools for Multiscale Modeling." Applications from Engineering with MATLAB Concepts. IntechOpen, 2016. - ↑ P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B, vol. 50, no. 24, pp. 17953–17979, Dec. 1994.
- ↑ J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized Gradient Approximation Made Simple,” Phys. Rev. Lett., vol. 77, no. 18, pp. 3865–3868, Oct. 1996.
- ↑ M. F. Horstemeyer, Integrated Computational Materials Engineering (ICME) for Metals: Using Multiscale Modeling to Invigorate Engineering Design with Science. John Wiley & Sons, 2012.
- ↑ Alers, G. A. "Elastic moduli of vanadium." Physical Review 119, no. 5 (1960): 1532.
- ↑ Zbib, Hussein M., Tomas Diaz de la Rubia, "A multiscale model of plasticity," International Journal of Plasticity, vol. 18, no. 9, pp. 1133-1163, 2002
- ↑ Mississippi State University, "MSU ICME," Mississippi State University, [Online]. Available: https://icme.hpc.msstate.edu/mediawiki/index.php/DMGfit_55p_v1p1. [Accessed Monday April 2019].
- ↑ Molecular dynamics simulation of vanadium using an interatomic potential fitted to finite temperature properties - ScienceDirect.” [Online]. Available: https://www.sciencedirect.com/science/article/pii/S002231150201019X. [Accessed: 09-Mar-2019]
- ↑ Gui, Li-Jiang, and Yue-Lin Liu, "First-principles studying the properties of oxygen in vanadium: Thermodynamics and tensile/shear behavior," Computational Condensed Matter, vol. 7, pp. 7-13, 2016
- ↑ Zhang, Xingming, et al., "The effects of interstitial impurities on the mechanical properties of vanadium alloys: A first-principles study," Journal of Alloys and Compounds, vol. 701, pp. 975-980, 2017
- ↑ Alexander Stukowski, "Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool.," Modelling and Simulation in Materials Science and Engineering, vol. 18, no. 1, 2009.
- ↑ Horstemeyer, M. F. (2012). Integrated Computational Materials Engineering (ICME) for metals: using multiscale modeling to invigorate engineering design with science. John Wiley & Sons.
- ↑ Kang, Keonwook, Vasily V. Bulatov, and Wei Cai, "Singular orientations and faceted motion of dislocations in body-centered cubic crystals," Proceedings of the National Academy of Sciences, vol. 109, no. 38, pp. 15174-15178, 2012.

## Appendix