I. INTRODUCTION
The conventional alloy design strategy is commonly based on one element as the principal constituent with other alloying elements added in minor concentrations for the optimization of various properties. In recent years, a new alloy design concept has been introduced, based on several principal elements in equiatomic or near-equiatomic compositions.Reference Murty, Yeh and Ranganathan1,Reference Zhang, Zuo, Tang, Gao, Dahmen, Liaw and Lu2 Due to the effects of high entropy of mixing and low heat of mixing, these new alloys can form a random solid solution in a single-phase structure (fcc or bcc) instead of ordered intermetallic compounds or decomposition into different phases.Reference Miracle and Senkov3,Reference Macdonald, Fu, Zheng, Chen, Lin, Chen, Zhang, Ivanisenko, Zhou, Hahn and Lavernia4 The design of high-entropy alloys (HEAs) requires a radical departure from conventional notions of alloy design. The high mixing entropy enhances the formation of solution-type phases by lowering the free energy of the disordered phase, and in general leads to a simple microstructure.Reference Cantor, Chang, Knight and Vincent5,Reference Otto, Yang, Bei and George6 “HEAs” were named as such because of this effect in 2004.Reference Cantor, Chang, Knight and Vincent5,Reference Yeh, Chen, Lin, Gan, Chin, Shun, Tsau and Chang7
As pointed out in recent reviews,Reference Zhang, Zuo, Tang, Gao, Dahmen, Liaw and Lu2–Reference Macdonald, Fu, Zheng, Chen, Lin, Chen, Zhang, Ivanisenko, Zhou, Hahn and Lavernia4 this concept is beyond the conventional scope of materials design and introduces a new path of developing advanced materials with unique properties, which cannot be achieved by the conventional micro-alloying approach based on only one dominant element. Many HEAs with promising properties have been reported in the literature.Reference Zhang, Zuo, Tang, Gao, Dahmen, Liaw and Lu2 Mechanical properties are particularly interesting, with attractive combinations of strength and ductility. For example,Reference Liu, Wu, He, Nieh and Lu8 a high-entropy Fe–Co–Ni–Cr–Mn alloy with a single face-centered cubic (fcc) phase was synthesized and subsequently annealed at different temperatures to systematically investigate the grain growth behavior. The hardness of the alloys was found to follow the classical Hall–Petch strengthening, with a relatively high hardening coefficient.Reference Liu, Wu, He, Nieh and Lu8,Reference Zaddach, Niu, Koch and Irving9 Glutovitz and co-workersReference Gludovatz, Hohenwarter, Catoor, Chang, George and Ritchie10 recently examined a five-element HEA, Cr–Mn–Fe–Co–Ni, which forms a single-phase fcc solid solution. They found it to have exceptional damage tolerance, tensile strength, and fracture resistance. This alloy’s mechanical properties actually improved at cryogenic temperatures. The authors attribute this to a transition from planar-slip dislocation activity at room temperature to deformation by mechanical nanotwinning with decreasing temperature, which results in continuous steady strain hardening.
Some basic guidelines for the formation of single-phase structures have been investigated recently.Reference Otto, Yang, Bei and George6,Reference Zhang, Lu, Ma, Liaw, Tang, Cheng and Gao11–Reference Poletti and Battezzati16 However, guidelines for searching alloys that exhibit optimum properties are not available because of particular challenges, difficulties, and complications in performing either controlled systematic experiments or accurate theoretical analysis. In particular, lattice distortion in these alloys is important because it can control solid solution strengthening and overall mechanical behavior.Reference Diao, Feng, Dahmen and Liaw17 Lattice distortions can also significantly affect dislocation core structure and motion. Lattice distortion is induced by atomic size mismatch and has been investigated recently as a possible explanation for the particular properties of HEA alloys.Reference Wang, Wu, Liu and Wang18–Reference Oh, Ma, Leyson, Grabowski, Park, Kormann and Raabe20 Lattice distortion can also be influenced by short-range order, and these effects have been studied in a variety of alloys.Reference Zhang, Jin, Xue, Lu, Olsen, Beland, Ullah, Zhao, Bei, Aidhy, Samolyuk, Wang, Caro, Caro, Stocks, Larson, Robertson, Correa and Weber21–Reference Zhang, Stocks, Jin, Lu, Bei, Sales, Wang, Beland, Stoller, Samolyuk, Caro, Caro and Weber23
Computational work has been performed using a variety of techniques to address various aspects of the stability and properties of these alloysReference Gao, Gao, Hawk, Ouyang, Alman and Widom24 As in many other classes of materials, to realize the potential of HEA alloys, simulation work is criticalReference Pollock and LeSar25 in providing insight into the mechanisms of deformation and fracture in these new alloys. However, this poses great challenges.Reference Zhao, Weber and Zhang26 For example, dealing with disordered solid solution using first principles density functional theory (DFT) calculations is a major challenge when the system contains many components due to the large number of possible atomic configurations. To treat disordered solid solutions using DFT methods, special methods are required, such as the special quasirandom structure method.Reference Wei, Ferreira, Bernard and Zunger27 In addition, to achieve our objective of understanding trends in deformation and fracture behavior, one requires a large number of atoms in the simulation as well as a large number of simulations, which have to be performed in a high throughput manner. It appears that high throughput simulations of these disordered single-phase multicomponent alloys are not within reach using first principle techniques.Reference Toda-Caraballo, Wrobel, Nguyen-Manh, Perez and Rivera-Diaz-Del-Castillo28 The study of mechanical behavior requires simulations with large numbers of atoms that are only accessible through empirical interatomic potentials. The challenge then is that the development of accurate empirical interatomic potentials that can be said to represent a particular multicomponent alloy is also a major task and the results would be of questionable accuracy. Even for binary alloys, this is challenging and questions about the actual accuracy of the empirical force laws to represent a specific alloy system always remain. These questions would make it impossible to trust the actual comparison of simulation results with experiments for a specific multicomponent HEA. Historically, the value of empirical interatomic potentials has been in understanding trends in the results. For example, the understanding of the inverse Hall–Petch phenomenon that came from simulations was the result of studying trends with grain size and observing new deformation mechanisms that appear as the grain size is decreased.Reference Van Swygenhoven, Spaczer, Caro and Farkas29 Other examples investigated include trends in dislocation core structure with ordering energy in B2 alloys.Reference Vailhe and Farkas30 The technique yielded an understanding of how deformation was influenced by the alloy ordering energy, through changing the details of the dislocation core structure.Reference Vailhe and Farkas30,Reference Vailhe and Farkas31 More recently, we have used the technique to show how the yield and flow stress observed in an fcc alloy is clearly related to the unstable stacking fault energy characteristic of the potential.Reference Smith and Farkas32
For the case of HEAs, one possible path is to utilize “toy model” interaction potentials that can be developed to depict a trend in the values of significant alloy properties. While these model interactions will not accurately represent a particular alloy, they are the ideal way of studying trends indicating how certain material parameters influence final properties, and therefore guide the design process. This is the approach taken in the present work. In this paper, we report the development of one such set of embedded atom method (EAM) model interaction potentials with core characteristics similar enough to HEAs that are likely to be relevant for the study of trends in alloy behavior. We present a potential set for a five-component alloy in an fcc structure. The set is characterized by atomic sizes that differ in no more than 3% and is developed to have binary heats of mixing of no more than 0.7 kJ/mol, well within the range of existence of these alloys (5 < ΔH mix < −15 kJ/mol).Reference Guo, Hu, Ng and Liu33 Parameters used here are such that the set resembles (in the sense of basic elastic and thermodynamic properties and size mismatches of binaries) fcc mixtures of Fe–Ni–Cr–Co–Cu, but does not necessarily model other properties that were not used in the fitting procedure (e.g., defect energetics). We then report a detailed study of the lattice distortions observed in this model quinary alloy. We study these distortions at a temperature of 1 K to separate them from the thermally induced distortion. Simulations at 300 K are also performed and the two are compared, to separate lattice distortions that are intrinsic of the near-equiatomic mixture from the dynamic lattice distortions introduced by temperature.
II. METHODS
Even though the present work centers on developing “toy model” interaction potentials, the approach utilized is similar to that used in creating reliable interatomic potentials for monoatomic metals.Reference Mishin, Farkas, Mehl and Papaconstantopoulos34 We use a basis set of cubic splines and a fitting procedure designed for developing the potential functions to reproduce desired values a 0, E 0, B, c 11, and c 12, unrelaxed values of the vacancy formation energy, stacking fault energy, relative phase stability, and surface energies for the pure elements. The procedure ensured that the fcc structure was the stable structure for all five components. The electronic density and embedding functions for the pure elements were set in the effective pair representation, i.e., the electronic density at equilibrium is equal to 1 and the derivative of the embedding function equals zero, to ensure that alloy energetics could be fitted via the cross pair term of the potential alone.
The main requirement for the interatomic potential set to be useful in understanding the properties of HEA alloys, we are interested on, is that the single-phase fcc structure for the equiatomic random mixture be the equilibrium phase. To accomplish this, we developed the mixed pair interactions as averages of the pair interactions of the two components of the pair. The idea was that this procedure results in very small heats of mixing and consequently nearly ideal solution behavior, maximizing the entropy. The fitting procedure was automated, building on previous work in the development of interatomic potentials.Reference Mishin, Farkas, Mehl and Papaconstantopoulos34–Reference Farkas, Schon, DeLima and Goldenstein38 We stress that we are not trying to model a particular alloy exactly but rather establish toy models to study the effects of the critical potential characteristics.Reference Mishin, Farkas, Mehl and Papaconstantopoulos34
Atomistic simulations to test the potentials were performed using LAMMPS molecular dynamics software.Reference Plimpton39 Simulations were carried out at 1 K and using random configurations of 4000 atoms for the binary mixtures and 16,384 atoms for the quinary alloy. Simulations at 300 K were also performed for comparison. Model structure visualization and statistical analysis of lattice distortions was performed using OVITO software.Reference Stukowski40
III. RESULTS AND DISCUSSION
A. Component potentials
The five-component potentials developed are meant to resemble Fe, Ni, Cr, Co, and Cu. While we are not claiming to describe a particular element or alloy exactly, our goal is to make a model that is as similar as possible to a real HEA alloy, within the limitations of the method used here, that is the EAM.
The data for fitting the potentials correspond to stable fcc structures. For the case of components 2 and 5, these data are the actual experimental properties of fcc Ni and fcc Cu, respectively. Component 1 was developed to resemble fcc Fe and the data were obtained from the values known for austenite. We note that the potential predicts the fcc phase to be stable even at 1 K. For components 3 and 4, the data were taken from available first principles calculations for the fcc phases of Cr and Co.Reference Hammerschmidt, Bialon, Pettifor and Drautz41
Table I lists the basic properties given by the potentials for the five-component elements in the potential set developed. For each property, Table I also lists the average values of the various properties for all five components and the standard deviation within the set. The largest variation among the properties is in the values of the stable and unstable stacking faults. The average values for the elastic constants that our potentials predict are within the range of predictions for fcc HEAs using first principles calculations.Reference Tian42 Figure 1 shows the pair components of the potentials developed.
Basic characteristics of the five-component potentials developed.
Pair components of the five-component potentials in the set.
B. Mixed interaction potentials
The objective here is to develop a potential set that would represent nearly ideal mixtures of the component. The approach to achieve this was to develop the mixed pair functions as a simple average of the component pairs to obtain very small heats of mixing. Table II lists the properties obtained for the ten random fcc equiatomic binary compounds. As expected from the construction of the mixed pair potentials, all binary random mixtures indeed show very small heats of mixing and very small deviations from Vegard’s law. The numbers in Table II correspond to random solid solutions since no Monte Carlo simulations were performed. However, because the heats of mixing that the potentials predict are small, we believe that this is a reasonable assumption.
Properties of the binary random equiatomic mixtures.
C. Quinary random fcc alloy
The quinary fcc alloy was equilibrated at 1 K using a block of 16,384 atoms with random configurations. The equilibration was performed for 10,000 molecular dynamics steps of 1 fs. Table III shows the 0 K properties of the relaxed equimolar five-component random alloy using the potentials developed. The fcc phase is indeed the equilibrium phase. The lattice parameter obtained is 0.3552 nm, very similar to the average of the five components of 0.3556. It also compares well with the experimentally measured lattice parameter of an Fe–Ni–Cr–Co–Cu HEA alloy of approximately 0.357 nm. Lattice parameters for similar alloys have been measured as 0.3551.Reference Owen, Pickering, Playford, Stone, Tucker and Jones43 The overall heat of mixing of the alloy is very small, 0.002 eV. To reach this value, we calculated the heat of mixing of the quinary alloy using different seeds for the random distribution of the five elements. This resulted in a more accurate value of the heat of mixing and confirmed that it is indeed a low value suggesting near random structure for the alloy. Using the different ransom distributions also allowed us to estimate the effect of different elemental distributions on the heat of mixing as the standard deviation of the heats of mixing obtained for the various distributions. This value is 0.003 eV.
Integral properties of the random equiatomic quinary alloy.
The overall cohesive energy of the fcc random solution was lower than random alloys in the bcc or hcp structures. This confirms that the potentials predict a stable fcc phase with nearly ideal solution behavior. The energy difference that our potentials predict between the fcc phase and the hcp phase is 0.017 eV/atom (1.7 kJ/mol), which compares well with thermodynamic calculations using the Calphad approach, ranging from 0.4 to 2 kJ/mol for similar HEA alloys.Reference Li and Raabe44 The nearly ideal behavior suggests that the structure of the alloy will indeed be close to random, and therefore, in the present work, no Monte Carlo studies of possible ordering were performed. A melting and quenching simulation was also performed using the potentials developed here, resulting in a metastable glassy structure. This metastable structure had an energy at 1 K that was 0.14 eV higher than the fcc phase.
D. Lattice structure and thermodynamics
Figure 2 shows the displacements from perfect lattice positions at 1 K that these potentials predict. The observed lattice displacements average about 0.0027 nm, representing distortions of about 0.7%. These distortions are of similar magnitude to those found using first principles calculations in a Fe–Ni–Cr–Co–Mn alloy.Reference Oh, Ma, Leyson, Grabowski, Park, Kormann and Raabe20
Lattice distortions visualized as atomic displacements from perfect lattice fcc positions. The displacements have been magnified 30 times. Red, pink, green, gray, and blue represent atomic types 1 through 5 ( ∼Fe, ∼Ni, ∼Cr, ∼Co, ∼Cu).
Figure 3 shows the potential energy distribution of the individual atoms in the alloy. For each component, the energy distribution shows standard deviations of 0.2 eV. The peaks corresponding to components 1 through 4 (Fe, Ni, Cr, Co) are very close, corresponding to close values of their energies in the pure state. Component 5 (Cu) is characterized by lower potential energies, also corresponding to the lower cohesive energy of the pure component. The average potential energies for each component in the quinary alloy were computed from the individual atom potential energy data and allow estimation of the deviations from ideality for each individual component. The calculated partial molar heats of mixing for each component are shown in Table IV. Components 1 and 2 show positive partial molar heats of mixing, whereas components 4 and 5 show negative values.
Distribution of potential energies in the alloy.
Partial molar properties of the random equiatomic quinary alloy.
The distribution of bond lengths was studied by calculating the radial distribution function. Figure 4 shows the first peak of the radial distribution function in the alloy, corresponding to nearest neighbor bonds. At the 1 K temperature, the width of each peak in the radial distribution function in a perfect fcc structure should be nearly zero. The peak shown in Fig. 4 calculated at 1 K for the quinary alloy presents a nonzero width, which is entirely due to lattice distortion and not due to thermal vibrations. The width at half maximum is 0.005 nm, attributable to the complex composition of the alloy.
First peak of the radial distribution function, showing a full width at half maximum of 0.005 nm.
E. Lattice strains and individual atomic stresses
These lattice distortions imply a strain on the lattice that can be calculated and analyzed statistically. Figure 5(a) is the calculated volumetric (hydrostatic) part of the strain tensor for the individual atom types in the alloy. As expected from the atomic sizes, the element 5 atoms (Cu) are larger and therefore present a positive volumetric strain of about 0.5%. This is compensated by volumetric strains of the opposite sign in the other four components. The strains of the other four components are similar, as would be expected from their respective atomic sizes. In Fig. 5(b), the averages of the atomic volumetric strains are shown resolved by the component.
Volumetric strain of the components in the quinary alloy. (a) Histogram of the distribution. (b) Average values for each component.
Using a Voronoi construction, the atomic volume occupied by each atom can be estimated. Combining this estimate with the force calculated on each atom, the actual stress state of each individual atom can be obtained. We note here that while each atom is under stress, the overall stress on the sample is zero after the relaxation performed. Figure 6 shows the resulting calculated distribution of atomic stresses on the individual atoms in the alloy. The finding is that individual atomic stresses due to intrinsic lattice distortions can be in the range of several GPa, either in tension or compression. The standard deviation of the stress distribution is of the order of 2 GPa. As expected, the atoms of the slightly larger size component 5 (∼Cu) are under compression. The same is true for the element 4 atoms (∼Co) and element 1 (∼Fe), even though it would not be expected from the atomic sizes. These compressive stresses are compensated by tensile stresses in the other two component elements (∼Ni and ∼Cr). These local stresses may have a significant effect on the structure and motion of dislocation cores in these materials.
Individual atomic stress of the components in the quinary alloy. (a) Histogram of the distribution. (b) Average values for each component.
F. Dynamical behavior at room temperature
Simulations were also performed at 300 K to study the lattice dynamics. Figure 7 shows the radial distribution function obtained at that temperature. The peaks are much broader than the ones at 1 K. The full width at half maximum is 0.02 nm as opposed to the 0.005 nm found at 1 K. Figure 7 also shows the radial distribution function obtained for pure Ni at 300 K. The full width at half maximum obtained here at 300 K compares well with that measured for a similar alloy,Reference Owen, Pickering, Playford, Stone, Tucker and Jones43 which is also 0.02 nm. The differences between the radial distribution in the quinary alloy and that in pure Ni at 300 K are only due to the slightly different lattice parameters, but the dynamical behavior is essentially the same. Overall the dynamic distortions due to thermal vibrations are about four times larger than the intrinsic lattice distortions due to the multicomponent nature of the alloy. This result indicates that the intrinsic lattice distortions due to the complexity in composition of this alloy may be less important than those introduced by the thermal vibrations. The comparison in Fig. 7 suggests that the dynamic behavior of lattice distortion is not significantly different in this alloy than in a pure component. This is consistent with the fact that the intrinsic lattice distortions are much smaller in this alloy than the thermal distortions at 300 K. This is true in the system studied here, which is characterized by small size differences. HEA alloys containing components with more significant size differences are expected to present larger intrinsic lattice distortions that may be of the same order as the thermal fluctuations at room temperature.
Radial distribution function calculated at 300 K for the quinary alloy (black). The corresponding data for pure Ni (red) are shown for comparison.
IV. CONCLUSIONS
The important goal of modeling work in HEA alloys is to help understand how the remarkable properties of these alloys emerge from correlations of the atomic constituents, and how we can control these properties. This would break new ground by understanding basic mechanisms of damage and failure in this new and potentially revolutionary class of materials. The potentials presented here are adequate for simulations of mechanical behavior of a solid solution quinary alloy with nearly ideal behavior and atomic size differences of around 3%. The potentials present much more significant variation in stacking fault energies, enabling investigation of these effects. This will be most likely to be of use in alloy design. The relationship between the atomic-level chemical complexity and lattice distortion can be different from traditional alloys due to the presence of the intrinsic lattice distortion and local stress fluctuations that are shown here. This can have profound effects on the mechanical properties as well as radiation resistance.Reference Zhang, Zhao, Weber, Nordlund, Granberg and Djurabekova45
We have used the potentials developed to study lattice distortion. Molecular dynamics simulation techniques allow the separation of the static lattice distortion from the dynamic effects of thermal vibration. The detailed analysis shows that the static distortions are of the order that is expected from the size difference in the component alloys. In the case presented here, the relatively small size differences among the alloy components results in static distortions that are smaller than the dynamic distortions observed at room temperature.
The potential set that we developed here is expected to be adequate for general investigation of mechanical response in a model system. Because the individual components were found to be under significant individual stresses, these effects are likely to be significant. Another important point is that the component potentials of the set presented here present a large variation in stacking fault energies, both stable and unstable. This makes this particular set useful to investigate the effects of a significant variation in these properties on mechanical behavior. The stacking fault value will be extremely dependent on the local atomic configuration in the stacking fault region. The specific local stacking fault values are very important for the study of dislocation behavior and therefore mechanical response.
The results presented here indicate that the approach of developing “toy model” interatomic potentials can indeed be valuable. Using this approach, sets could be developed that address more significant size differences, or more significant deviations from ideality. In the latter case, Monte Carlo simulations would be required to include the effects of short-range order.
ACKNOWLEDGMENT
This work was supported by the National Science Foundation, Division of Materials Research, under Grant No. 1507846.
