First-principles calculation of the effect of Ti content on the structure and properties of TiVNbMo refractory high-entropy alloy

Due to the high mixing entropy, high entropy alloys usually form simple disordered solid solutions, such as BCC, FCC, and HCP. In order to predict the phase structure of HEAs, it is necessary to focus on thermodynamic parameters such as mixing entropy ΔS_mix, mixing enthalpy ΔH_mix, atomic size difference δ, electronegativity difference Δχ and valence electron concentration (VEC), and they all follow the classical Hume Rothery rule [20]. Their values can be calculated by a series of equation (1)

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\rm{\Delta }}{S}_{mix}=-R\displaystyle \sum {c}_{i}\,\mathrm{ln}\,{c}_{i}\\ {\rm{\Delta }}{H}_{mix}=\displaystyle \sum _{i=1,i\ne j}^{n}4{\rm{\Delta }}{H}_{ij}^{mix}{c}_{i}{c}_{j}\\ \delta =\sqrt{\displaystyle \sum _{i=1}^{n}{c}_{i}{\left(1-\displaystyle \frac{{r}_{i}}{\bar{r}}\right)}^{2}};\bar{r}=\displaystyle \sum _{i=1}^{n}{c}_{i}{r}_{i}\\ {\rm{\Delta }}\chi =\sqrt{\displaystyle \sum _{i=1}^{n}{c}_{i}{\left({\chi }_{i}-\bar{\chi }\right)}^{2}};\bar{\chi }=\displaystyle \sum _{i=1}^{n}{c}_{i}{\chi }_{i}\\ VEC=\displaystyle \sum _{i=1}^{n}{c}_{i}{\left(VEC\right)}_{i}\end{array}\right.\end{eqnarray} \tag{ 1 }$$

Where ΔS_mix is the mixing entropy of the HEA, c_i is the atomic ratio of the i-th alloy component, and R is the molar gas constant. ΔH_mix is the mixing enthalpy of the HEA, ${\rm{\Delta }}{H}_{ij}^{mix}$ is the enthalpy of mixing of binary liquid alloys, c_i and c_j are the atomic ratios of the i-th and j-th alloy components, respectively. δ is the atomic radius difference of the HEA, r_i is the atomic radius of the i-th alloy component, $\bar{r}$ is the average atomic radius. ${\rm{\Delta }}\chi$ is the electronegativity difference of the HEA, $\bar{\chi }$ is the average electronegativity difference of the constituent elements calculated by the mixing rule, ${\chi }_{i}$ is the Pauling electronegativity of the i-th element; VEC is the valence electron concentration of HEA, ${(VEC)}_{i}$ is the valence electron concentration of the i-th element. At the same time, some extended Hume-Rothery rule parameters will also affect the phase formation results of HEA, such as: Ω, Λ and γ [20], calculated according to equation (2)

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Omega }}=\displaystyle \frac{{T}_{m}{\rm{\Delta }}{S}_{mix}}{\left|{\rm{\Delta }}{H}_{mix}\right|};{T}_{m}=\displaystyle \sum _{i=1}^{n}{c}_{i}{\left({T}_{m}\right)}_{i}\\ \gamma =\displaystyle \frac{{\omega }_{S}}{{\omega }_{L}};{\omega }_{S}=1-\sqrt{\displaystyle \frac{{\left({r}_{S}+\bar{r}\right)}^{2}-{\bar{r}}^{2}}{{\left({r}_{S}+\bar{r}\right)}^{2}}};\\ {\omega }_{L}=1-\sqrt{\displaystyle \frac{{\left({r}_{L}+\bar{r}\right)}^{2}-{\bar{r}}^{2}}{{\left({r}_{L}+\bar{r}\right)}^{2}}}\end{array}\end{eqnarray} \tag{ 2 }$$

Where ${({T}_{m})}_{i}$ is the melting temperature of the i-th element, and Ω takes into account both effect of the entropy and the enthalpy. γ describes the instability of atomic packing and is defined as the ratio of the smallest and largest solid angles, ${\omega }_{S}$ is the solid angle of the smallest atom, ${\omega }_{L}$ is the solid angle of the largest atom, and ${r}_{S}$ and ${r}_{L}$ are the radii of the smallest and largest atoms. When Ω ≥ 1.1, δ ≤ 6.6%, and γ < 1.175, The parameter values listed in table 1 fully meet the limitation of obtaining a single solid solution structure [23]. TiVNbMo as a potential aerospace application material, low density is a significant advantage, the theoretical alloy density ρ_theor determined by equation (2) [1]

$$\begin{eqnarray}&&{\rho }_{{\rm{theor}}}=\displaystyle \frac{\displaystyle {\sum }_{i=1}^{n}{c}_{i}{A}_{i}}{\displaystyle {\sum }_{i=1}^{n}\tfrac{{c}_{i}{A}_{i}}{{\rho }_{i}}}\end{eqnarray} \tag{ 3 }$$

Where ${A}_{i}$ and ${\rho }_{i}$ are the atomic weight and density of element i.

In this work, the energy and elasticity of Ti_xVNbMo (x = 1.00, 1.25, 1.75, 2.00) RHEAs are calculated by the VCA model. All the first-principles calculations are based on the Cambridge sequential total energy package (CASTEP) [24]. Exchange-Correlation functional was selected the Perdew–Burke–Ernzerhof (PBE) in generalized gradient approximation (GGA), choose the TPSD algorithm when geometrically optimizing the structure, pseudopotential selects OTFG norm serving, custom cutoff energy was set to 950 eV, and a k-points was setting to 20 × 20 × 20. In this MonkhorstPack scheme, the Brillouin zone is sampled to ensure that the total energy error is less than 1 × 10⁻⁵ eV atom⁻¹. Besides, in the subsequent geometric optimization, all the forces on a single atom converged to less than 0.03 eV Å⁻¹, and the total stress tensor was reduced to the order of 0.05 GPa, max placement of an atom is less than 0.001 Å.

2.2.1. Lattice constant

A more efficient TPSD structure optimization method is selected for lattice constant calculation, and the optimization results match well with Birch–Murnaghan equation of state (EOS) [25, 26]. The EOS method calculates the energy of crystals of different volumes (here the shape of the crystal lattice is retained as the initial structure, namely BCC), and then fits the E-V curve through the third-order Birch-Murnaghan equation (4) of the EOS method, as shown in figure 1. The volume with the minimum energy corresponds to the equilibrium volume. The initial VCA-BCC crystal structure is based on the original structure of Nb, and its lattice parameters a = 3.3006 Å$\unicode{x002F3}$

$$\begin{eqnarray}&&E\left(V\right)={E}_{0}+\displaystyle \frac{9{V}_{0}{B}_{0}}{16}{\left\{\left[{\left(\displaystyle \frac{{V}_{0}}{V}\right)}^{2/3}-1\right]\right.}^{3}\left.{B}_{0}^{{\prime} }+{\left[{\left(\displaystyle \frac{{V}_{0}}{V}\right)}^{2/3}-1\right]}^{2}\left[6-4{\left(\displaystyle \frac{{V}_{0}}{V}\right)}^{2/3}\right]\right\}\end{eqnarray} \tag{ 4 }$$

Where V₀ and E₀ are the equilibrium volume and the corresponding energy, respectively. V is the volume of change, B₀ is the bulk modulus, and ${B}_{0}^{{\prime} }$ is a constant. The deviation between the EOS fitting value of lattice constant and that of the vacuum arc melting experiment is only 0.13%–0.91%, their values are listed in table 2, which indicates that the parameter setting of the calculation method is reasonable.

Zoom In Zoom Out Reset image size

Table 2. The lattice constants of Ti_xVNbMo (x = 1.00, 1.25, 1.75, 2.00) refractory high-entropy alloy.

			VCA Calculations (this work)
Lattice constants	TiVNbMo [11]	TiVNbMo [27]	Ti1.00	Ti1.25	Ti1.50	Ti1.75	Ti2.00
a Å−1	3.211	3.186	3.182	3.187	3.193	3.198	3.202

2.2.2. Elastic constant

For cubic crystals, there are three independent elastic constants: C_11, C_12, and C_44, which should meet the requirements of Born-Huang mechanical stability criterion [28], as shown in equation (5)

$$\begin{eqnarray}&&{C}_{11}-{C}_{12}\gt 0,{C}_{11}+2{C}_{12}\gt 0,{C}_{44}\gt 0\end{eqnarray} \tag{ 5 }$$

When the Cauchy pressure C_p is greater than 0, it indicates that the alloy is completely metal bond, which is defined as equation (6) [28]

$$\begin{eqnarray}&&{C}_{p}={C}_{12}-{C}_{44}\end{eqnarray} \tag{ 6 }$$

The calculation method of modulus for cubic crystal is shown in equation (7) [16]

$$\begin{eqnarray}&&\left\{\begin{array}{l}{B}_{V}={B}_{R}=\left({C}_{11}+2{C}_{12}\right)/3,\\ {G}_{V}=\left({C}_{11}-{C}_{12}+3{C}_{44}\right)/5,{G}_{R}=5\left({C}_{11}-{C}_{12}\right){C}_{44}/\left[4{C}_{44}+3\left({C}_{11}-{C}_{12}\right)\right]\end{array}\right.\end{eqnarray} \tag{ 7 }$$

According to the Voigt-Reuss-Hill modulus principle [16], the calculation equations for the modulus B, G, E and Poisson's ratio σ of cubic crystals can be obtained (8)

$$\begin{eqnarray}&&\left\{\begin{array}{l}B=\left({B}_{V}+{B}_{R}\right)/2;G=\left({G}_{V}+{G}_{R}\right)/2\\ E=9BG/\left(3B+G\right);\sigma =\left(3B-2G\right)/\left(6B+2G\right)\end{array}\right.\end{eqnarray} \tag{ 8 }$$

The general anisotropy index A^U is calculated by equation (9) [29]

$$\begin{eqnarray}&&{A}^{U}=5\displaystyle \frac{{G}_{V}}{{G}_{R}}+\displaystyle \frac{{B}_{V}}{{B}_{R}}-6\geqslant 0\end{eqnarray} \tag{ 9 }$$

The hardness value of the alloy is calculated according to the equation (10) proposed by Chen [30]

$$\begin{eqnarray}&&{H}_{V}=0.151G\end{eqnarray} \tag{ 10 }$$

2.2.3. Thermodynamic properties

Thermodynamic properties are very important, and they play an important role in understanding the thermal response of solids. In order to study the thermodynamic properties of Ti_xVNbMo refractory high-entropy alloys, the quasi-harmonic Debye model using the Gibbs1.0 program is applied here [31]. This tool can be used to calculate the thermodynamic properties of the phases in the alloy, such as pressure-volume-temperature relationship, isothermal volume Modulus, volume thermal expansion coefficient, heat capacity, entropy, and Debye temperature. Since the detailed introduction of this method and its application has been described in other literature reports [32, 33], we only briefly summarize the method here. The unbalanced Gibbs function G*(V; P, T) has the form

$$\begin{eqnarray}&&G* \left(V;P,T\right)=E\left(V\right)+PV+{A}_{{\rm{vib}}}\left({\rm{\Theta }}/T\right)\end{eqnarray} \tag{ 11 }$$

Where E(V) is the total energy of each unit cell of Ti_xVNbMo at 0 K, which can be calculated by EOS, PV corresponds to a constant hydrostatic pressure condition, A_vib is the Helmholtz free energy of vibration, and Θ is Debye temperature. Using the complete phonon density of states (PDOS) requires many calculations. The GIBBS method uses the Debye model and considers the quasi-harmonic approximation. The A_vib Helmholtz free energy is expressed as Debye temperature

$$\begin{eqnarray}&&{A}_{{\rm{vib}}}\left({\rm{\Theta }},T\right)=nkT\left[\displaystyle \frac{9}{8}\displaystyle \frac{{\rm{\Theta }}}{T}+3\,\mathrm{ln}\left(1-{e}^{-{\rm{\Theta }}/T}\right)-D\left({\rm{\Theta }}/T\right)\right]\end{eqnarray} \tag{ 12 }$$

Where D(Θ/T ) is the Debye integral, defined as equation (13), n is the sum of the number of atoms in the chemical formula, and Θ is the Debye temperature, defined as equation (14).

$$\begin{eqnarray}&&D\left({\rm{\Theta }}/T\right)=\displaystyle \frac{3}{{y}^{3}}\displaystyle {\int }_{0}^{y}\displaystyle \frac{{x}^{3}}{{e}^{x}-1}dx\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\rm{\Theta }}=\displaystyle \frac{\hslash }{k}{\left[6{\pi }^{2}{V}^{1/2}n\right]}^{1/3}f\left(\sigma \right)\sqrt{\displaystyle \frac{{B}_{s}}{M}}\end{eqnarray} \tag{ 14 }$$

Where M is the molecular weight of each molecular formula, B_s is the adiabatic bulk modulus, which can be approximated as equation (15), σ is the Poisson's ratio calculated by the elastic constant, f(σ) is given by equation (16)

$$\begin{eqnarray}&&{B}_{s}\approx B\left(V\right)=V\left(\displaystyle \frac{{d}^{2}E\left(V\right)}{d{V}^{2}}\right)\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&f\left(\sigma \right)={\left\{\left.3{\left[2{\left(\displaystyle \frac{2}{3}\displaystyle \frac{1+\sigma }{1-2\sigma }\right)}^{3/2}+{\left(\displaystyle \frac{1}{3}\displaystyle \frac{1+\sigma }{1-\sigma }\right)}^{3/2}\right]}^{-1}\right\}\right.}^{1/3}\end{eqnarray} \tag{ 16 }$$

When the Gibbs free energy is determined, other thermodynamic properties can be calculated. For example, vibration internal energy U_v equation (17), heat capacity C_v equation (18), thermal expansion is described by the Debye-Grüneisen model [26]. In this model, the volume expansion coefficient can be calculated by equation (19)

$$\begin{eqnarray}&&{U}_{{\rm{v}}}=nkT\left[\displaystyle \frac{9}{8}\displaystyle \frac{{\rm{\Theta }}}{T}+3D\left({\rm{\Theta }}/T\right)\right]\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{C}_{{\rm{v}}}=3nk\left[4D\left({\rm{\Theta }}/T\right)-\displaystyle \frac{3{\rm{\Theta }}/T}{{e}^{{\rm{\Theta }}/T}-1}\right]\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\alpha }_{{\rm{v}}}=\displaystyle \frac{\gamma {C}_{{\rm{v}}}}{{B}_{T}V}\end{eqnarray} \tag{ 19 }$$

Where γ is the Grüneisen parameter, defined as equation (20), and B_T is the bulk modulus at temperature T, which can be calculated using equation (21)

$$\begin{eqnarray}&&\gamma =-\displaystyle \frac{{\rm{d}}\,\mathrm{ln}\,{\rm{\Theta }}\left(V\right)}{{\rm{d}}\,\mathrm{ln}\,V}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{B}_{T}=V\left(\displaystyle \frac{{\partial }^{2}G\left(V;P,T\right)}{\partial {V}^{2}}\right)\end{eqnarray} \tag{ 21 }$$

The elastic stiffness, polycrystalline elasticity, and hardness of Ti_xVNbMo were calculated with the VCA model. The influence of Ti content on the lattice constant a and theoretical density ρ_theor of the alloys are shown in figure 2. Because Ti has the largest atomic radius and the smallest density, the lattice constant a of Ti_xVNbMo alloy gradually increases from 3.182 Å to 3.202 Å, and the theoretical density ρ_theor gradually decreases from 4.28 g cm⁻³ to 3.85 g cm⁻³.

Zoom In Zoom Out Reset image size

The elastic constants determine the elastic properties of materials, which can directly reflect the mechanical properties of materials. Figures 3(a) and (b) show that with the change of Ti content, the elastic constant C_ij of Ti_xVNbMo alloy fully meets the Born-Huang mechanical stability criterion, and they are both stable BCC structure and all of them are metal bond binding (C_p > 0). The elastic properties of materials are mainly reflected by bulk modulus B, shear modulus G, Young's modulus E, and Poisson's ratio σ. The bulk modulus reflects the ability of the alloy to resist volume changes and is a manifestation of the average bonding strength between atoms in the material. The shear modulus expresses the ability of the material to resist shear deformation. Young's modulus reflects the ability of the alloy to resist elastic deformation and is a manifestation of alloy rigidity. Poisson's ratio σ is a physical quantity used to quantify the stability of crystal materials against shear deformation. Higher Poisson's ratio σ, better plasticity. B/G the ratio of shear modulus to bulk modulus can also reflect the brittleness and toughness of metal materials. The higher the ratio, the better the toughness of metal [28]. Figure 3(c) is the first-principle calculation result of the elastic properties of Ti_xVNbMo RHEA. With the increase of Ti content, B, G, and E all gradually decrease. That is, as the Ti content increases, the average bonding ability between alloy atoms and the ability to resist deformation simultaneously become weaker, the rigidity of the material decreases, and the material are more inclined to deform elastically. Figure 3(d) shows the effect of Ti content on Poisson's ratio σ and B/G. They all increase gradually with the increase of Ti content, indicating that the plasticity and toughness of Ti_xVNbMo RHEAs are positively correlated with the Ti content, and this trend is more significant when x ≥ 1.50.

Zoom In Zoom Out Reset image size

Figure 4 is a dotted line graph of the predicted values of the general anisotropy index and the hardness of Ti_xVNbMo RHEA with different Ti content calculated from equations (9) and (10). The smaller the value of A^U, the closer to 0, the more isotropy of the alloy. When A^U is equal to zero, the alloy is isotropic. It can be seen from figure 4 that as the Ti content increases, the alloy gradually tends to be isotropic, and the hardness value gradually decreases.

Zoom In Zoom Out Reset image size

In this work, in order to observe the mechanical anisotropy of the material more clearly, figure 5 shows the three-dimensional (3D) surface of Young's modulus of Ti_xVNbMo RHEA BCC phases with different Ti content in spherical coordinates. The calculation process is consistent with that of [34]. For cubic structures, the 3D representation of Young's modulus is given by equation (22) [35]. As shown in figures 5(a)–(e), with the Ti content increases, Young's modulus of Ti_xVNbMo RHEA gradually tends to form a completely spherical shape, indicating that they gradually tend to be isotropic. Projection of Young's modulus of Ti_xVNbMo alloy on crystal plane (001) and crystal plane (110) are shown in figures 6(a) and (b), respectively. The Young's modulus anisotropy of the (001) crystal plane is weaker than that of the (110) plane. When x = 1.00, that is, TiVNbMo alloy has the strongest Young's modulus anisotropy, because its shape deviates from a perfect sphere.

$$\begin{eqnarray}&&\displaystyle \frac{1}{E}={l}_{1}^{4}{S}_{44}+{l}_{2}^{4}{S}_{22}+{l}_{3}^{4}{S}_{33}+2{l}_{1}^{2}{l}_{2}^{2}{S}_{12}+2{l}_{1}^{2}{l}_{3}^{2}{S}_{13}+2{l}_{2}^{2}{l}_{3}^{2}{S}_{23}+{l}_{1}^{2}{l}_{2}^{2}{S}_{66}+{l}_{1}^{2}{l}_{3}^{2}{S}_{55}+{l}_{2}^{2}{l}_{3}^{2}{S}_{44}\end{eqnarray} \tag{ 22 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the previous calculation, the total energy and equilibrium volume of the system in the 0 K state conforms to the EOS equation of state. Taking the E–V data in figure 1 as the input data of the Gibbs1.0 program, the heat capacity C_v, the volume thermal expansion coefficient α_v, the isothermal body modulus B_T and the Grüneisen parameter γ of the Ti_xVNbMo alloy are calculated. The curve of constant heat capacity C_v with the temperature of the Ti_xVNbMo alloy is drawn as figure 7. The results show that the heat capacity C_v value strongly depends on the temperature at low temperatures, which is due to the limitation of the Debye model at low temperatures. However, at a higher temperature, the C_v value is gradually close to the Dulong-Petit limit [33] and increases gradually with the increase of Ti content.

Zoom In Zoom Out Reset image size

The thermal expansion coefficient is a key parameter to characterize the thermodynamic and thermoelastic deformation behavior of alloys at high temperatures [36]. Applying the quasi-harmonic Debye model to calculate the thermal expansion rate is a feasible strategy, which has been verified in many papers [19, 31, 33, 37, 38]. Figure 8 shows the predicted value curve of the volumetric thermal expansion coefficient α_v of Ti_xVNbMo alloy with temperature. It can be seen from figure 8 that at low T, the volume thermal expansion coefficient αv increases with the increase of T3, and gradually increases rapidly with the increase of T, and then the increasing trend gradually becomes flat, and with the Ti content, this trend is more obvious.

Zoom In Zoom Out Reset image size

Not only the thermal expansion coefficient but also bulk modulus is a very important parameter to understand the thermal properties of solid materials at high temperature [39]. Figure 9 shows the predicted value of the isothermal body modulus B_T of Ti_xVNbMo RHEA with different Ti content between 0 ∼ 1150 K calculated by the Quasi-Harmonic Debye model. It can be easily seen from figure 9 that the isothermal modulus B_T gradually decreases with increasing temperature and Ti content, and the bulk modulus value at 0 K is equivalent to the value obtained from the previous DFT elastic constant calculation.

Zoom In Zoom Out Reset image size

Figure 10 shows the curve of the Grüneisen parameter γ with Ti content and temperature. The parameter γ increases gradually with the increase of temperature. At the same time, it can be seen from figure 10(a) that in the lower temperature stage, the Grüneisen parameter γ is negatively correlated with the increase of Ti content. However, this trend gradually changes after the temperature reaches 775 K; this negative correlation gradually changes to positive correlation, as shown in figure 10(b).

Zoom In Zoom Out Reset image size