A systematic investigation on quaternary NbTiZr-based refractory high entropy alloys using empirical parameters and first principles calculations

The overall research route of present work is shown in figure 1. Firstly, the X element in the quaternary NbTiZrX alloy system is determined according to some simple rules based on the properties of the element, such as the boiling point, the ratio of melting point to the density. Secondly, some commonly used empirical parameters, such as the atomic size mismatch, the average mixing enthalpy and Ω, are employed to predict the phase structure and stability of each alloy system. The alloy systems with single-phase solid solutions are left for the follow-up study. Thirdly, first principles calculations based on SQS method are used to research the elastic properties of all NbTiZrX alloy systems screened out from previous step. By analyzing these elastic properties, the stability and mechanical performance of each NbTiZrX alloy can be evaluated. Finally, the quaternary NbTiZrX alloy systems with stable single-phase solid solution structures and excellent mechanical properties are obtained, which can be regard as potential RHEAs for further investigation. A more detailed introduction to the theories or methods used in present work is given below.

Zoom In Zoom Out Reset image size

According to the definition of HEA, there is no provision for the phase structure of the alloy, which means that the corresponding phase structure can be either a single-phase solid solution, or an intermetallic compound, or a combination of the two. Nevertheless, most HEAs reported in the literatures are single-phase solid solutions, or it can be said that researchers expect to obtain single-phase solid solution alloys. The main reason for such phenomenon is that HEAs with a single-phase solid solution generally have better comprehensive mechanical properties. In addition, single-phase solid solution HEAs are more convenient to study. Hence, it's particularly important to predict the phase structure of HEA based on its composition. In the past few decades, many methods have been proposed to solve this problem, such as CALPHAD method [12–16], empirical parameter method [17, 18] and machine learning method [19]. Among them, the simplest and most commonly used method is the empirical parameter method, which is also used in present work.

Researchers have proposed many different empirical parameters according to different properties such as atomic radius, melting point, electronegativity, mixing entropy, etc. In present work, the most commonly used empirical parameters are chosen as the criteria for predicting the phase structure or stability of the quaternary NbTiZr-based RHEAs, including the atomic size mismatch (δ), the mixing entropy (ΔS_mix), the average mixing enthalpy (ΔH_mix), the average melting point (T_m), Ω and the average valence electron concentration (VEC) [20]. To maintain the integrity of present paper, the definitions of these parameters are given directly below:

$$\begin{equation}\delta =\sqrt{\sum\limits _{i=1}^{n}{x}_{i}{\left(1-\frac{{r}_{i}}{a}\right)}^{2}}\end{equation} \tag{ 1.1 }$$

$$\begin{equation}{\Delta}{\mathit{\text{S}}}_{\mathrm{m}\mathrm{i}\mathrm{x}}=-R\sum\limits _{i=1}^{n}{x}_{i}\enspace \mathrm{ln}\enspace {x}_{i}\end{equation} \tag{ 1.2 }$$

$$\begin{equation}{\Delta}{\mathit{\text{H}}}_{\mathrm{m}\mathrm{i}\mathrm{x}}=4\sum\limits _{i=1,j\ne i}^{n}{x}_{i}{x}_{j}{H}_{ij}\end{equation} \tag{ 1.3 }$$

$$\begin{equation}{\mathit{\text{T}}}_{\text{m}}=\sum\limits _{i=1}^{n}{x}_{i}{T}_{\text{m}i}\end{equation} \tag{ 1.4 }$$

$$\begin{equation}{\Omega}=\frac{{T}_{\text{m}}{\Delta}{S}_{\text{mix}}}{\vert {\Delta}{H}_{\text{mix}}\vert }\end{equation} \tag{ 1.5 }$$

$$\begin{equation}\text{VEC}=\sum\limits _{i=1}^{n}{x}_{i}{\text{VEC}}_{i},\end{equation} \tag{ 1.6 }$$

where x_i , r_i , T_mi and VEC_i are the atomic percent, the atomic radius, the melting point and the valence electron number of the ith component, respectively, and a and H_ij are the average atomic radius ($a=\sum _{i=1}^{n}{x}_{i}{r}_{i}$) and the mixing enthalpy of the corresponding binary alloy, respectively. In present work, the average mixing enthalpies of all alloys were calculated using Miedema model [21, 22], and the data of binary enthalpy was acquired from the work of Takeuchi et al [23].

The atomic size mismatch δ, describing the distortion degree of the solid solution alloy structure, has a significant effect on the phase structure and stability of the corresponding alloy. Empirically, δ should be less than 6.6% if a single-phase solid solution alloy is to be formed [24]. According to the second law of thermodynamics, the Gibbs free energy, which can reflect the phase stability of the alloy, is determined by its mixing entropy and mixing enthalpy at a specific temperature. Therefore, both the entropy and mixing enthalpy can be used to characterize the stability of the alloy structure. In terms of mixing entropy, as the temperature increases, it will play an increasingly important role in stabilizing the alloy. As for the mixing enthalpy, it is suggested that ΔH_mix should between −15 KJ mol⁻¹ and 5 KJ mol⁻¹ so as to form a single-phase solid solution HEA [24]. Besides, the parameter Ω, which is a combination of ΔS_mix, ΔH_mix, and T_m, was proposed to predict the phase structure of HEAs [20, 25, 26]. To form single-phase solid solution HEAs, Ω must be greater than 1. It is found that these empirical parameters are not always consistent in predicting the phase structure of HEAs. Therefore, it should be more reliable to apply these parameters comprehensively.

For single-phase solid solution HEAs, the crystal structure has a crucial influence on its mechanical properties. Generally speaking, alloys with bcc structure have higher strength, while alloys with fcc or hexagonal close-packaged (hcp) structure have better ductility. It is shown that the parameter VEC plays a vital role in predicting the crystal structure of single-phase solid solution alloys. According to experience, when VEC is less than 3.0, the structure of the HEA will be hcp; when VEC is greater than 3.0 and less than 7.5, the structure will be bcc; when VEC is greater than 7.5 and less than 7.8, the structure will be a combination of bcc and fcc; when VEC is greater than 7.8, the structure will be fcc.

Elastic properties are the basic properties of alloys, which cannot only characterize the stability of alloys, but also reflect the corresponding mechanical performance to a certain extent. In present work, various elastic properties of quaternary NbTiZr-based RHEAs were investigated, including elastic constants, Born stability criteria, elastic moduli, Poisson's ratio and Pugh ratio, elastic anisotropy indices, etc. The relevant calculation methods of these properties are as follows.

Elastic constants, which are the most basic elastic properties and can be used to derive many elastic properties, were obtained by first principles calculations in present work. Considering that the atomic distribution in the relatively small unit cell may lead to an anisotropic chemical environment, an averaging scheme [2] was adopted to obtain the final elastic constants:

$$\begin{equation}{C}_{11}=\frac{{C}_{11}+{C}_{22}+{C}_{33}}{3}\end{equation} \tag{ 2.1 }$$

$$\begin{equation}{C}_{12}=\frac{{C}_{12}+{C}_{23}+{C}_{13}}{3}\end{equation} \tag{ 2.2 }$$

$$\begin{equation}{C}_{44}=\frac{{C}_{44}+{C}_{55}+{C}_{66}}{3},\end{equation} \tag{ 2.3 }$$

where C_ij are single crystal elastic constants.

Based on the elastic constants, Born stability parameters [27], which include C₄₄ , C₁₁ − C₁₂ , C₁₂ + 2C₁₂ , and can be used to measure the mechanical stability of a crystal lattice, can be calculated accordingly. For a cubic crystalline lattice to be stable, the following inequalities must be satisfied: C₄₄ > 0, C₁₁ − C₁₂ > 0, C₁₂ + 2C₁₂ > 0, which can be called Born stability criteria. Generally speaking, the HEAs satisfying these criteria can be considered mechanically stable.

An elastic modulus is a quantity that measures an object or substance's resistance to being deformed elastically when a stress is applied to it. Specifying how stress and strain are to be measured, allows for many types of elastic moduli to be defined. The three primary ones are Young's modulus (E), shear modulus (G) and bulk modulus (B). Based on the obtained elastic constants, elastic moduli can be calculated by Voigt–Reuss–Hill method [5]. For cubic structure, the bulk and shear moduli in Voigt and Reuss approximation are as follows:

$$\begin{equation}B=\frac{{B}_{\text{V}}+{B}_{\text{R}}}{2}\end{equation} \tag{ 3.1 }$$

$$\begin{equation}G=\frac{{G}_{\text{V}}+{G}_{\text{R}}}{2},\end{equation} \tag{ 3.2 }$$

where B_V and G_V are the bulk and shear moduli calculated using Voigt model, while B_R, G_R are those calculated using Reuss model. They can be obtained using the following formulas:

$$\begin{equation}{B}_{\text{V}}=\frac{{C}_{11}+2{C}_{12}}{3}\end{equation} \tag{ 4.1 }$$

$$\begin{equation}{G}_{\text{V}}=\frac{{C}_{11}-{C}_{12}+3{C}_{44}}{5}\end{equation} \tag{ 4.2 }$$

$$\begin{equation}{B}_{\text{R}}=\frac{1}{\left({S}_{11}+{S}_{22}+{S}_{33}\right)+2({S}_{12}+{S}_{13}+{S}_{23})}\end{equation} \tag{ 4.3 }$$

$$\begin{equation}{G}_{\text{R}}=\frac{15}{4\left({S}_{11}+{S}_{22}+{S}_{33}\right)-4\left({S}_{12}+{S}_{13}+{S}_{23}\right)+3({S}_{44}+{S}_{55}+{S}_{66})},\end{equation} \tag{ 4.4 }$$

where S_ij are the elastic compliances, and the matrix composed by S_ij is the inverse of the elastic matrix composed by C_ij . Then the Young's modulus E is given by:

$$\begin{equation}E=\frac{9BG}{3B+G}.\end{equation} \tag{ 5 }$$

To estimate the ductility or brittleness of materials, Pugh introduces the B/G ratio (sometimes G/B ratio is used) [27]. Generally speaking, the larger the B/G ratio is or the smaller the G/B ratio is, the better the ductility of a material is. Quantitatively, if the B/G ratio is larger than 1.75 or G/B is less than 0.57, the ductile behavior can be predicted, and otherwise the material behaves in brittle manner. Besides, Poisson's ratio (ν) [27] can also be used to measure the ductility or brittleness of a material, which is given by:

$$\begin{equation}\nu =\frac{3B-2G}{2(3B+G)}.\end{equation} \tag{ 6 }$$

According to experience, the material with a Poisson's ratio greater than 0.31 are ductile and the larger the value is, the more ductile the system is. In fact, there is an inherent correlation between the Poisson's ratio and the Pugh ratio, which can be seen from their respective formulas clearly. In addition to these two ratios, Cauchy pressure (C₁₂ − C₄₄ ) [28] can also be applied to evaluate the ductility of a material. For a material with good ductility, positive Cauchy pressure can usually be observed, while negative Cauchy pressure usually indicates brittleness.

Elastic anisotropy of alloys is important since it correlates with the possibility to induce micro-cracks in materials. In present work, two different indices [29, 30] Zener anisotropy index A_Z and Chung–Buessem anisotropy index A_VR are used to assess the elastic anisotropy of RHEAs, for isotropic materials, the values of A_Z and A_VR will be close to 1 and 0, respectively and can be calculated using the following formulas:

$$\begin{equation}{A}_{\text{Z}}=\frac{2{C}_{44}}{{C}_{11}-{C}_{12}}\end{equation} \tag{ 6.1 }$$

$$\begin{equation}{A}_{\text{VR}}=\frac{{G}_{\text{V}}-{G}_{\text{R}}}{{G}_{\text{V}}+{G}_{\text{R}}}.\end{equation} \tag{ 6.2 }$$

As mentioned above, first principles method was employed to calculate the elastic constants of each quaternary NbTiZr-based RHEA. To compare with NbTiZrX alloys, the elastic constants of ternary NbTiZr alloy were also calculated. In the process of calculation, a reasonable crystal model is crucial to the reliability of the final results. In present work, the SQS method based on the cluster expansion formalism [31] was chosen to construct the crystal structures of both NbTiZr and NbTiZrX alloys. According to the literatures and the predicted results by empirical parameters (see below), if the NbTiZrX alloy can form a single-phase solid solution, it tends to be a bcc structure. Hence, bcc was used when constructing the crystal structure of NbTiZrX alloy. Besides, supercells were constructed for both NbTiZr alloy and NbTiZrX alloy to simulate the bulk phase. Specifically, a 3 × 3 × 3 supercell with 54 atoms was constructed for NbTiZr alloy, while a 4 × 2 × 2 supercell with 32 atoms was constructed for NbTiZrX alloy. Besides, the lattice constant for each quaternary NbTiZrX alloy was set according to Vegard's Law [32], which can be expressed as:

$$\begin{equation}{a}_{\text{RHEA}}=\frac{1}{3+x}\left({a}_{\text{Nb}}+{a}_{\text{Ti}}+{a}_{\text{Zr}}\right)+\frac{x}{3+x}{a}_{\text{X}},\end{equation} \tag{ 7 }$$

where x is the molar ratio of undetermined element X, and a_Nb, a_Ti, a_Zr, and a_X are the equilibrium lattice constants of the elements Nb, Ti, Zr and X in bcc structure, respectively. Figure 2 shows the unrelaxed supercell models of the ternary NbTiZr alloy and the quaternary NbTiZrX alloy constructed using the SQS method.

Zoom In Zoom Out Reset image size

In terms of computational detail, Vienna ab initio simulation package (VASP) [33] was adopted to perform all the first principles calculations. Exchange and correlation were treated at the Perdew–Burke–Erzernhof functional level [34]. The plane ware energy cutoff was set 400 eV. Brillouin-zone sampling was performed using the Monkhorst–Pack scheme [35] with a 3 × 5 × 5 k-point mesh. Before calculating the elastic constants, the crystal structures were relaxed to minimize the atomic forces until the force on each atom is less than 10⁻² eV Å⁻¹.

As mentioned above, the first step in present work is to determine the elemental composition of the NbTiZrX alloy, that is, to determine the fourth principal alloying element X in the alloy in addition to Nb, Ti and Zr. In order to reduce the subsequent calculations, some basic rules were applied to screen all the 118 elements in the periodic table of chemical elements.

Firstly, the three elements (Nb, Ti, Zr) that already exist in the NbTiZrX alloy can be excluded, leaving 115 elements. Although non-metallic elements can appear in HEAs, they are mainly in the form of minor elements rather than major elements. Therefore, the second rule is that the element X in the NbTiZrX alloy must be a metal element. According to this rule, the range of candidate elements is reduced to 91 metal elements. The third rule is that the element X must be non-radioactive, which narrows the candidate range to 59 non-radioactive metal elements. The fourth rule is that the boiling point of the element X mush be greater than the maximum melting point of the three elements Nb, Ti and Zr, so as to avoid the element X from burning out during the alloy forming process. The melting point of Nb is the largest among the three elements of Nb, Ti and Zr, reaching 2750 K, so the boiling point of the fourth element X in the NbTiZrX alloy should be greater than this value. According to this rule, 25 elements can be excluded, leaving 34 elements. The fifth and final rule is based on the melting point and density of the element. For RHEAs, high melting point and low density are two requirements that need to be considered comprehensively to meet the application in aerospace field. Therefore, the elements in RHEAs should meet the two requirements as much as possible. For this reason, the ratio of melting point to density, or called specific melting point, was created as a criterion to screen the remaining elements. Specifically, the specific melting point of Hf was employed as the threshold value, and the elements whose specific melting point is greater than the threshold value were retained. Based on this rule, 13 more elements were excluded, leaving 21 elements, which are Al, Co, Cr, Dy, Er, Fe, Gd, Ge, Hf, Ho, La, Lu, Mo, Ni, Ru, Sc, Ta, Tb, V, W, Y. Table 1 lists all the elements that were screened out and the corresponding melting points, densities and the specific melting points. By observing these elements, it can be found that the most frequently researched elements in quaternary NbTiZrX RHEAs are all included, such as Cr, Hf, Mo, Ta, V and W, which reflects the rationality of the screening rules used in present work. These screened elements would be regarded as element X and combined with Nb, Ti, Zr to form quaternary NbTiZrX alloy systems, and subsequent investigations would be carried out on these alloy systems.

To predict the phase structure and stability of all involved NbTiZrX RHEAs, the empirical parameters introduced in section 2.2 were calculated for each alloy system, which includes the atomic size mismatch (δ), the mixing entropy (ΔS_mix), the average mixing enthalpy (ΔH_mix), the average melting point (T_m), the parameter Ω and the average VEC. Specifically, four alloys with different contents of element X, including NbTiZrX_0.5, NbTiZrX_1.0, NbTiZrX_1.5, NbTiZrX_2.0, were considered for each alloy system. Besides, the empirical parameters of the ternary NbTiZr alloy were also calculated as a comparison. According to the calculated value of mixing enthalpy ΔH_mix, all 21 NbTiZrX alloy systems screened out in the previous step can be divided into two categories, one of which meets the requirements of forming a single-phase solid solution, including Al, Cr, Hf, Mo, Ta, V, W, and the other does not meet the requirements, including Co, Dy, Er, Fe, Gd, Ge, Ho, La, Lu, Ni, Ru, Sc, Tb, Y. The calculated empirical parameters of the former category are listed in table 2.

For Al, as can be seen from table 2, with the increase of Al in NbTiZrAl_x , the atomic size mismatch δ and the average mixing enthalpy ΔH_mix both decrease. The decrease of δ will make it easier for the alloy to form a single-phase solid solution, while the decrease of ΔH_mix will make it easier to form an intermetallic compound. Combined with the parameter Ω, it can be determined that when the molar ratio of Al is relatively low, the NbTiZrAl_x alloy will form a single-phase solid solution, and when the molar ratio of Al is greater than 1.5, the value of ΔH_mix no longer satisfies the requirement of forming a single-phase solid solution for HEAs. This result is consistent with the conclusions reported by Stepanov et al [9, 36] and Senkov et al [37, 38], that is, when the content of Al is very low, it serves as a stabilizer for the bcc/B2 phase, and when the content of Al gets higher, it tends to form short-range ordered structures with Hf and Zr in the alloy, which promotes the formation of hcp phase.

For Cr, referring to the calculated values of ΔH_mix and Ω, it can be determined that the NbTiZrCr_x should be able to form a single-phase solid solution. However, the values of δ have exceeded the corresponding threshold of 6.6%. Therefore, the NbTiZrCr_x alloys are not necessarily able to form a single-phase solid solution. Moreover, with the increase of Cr, the possibility of forming a single-phase solid solution decreases. Besides, since the ratio of the atomic radius of Cr to that of Zr is between 1.1 and 1.6, it is likely to form Laves phases by combining with Zr in NbTiZrCr_x alloy. Miracle and Senkov [39] used Al to replace Cr when studying RHEAs so as to avoid the formation of brittle, topologically tightly packed Laves phases, thereby improving the ductility at room temperature, which is also consistent with the conclusions on Cr in present study.

For Hf, Ta and W, the corresponding NbTiZrX alloys are predicted to be single-phase solid solutions no matter which empirical parameter is used. Furthermore, it can be observed that, for NbTiZrHf_x and NbTiZrTa_x alloy systems, the calculated empirical parameters differ greatly with their respective thresholds, and the values of the parameter Ω increase with the increase of Hf and Ta, indicating that a single-phase solid solution structure can still be maintained for alloys with higher content of Hf or Ta. As for NbTiZrW_x alloy systems, the values of the parameter Ω are only slightly larger than the corresponding threshold 1.0, and decrease with the increase of W, indicating that the content of W is approaching the limit value of being able to form a single-phase solid solution.

For Mo and V, all involved NbTiZrMo_x and NbTiZrV_x alloys are predicted to be single-phase solid solutions according to the calculated empirical parameters listed in table 2. However, when the molar ratio of V increases from 1.0 to 2.0, the value of the parameter Ω for corresponding alloy decreases dramatically from 103.397 to 19.289, which suggests that when the molar ratio of V is larger than 1.0, the NbTiZrV_x alloy tends to become unstable with the increase of V. For NbTiZrMo_x , similar tendency can be observed, indicating that alloys with higher content of Mo are likely to be unstable. Such conclusions agree well with those reported by Zhang et al [40] and Hui et al [41], that is, the addition of V can promote the phase separation of bcc, and the addition of Mo can enhance the tendency of phase decomposition. Nevertheless, for NbTiZrX alloys with a molar ratio of Mo or V less than 2.0, their stability of single-phase solid solution is very good.

As for the remaining 14 elements, they can be further divided into two sub-categories, one is transition metals (including post-transition metals) and the other is lanthanide metals. The former sub-category includes 6 transition metals (Co, Fe, Ni, Ru, Sc, Y) and one post-transition metal (Ge), and the latter sub-category includes 7 lanthanide metals (Dy, Er, Gd, Ho, La, Lu, Tb). For NbTiZrX alloy systems where the X element belongs to these two sub-categories, the calculated empirical parameters are listed in tables S1 and S2 (https://stacks.iop.org/MSMS/29/075002/mmedia) of the supplementary information file, respectively. For transition metals, it can be seen from table S1, although the atomic size mismatch δ of NbTiZrX alloys for some X elements, such as Sc and Ge, can meet the requirements of forming a single-phase solid solution, the average mixing enthalpy ΔH_mix of all NbTiZrX alloys do not meet the requirements. As for lanthanide elements, as can be seen from table S2, none of their NbTiZrX alloys can form a single-phase solid solution whether it is based on the atomic size mismatch δ, or the average mixing enthalpy ΔH_mix, or the parameter Ω.

In summary, according to the calculated empirical parameters, among the 21 elements screened out in the previous step, only 7 elements, namely Al, Cr, Hf, Mo, Ta, V, W, can form a quaternary NbTiZr-based RHEA with a single-phase solid solution structure. Furthermore, according to the value of VEC, it can be easily determined that all the NbTiZr-based alloys listed in table 2 are predicted to be bcc structure.

By observing these remaining elements screened out by empirical parameters, it can be found that all elements except Al are around the three elements Nb, Ti, and Zr in the periodic table of chemical elements, that is, transition metal elements from groups IVB to VIB. Generally, the closeness of elements in the periodic table means that some of their basic properties, such as atomic size, electronegativity, etc, are relatively close, which can ensure their chemical compatibility with each other [42]. In addition, from the perspective of crystal structure, all the remaining elements except Al (fcc) and Hf (hcp) have a bcc structure, which are the same as the final crystal structure of the single-phase solid solution of NbTiZr-based RHEA. Conversely, the crystal structures of those elements excluded during the screening process are mostly hcp, except Ge (diamond), Ni (fcc) and Fe (bcc). Therefore, a simple, slightly rough rule to determine whether an NbTiZr-based RHEA can form a single-phase solid solution can be concluded, that is, the added element should be around the basic elements of the alloy in the periodic table and the corresponding crystal structure should be the same as that of the final single-phase solid solution.

Based on the assumption of ideal mixing, MacDonald et al [43] proposed that the atomic composition with equimolar ratio can maximize the configurational entropy of corresponding HEA. Therefore, the elastic properties of equimolar NbTiZrX RHEAs were investigated, where the X element is among the 7 elements screened out by the empirical parameters in the previous step.

The elastic constants for both NbTiZr and NbTiZrX alloys (X = Al, Cr, Hf, Mo, Ta, V, W) with bcc structure were calculated using first principles method. Table 3 lists the calculated results of the elastic constants C₁₁ , C₁₂ and C₄₄ for all involved alloys and more detailed results of C₁₁ ∼ C₆₆ can be found in table S3 of the supplementary information file. To facilitate the comparison of the calculated elastic constants for all involved alloys, figure 3 is drawn according to the relevant data. In the figure, the label none denotes the ternary NbTiZr alloy (same below). It can be seen that C₁₁ of NbTiZrW alloy is the largest, followed by NbTiZrMo alloy, and that of NbTiZrAl alloy is the lowest. C₁₂ of all alloys are almost the same. C₄₄ of NbTiZrAl alloy is the largest, followed by NbTiZrHf alloy and NbTiZrTa alloy, and the rest are almost the same.

Zoom In Zoom Out Reset image size

Based on the obtained elastic constants of these NbTiZrX alloys, other elastic properties, which includes Born stability parameters (C₄₄ , C₁₁ − C₁₂ , C₁₂ + 2C₁₂ ), elastic moduli (shear modulus G, bulk modulus B and Young's modulus E), Pugh ratio (B/G), Poisson's ratio (ν), Cauchy pressure (C₁₂ − C₄₄ ), elastic anisotropy indices (Zener index A_Z and Chung–Buessem index A_VR), can be calculated by referring to the formulas introduced in section 2.3. The corresponding results are listed in table 3. To show these results clearly, figure 4 with four sub-figures are plotted.

Zoom In Zoom Out Reset image size

Figure 4(a) shows the calculated Born stability parameters for both the ternary NbTiZr alloy and the 7 quaternary NbTiZrX alloys. It can be seen easily that all these values are larger than 0, which means that all these involved alloys are predicted to be mechanically stable according to the Born stability criteria. What should be mentioned, although the C₁₁ − C₁₂ values of NbTiZrAl alloy and NbTiZrHf alloy are greater than 0, they are closer to 0 compared with other alloys, indicating that the mechanical stability of these two alloys is not as good as other alloys. The main reason for such result is that the crystal structure of Al (fcc) and Hf (hcp) is different from the corresponding NbTiZr-based alloy (bcc), which makes the structure of the alloy more distorted and less stable.

Figure 4(b) displays the calculated elastic moduli, including shear modulus, bulk modulus and Young's modulus, for ternary NbTiZr alloy and quaternary NbTiZrX alloys. Combined with table 3, it can be seen that each modulus of any quaternary NbTiZrX alloy is greater than that of the ternary NbTiZr alloy, except for the bulk modulus of NbTiZrHf alloy, which indicates that the addition of the fourth element to NbTiZr alloy can strengthen the alloy. Secondly, although the three moduli of different alloys can be quite different, the order of each modulus among these involved alloys is almost the same, that is, the alloy with a larger shear modulus will also have a larger bulk modulus or Young's modulus, which reflects the intrinsic relationship among these three moduli to some extent. Quantitatively, the three moduli of NbTiZrW alloy are the largest among all alloys, which indicates that the W element has the best strengthening effect on the ternary NbTiZr alloy. NbTiZrMo alloy and NbTiZrTa alloy also have good performance in these three moduli. In contrast, the moduli of NbTiZrAl and NbTiZrHf alloys are relatively low.

From the perspective of Cauchy pressure, it can be seen from table 3 that all the calculated results of NbTiZr and NbTiZrX alloys are positive, which demonstrates that all alloys involved are ductile. Nevertheless, the ductility of these alloys remains different, which can be observed from the calculated values of Poisson's ratio and Pugh ratio, as drawn in figure 4(c). From the figure, it can be seen that the Poisson's ratio and Pugh ratio of quaternary NbTiZrX alloys are smaller than those of ternary NbTiZr alloy, which indicates that the addition of the fourth alloying element to the NbTiZr alloy will reduce the ductility of the resulting alloy. This is exactly the opposite of the influence of alloying elements on the strength of the alloy mentioned above, which demonstrates again that strength and ductility are a pair of mutually exclusive properties to some extent. What can be found next is that the Poisson's ratio and the Pugh ratio are consistent in changing trends for all alloys. Specifically, the NbTiZrX alloys with Ta, Mo and W as the X element has relatively small Poisson's ratio and Pugh ratio, which means that the ductility of these three alloys is not as good as other alloys. On the contrary, NbTiZrAl alloy has relatively good ductility, followed by NbTiZrV and NbTiZrHf.

The calculated results of anisotropic indices A_Z and A_VR for NbTiZr and NbTiZrX alloys are drawn in figure 4(d). It can be seen clearly that the A_Z value and A_VR value of most alloys are close to 1 and 0, respectively, except those of NbTiZrAl alloy. Therefore, it can be inferred that most NbTiZrX alloys here are elastically isotropic. As for NbTiZrAl alloy, from table 3, the A_Z value and A_VR value of NbTiZrAl alloy have reached 17.749 and 0.655, exceeding their respective thresholds, which indicates that the alloy has obvious anisotropy. In addition, the NbTiZrHf alloy also shows a certain degree of anisotropy. By analyzing the crystal structure of these elements and the corresponding NbTiZr-based RHEA, it can be inferred that the anisotropy of the alloy is closely related to the crystal structure of the added element. If the crystal structure of the added element is consistent with that of the final alloy, the anisotropy of the alloy will be weaker, and vice versa.

From the introduction of the above sub-sections, it can be seen that the quaternary NbTiZrX alloy, with the X element in Al, Cr, Hf, Mo, Ta, V, W, has relatively outstanding performance compared with other alloys. For a clear comparison, the various properties of these alloys in terms of relative degrees, such as good, average, fair, etc., are listed in table 4.

As a whole, the NbTiZrX alloys listed in table 4 can form a stable single-phase solid solution structure, although the stability can be different for different alloys. From the perspective of strength and ductility, no alloy can have the best performance in both properties. That's to say, the alloys with higher strength are generally less ductile than other alloys. Specifically, the alloys containing W, Ta, and Mo have higher strength, while the alloys containing Al, Hf, and V have better ductility. It should be noted that although the strength of some alloys is not as good as other alloys, it does not mean that these alloys are not competent for the environments where high strength is required. On the one hand, the properties shown here are intrinsic properties, while the actual properties are determined by many factors, such as the fabrication process. On the other hand, the properties listed in table 3 are relative, not absolute. Therefore, NbTiZrHf alloy may also be able to meet the requirements of high strength. In addition to the above, it is necessary to consider the density or melting point sometimes. For example, for those occasions that require high temperature but not too much density, NbTiZrW alloy should be a good choice, and if there is a requirement for low density, NbTiZrMo or NbTiZrV should be a better choice.