Molecular Dynamics on Hf-Nb-Ta-Ti-Zr High Entropy Alloy

Abstract

Classical molecular dynamics simulations were used to investigate the structure and mechanical properties in the equiatomic Hf-Nb-Ta-Ti-Zr high entropy alloy. The open-source code LAMMPS was used to generate alloys with different crystalline lattices to determine the stable structure at 300 K. Alloying elements interacted under the action of the MEAM interatomic potential. The result showed that the alloy stabilizes in body-centered cubic (BCC) structure at 300 K. However, a wide dispersion of potential energy data as a function of atomic separation suggests the coexistence of another crystalline phase. Heating tests indicated a polymorphic phase transformation from BCC to hexagonal close-packed (HCP) at temperatures around 1100 K. Uniaxial tensile tests at a rate of 1×10¹⁰ s⁻¹ along the [001], [110], and [111] crystallographic directions in cylindrical monocrystalline bars at 300 K were conducted. The results revealed a strong anisotropy of mechanical properties. This work provides a microscopic understanding of the mechanical behavior of the multicomponent alloy and aligns with the macroscopic theory of plastic deformation of single crystals.

1. Introduction

For a long period, the development of technological alloys has been grounded on the traditional concept of solvents and solutes. The main goal of this approach was to improve the properties of the solvent element, also known as matrix, by incorporating certain quantities of solutes [1, 2]. In various situations, this addition can lead to the formation of a second phase that modifies several properties. However, this paradigm changed in 2004 when the concept of high entropy alloys (HEAs) was postulated as single-phase alloys [3, 4], where all the elements have similar importance. HEAs are based on the metallurgical concept that several chemical elements form mainly a disordered solid solution phase [5]. The content of each multiple principal element ranges from 5 to 35 atom percent [6, 7]. Due to its chemical homogeneity, HEAs crystallized mainly in one of the simple metallic lattice structures, such as face-centered cubic (FCC), body-centered cubic (BCC), or hexagonal close-packed (HCP). HEAs present a set of superior mechanical properties such as high ductility, toughness, creep, hardness and fatigue resistance when compared to conventional metallic alloys [8]. This outstanding array of physical and mechanical properties [9] opens up possibilities for applications in various fields such as structures [10], refractory [11], cryogenic [12], energy storage [13], biomedical [14], aerospace [15], electrocatalysis [16], among others [17].

Based on the concept of entropy, researchers have identified two types of alloy systems [18]. The first type is named as HEAs due to its high configurational entropy. Usually, these alloys are composed by at least five elements mixed in equiatomic amounts [19], or in compositions close to that. Due to their high configurational entropy of mixing, the most probable microstructure is composed of a single stabilized solid solution phase [20]. However, in some cases, especially when compositions are out of the equiatomic fraction, HEAs do not exhibit single-phase microstructure [21]. Moreover, HEAs can be formed in relatively large compositional fractions. The second type of alloy system is termed as medium entropy alloys (MEAs) [22]. In this case, the alloys are composed by three or four multi-principal elements and similarly to the HEAs, the MEAs show improved mechanical and physical properties in comparison with ordinary alloys. Although the criterion based on the number of multielements is valid for most compositions, there are exceptions to this classification [23].

While HEAs or multi-principal element alloys have been designed with a focus on improving mechanical and chemical properties, it is notable that among the extensive array of alloy systems, the development of lightweight alloys is crucial for industries such as aeronautics and transportation. Furthermore, because over the years various HEAs have been developed, different reviews have been written paying attention to the processing routes [24], microstructure [25], mechanical [9, 26], and chemical properties [27]. Nonetheless, during the last years, diverse computer simulation techniques have been used to study and design new entropic alloys [28, 29]. Even for equiatomic composition, the number of combinations of alloys with four or five elements, considering all metallic elements of periodic table, is enormous. Therefore, the use of previous computational modeling is crucial in the design of entropic alloys.

Refractory HEAs have been identified as promising candidate materials for aerospace propulsion systems, land-based gas turbines, nuclear reactors, heat exchanger tubing [30], and other applications in the chemical process industry. Among those systems, the equiatomic Hf-Nb-Ta-Ti-Zr HEA stands out by possessing outstandingly high strength combined with good ductility at room temperature [31]. Also, that alloy has attracted attention in the fields of hydrogen storage [32].

In this sense, the present work aims to present the application of molecular dynamics simulation on the study of structure and mechanical behavior at room temperature of the equiatomic Hf-Nb-Ta-Ti-Zr HEA due to its several possibilities of applications.

1.1 Classical molecular dynamics

Several semi-empirical methodologies developed for predicting HEAs provide partial success, and their predictions are not entirely satisfactory. Although the resulting parameters offer valuable insights into the phase formation of HEAs, it is not possible to generalize to all alloy systems and find a comprehensive model remains as a significant challenge to the scientific community.

On the other hand, with the advances in the development of powerful computer devices and robust algorithms, various tools now exist to assist the search for the most promising HEA compositions, whether they exhibit a single or multiphase structure. Among those different approaches, we can find the calculation of phase diagrams (CALPHAD) [33], ab-initio molecular dynamics (AIMD) simulation [34], Monte-Carlo (MC) simulations [35], classical molecular dynamics (MD) [36, 37, 38] simulations, and the more recently developed machine learning approach [39, 40, 41]. Each approach has its own advantages compared to experimental research. Notably, cost savings can be achieved as all experiments are conducted within a computer, which is more economical than traditional experimental setups. Nevertheless, it is important to highlight that to obtain satisfactory results is necessary to have access to extensive base data [42, 43] or, regardless the study type, to count with powerful computer capacity.

The MD methodology allows for determining both the equilibria as well as out-of-equilibrium thermodynamic properties. In addition, the dynamical properties of a system at finite temperature are also simulated by MD. MD is an efficient tool for exploring the type of structure, thermodynamics, physical, and mechanical properties, offering a nanoscale understanding of the structure and deformation behaviors of materials. However, the quality of MD simulations largely depends on the method used to specify forces, specifically, the interatomic potential among the atoms or force field function. Thus, the interatomic potential is considered as the heart of the simulations, and the accuracy of the obtained results strongly depends on it.

It is widely recognized that there are currently numerous open-source as well as commercial packages available for conducting MD simulations. A simple digital search on MD open-source packages reveals various options each one with its own characteristics. Nonetheless, it is worth noting that most of scientific papers in the field of materials science use the large-scale atomic/molecular massively parallel simulator (LAMMPS) package [44, 45].

LAMMPS is a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, which can simulate both the MD and MC approaches. Additionally, it is an open-source code that operates on various operational systems, including Windows, Linux, or Mac. Furthermore, although LAMMPS does not have a user-friendly interface its programming is not too complicated. This tool is accessible for both trained programmers and novice researchers.

1.1.1 Interatomic potential

Molecular dynamics simulation is divided into two categories. The first category corresponds to AIMD or first principles simulation, and the second category corresponds to the classical MD. AIMD is a powerful tool with a great accuracy in predicting the parameters and behaviors of different materials. Its principal drawback is the significant amount of time required to obtain the results, in addition to supporting a small number of atoms. Nonetheless, it is strongly recommended to study the physical and chemical properties of materials. On the other hand, classical MD is faster than AIMD, although at the expense of a slight reduction in prediction accuracy. MD is based on Newton’s second law, where its primary control mechanism is the force field or interatomic potential (IP) that governs the interaction of atoms. The IP plays a fundamental role in the MD simulations, and the predictive accuracy depends on how well it describes the interaction among atoms.

All the IPs enable the rapid calculations of the total energy of the atomic systems, facilitating the behavior simulation of metals, ceramics, polymers, and composites. The development of any IP should consider a compromise between computational speed and the accuracy of modeling. For each type of material, it is necessary to consider specific aspects to adequately represent atomic bonding. Hence, for instance, an IP for ceramics should differ from that of the metallic alloy.

IPs also known as force fields, are mathematical functions dependent on atomic coordinates, approximating either the electronic ground state energy or, in the case of metals, the electronic free energy within a system [46]. In MD simulations, the forces acting on atoms, essential for their motion, can be determined by computing the gradient of these functions with respect to the atomic coordinates. Understanding the importance of IP is essential for accurately predicting and interpreting the dynamic evolution of a molecular system.

There are several types of IPs. However, here we are focusing just on the IP used to simulate the equiatomic HEA that corresponds to one embedded atom method (EAM) type. Daw and Baskes [47, 48] proposed an alternative to the pair potential description by employing density functional ideas. This approach, known as the EAM potential, has been applied in MD to characterize diverse and intricate physical scenarios in metallic systems [49]. EAM potentials distinguish themselves from two-body potentials by incorporating a many-body term, dependent on an effective electron density, alongside the pair potential. This many-body term is employed to encapsulate the metallic bonding of the system [50].

In the EAM, each atom is considered to be embedded within the surrounding electron density contributed by neighboring atoms. As a result, the potential energy of a group of atoms is determined by the sum of the pair interaction energy between nuclei of atoms i and j and the embedding energy. The embedding energy is a function of the local background electron density around the i-th atom [51]. The EAM is particularly well-suited for investigating systems featuring crystalline defects such as dislocations and grain boundaries, owing to its capacity to handle local background electron density. However, a crucial assumption in the derivation of the EAM expression is that the electronic cloud around each atom has a spherical shape. While this approximation holds true for FCC crystal structures, the EAM faces challenges in accurately representing systems where directional bonding is crucial, such as BCC, HCP, and nonmetallic materials [51]. The total energy of a system with N atoms within the EAM formalism is given by:

$$U_{\mathrm{EAM}}=\sum _{i=1}^N\left[F_i\left({\rho }_i\right)+\frac{1}{2}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{u}_{\mathit{ij}}\left(r_{\mathit{ij}}\right)\right]$$E1

where $F_i\left({\mathrm{\rho }}_i\right)$ is the energy required to insert atom i into the electron density ${\mathrm{\rho }}_i$ at position i, and $u_{\mathit{ij}}\left(r_{\mathit{ij}}\right)$ is the pair-interaction potential between atom i and j with distance $r_{\mathit{ij}}$. The electron density can be expressed as:

$${\rho }_i=\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N\varphi \left(r_{\mathit{ij}}\right)$$E2

where $\varphi \left(r_{\mathit{ij}}\right)$ is the electron density of atom j at a distance $r_{\mathit{ij}}$ from the nucleus of atom i. Used either independently or in a hybrid configuration, the EAM potential is extensively applied in the study of HEA [30, 52, 53, 54, 55].

EAM potentials demonstrate a limitation in accurately representing non-close-packed structures, specifically those differing from FCC. To enable a single formalism for a broad range of structures and elements, including FCC, BCC, HCP, diamond-structured, and even gaseous elements, Baskes [56] developed the modified embedded-atom method (MEAM). Similar to the EAM, the MEAM formulation includes pairwise repulsions and an embedding function. However, unlike the EAM, the MEAM potential introduces angular-dependent interactions through the electron density term. This feature allows MEAM to model directional bonding [57]. Therefore, the overall formalism in MEAM closely resembles that of EAM (Eq. (1)), with a notable distinction lying in MEA’s effective electron density, which includes angular contributions [57]:

$${\rho }_i^{\left(0\right)}=\sum _j{\varphi }^{\left(0\right)}\left(r_{\mathit{ij}}\right)$$E3

$${\left({\rho }_i^{\left(1\right)}\right)}^2=\sum _{\alpha }{\left[\sum _j{\varphi }^{\left(1\right)}\left(r_{\mathit{ij}}\right)\frac{{\alpha }_{\mathit{ij}}}{r_{\mathit{ij}}}\right]}^2$$E4

$${\left({\rho }_i^{\left(2\right)}\right)}^2=\sum _{\alpha ,\beta }{\left[\sum _j{\varphi }^{\left(2\right)}\left(r_{\mathit{ij}}\right)\frac{{\alpha }_{\mathit{ij}}{\beta }_{\mathit{ij}}}{r_{\mathit{ij}}^2}\right]}^2-\frac{1}{3}{\left[\sum _j{\varphi }^{\left(2\right)}\left(r_{\mathit{ij}}\right)\right]}^2$$E5

$${\left({\rho }_i^{\left(3\right)}\right)}^2=\sum _{\alpha ,\beta ,\gamma }{\left[\sum _j{\varphi }^{\left(3\right)}\left(r_{\mathit{ij}}\right)\frac{{\alpha }_{\mathit{ij}}{\beta }_{\mathit{ij}}{\gamma }_{\mathit{ij}}}{r_{\mathit{ij}}^3}\right]}^2-\frac{1}{3}{\left[\sum _j{\varphi }^{\left(2\right)}\left(r_{\mathit{ij}}\right)\right]}^2$$E6

where ${\varphi }^{\left(l\right)}$ represents the atomic electron densities weighted by the x, y and z components of the distances between atoms (identified by α, β, γ). Note that the spherically symmetric partial electron density ${\rho }_i^{\left(0\right)}$ is the same as the electron density in the EAM (Eq. (2)). Finally, it is worth mentioning that the MEAM, in contrast to EAM, is an empirical extension and has not been justified by strong physical arguments [58]. Though this potential has been used for different simulation systems, recent works usually avoid using it because it considers only first nearest-neighbor interactions using instead a modified version known as second nearest-neighbor modified embedded-atom method potential (2NN-MEAN) that rectifies this by incorporating second nearest-neighbor interactions [55].

2. Computational methodology

2.1 Processing

We conducted MD simulations using the open-source code LAMMPS. The atomic interactions in all the simulations were modeled using the MEAM interatomic potential, as parameterized by Huang et al. [59].

To accurately determine the lattice structure of the Hf-Nb-Ta-Ti-Zr equiatomic alloy, we created single-crystalline samples with BCC, FCC, and HCP structures. Initially, these samples had distinct solid solution configurations composed by a total of 24,000 atoms, with 4800 atoms of each element. Periodic boundary conditions were applied in all simulations along the three orthogonal directions to eliminate surface effects. To obtain the equilibrated structure of the HEA, samples underwent energy minimization at 0 K using the conjugate gradient method. Figure 1 depicts snapshots of the solid solution sample under a BCC structure captured in the stages before energy minimization (a), and after heating and stabilization at 300 K (b). The coloring follows the CNA [60] algorithm, enabling a clear observation of the structural differences of surface-defect atoms in gray color and volumetric ones in blue.

Figure 1.

Snapshot of BCC Hf-Nb-Ta-Ti-Zr HEA (a) before energy minimization step and, (b) after heating and stabilization at 300 K.

After the minimization stage to ensure the stable structure of HEA, the sample was submitted to heating processes until reaching room temperature, set as 300 K, with a heating rate of 1.0×10¹⁰ K/s. In this process, the isothermal-isobaric (NPT) ensemble, where temperature and pressure are controlled by the Nose-Hoover thermostat and barostat, was used. In a second step, another heating was carried out until the sample reached a temperature of 1500 K, while maintaining the same thermodynamic conditions. To achieve a stable structure at that temperature, a thermal stabilization of 0.1 ns was conducted under the control of the canonical (NVT) ensemble. After the heating and stabilization stages, the system underwent controlled cooling at a rate of 1×10¹⁰ K/s under the same conditions as during the heating, until reaching a temperature of 300 K. In addition, the sample was stabilized at this temperature for a period of 0.1 ns to remove possible residual stresses.

Once the equilibrium structure at room temperature was determined, single-crystal samples of the equiatomic Hf-Nb-Ta-Ti-Zr HEA, comprising a total of 113,800 atoms and arranged in BCC crystalline structure, were generated in a cylindrical shape with a diameter and height around of 10.1 and 27.2 nm. The effects of anisotropy on the mechanical behavior of the alloy were studied by adjusting the tensile efforts parallel to the z-axis along the [001], [110], and [111] crystallographic directions. Figure 2 displays a snapshot of the samples in the initial stage of the simulation, before the energy minimization process.

Figure 2.

Snapshots of cylindrical sample with the z-axis parallel to the displacement application. The atoms of Hf, Nb, Ta, Ti, and Zr are represented as spheres with colors red, blue, yellow, pink, and green, respectively.

Subsequently, the samples were thermalized and stabilized at 300 K for 0.1 ns under the control of the NPT ensemble. Samples underwent uniaxial tensile tests in the [001], [110], and [111] crystallographic directions at a strain rate of 1.0 × 10⁹ s⁻¹ and at 300 K. The stresses in the orthogonal directions to the tensile axis were set to zero during deformation.

2.2 Post-processing

The post-processing stage is extremely important in MD studies, producing results that correspond to the properties and characteristics of the materials. The structural characteristics of the material were analyzed using the algorithms: total pair distribution functions (PDF), common neighbor analysis (CNA) [60], and Warren-Cowlay (WC) [61]. In addition, the X-ray diffraction technique, widely used in experimental studies, was utilized in the simulations. The XRD parameters were set to consider radiation from copper Kα with a wavelength of 0.15418 nm and an angular 2θ range from 30 to 100 degrees. The post-processing analysis corresponding to mechanical deformation samples was carried out by applying the displacement extraction algorithm (DXA) [62], as implemented in the OVITO package [63].

3. Results and discussion

3.1 Analysis of the structure of Hf-Nb-Ta-Ti-Zr alloy

In the initial steps of DM simulation, the natural choice is to achieve thermodynamically and structurally stable systems. The structural stability strongly depends on the interatomic potential used, which is a key factor in producing reliable results. The stable structure of high entropy alloys is determined by performing energy minimization in the simulations, while also considering ideal lattices as BCC, FCC, and HCP. The simulations allow the determination of two factors of fundamental importance for the accuracy of the results: the cohesion energy, which depends on potential energy, and the lattice parameter a. As an example, Figure 3(a) shows the curves of potential energy (PE) as a function of the atomic separation distance (r) for the BCC crystalline sample. Figure 3(b) depicts curves fitted to different sets of data obtained from energy minimization stages. It should be noted that the BCC structure produced the lowest PE value among the other structures, then this structure was considered as the stable equilibrium of the equiatomic Hf-Nb-Ta-Ti-Zr HEA.

Figure 3.

(a) Potential energy (PE) curve as a function of atomic separation distance (r) for HEA sample with BCC crystalline lattice. (b) Fitted curves to different sets of data obtained from energy minimization stages.

During the energy minimization stage, the PE vs. r curves generally show data with minimal dispersion in the energy, with increasing or decreasing in the separation distances. However, in the present case, a large dispersion in the PE values was found, which is unexpected, despite the alloy having several components that typically lead to strong distortion of the crystalline lattice. A non-linear fit of the data using a third-order polynomial equation, typically used in similar cases, yields a curve described by the eq. PE = 360.27–302.8r + 82.95r²–7.54r³. The fitting presents a correlation factor (R²) of 0.93, with energy and lattice parameter of −7.010 eV and 3.425 Å, respectively. The estimated standard deviation was around 3%. Furthermore, it can be observed that the fitting curve exhibits two distinct envelopes of data: one above the average fitting curve and the other below it. Taking into account the observable scattering, we considered a possible mixture of phases with BCC lattices. Subsequently, the data were separated into two groups, and new fittings were made, under the same conditions as those for the initial fitting. It is worth noting that this situation is not typical, and great care is required for its use. Figure 3(b) shows both curves, labeled as BCC-1 for the data located above the average fitting curve and BCC-2 for the data located below that curve. The fittings of the two set of data result in the relationship: PE = 363.85–306.0r + 83.92r²–7.63r³ for BCC-1 and PE = 290.66–242.5r + 65.51r²–5.86r³ for BCC-2. The potential energy and lattice parameter values are shown in Table 1. Both fittings exhibit a correlation factor of 0.98, with an estimated standard deviation of less than 2%.

Lattice structure	BCC	BCC-1	BCC-2
Potential energy (eV)	−7.010	−6.769	−6.813
Lattice parameter (Å)	3.425	3.413	3.418

Table 1.

Potential energy and lattice parameters for BCC lattices in thermodynamic equilibrium state.

The results of Table 1 align well with those obtained by Huang et al. [60], who reported a BCC structure for the equiatomic alloy in the same system. However, the simulated values deviate from those obtained by experimental XRD of samples produced by powder metallurgy [64] and subjected to Rietveld refinement. The measured lattice parameter of real samples is 3.4024 Å for the same equiatomic HEA. This slight difference is expected due to variations in methods, such as the system size [65, 66, 67] and the inherent variability in the typical characterization of real samples. Thus, the values obtained in this work by MD are considered acceptable when compared with the experimental data. Furthermore, the good correlation in the fittings of the separated curves indicates a high probability of the material presenting BCC phase separation. This hypothesis was corroborated by analyzing of the simulated XRD and comparing it to experimental data reported by Málek et al. [64].

Figure 4 shows the X-rays diffractograms (XRD) obtained from experiments and simulations under two temperatures. The experimental data were obtained by capturing the XRD curves from Málek et al. [64] using a custom Python code. At 300 K, the simulated diffractogram shows diffraction peaks corresponding to both BCC phases, as determined by Lukac et al. [68], with lattice parameters of a = 3.3197 Å for the BCC-1 phase, and a = 3.4181 Å for the BCC-2. However, the simulated curves do not show the corresponding peaks for the experimentally determined HCP phase.

Figure 4.

X-ray diffractograms at temperatures of 300, and 1200 K together with the curves obtained experimentally by Málek et al. [366] for the Hf-Nb-Ta-Ti-Zr HEA.

In Figure 5(a), the curve of PE vs. temperature is depicted during the heating stage of the equiatomic Hf-Nb-Ta-Ti-Zr HEA. Two approximately linear regions emerge, each characterized by distinct slopes. The inflection point occurs at a temperature around 1100 K, suggesting that a solid-state phase transformation takes place. As the alloy was only heated, it is supposed that the phase transition is likely thermally activated. However, it occurs in a short range of temperatures and very quickly, which allows us to infer that this phase transition is similar to a martensitic transformation.

Figure 5.

(a) Potential energy vs. temperature during the heating stage at a rate of 1×10¹⁰ K/s of the Hf-Nb-Ta-Ti-Zr alloy. (b) Evolution of BCC, FCC, HCP, and other atomic configuration according to CNA algorithm. Both set of data show the fast phase transition from BCC to HCP structures. T_PT corresponds to the temperature of phase transition.

A clearer observation can be made in the curves of Figure 5(b), which display the evolution of BCC, FCC, and HCP structures, along with other atomic configurations following the CNA algorithm. The CNA algorithm can detect various atomic structures, including BCC, FCC, HCP, icosahedral, and other structural configurations. These additional configurations mainly correspond to certain crystalline defects, such as dislocations, twins, and surface atoms, as illustrated in Figure 1, where blue-colored atoms represent volumetric atoms and gray atoms represent surface atomic configurations, among other defects. In Figure 5(b), as the temperature increases, it is possible to observe a rapid decrease in the volume of BCC structures, accompanied by a fast increase in other types of structures. Additionally, at temperatures around 1100 K, the population of the BCC structure decreases to zero, while the HCP structure increases its population to about 37%.

On the other hand, the application of the DXA algorithm at the evolution of the atomic structure in function of temperature shows that initially at 300 K, there are a small population of ½<100 > dislocations. When the sample is stressed, the dislocation density increases rapidly, however, the atom population involved in the other type of structures is higher than that of the involved in the dislocations. Thus, we can infer that many atom fractions are involved in defects such as stacking faults and twins. In temperatures close to the T_PT the density of ½<100 > dislocations becomes to zero, nevertheless, a small fraction of the HCP phase starts to form and its fraction increase rapidly giving rise to the formation of $1/3<1\overline{1}00>$ dislocation.

Figure 6 depicts the analysis of small range order (SRO) by using the Warren-Cowley (WC) parameters (α_ij) in the quinary Hf-Nb-Ta-Ti-Zr HEA after the stabilization stage at 300 K. The α_ij parameters are determined around an atom i within the first nearest-neighbor shell [59]. A positive value indicates that the atom pair is unfavored while for a negative value the atom pair is preferred. The formation of Hf-Nb, Nb-Ta, and Hf- Zr atom pairs are favored. Other atom pairs are not observed, which means that these pairs are randomly distributed inside the material. It is worth noting that in spite of the little difference in fraction of Hf-Nb pairs this result is in agreement with that obtained by Huang et al. [59] using a hybrid Monte Carlo MC/MD simulation. The authors found that in the quinary HEA the main atom pairs are Hf-Nb, Nb-Ta, and Hf-Zr.

Figure 6.

Warren-Cowley parameters for determining the atom pairs formation.

The population of surface atoms on an ideal defect-free sample corresponds to 18.4%. However, when stabilized at 300 K this fraction increases to 22.4%. Consequently, the actual change in the atom population is approximately 15%. On the other hand, the DXA analysis shows that the increase in the population of other atomic configurations is not associated with dislocation formation. In fact, at 300 K, the sample presents just a few <100> dislocations, and this number decreases to zero as the temperature increase. Thus, the increase in atom population associated with other types of structures may be attributed to the formation of twins, which are not clearly identified by the DXA algorithm.

Furthermore, as the material undergoes a thermally induced transformation, it is crucial to determine the crystalline structure at high temperatures. To identify that structure, the sample obtained after the phase transition was subjected to an XRD test, specifically at a temperature of 1200 K. The diffractogram indicates a mixture of crystalline phases, with the matrix corresponding to an HCP structure, while a fraction of the remaining material shows a BCC-1 as the secondary phase (see Figure 4).

In general, metallic materials tend to expand when subjected to heat. Thermal expansion is usually observed as a straight line with a certain slope that indicates the volumetric expansion modulus. In the case of PE vs. T curves, at temperatures below as well as above 1100 K, there is a slight curvature instead of the expected linearity. The curvatures around 1100 K occur due to the mixing of two different phases. At low temperatures, starting from room temperature, the material is composed by a mixture of the two BCC phases (BCC-1 and BCC-2) with very close lattice parameters. Above 1100 K, the phases present are HCP and BCC-2. Although the alloy processed by Lukac et al. [68] using spray technique has both cubic phases, when the powder is consolidated and sintered at temperatures above 1400 K, the alloy no longer exhibits the HCP phase found by MD.

The structural analysis by using PDFs enables us to infer interesting characteristics of the alloy. Figure 7 shows the PDF curves at 0 K, obtained after energy minimization, 300 and 1200 K. The PDF curves show peaks corresponding to the first, second, third, etc., nearest neighbors when considering a central atom as a reference. The curves have the well-established shape characteristic of the BCC structure. Moreover, there are small peaks that are not typically observed in the PDF patterns of BCC structures in pure metals. These peaks may be related to the surface atoms, or possibly linked to other structures present in the material at a small volume fraction. At the temperature of 300 K, the PDF curves show low-intensity peaks with an elongated base, indicating a change in atomic separation distances from the initial structure at 0 K. This change is associated with crystalline lattice distortion and possibly the formation of the second cubic phase. Indeed, the formation of approximately 0.2% of the HCP phase which was detected by the CNA algorithm. At 1200 K, there is a decrease in the intensity of the peaks and a wider base, particularly noticeable in the first peak where there is a coalescence of the two peaks. In addition, the secondary peaks shift to other atomic separation distances, suggesting the formation of a second phase.

Figure 7.

Pair distribution functions of the Hf-Nb-Ta-Ti-Zr HEA, at temperatures of 0, 300, and 1200 K.

The phase transformation from BCC to HCP structures was observed during heating, considering that atomic diffusion in HEAs is relatively slow. Hence, it is thought that the nucleation of the phase transformation caused lattice distortions that influenced the elastic deformation. It should be noted that there is a high percentage of unidentified structures at both temperatures of 300 and 1200 K, constituting 22.4 and 37.0% of the total atoms, respectively. At 0 K, the fraction was of 18.4%, primarily located on the surfaces of the sample, possibly associated with surface defects as these atoms are not fully bonded. Considering that the 18.4% fraction is attributable to the imperfections in bonding, the atomic fraction immersed in the unidentified structures was of 4.0 and 18.6% at 300 and 1200 K.

3.2 Analysis of the mechanical behavior of Hf-Nb-Ta-Ti-Zr alloy

The technological applications of any material strongly depend on its physical, chemical and mechanical properties. For instance, in structural materials, mechanical properties are crucial. Therefore, it is very important to determine their responses to various forms of external stress or forces, such as tensile, compression, shear, or flexural loads. To determine the mechanical properties in metallic materials, it is customary to characterize by tensile tests, then this technique was adopted here.

In the present work, we investigated the effect of single-crystalline directions on the mechanical properties of the Hf-Nb-Ta-Ti-Zr HEA. The elongation was applied in three crystallographic directions, namely [001], [110], and [111] of a BCC structure. Figure 8 displays the stress (σ)–strain (ε) curves obtained from MD simulations in the conditions where the anisotropic effects are clearly established. Thus, the values of ultimate stress (σ_u), elastic modulus (E), and elongation of the same HEA alloy vary with the direction of elongation. The inset of Figure 8 shows the structural evolution of sample in [0 0 1] direction.

Figure 8.

σ–ε curves in single-crystal samples of the Hf-Nb-Ta-Ti-Zr HEA. Displacement applied in the [001], [110], and [111] crystallographic directions. T = 300 K. The inset display the tensile curve of [001] load direction with different stages of structures and mechanical phenomena.

For the BCC Hf-Nb-Ta-Ti-Zr HEA, the [111] direction represents higher compactness, hence the resistance of the sample to deformation exhibits higher resistance. The lowest mechanical strength is observed when the deformation is applied in the [001] direction, reaching an ultimate stress of approximately 300 MPa. Subsequently, the material undergoes a softening, with the stress stabilizing at around 130 MPa. However, the material exhibits a pseudo-elastic behavior with ductility close to 43%. A high strain induces a phase transformation in which the material hardens slightly, reaching an σ_u of 220 MPa, followed by pseudo-elastic behavior up to approximately 70% strain.

Figure 9 depicts the evolution in atomic configurations during the tensile test with a load applied in [001] crystallographic direction. Figure 9(a) shows a snapshot of the unloaded cylindrical sample, i.e., before the tensile test. Figure 9(b) depicts the sample at the instant when the σ_u is achieved. According to the DXA algorithm, at that point, there are some 1/6 ⟨112⟩ partial Shockley dislocations [69], which can give rise to the formation of deformation twins. As the deformation increases and reaches the σ_u, there is a little increase in the density of partial Shockley dislocations, and this phenomenon is accompanied by the appearance of some Hirth 1/3 ⟨100⟩ [70], and a great number of the screw ½<111 > dislocations. In single-component BCC systems, the shortest Burgers vector correspond to the perfect ½<111 > dislocation. In the present work, screw dislocations are present in the sample after the initial stage of homogeneous deformation. With the increase in deformation, the partial dislocations are no longer detected, also, the screw dislocation population decreases. However, according to the DXA algorithm the structure is predominantly BCC until 43% of strain. At that point, the HCP phase starts to appear and the structure increases rapidly, however, there is a remanent fraction of the BCC structure. Furthermore, the sample undergoes a change in shape from circular to elliptical, as shown in Figures 9(c) and (d). The elliptical shape persists during homogeneous pseudo-elastic deformation and is sustained during the fracture process of the material, as depicted in Figure 9. Also, a significant level of deformation leads to the formation of a very thin neck, indicative of high ductility.

Figure 9.

Shapes of Hf-Nb-Ta-Ti-Zr sample during uniaxial tensile test along the [001] direction: (a) before the test, (b) at ultimate stress, (c) onset of elliptical shape, (d) full elliptical shape, and (e) close to fracture.

The analysis of the atomic structure under deformation by using the CNA tool indicates a continuous variation in the atomic fraction belonging to the BCC or HCP structures. Initially, the entire sample is composed mainly of a BCC matrix phase, and with deformation progress, part of it is continuously transformed into the HCP structure. This result is corroborated by the DXA results that show that during the first stage of strain the structure is predominantly BCC, but the phase transformation to HCP structure occurs at 43% of strain, as can be seen on the inset of Figure 8. During a short range of strain, the HCP phase reaches a high fraction, even though the BCC phase still remains as a minority phase. Its evolution persists until the formation of the neck. This phase transformation significantly influences the mechanical behavior of the HEA.

It should be noted that the main dislocation type formed during the deformation process of the BCC phases is the ½ < 111>. After the phase transformation, two distinct types of HCP dislocation emerge: $1/3\mathit{<}\mathbf{1}\overline{\mathbf{2}}\mathbf{10}\mathit{>}$ and $1/3\mathit{<}\mathbf{1}\overline{\mathbf{1}}\mathbf{00}\mathit{>}$, which are not observed in the case of the thermally induced phase transformation.

When the deformation is applied in the [110] direction, there is an increase in the σ_u, reaching close to 540 MPa. However, after a slight drop in the stress, there is continuous hardening until it reaches 980 MPa. This fact is related to the formation of a large number of dislocations, leading to a new phase induced by deformation. The material then undergoes continuous softening until fracture, initiating at approximately 30% of deformation.

The application of the tensile efforts in the [111] crystallographic direction produces a higher ultimate stress, exceeding 1400 MPa. This stress is the highest value reached among the three studied directions. Nevertheless, this elevated stress level is associated with a significant reduction in ductility. In addition, the stress drop occurs relatively fast, when the phase transformation is not well defined, with the stress-strain halo corresponding to the second phase being almost imperceptible. In this configuration, there was no observed presence of necking, and the sample failed with a low ductility.

The ultimate stress, relative ductility, and elastic modulus values corresponding to each crystallographic direction are summarized in Table 2. The anisotropic effect is clearly detected by the MD simulation for the monocrystal. Nevertheless, the elastic modulus is very low, with values lower than those of pure aluminum (E = 61 GPa) and magnesium (E = 44 GPa) [71]. However, it should be noted that these moduli correspond to polycrystalline materials, while in the current work, values are reported for three crystallographic directions of a single crystalline sample. Moreover, the simulated results are always influenced by the strong impact of the chosen interatomic potential. Therefore, any results must be interpreted with a certain level of confidence.

Direction	[001]		[110]		[111]
Phase	BCC	HCP	BCC	HCP	BCC
σu (MPa)	303	219	538	981	1420
Elongation (%)	42	28	14	21	15
E (GPa)	11.0	3.4	16.4	7.8	28.9

Table 2.

Mechanical properties obtained from the tensile curves at 300 K.

The tensile test applied to different crystallographic directions at room temperature shows distinct mechanical behavior. The BCC to HCP phase transformation induced by strain is strongly dependent on the crystal orientation. The [001] uniaxial tensile test produces a lower density of dislocation leading to a significant BCC to HCP transformation, probably due to the emissions of dislocations and twins. Uniaxial test applied in [110] direction allows a partial phase transformation; however, a density of dislocation increases in respect to the [001] direction. The tensile test in the [111] direction reveals a higher dislocation density, with no noted phase transformation. Those results are in good agreement with the study of Hsieh et al. [72], who reported phase transformation induced by stress on Co-Cr-Fe-Mn-Ni HEA with crystal orientation dependence.

4. Conclusions

The equiatomic Hf-Nb-Ta-Ti-Zr high entropy alloy was studied by classic molecular dynamics simulation in a single crystal. The structural features and the mechanical properties were simulated using the LAMMPS code and partial data analysis was carried out in Ovito code. At 300 K, the structure is composed by two BCC structures with very similar lattice parameters. Heating from 300 to 1200 K induces a phase transformation from BCC to HCP structures. The tensile mechanical test used a round sample with a diameter of 10.1 nm and a length of 27.2 nm. The samples were elongated at a strain rate of 1.0 × 10⁹ s⁻¹ in three crystallographic directions: [001], [110], and [111]. The results revealed a strong anisotropy in mechanical properties, with the [111] direction exhibiting greater resistance. The tensile test in [111] direction showed a higher dislocation density; however, no phase transformation was observed during the tensile test in this direction. Moreover, the phase transformation is strongly dependent on the crystallographic direction.

Acknowledgments

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Finance Code 001. The authors also are grateful to FAPERJ – Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (Grants: SEI-260003/014215/2021 and SEI-260003/001582/2022) and CNPq (Grant: 303023/2019-8).