Temperature dependent mechanical properties of Mo–Si–B compounds via ab initio molecular dynamics

Rapid development of technology is in constant pursuit of advanced materials with higher efficiency and durability. In the areas such as power generation and aerospace, the traditional structural materials such as Ni-based super-alloys have reached their physical limits in terms of their maximum operational condition temperature (∼1200 °C), which is already approaching 90% of its melting point.¹ By utilizing the high temperature coating and sophisticated cooling systems, the operating temperature can be increased only up to 1500 °C at most.¹ The circumstance requires materials not only having a higher melting point but also considerable strength, stiffness, and oxidation resistance at high temperatures as well as adequate mechanical behavior at room temperature.^2,3 Hence, the vigilant search for alternative materials has brought many researchers' attention to Mo–Si–B class of alloys.^2,4 Within this class of multiphase materials, the T1 (Mo₅Si₃) and T2 (Mo₅SiB₂) phases have shown high potential for the aforementioned applications by many recent studies.^2,5,6 The T1 phase has very high melting point (2180 °C) and excellent high temperature mechanical strength. The T2 phase also has excellent high temperature mechanical strength and good oxidation resistance.⁷ Further, T2 phase has a higher creep resistance⁸ comparable to high-temperature ceramics.⁵ The main drawback for the T1 phase, however, is its high thermal expansion anisotropy (TEA) with the α_a/α_c ratio measured to be as high as 2.0, resulting in high residual thermal stresses, which consequently limit its structural applications.^8,9

Previous works on the mechanical properties at high temperatures for these two compounds have been limited to a modest (∼1200 °C) temperature range. These properties, especially the elastic tensor, are not easily accessible via experimental methods due to the difficulty in synthesizing and performing high temperature experiments on single crystal T1 and T2 phases.⁵ In this regard, ab initio computational methods offer a viable alternative to obtain these critical properties. Numerous Density functional theory (DFT)-based ab initio methods have been developed for this purpose using total energy approaches,¹⁰ stress-strain approaches,¹¹ and density functional perturbation theory approaches.¹² Most of these methodologies work remarkably well for 0 K, but inadequate for higher temperature.¹³ The main objective of our study is to apply ab initio molecular dynamics (AIMD) to access the properties at high temperatures. AIMD method has become very popular as it encompasses, up to a greater extent, the accuracy of the DFT and ability to address high temperature^14–16 properties. However, due to the demand for computational resources, typical AIMD calculations are only applied to systems in the order of a few hundreds of atoms at most.^13,15,17 For Mo–Si–B system, such a size is adequate owing to the relatively small number of atoms in the unit cell (32 atoms) and thus a 2 × 2 × 2 supercell of 256 atoms would be sufficient. Since the study does not consider off-stoichiometric models, the demand for the size of the supercell is comparatively less. In our study, the crystal structures were optimized with AIMD method for each temperature and used to extract the high-temperature properties.

Both the T1 and T2 phases have body-centered tetragonal structure. T1 forms in a W₅Si₃ type crystal with D8_m space group symmetry. This crystal is stacked in a repetition of ABAC… layers along the c direction (Fig. 1(a)). The A-layers appear to have lower atom packing density than the B- and C-layers, which have Mo–Si bonds. However, the A-layers are close enough to have Mo–Mo bonds through B- and C-layers. The B-layer is 180° rotated about the c axis. The T2 unit cell consists of 20 Mo (A type layers), 4 Si (B type layers), and 8 B atoms per cell. The C-layer, however, consists both of Mo and B atoms. The A_{1/2 1/2} layer is similar to A-layer but shifted half of the diagonal of the unit cell. Figure 1(b) clearly shows the four layer types.

For all calculations we used the Vienna ab initio simulation package (VASP).^14,18 AIMD calculations in this study within NVT ensemble used following parameters: (1) The projector augmented wave method (PAW-PBE) potentials¹⁹ with the generalized gradient approximation (GGA); (2) electronic convergence criterion is set at 10⁻⁴ eV; (3) MD time step of 2 fs; (4) a high energy cutoff of 500 eV; an energy cutoff of 600 eV shows no discernible difference; (5) Γ point sampling was used.

Our approach is to first use AIMD^14,15 to get the equilibrium structure and lattice parameters at finite temperature. We assume that the structures obtained from AIMD embody all physical properties of the material at that temperature. Using the lattice parameters and atomic positions at respective temperatures, first-principles stress-strain method²⁰ was used to calculate the elastic tensor and the bulk elastic properties. It is essential in AIMD that the periodic supercell is sufficiently large to manage the fluctuations of thermodynamic variables, yet small enough to be computationally feasible.¹³ For the two crystals (T1, T2), a 2 × 2 × 2 supercell (256 atoms) was generated for AIMD to extract the T dependence of the volume. The procedure was to systematically perform many NVT ensembles by varying (isotropic) the volume at each temperature so as to obtain the equilibrium volume (V₀) at 1 atm (0.0001 GPa). Since T1 and T2 phases have a tetragonal crystal structure, we initially assume a volumetric contraction/expansion with a constant c/a ratio to estimate V₀ close to targeted pressure/temperature. We ensure the mean pressure and standard deviation of the mean is close to the targeted values. The data was then fitted to third-order Birch-Murnaghan isothermal equation of state²¹ (BM-EOS) to extract an initial estimate of the bulk modulus (K₀) and V₀. Since the crystals are not isotropic, the c/a ratio for a given temperature must be further modified. The above procedure would only enable us to provide a structure close to true V₀ at the targeted pressure. While this procedure gives the K₀, it clearly deviates from the experimental trend as shown in Figs. 2(b) and 2(c). This suggests that the estimated V₀ from the isotropic expansion is not close to the required equilibrium volume. Further analysis shows that although the total pressure on the supercell is close to the required value, the P_xx, P_yy, and P_zz components vary considerably. This implies that the anisotropy ratio (c/a) has to be accounted for as well as the volume. Thus, we subsequently evaluated the tensor components of the pressure varying the c/a ratio for the given V₀ to ensure that both the total pressure and its components are close to targeted value. After getting the equilibrium volume and c/a ratio for each temperature, a longer AIMD run of 4 ps is sufficient for convergence in pressure and temperature. The mean pressure ranged from 0.0005 to 0.2542 GPa and the standard deviation of the mean ranged from 0.0028 to 0.0138. The pressure fluctuations at each temperature for both crystals are shown in Fig. S1 in the supplementary material.²² Test runs with 10 ps show no change in mean pressure. We next calculate the thermal expansion by plotting a and c as a function of averaged temperature. Assuming a linear thermal expansion for a and c, the coefficient of thermal expansion (CTE) (α) for each axis can be calculated using |$\alpha _a = a_0 ^{ - 1} (da/dt)$| $α a = a 0 − 1 ( d a / d t )$ ⁠, where a represents the axis of the volume and a₀ is the respective value of the axis at the ambient temperature. By calculating the ratio α_c/α_a, we get the TEA. The reliability of the CTE in comparison with the experiment²³ is an indicator of the reliability of the AIMD configurations, which greatly encourages the following C_ij calculation (Fig. 2).

Because of the fluctuations, final output from the AIMD optimization is not representative of the required because of the fluctuations as mentioned earlier. Therefore, a sampling procedure was adopted prior to stress-strain analysis. Therefore, we took the last four structures of the most relaxed atomic configurations at a given NVT simulation that have the desired instantaneous temperature and pressure. The structures were then used to calculate the elastic tensors using the stress-strain method^20,24 as in the usual calculation for crystals at zero temperature. We assume that all the temperature effect is reflected in the structure obtained from filtered AIMD structures at that temperature. For this purpose, a strain of +1% and −1% is applied to the supercell obtained from AIMD and σ_j's data are obtained by taking the difference (σ_j⁺ − σ_j⁻)/2. From the calculated σ_j’s, C_ij (i, j = 1, 2, 3, 4, 5, 6) were evaluated by solving the linear equation |$\sigma _{\rm j} = \sum\nolimits_{j = 1}^6 {C_{ij} \varepsilon _j }$| $σ j = ∑ j = 1 6 C i j ɛ j$ ⁠. Thus, a set of 6 linear equations was obtained with 6 stress components and 21 elastic constants (C_ij is symmetric). The resulting C_ij are then averaged over the symmetry equivalent elements. From the C_ij values, we obtained the averaged bulk mechanical properties: K (bulk modulus), G (shear modulus), E (Young's modulus), and γ (Poisson's ratio) based on the Voigt-Reuss-Hill (VRH) approximation for polycrystals.²⁵ 

Thermal expansion is a good indicator to validate our method. Many previous works^6,8,26 had assumed linear thermal expansion within the temperature range they have considered. The expansion of a and c from this study is shown in Figs. 3(a)–3(d) together with the experimental data. The agreement is excellent for T1-a with Ref. 8, whereas the c-axis is overestimated by only 1.5% at low temperatures, but improves at higher temperatures. The alloy sample used in Ref. 8 is slightly off stoichiometric which could contribute to the discrepancy. A similar comparison for T2 phase shows the agreement improves with increasing temperature. Following the method used in the experimental work, we consider a linear fit for the thermal expansion of the a-axis and c-axis to calculate the CTE. The calculated CTE and TEA is tabulated in Table I along with the experimental data for the T1 alloys in Ref. 8, and the T2 single crystal data in Ref. 26. Experimental data in Ref. 8 show that the CTE is quite sensitive to the Si concentration and our results are closest to Mo_4.81Si₃, whose composition is the closest to Mo₅Si₃. We do however note that the measured CTE and TEA for single crystal Mo₅Si₃ in Ref. 6 are higher (22%) than our calculated values as well as those measured in Ref. 26. Overall, the AIMD simulation result shown in Table I was able to depict the high anisotropy of CTE in T1 phase, and lack of it in T2 phase is consistent with the experimental observations.

The T dependence of the C_ij for T2 using the above scheme is shown in Fig. 2(a), and compared with the experimental results up to 1400 K from Ref. 5. The agreement is reasonably well with a maximum deviation with the experiments on the C₆₆ of about 11%. It should be noted here that even though the AIMD simulations do not assume any symmetry, the temperature dependence of C₁₁ equals to C₂₂ as expected from a body centered tetragonal symmetry. This is also another confirmation that our 2 × 2 × 2 supercell is adequate to represent the elastic property of these materials. The stiffness in a and c directions that decreases with temperature is consistent with experimental observations. The relative strength of C₁₁ and C₂₂ is less compared with the experiment. This could be due to the fact that the c/a ratio is still not exact enough, so the stiffness C₃₃ is more than it should be. We found no experimental data for the T dependence of C_ij for the T1 phase but the calculated results are tabulated in Table II. Elastic constants of T1 are less compared to T2, indicating that T2 possesses mechanical properties better suited for high temperature applications. However, we found experimental C_ij data for T1 at room temperature,⁶ which is in excellent agreement with our results. Table II lists the averaged values of the 4 configurations filtered from the 4 ps AIMD run. The accuracy of the method depends on the standard deviation of the sampled models. Table ST I in the supplementary material²² listed each value for all the models considered, which shows very small deviation for the four sample models for the temperatures considered. This is yet another indication that a 2 × 2 × 2 supercell is sufficient for this study.

To compare our results with the experimental data, first we applied the VRH approximation²⁵ on the full elastic tensors (C_ij) from Ref. 5 to obtain the bulk elastic properties as shown in Fig. 2(c). This procedure is identical to the method used in our calculation to obtain bulk elastic properties. In addition, we also plot the K₀ by the fitting procedure as discussed in the methodologies, which has been a popular method in the previous DFT studies.^17,27 As clearly shown in Figs. 2(b) and 2(c), with the assumption of hydrostatic volume expansion (a constant c/a), the K₀ values derived from BM-EOS²¹ as a function of temperature deviate quite significantly from both the experimental data as well as the AIMD results even though they do follow the same trend. By contrast, our calculated results in Fig. 2(c) are much closer to those derived from the experiments.

We are further able to derive three other important elastic parameters; shear modulus (G), Young's modulus (E), and Poisson's ratio (γ) based on the analysis of the full elastic tensor results from stress-strain method. As shown in Fig. 2(c), these bulk elastic properties for the T2 phase agree with the results derived from Ref. 5. To the best of our knowledge, there are no experimental data available on the temperature dependent bulk elastic properties of the T1 phase. Nevertheless, the available values at room temperature⁵ are in excellent agreement with calculated ones. Our results show the bulk, shear, and Young's modulus decrease monotonically with increasing T, whereas γ increases with T and has a higher slope for T1 than for T2. At low T, γ is close to those typically obtained ceramics.²⁸ As the T increases γ increases to about 0.3, a typical value of a metal.²⁸ 

The analysis is limited to stoichiometric composition of T1 and T2 phases. Further, it should be mentioned that with regard to the high temperature behavior of intermetallic crystals, these present additional factors such as the effect of off-stoichiometry and environments as well as their potential impact on the dislocation mobility and the slip planes at very high temperatures are not considered. The effect of dislocations, we believe, should be minor at least up to 1500 K. Reference 5 noted that below 1500 °C (1773 K) only limited slip planes are active in the T2 phase. Similarly, the T1 phase has limited dislocation mobility at least below 1500 K.²⁹ The experimental data were conducted either under vacuum or inert gases,^5,6,8,26 and thus environment effects on the experimental results that we used for comparison should be minimal. Overall, above results suggest that the AIMD simulation in combination with the stress-strain method can be applied to obtain accurate C_ij values and mechanical properties at high temperature.

We have implemented a simple but very effective method to obtain the temperature dependent elastic and thermal properties of intermetallics in the ternary Mo–Si–B system via AIMD and the stress-strain response method. The combination of AIMD results and subsequent DFT's stress-strain method analysis shows that reliable assessment of elastic and thermal properties is feasible without resorting to explicit T corrections. The results are in excellent agreement with available experimental data. The calculations confirm the abnormal TEA of the T1 phase compared to T2 phase. Furthermore, this effect is shown by the Poisson's ratio as well. We observe a monotonic decrease in principle elastic constants with increasing T and the mechanical strength is higher in T2 compared to T1 and continues to follow the trend at higher T. This method with its limitations is fairly straightforward and has no inherent assumptions or extrapolations, which make it plausible to be applied to T1 and T2 phases in other TM–Si–B (TM = Transition Metal) systems so that the critical high temperature thermal/mechanical properties can be extracted.

This work was supported by the US Department of Energy (DOE), National Energy Technology Laboratory (NETL) under Grant No. DE-FE0004007. R.S. would also like to acknowledge a partial support from the NETL (DOE) under Grant No. DE-FE0007377 (PI: J. H. Perepezko). This research used the resources of the National Energy Research Scientific Computing Center supported by the Office of Basic Science of DOE under Contract No. DE-AC03-76SF00098.