Finite-temperature Mechanical Properties of fcc-disordered PtRh Alloys from First-principles Calculations

1. Introduction

Only a limited numbers of alloy systems are suitable for high-temperature and high-corrosion environments. Platinum (Pt) is one of the few element that exhibited a superior corrosion properties, owing to its non-reactivity with other elements. Pt maintains a single fcc crystal structure from room temperature up to the melting point. Rhodium (Rh) is one of the common alloying element for Pt. Rh has higher melting point than Pt and, similar to Pt, also maintains a single fcc structure until melting. The PtRh alloys formed a completely fcc solid-solution phase across the composition range up to very high temperature[1]. The strengthening mechanism in the Pt-Rh alloys system is mainly come from solid-solution strengthening. The PtRh alloys do not suffer from the degradation of mechanical strength on prolonged high-temperature exposure like other high-temperature alloys systems (Ni or Co superalloys) that depend on the precipitation strengthening.

For the development of high temperature applications, apart from phase stability of the alloys system, mechanical properties such as bulk, shear, Young’s modulus, and thermal expansion are very important. From the micro-scale microstructure prediction method such as phase field simulation (PF) to macro-scale stress analysis tool such as finite element method (FEM), it requires variety of input data such as thermodynamic stability, interfacial energy, elastic properties, and thermal expansion. Unfortunately, there are limited experimental data, especially, on both elastic properties and thermal expansion of the PtRh alloys. Most of available experimental data are related to the pure elements[2,3]. However, there are numbers of theoretical studies on the phase stability of the Pt-Rh binary system that determined the present of ordered phases (PtRh-I4₁/amd and PtRh₃-I4/mmm) at low temperature[4,5,6,7,8,9].

In this work, we generated the special quasirandom structure that represent the mimic the disordered properties of the fcc structure and performed the first-principles calculations to obtain the ground state energy of the disordered configurations. Finite-temperature free energy and lattice thermal expansion are calculated via the vibrational contribution using phonon quasi-harmonic approximation. The finite-temperature elastic properties of the fcc-disordered PtRh alloys are calculated using stress-strain method.

2. Methods

2.1. First-Principles Calculations

First-principles calculations, based on density functional theory (DFT), are performed to obtain the necessary data for ground state structures and elastic properties calculations. The Helmholtz free energy as a function of volume and temperature of a single structure ($F(V,T)$) can be partitioned into the different contributions:

$$F\left(V,T\right)=E_{0K}\left(V\right)+F_{vib}\left(V,T\right)$$

(1)

where $E_{0K}\left(V\right)$ and $F_{vib}(V,T)$ are ground state energy and lattice vibrational contributions, respectively.

Lattice vibrational contribution or phonon contribution is calculated through the phonon density of states (PDOS). Supercell method, as implemented in the Alloy Theoretical Automated Toolkit (ATAT)[10], are used for the PDOS calculations. The minimum separation between perturbed atoms are set to be about 10 - 12 Å to ensure that the interference from perturbed atom are isolated from periodic boundary condition. The 0.1 Å displacement for each perturbed atom is selected for the phonon calculations. For better accuracy in force constant matrix calculation during the phonon calculations, the 1^st order Methfessel-Paxton[11] is used. The lattice vibrational contribution is calculated using quasi-harmonic approximation[12] and obtained from the following equation[13]:

$$F_{vib}\left(V,T\right)=k_{B}T{\int }_0^{\infty }ln\left[2sinh{\displaystyle \frac{\hslash \omega \left(V\right)}{2k_{B}T}}\right]g\left(\omega \left(V\right)\right)d\omega$$

(2)

where ℏ, $\omega \left(V\right)$, and $g\left(\omega \left(V\right)\right)$ are the reduced Planck constant, volume-dependent phonon frequencies, and volume-dependent PDOS, respectively. The equilibrium volume at selected temperature can be obtained by minimizeing the $F_{vib}\left(V,T\right)$ term with respect to volume. The volumetric coefficient of thermal expansion (CTE) (${\alpha }_V$) at constant pressure can be expressed as:

$${\alpha }_V\left(T\right)={\displaystyle \frac{1}{V\left(T\right)}}{\left({\displaystyle \frac{\partial V\left(T\right)}{\partial T}}\right)}_P$$

(3)

Assuming the isotropic properties of cubic symmetry, the linear CTE can be calculated.

Vienna Ab initio Simulation Package (VASP)[14,15] is used for all the calculations to obtain the ground state energy of any structures. The projector augmented wave (PAW) method[16] is used for the electron-ion interactions descriptions. The exchange and correlations energy are described by the generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE)[17]. The PAW potential sets are comprised of a 5s¹5d⁹ configurations of valance electrons for Pt and 4s¹4d⁸ for Rh. All the calculations are performed with spin polarized enable. The k-point mesh of $\Gamma$-centered Monkhorst-Pack[18] is used for integration of the Brillouin zone during both structure relaxations and final static calculations during the total energy calculations. The 1^st order Methfessel-Paxton smearing method[11] is used with a 0.2 eV temperature broadening parameter for the electronic occupation. In order to obtain highly accurate energy and density of states (DOS) during the final static calculation, the Blöchl correction linear tetrahedron method[19] is performed. Table 3 lists k-point meshes for each structure. Each k-point mesh is selected to ensure that convergence of the electronic energies is lower than 0.1 meV with the cutoff energy of 500 eV.

2.2. Elastic Properties Calculations

The elastic stiffness constant ($c_{ij}$) is calculated using the stress-strain approach[20]. The fully relaxed structure ($\mathbf{R}$) is homogeneous deformed:

$${\mathbf{R}}^{\prime }=\mathbf{R}\left(\mathbf{I}+\mathbf{D}\right)$$

(4)

$$\mathbf{D}=\left(\begin{array}{ccc}{d}_{xx}& d_{xy}& d_{xz}\\ d_{yx}& d_{yy}& d_{yz}\\ d_{zx}& d_{zy}& d_{zz}\end{array}\right)$$

(5)

where ${\mathbf{R}}^{\prime }$, $\mathbf{I}$, and $\mathbf{D}$ are the deformed structure, 3x3 identity matrix, deformation matrix, respectively. By the definition, the $\mathbf{D}$ matrix possesses 9 degree of freedoms to deform a crystal. However, 3 out of 9 degree of freedoms are pure rotations which do not contribute the the change in the energy. $\mathbf{D}$ can be rewritten to the following:

$$\mathbf{D}=\mathbf{E}+\mathbf{U}$$

(6)

$$\mathbf{E}=\left(\begin{array}{ccc}{e}_1& e_6/2& e_5/2\\ e_6/2& e_2& e_4/2\\ e_5/2& e_4/2& e_3\end{array}\right)$$

(7)

$$\mathbf{U}=\left(\begin{array}{ccc}0& u_6& u_5\\ u_6& 0& u_4\\ u_5& u_4& 0\end{array}\right)$$

(8)

where $e_1=d_{xx}$, $e_2=d_{yy}$, $e_3=d_{zz}$, $e_4=d_{yz}+d_{zy}$, $e_5=d_{xz}+d_{zx}$, $e_6=d_{xy}+d_{yx}$, $u_4=\left(d_{yz}-d_{zy}\right)/2$, $u_5=\left(d_{xz}-d_{zx}\right)/2$, and $u_6=\left(d_{xy}-d_{yx}\right)/2$. For a small deformation, transformation matrix can be rewritten as:

$$\mathbf{I}+\mathbf{E}+\mathbf{U}\cong \left(\mathbf{I}+\mathbf{E}\right)\left(\mathbf{I}+\mathbf{U}\right)$$

(9)

It is evident that $\mathbf{I}+\mathbf{U}$ is a pure rotation, so the Equation 4 can be rewritten as:

$${\mathbf{R}}^{\prime }=\mathbf{R}\left(\mathbf{I}+\mathbf{E}\right)$$

(10)

In the elastic theory, $e_i$ are called strains. The energy change ($\Delta E$) due to elastic deformation can be expressed as:

$$\Delta E=V_0\left({\displaystyle \sum _{i=1}^6}{\sigma }_i^o{e}_i+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i,j=1}^6}{e}_i{c}_{ij}{e}_j\right)$$

(11)

and the correspondent stress ${\sigma }_i$ are:

$${\sigma }_i={\sigma }_i^o+{\displaystyle \sum _{j=1}^6}{c}_{ij}{e}_j$$

(12)

where $c_{ij}$ is elastic constants and ${\sigma }_i^o$ is the stresses at the reference point. The temperature-dependent elastic constant can be calculated by performed a series of calculations at different volumes. With the thermal expansion coefficient obtained from lattice vibration calculations, temperature-dependent elastic constant can be derived. The bulk modulus (B), shear modulus (G), and Young’s modulus (E) are calculated using Voigt’s method[21]:

$$B={\displaystyle \frac{{\overline{c}}_{11}+2{\overline{c}}_{22}}{3}}$$

(13)

$$G={\displaystyle \frac{{\overline{c}}_{11}-{\overline{c}}_{12}+3{\overline{c}}_{44}}{5}}$$

(14)

$$E={\displaystyle \frac{9BG}{G+3B}}$$

(15)

where

$${\overline{c}}_{11}={\displaystyle \frac{c_{11}+c_{22}+c_{33}}{3}}$$

(16)

$${\overline{c}}_{12}={\displaystyle \frac{c_{12}+c_{13}+c_{23}}{3}}$$

(17)

$${\overline{c}}_{44}={\displaystyle \frac{c_{44}+c_{55}+c_{66}}{3}}$$

(18)

The Poisson’s ratio ($\nu$) can be calculated from the following expression:

$$G={\displaystyle \frac{E}{2\left(1+\nu \right)}}$$

(19)

2.3. Disordered Structure

There are several methods to represent the disordered properties of the alloys structure. The most widely used techniques are the single-site coherent potential approximation (CPA)[22], the "coarse-graining" cluster expansion method (CEM)[23], and the special quasirandom structure (SQS) approach[24]. While CPA can successfully treats chemical and magnetic disordering for random alloys, owing to the nature of mean-field approximation, it fails to capture the locally dependent effects such as local displacement of atoms from the ideal lattice positions. CEM is very versatile and powerful approach for capturing the short-range order effects. However, CEM can be computationally expensive when consider the finite-temperature properties. The SQS represents the best possible periodic supercell structure that matches the local multisite correlation functions of a theoretical random alloys under a given N atoms in unit cell. It is very suitable for elastic properties calculations, since it provides a structure that represents the disordered properties of specific concentration of alloys. The SQS is proved to successfully use in the calculation of the elastic properties of multi-component alloys[25]. SQS is closely related to the CEM, since it uses the same correlation functions in CEM to determine the randomness of the structure. The correlation function for $\sigma$ configuration associated with a certain cluster $\alpha$ (${\rho }_{\alpha }\left(\sigma \right)$) can be written as:

$${\rho }_{\alpha }\left(\sigma \right)\equiv {〈{\Gamma }_{{\alpha }^{\prime }}\left(\sigma \right)〉}_{\alpha }$$

(20)

where ${\Gamma }_{{\alpha }^{\prime }}$ is cluster functions for cluster ${\alpha }^{\prime }$ and ${〈\dots 〉}_{\alpha }$ is the average taken over all cluster ${\alpha }^{\prime }$ that are equivalent by symmetry to cluster $\alpha$. The correlation function for the fully disordered state is site-independent and can be directly calculate from the average composition. SQS is determined by comparing the correlation functions of its structure with the theoretical disordered state. The SQS is generated using program "mcsqs" from ATAT package[10]. Table 1 shows the correlation functions of 16 atoms equiatomic-fcc-SQS against theoretical random configuration. It can be seen that the calculated pair (2-body) correlation functions of SQS are matched with the random structure up to 7^th pair and all the triplet (3-body) clusters up to 1.7321 times of smallest cluster are matched.

3. Results and Discussions

3.1. Disordered Structures

Three SQSs are generated for fcc-disordered PtRh alloys at 12.5, 25, and 50 at.% Rh. Two additional concentrations, 75 and 87.5 at.% Rh, are inverted elements from 25 and 12.5 at.% Rh, respectively. Figure 1 illustrates the unit cell at 12.5, 25, and 50 at.% Rh. The unit cell shape of generated SQS are constrained to make it easier for the elastic constant calculations. Thus, all the length of lattice vectors are equal. As an example, Table 2 shows the calculated correlation functions of the 50 at.% Rh SQS. The correlation functions of generated SQS are matched with theoretical random structure up to the 5^th pair cluster and all the triplet cluster that smaller than 2 times of smallest cluster. The 64 atoms unit cell are selected for the other composition owing to the requirement of the larger number of atoms to ensure that the possible configurations space are large enough.

Figure 2 shows the calculated 0 K enthapy of formation from DFT of SQS at different Rh compositions with three different ordered structures (Pt₃Rh-I4/mmm, PtRh-I4₁/amd, and PtRh₃-I4/mmm). It can be seen that, at 0 K, fcc-disordered at 50 at.% Rh is not stable compared to the PtRh-I4₁/amd ordered phase. Similar to PtRh-I4₁/amd, both Pt₃Rh-I4/mmm and PtRh₃-I4/mmm are also more stable than the fcc-disordered structure with the same composition at 0 K. However, if the convex hull is connected, Pt₃Rh-I4/mmm will not be part of the convex hull, since the coexistence between pure fcc-Pt and PtRh-I4₁/amd will have lower energy than Pt₃Rh-I4/mmm at the same composition. That is, Pt₃Rh-I4/mmm enthalpy of formation sits just slightly above the convex hull line that connected between fcc-Pt and PtRh-I4₁/amd. The difference between the energy of Pt₃Rh-I4/mmm and the coexistence between two phases is only 0.7 meV/atom (68.27 J/mol-atom). Pohl, J.[8] investigated the phase equilibrium of Pt-Rh binary system using semi-grand canonical lattice Monte Carlo simulation and shows that only PtRh-I4₁/amd and PtRh₃-I4/mmm are stable at very low temperature (below 250 K). Table 3 shows the ground state properties and k-mesh used for the structures.

3.2. Coefficient of Thermal Expansion

Figure 3 shows the calculated CTEs of fcc Pt and fcc Rh from DFT compared with two set of experimental measurements of Pt from literature[27,28]. The calculated results agree with low temperature data while high temperature measurement data start to deviate beyond 1500 K. The deviation from the experimental data at very high temperature might originate from the limitation of quasi-harmonic approximation. At high temperature, the contribution from the anharmonic vibration mode might be large enough so that it can affect the CTE. The calculated results for all three structures behaved similarly. The CTE of all three structures behaved similarly. The CTE values of 25 at.% Rh SQS fall between the CTE values of both Pt and Rh from low temperature until above 1500 K, where the CTE of SQS become higher than PT. The calculated CTE at 300 K for Pt and Rh are 8.65 x ${10}^{-6}{K}^{-1}$ and 8.07 x ${10}^{-6}{K}^{-1}$, respectively. It is very similar to the experimental measurement data[29] at 8.9 x ${10}^{-6}{K}^{-1}$ for Pt and 8.4 x ${10}^{-6}{K}^{-1}$ for Rh. It indicates that the calculated CTE results from DFT are reliable. Table 4 lists the corresponding CTE equation for Pt, Rh, and SQS at different concentrations. With the CTE results, lattice volume as a function of temperature can be obtained and correlated with the elastic properties calculations at different lattice volume. The elastic properties as a function of temperature can be estimated.

3.3. Elastic Properties

Figure 4 shows the calculated 0 K elastic properties ($c_{11}$, $c_{12}$, $c_{44}$, B, G, E, and $\nu$) of fcc-disordered at different Rh concentrations. It shows the monotonic change of all elastic properties across the composition range. There is no noticeable sudden change of in any elastic properties. The bulk modulus stays almost flat as the Rh concentration increase, but the Young’s modulus drastically increases as the Rh concentration increases. This is due to the increasing of the shear modulus with the Rh concentration. The enhancement of Young’s modulus is observed in the experiment and is also one of the reasons that Rh is used for alloying element of Pt.

Figure 5(a) shows the calculated temperature dependent Young’s and shear modulus for pure fcc Pt compared with the experimental data from Hamada et al.[30]. It can be seen that the calculated elastic properties are in agreement with the experimental data. The calculated results are slightly underestimated comparing to the experiment. The underestimated properties are possibly come from the selection of exchange and correlation functions used in the calculations. GGA exchange and correlation function is known to underestimate the bonding strength between atoms. This results is a slightly overestimation of lattice parameters and softening elastic properties. On the other hand, if local density approximation (LDA) is used for the exchange and correlation energy, underestimation of lattice parameters is expected. The calculated elastic properties of pure Rh in Figure 5(b) also shows the similar trend as appear in the pure fcc Pt. The extrapolation of elastic stiffness from experimental measurements[3] are slightly lower than the DFT calculations. The behaviors of Poisson ratio are clearly different between two elements. As the temperature increases, the Poisson’s ratio of fcc Pt also increase, on contrary, the values of fcc Rh decrease.

From room temperature up to high temperature applications, fcc-disordered phase is the stable structure across all the compositional range[1]. To calculate the elastic properties of fcc-disordered PtRh alloys, the SQSs are used for representing the disordered structures. Owing to the limited computational resources, only the temperature dependent elastic properties of SQS at 25, 50, 75 at.% Rh are performed. Figure 6 shows the calculated finite temperature elastic properties of fcc-disordered PtRh at three different Rh concentrations, (a) 25, (b) 50, and 75 at.% Rh. Almost all of the elastic properties are monotonically decreased as the temperature increased, except for the Poisson’s ratio of fcc-disorderd PtRh alloys at 25 at.% Rh. Figure 6(a) shows the increases of Poisson’s ratio as the temperature rising with the maximum around 400 K before trending downward. It seems that the raising of Poisson’s ratio in the low temperature region at 25 at.% Rh is inherited from the pure fcc Pt. As the temperature increase, the elastic properties is shifted toward the fcc Rh.

For typical Pt-Rh alloys, up to 7 wt.% Rh (∼ 12.5 at.% Rh) is added to improve both mechanical and corrosion properties. To calculate the elastic properties at this concentration, a larger SQS is needed. The computational resources required for the calculation will be too expensive. Using the interpolation from nearest calculated data point is one of the reasonable option. Figure 7 shows the calculated finite temperature elastic stiffness $c_{11}$ between direct SQS and linear interpolation from pure elements at equiatomic composition (50 at.% Rh). It can be seen that the difference between direct SQS calculation and interpolation is small. Even at 2000 K, the difference of $c_{11}$ is less than 20 GPa, which is less than 8 % error. The possible explanation that the linear interpolation works fairly well is possibly comes from the mixing energy behavior of fcc-disordered Pt alloys. According the CALPHAD methodology[31], Gibbs free energy of fcc-disordered phase can be expressed from the following equation:

$$G_m^{fcc}=x_{Pt}^0{G}_{Pt}^{fcc}+x_{Rh}^0{G}_{Pt}^{fcc}+RT\left(x_{Pt}ln{x}_{Pt}+x_{Rh}ln{x}_{Rh}^{ex}\right)+G_{mix}^{fcc}$$

(21)

where $G_m^{fcc}$, ${}^{0}G_i^{fcc}$, $x_i$, R, T, and ${}^{ex}G_{mix}^{fcc}$ are molar Gibbs free energy, Gibbs free energy of i element in $fcc$ structure, mole fraction of i element, gas constant, temperature, and excess mixing Gibbs free energy, respectively. Excluding the excess mixing term, the Gibbs energy expression represented the "ideal solution model", where the Gibbs free energy is the linear combination of the Gibbs energy of pure Pt and Rh plus the Gibbs free energy arising from fully random configuration of Pt and Rh. The excess term indicates the deviation of Gibbs free energy from the ideal solution model. Figure 8 shows the calculated excess mixing Gibbs free energy of fcc-disordered derived from DFT calculations against the Gibbs free energy of fully random configuration at 300 and 500 K. It can be seen that the excess mixing Gibbs free energy of fcc-disordered is smaller than the random configuration at 300 K. At 500 K, the excess term becomes only a fraction of random configuration. The values of excess mixing Gibbs free energy on the Pt rich are even smaller compared to the random configuration Gibbs free energy, which should make the mixing of Pt rich PtRh alloys behaved closer to the ideal solution. The same argument should be able to apply to the elastic properties of the fcc-disordered. The mixing energy of fcc-disordered PtRh is slightly negative. This should result in the strengthening of the bonding and increasing of elastic properties. As can be seen in Figure 7, the $c_{11}$ from SQS is actually higher than the interpolated $c_{11}$. Figure 9 shows the interpolation of elastic properties for fcc-disordered at Pt-12.5 at.%Rh. The elastic properties are monotonically changed across the temperature range. The Poisson’s ratio remains almost constant across the temperature range. According to Equation 19, it suggests that the ratio between Young’s and shear modulus remains almost constant across the temperature range.

4. Conclusions

Special quasirandom structures of fcc-disordered PtRh alloys are generated at different Rh concentrations with special constraints to accommodate the elastic constant calculations. Using the generated SQSs, the CTE are calculated via quasi-harmonic phonon approximation. The 0 K elastic properties of fcc-disordered PtRh alloys at different Rh concentrations and volumes are calculated using first-principles calculations. As the Rh concentration increase, the increasing of Young’s modulus is observed and is in agreement with the experimental observation. By coupled the elastic constant results with the CTE, the finite-temperature elastic properties can be obtained. Owing the the near ideal solution of the fcc-disordered PtRh alloys, the interpolation of elastic properties compared with direct calculation from SQS show very small differences. Thus, the finite-temperature elastic properties for fcc-disordered at PtRh alloys at 12.5 at.%Rh concentration are calculated using the interpolation from elastic properties of the pure Pt and SQS at 25 at.%Rh.