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

The virtual crystal approximation (VCA) method based on the Cambridge Sequential Total Energy Package (CASTEP) was used to establish the TiVNbMo refractory high-entropy alloy structure model. The effects of different Ti contents on the elastic and thermodynamic properties of TixVNbMo (x = 1.00, 1.25, 1.50, 2.00) high entropy alloys were calculated. The lattice constants calculation results of TiVNbMo with equal atomic ratio match well with the experimental values of vacuum arc melting, indicating that the VCA method is suitable for the first-principles calculation of TixVNbMo random solid solution. The EOS equation of state is used to determine the energy and volume of the equilibrium structure of the alloy. The elastic constants of TixVNbMo (x = 1.00, 1.25, 1.50, 2.00) high entropy alloys are calculated based on the body-centered cubic structure, and their Young’s modulus anisotropic three-dimensional contour stereograms are drawn. Moreover, the quasi-harmonic Debyeg-Grüneisen model is used to calculate the thermodynamic properties, such as thermal capacity, isothermal body modulus, volumetric thermal expansion coefficient, and Grüneisen parameter with Ti content and temperature.

For example, TiVNbMo has a single BCC solid solution phase and with a typical dendritic structure. Its lattice constant a=0.3211 nm, the yield strength σ y is up to 1200 MPa, the tensile strain at room temperature can reach 25.62%, and the hardness of TiVNbMo at room temperature is about 440.7 Hv [11]. These parameters indicate that TiVNbMo refractory high entropy alloy is a high entropy alloy with good ductility and high strength. Besides, first-principles calculations also show that TiVNbMo high entropy alloy has good ductility [12].
First-principles calculation is a popular theoretical research method in the field of materials research. The components in the HEA are randomly distributed in the structure, which leads to severe lattice distortion. These factors make modeling and computing elastic properties of the HEA difficult and expensive. At present, there is no uniform opinion on the establishment of high-entropy alloy models. Researchers are still in the exploratory stage. There are four common methods for establishing random models of high-entropy alloys: coherent potential approximation (CPA) [13], virtual crystal approximation (VCA) [14], special quasi-random structure (SQS) [15] and similar local atomic environment ( SLAE). Nevertheless, in the Exact Muffin-Tin Orbitals combined with the CPA (EMTO-CPA) method [16], the finite deformation of the closed structure is available [17,18]. SQS and SLAE require the construction of the supercell. The symmetry of the supercell is often broken to satisfy the correlation of different elements, which requires many calculations. On the other hand, although VCA is considered an oversimplified method to replace solid solution, many examples show that VCA can be reliably used to study high-entropy alloys composed of refractory elements [19][20][21][22]. Therefore, it is reasonable to explore the mechanical properties and thermodynamic stability of refractory HEA through the VCA method. It is reported that there is no literature report using ab initio method to study the effect of different Ti content on the mechanical properties and thermodynamic stability of TiVNbMo refractory HEA.
Therefore, in this work, based on DFT and plane wave pseudopotentials, the structure, elasticity, and thermodynamic properties of Ti x VNbMo (x=1.00, 1.25, 1.50, 2.00) RHEA with different Ti content were systematically studied through VCA. First of all, the accuracy of the calculation results is verified by comparing the calculation results when the equal atomic ratio, that is, x=1.00, with the known experimental results (table 2). Then the elastic and thermodynamic properties of Ti x VNbMo (x=1.00, 1.25, 1.50, 2.00) are studied through the VCA model, equation of state (EOS), and quasi-harmonic Debye-Grüneisen model. For the convenience of data recording, these five alloys will be marked as Ti 1.00 , Ti 1.25 , Ti 1.50 , Ti 1.75 , and Ti 2.00 .
2. Alloy composition design and calculation method 2.1. Composition design of Ti x VNbM RHEA 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 ( 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, DH 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, r is the average atomic radius. c D is the electronegativity difference of the HEA, c is the average electronegativity difference of the constituent elements calculated by the mixing rule, c 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)  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, w S is the solid angle of the smallest atom, w 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) Where A i and r i are the atomic weight and density of element i.

Calculation details
In this work, the energy and elasticity of Ti x VNbMo (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 −5 eVatom −1 . Besides, in the subsequent geometric optimization, all the forces on a single atom converged to less than 0.03 eV Å −1 , and the total stress tensor was reduced to the order of 0.05 GPa, max placement of an atom is less than 0.001 Å.

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 Å ˳ Where V 0 and E 0 are the equilibrium volume and the corresponding energy, respectively. V is the volume of change, B 0 is the bulk modulus, and ¢ B 0 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.

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) - ( ) 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] ( ) The calculation method of modulus for cubic crystal is shown in equation (7) [16] 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) The general anisotropy index A U is calculated by equation (9) [29] The hardness value of the alloy is calculated according to the equation

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 x VNbMo 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 is the total energy of each unit cell of Ti x VNbMo 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 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).
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) 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 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

Elastic constants of TixVNbMo
The elastic stiffness, polycrystalline elasticity, and hardness of Ti x VNbMo 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 x VNbMo alloy gradually increases from 3.182 Å to 3.202 Å, and the theoretical density ρ theor gradually decreases from 4.28 g cm −3 to 3.85 g cm −3 . 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 x VNbMo 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 firstprinciple calculation result of the elastic properties of Ti x VNbMo 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 x VNbMo RHEAs are positively correlated with the Ti content, and this trend is more significant when x 1.50. Figure 4 is a dotted line graph of the predicted values of the general anisotropy index and the hardness of Ti x VNbMo 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.
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 x VNbMo 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 x VNbMo RHEA gradually tends to form a completely spherical shape, indicating that they gradually tend to be isotropic. Projection of Young's modulus of Ti x VNbMo 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.

Thermodynamic properties of TixVNbMo
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 x VNbMo alloy are calculated. The curve of constant heat capacity C v with the temperature of the Ti x VNbMo 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. 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 x VNbMo 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.
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 x VNbMo 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. 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).

Conclusion
Based on the CASTEP computing platform, the crystal structure is approximated by VCA, combined with the EOS equilibrium equation of state and the quasi-harmonic Deby-Grüneisen model, the phase structure, elastic constants and thermodynamic properties of Ti x VNbMo RHEA are systematically studied. Through comparison  with other experiments, it is confirmed that VCA modeling is suitable for the calculation of lattice constants, elastic properties, and thermodynamic properties of Ti x VNbMo RHEA.
The calculation results show that the Ti x VNbMo RHEAs has a stable metal bond bonding BCC structure. With the increase of Ti content, the B, G, and E of Ti x VNbMo RHEA gradually decrease, indicating that the average bonding capacity and deformation resistance between atoms become weaker at the same time, and the alloy is more prone to elastic deformation. The values of Poisson's ratio σ and B/G ratio gradually increase with the increase of Ti content, which means that the plasticity and toughness of the material gradually increase, and when x 1.50, this trend becomes more and more significant. The Young's modulus of Ti x VNbMo RHEA is the most anisotropic when x=1.00. With the increase of Ti content, the alloy tends to be isotropic more and more, and the hardness H v decreases gradually.
The calculation of the Gibbs1.0 program proves that the heat capacity C v of Ti x VNbMo RHEAs is positively correlated with the increase of Ti content, and the isothermal bulk modulus B T and volume thermal expansion coefficient α v at different temperatures are negatively correlated with the increase of Ti content. When the temperature is higher than 775 K, the Grüneisen parameter γ and Ti content of Ti x VNbMo RHEA gradually change from a negative correlation to a positive correlation. In summary, the work in this paper provides good guidance for the design of refractory alloys for BCC structures.