Introduction

As the world pursues carbon-neutral power generation, the gas turbine is expected to play a vital role. Modern gas turbines can burn hydrogen, which can be produced, amongst other methods, via water electrolysis using wind power during off-peak hours1,2. Hydrogen burns hotter than natural gas3, which makes gas turbines more efficient assuming operational temperatures can be pushed higher without damaging the turbine blades. The current generation high-speed blade materials are Ni-based superalloys with melting points near 1350 °C, reduced below 1250 °C at the zone between the bond coat and single-crystal blade. They are operated at 1100 °C and require thermal-barrier coatings and internal cooling to withstand the working gas as hot as 1600 °C4,5. A next-generation material that can operate at/above 1300 °C, remain phase stable while retaining strength, and permit the use of hotter gas (while reducing the need for internal cooling) will improve fuel efficiency by as much as 7% and reduce environmental-damaging emissions6.

Refractory multi-principal-element alloys (RMPEA) are promising alloys for high-temperature (HT) applications7 if stability, strength, and oxidation resistance in harsh-service conditions can be addressed. Using only the high-entropy alloy concept presupposes that higher configurational (mixing) entropy \({S}_{{{\rm {conf}}}}\) that competes with formation energy \({E}_{{\rm {f}}}\) [not too positive (clustering) or too negative (strongly ordering)] and lattice strain lead to higher stability and strength8,9. The high entropy alloy (HEA) concept10 has resulted in a class of materials that promise improved strength or toughness11,12,13, especially at cryogenic temperatures8,14,15,16,17,18, and corrosion properties19,20,21,22,23. High entropy was identified for higher phase stability in the earlier studies, but a more careful consideration shows it is not a major factor24. Notably, with more alloying elements, N, introduces the more likely possibility (e.g., from electronegativity differences) that at least one pair of \(\tfrac{1}{2}\) N(N−1) pairs will be strongly ordered, so alloying effects and chemistry, i.e., non-equiatomic compositions, remain very important. Ideally, RMPEAs are single-phase systems with extremely high melting points, and high specific strength via solid-solution strengthening well above 1000 °C. The first two RMPEAs reported with HT strength well above other metal alloys are equiatomic solid-solutions MoNbTaW and MoNbTaVW9,25. They offer exceptional compressive yield strength of 400–700 MPa at 1600 °C, whereas in typical Ni-based superalloys, e.g., Inconel 718, yield strength drops to near zero near 1200 °C8. Such a significant difference in mechanical strength at high temperatures testifies to the potential of RMPEAs. Unfortunately, these two RMPEAs are brittle at room temperature (RT). Efforts were made to improve their ductility, but always resulted in a significant loss of HT strength. For example, the equiatomic HfNbTaTiZr alloy exhibits good RT ductility, up to 50% deformation can be reached under compression26, but the HT yield strength is low, only 92 MPa at 1200 °C27. Another alloy that is promising is HfMoNbTaTi alloy, which shows 27% RT compressive strain. Its compressive yield strength maintained at 699 MPa at 1200 °C then dropped to 367 MPa at 1400 °C28. When constrained by high configurational entropy (i.e., N > 4 with equiatomic composition), finding an RMPEA viable for HT application is quickly limited. Fortunately, numerous non-equiatomic RMPEAs do indeed retain single-phase solid solutions (and multiphase for strengthening, controlled by composition) that exhibit superior mechanical properties (e.g. refs. 8,29,30). By removing the equiatomic constraint, the MPEA composition space expands significantly yet demands an excessive amount of computational and experimental work. The high-throughput (HTP) scheme is one such approach31,32,33,34,35,36,37,38,39,40,41 which, when combined with numerically efficient density-functional theory (DFT) methods42, can enable accurate MPEA design29. Experimentally, methods like rapid alloy prototyping31, diffusion-multiples32, additive manufacturing43,44, and thin-film-related co-deposition33,45 have been successful. Recently, computation materials science has started to play a critical role in HEA development, including works like CALPHAD-based high-throughput computational method24,34, machine learning assisted design coupled with modeling/experimental data35,36,46, and density functional theory (DFT) method37. However, past effort focused on either computation or experiments with a limited synergy of the two. Moreover, there has been limited or no focus on developing numerically efficient ways to directly use DFT in MPEA design.

Here, based on our recent success42, we developed a numerically efficient screening strategy combining HTP-DFT with experimental validation along with high-fidelity testing (RT and HT) for efficiently exploring the vast RMPEA compositional space. The strategy focuses on assessing and down-selecting worthwhile alloys under a chosen set of constraints: (Tier 1) HTP-DFT calculations (with specific metrics) are used; (Tier 2) secondary down-selection from exhaustive DFT results are then validated and additionally screened by rapid bulk-alloy synthesis and characterization at RT (modulus); and (Tier 3) small-sample punch tests are performed at high temperatures for creep and tensile mechanical properties. By applying Tier 2 and 3 screening, this strategy led us to the discovery of a Mo-rich RMPEA with well-balanced mechanical properties within 12 months.

Results and discussion

Alloy design and screening

Tier 1 screening and down-selection employed the well-established, DFT-based Korringa–Kohn–Rostoker coherent potential approximation (KKR-CPA) calculations to average directly all possible alloy configurations and to predict the relative structural phase stability (i.e., formation energy, \({E}_{{\rm {f}}}\); Fig. 1a), structural parameters (lattice constants, {\({a}_{{i}}\)}), and bulk moduli (B; Fig. 1b) for the entire design space, i.e., for any arbitrary composition solid-solution29,47,48,49,50. Importantly, recently we have shown that the interstitial electron density \({\rho }_{\rm {{o}}}\) (obtained directly by KKR-CPA) provides a universal metric \({r}_{\rm {{s}}}\) for maximum strength42 for any metal (elemental, crystalline, or amorphous); \({r}_{{\rm {s}}}\) (Fig. 1c) is a DFT exchange-correlation parameter (with \({\rho }_{{\rm {o}}}\) = \(\tfrac{4\pi }{3}{r}_{{\rm {s}}}^{3}\)) for a homogeneous electron gas—which exhibits universal behavior for all elastic moduli and shear strength. As \({r}_{{\rm {s}}}\) involves only interstitial electrons (independent of any alloying element), a rule-of-mixtures (RoM) \({r}_{{\rm {s}}}^{{{\rm {RoM}}}}\) is valid and, notably, provides a reliable a priori estimate of relative strength, as we have already verified42 (i.e., \({r}_{{\rm {s}}}^{{{\rm {RoM}}}}\) down-selected regions agree with direct KKR-CPA results as we verified experimentally in a subset of alloys). Hence, \({r}_{{\rm {s}}}^{\rm {{{RoM}}}}\) for any RMPEA composition directly narrows the search ranges (Fig. 1c) for alloys with the highest strength criterion (largest \({r}_{\rm {{s}}}^{\rm {{{RoM}}}}\)’s—with \({r}_{{\rm {s}}}\) post-verified by KKR-CPA). Subsequently, those alloys with the most stable criterion (\({E}_{{\rm {f}}}\) from KKR-CPA) were down-selected for single-phase (bcc) stability. The smaller the \({r}_{\rm {{s}}}\) the higher the strength, where the smallest \({r}_{{\rm {s}}}\) for elemental metals is 1.5 as found for Iridium42. From the four criteria, \({r}_{{\rm {s}}}\) and formation energy \({E}_{\rm {{f}}}\) can show competing behavior, for example, the alloy is more brittle if \({r}_{\rm {{s}}}\) indicates higher strength but \({E}_{\rm {{f}}}\) is too negative (too stable or too strongly bonded), as is verified by the ductility criteria, as discussed in results. Notably, the high-strength alloys down-selected and characterized here were the same as given initially by \({r}_{\rm {{s}}}^{\rm {{{RoM}}}}\).

Fig. 1: Contour plots for (MoW)–Ta–(TiZr) system (Ti:Zr is fixed at 1).
figure 1

DFT \({E}_{{\rm {f}}}\) [meV] (a), hydrostatic, bulk compressive single crystal modulus B [GPa] (b), DFT universal strength metric \({r}_{{\rm {s}}}\) (c), and Pugh’s ratio (G/B) (d); DFT results are derived from 0 K calculations.

For the five-element MoWTaTiZr system, we directly evaluated 8000 compositions using KKR-CPA with Mo, W, and Ta varied at 5 at.% intervals and with Ti:Zr ratio fixed at 1:1 (see ref. 29, small amounts of Ti:Zr provide an electron-dispersion effect to stabilize the alloy). To expedite the down-selection of compositions with −150 ≤ \({E}_{{\rm {f}}}\) ≤ +70 meV and bulk moduli B > 200 GPa, we used alloy lattice constants estimated by Vegard’s rule51 (i.e., \(\bar{a}\approx {\sum }_{i=1}^{N}{c}_{i}{a}_{i}^{{\rm {o}}}\), with compositions \(\{{c}_{{i}}\}\) and \({a}_{{i}}^{{\rm {o}}}\) are the optimized elemental lattice constants); down-selected alloy composition was then fully optimized as needed. The \({E}_{\rm {{f}}}\) criteria arise29 from the observation that complex phases (e.g., sigma phase) form when \({E}_{\rm {{f}}}\) > +70 meV and metastable metallic glasses form or ordered phases appear when \({E}_{{\rm {f}}}\) < −150 meV52. The alloys that met the \({E}_{\rm {{f}}}\), \({r}_{\rm {{s}}}\), B0, and Pugh’s ratio for ductility criteria are depicted in Fig. 1. A further screening was conducted based on forming a solid solution using the Hume-Rothery rules, which describe that solute and host atoms must be similar in electronegativities, crystal structure, atomic radii (difference < 15%), and valence-electron concentration (4 < VEC < 6). The compositions that satisfied the Hume–Rothery rules were examined by calculating the Warren-Cowley short-range order parameters using KKR-CPA linear response, which predicts (if the alloy forms a disordered solid solution) what low-T ordering or spinodal decomposition is expected53. The Tier 1 screening generated four compositions of interest: (MoW)80Ta12.5(TiZr)7.5, (MoW)80Ta10(TiZr)10, (MoW)80Ta20(TiZr)0, and (MoW)85Ta10(TiZr)5.

The ductility was screened using first-principle DFT calculated shear G and bulk B moduli, by assessing inverse Pugh’s ratio G/B empirical rule for a ductile material (Fig. 1d), i.e., G/B < 0.57–0.6054,55. The four compositions from Tier-1 screening all have inverse Pugh’s ratios G/B < 0.42.

Tier-2 screening used more exhaustive calculations to examine the 4 alloys down-selected from the Tier-1 screening and also the neighboring space (±10 at.%) around each group (MoW), Ta, and (Ti,Zr). Four groups, each containing 8 compositions, were designed around the 4 alloys (Table 1 and Supplementary Table 1). Some compositions also included Al and Cr additions for enhanced oxidation resistance and density reduction considerations. The calculated properties for these 32 compositions include formation enthalpy (\({E}_{{\rm {f}}}\)), universal strength metric (\({r}_{{\rm {s}}}\)), shear modulus (G), bulk modulus (B), inverse Pugh’s ratio (G/B), Poisson’s ratio (ν), and Young’s modulus (E), as well as yield strength (σy) and its temperature dependence, through a model56. In Table 1, we have provided the universal metric \({r}_{{\rm {s}}}\) (for maximum strength42) and Young’s moduli for the produced alloys. Given high B (~264 GPa) and E (~274 GPa) along with optimal formation energy (~−65 meV per atom) and Pugh’s ratio (~0.42 < 0.57) range, the tier-2 down-selection process predicts AL-A6 and AL-A8 as best MPEAs. The AL-D6 and AL-D8 are two other MPEAs that fall in a similar category (based on B and E) as AL-A6 and AL-A8. However, they do not satisfy the optimal \({r}_{\rm {{s}}}\) and ductility criteria.

Table 1 The nominal composition of 32 Ames Lab combinatorial arc-melted samples (AL) A1–A8, B1–B8, C1–C8, and D1–D8, along with DFT universal strength metric rs and Young’s moduli E.

The DFT-based \({E}_{\rm {{f}}}\) suggests 32 compositions are all energetically stable in bcc (A2) phase (Fig. 2). It shows that adding Al enhances (Cr decreases) the MPEA stability. Singh et al.29,53,57. attributed the effect of Al addition on stabilizing the bcc phase to stronger electronic hybridization.

Fig. 2: DFT formation energy \(E_{{\rm {f}}}\) of the 32 samples.
figure 2

Ames Lab (AL) A1–A8 (a), B1–B8 (b), C1–C8 (c), and D1–D8 (d). RT equivalent energy is denoted (horizontal dashed line) at 25 meV atom−1. Aluminum additions enhance MPEAs stability (\({E}_{\rm {{f}}}\) is more negative), while Cr decreases stability (\({E}_{\rm {{f}}}\) is less negative).

The valence-electron concentration (VEC) of the 32 compositions varies between 5.6 and 5.8, consistent with the formation of an A2 phase29,58. Though a VEC < 4.5 was suggested by Sheikh et al.59 for intrinsic ductility for refractory high-entropy alloys, concentrated refractory alloys seem to deviate from this criterion60. Electronegativity difference on the Allen scale61 was calculated to be 2–3%, which is not significant enough to favor the formation of intermetallics. Nearly all criteria for solid-solution formation are satisfied for the four groups of compositions. Notably, the failure mode or ductile–brittle transition in crystalline alloys is often connected to inverse Pugh’s (G/B) and Poisson’s (v) ratio62 as G and B directly represent resistance to plastic deformation and fracture in polycrystalline materials, respectively. The inverse Pugh’s ratio of these compositions is 0.35 ≤ G/B ≤ 0.43, below the ductility criteria of <0.5763 (see Supplementary Fig. 1), which allows initial filtering of useful RMPEAs.

Although tensile strength is most relevant (more on that below), e.g., for machinability, compressive strength for expediency is often reported in initial work. Hence, compressive yield-strength (solid solution portion only throughout the paper) for these 32 compositions was evaluated using a yield-strength (\({\sigma }_{\rm {{y}}}\)) model described in the “Methods” section29. The model was validated using well-characterized MoNbTaW and MoNbTaVW alloys25. At RT, the DFT-calculated and experimentally measured values are (\({\sigma }_{\rm {{y}}}^{{\rm {{DFT}}}}\), \({\sigma }_{\rm {{y}}}^{{\rm {{EXP}}}}\)) = (970, 1058) MPa for MoNbTaW and (1170, 1246) MPa for MoNbTaVW25. Similarly, at 1300 °C, (\({\sigma }_{\rm {{y}}}^{{\rm {{DFT}}}}\), \({\sigma }_{\rm {{y}}}^{{\rm {{EXP}}}}\)) = (367, 450) MPa for MoNbTaW and (655, 700) MPa for MoNbTaVW25.

To understand the temperature dependence of elastic properties, we fitted the DFT-calculated Young’s modulus E (0 K) using a linear model proposed by Wachtman et al.64 at 300 K, which was later applied to Cij’s by Steneteg et al.65. These DFT results agree with our measurements at RT with an average error < 1% (Fig. 3a). The insignificant change of E modulus is reasonable as the changes in bonding and thermal expansion from 0 K to RT are small (as expected from Grüneisen theory)66,67,68. Extrapolating from RT to 1300 °C, E drops more noticeably, by an average of 7.5% (Fig. 3b). Those compositions with higher melting temperature drop less, whereas the ones with lower melting temperature, e.g., Al-doped samples, show the largest drop (15%). This suggests that maximizing the melting temperature for MPEAs may maximize mechanical properties (a strategy however that has an adverse impact on RT ductility), which is well explained by the modulus-temperature relationship \({E}({T})/{E}(0)=\,[1-{\rm{\alpha }}({T}/{{T}}_{{\rm{m}}})]\), where E is the modulus at temperature T in Kelvin (E(0) is at 0 K) and α is a constant of order 0.569.

Fig. 3: Young’s moduli, E, by experimental and DFT.
figure 3

Comparison between RT experimental and DFT Young’s moduli, E (a), showing good agreement; E at 1300 °C shows a minor drop from RT values; however, E decays much faster (as marked) on Al added samples MPEAs (b).

Rapid sample preparation and verification

These 32 compositions were synthesized and characterized using a set of HTP methods tailored for this study. Small samples (1.5 g each) were prepared using a combinatorial arc melting system capable of making 32 samples in one run. The obtained samples were characterized for compositions using scanning X-ray fluorescence (XRF), crystal structures and phases using X-ray diffraction, density using the Archimedes’ method, elastic moduli (E, G, \(\nu\)) using ultrasonic pulse-echo measurement, and \({\sigma }_{\rm {{y}}}\) using RT compression test.

The as-cast samples (quench rate ~103 °C/s) are all nearly a single solid solution A2 (bcc) phase based on X-ray diffraction analysis (Fig. 4a–d). The measured lattice parameters match well with DFT optimized results (Fig. 4e) for group AL-A and AL-D, while only a very small deviation was observed for group AL-B and AL-C, indicating that the as-cast samples form a continuous solid solution across the phase space29. The densities of all samples are over 10 g/cm3 (Fig. 4f). The experimentally measured densities agree well with a rule-of-mixture estimate from the individual elements, again indicating solid-solution formation rather than intermetallic compound formation. Though the densities of RMPEAs are higher than typical nickel superalloys, their specific strengths are higher, especially above 1000 °C as shown below.

Fig. 4: Phase information, lattice parameter, and density of the 32 samples.
figure 4

XRD pattern of the samples (ad). Lattice parameter from theory and from XRD (e). Density from the rule of mixtures and Archimedes measurement (f). Compositions of the 32 samples can be found in Table 1 and Supplementary Table 1.

Standard size sample preparation and characterization

Based on the measured density, melting point, strength metric, and RT moduli, samples AL-A3: (Mo95W5)85Ta10(TiZr)5, AL-A6: (Mo85W15)85Ta10(TiZr)5, and AL-A8: (Mo95W5)85Ta5(TiZr)5Cr5 were selected for a detailed Tier 3 screening, showing good a priori strength metric. The AL-A3 alloy was selected for its relatively low density while having a high melting point and modulus; AL-A6 was selected for its high melting point and high modulus; AL-A8 was selected for its low density, strength metric, and high modulus.

The as-cast AL-A3, A6, and A8 exhibit dendritic structures typical for arc-melted RMPEA due to constitutional segregation during solidification8,9,29. The samples remain solid solution going from tier 2 (1.5 g) to tier 3 selection (~50 g) though there were minor reductions in the cooling rate with the larger mass. It is not clear if the identified alloys will remain solid solutions after casting when the size of the cast material is scaled up to kg or higher. Scaling may be addressed by powder metallurgy approaches like gas atomization followed by hot-isostatic-press approach. The Z contrast in the backscattered SEM images (Fig. 5) shows that the dendrites are enriched while the interdendritic regions are depleted in higher Tm atoms. The EDS line scan across the dendrites further confirms that the heavier atoms W and Ta are concentrated in the dendritic region while the lighter elements Zr, Ti, and Cr are concentrated in the interdendritic region.

Fig. 5: Microstructure and compositional information of the as-cast samples.
figure 5

SEM backscattered images of ALA3 (a), ALA6 (b), and ALA8 (c) samples. EDS line scan showing the element distribution for ALA3 (d), ALA6 (e), and ALA8 (f) samples.

Mechanical properties by compression test

The yield strength \({\sigma }_{0.2}\) and compressive strain before fracture for AL-A3, A6, and A8 are (820 MPa, 10.0%), (760 MPa, 9.2%), and (750 MPa, 7.8%), respectively (Fig. 6). An equiatomic MoNbTaW that has high strength but limited ductility, based on the work by Senkov et al.25 was prepared and tested in the same manner. The results show that MoNbTaW has higher \({\sigma }_{0.2}\) (1010 MPa) but less compressive strains (2.8%) than A3, A6, and A8. AL-A6 was also tested at 1300 °C in compression. The compressive yield strength at 1300 °C for Al-A6 is 446 MPa, similar to MoNbTaW (460 MPa) and less than MoNbTaVW (655 MPa)25. However, AL-A6 is 3× more ductile than the equiatomic MPEAs and shows less temperature dependence of its yield strength. The high W content in MoNbTaW alloy likely contributed to its high yield strength, increased density, and reduced ductility. AL-A3, A6, and A8 have larger atomic size differences (\(\delta\) = 2.64, 2.62, and 3.22, respectively) than MoNbTaW (\(\delta\) = 2.28), which may boost yield strength due to increased solute-induced stress26, but not enough to be comparable to MoNbTaW. AL-A8 has the highest \(\delta\) due to the introduction of smaller and lighter Cr atoms, and it has the highest specific yield strength. The improved ductility can be partly attributed to the Ti and Zr added to the alloys, which includes filling a favorable electronic dispersion at the Fermi energy29, enhancing the atomic interstitial electron density of the metal and thus its strength42. Similar effects of Ti and Zr on improving ductility were observed for HfNbTaTiZr26, AlHfNbTaTiZr70, HfMoNbTiZr71, HfNbTiVZr72, and NbTiVZr 73alloys. Compared to systems reported by Tseng et al.28 (HfMoNbTaTi, HfMoNdTaZr, and HfMoTaTiZr), the yield strength of AL-A6 is much higher than that of HfMoTaTiZr and comparable to the other two at 1300 °C. AL-A6 also shows higher ductility than HfMoTaTiZr, although it is not as ductile at RT as the other two (see Table 2). Tseng et al.28 contribute higher ductility to the Nb addition. As Nb has the same VEC as Ta, we did not include Nb in our current HTP study. It should be noted that since Nb and Ta can have large differences when interacting with Ti and Zr, the phase formation and ductility can differ significantly.

Fig. 6: Compressive properties of the samples.
figure 6

Compression stress–strain curves at RT for AL-A3, A6, A8 samples (a); at 1300 °C for ALA6 (b); and comparison of compressive yield strength (YS) (c) and specific compressive YS (d) of ALA6 with other MPEAs11,14,36—open (filled) symbols show <5% RT (≥10% RT) fracture strain.

Table 2 Density and compressive (RT and 1300 °C) properties for AL-A3, A6, A8, Mo NbTaW, HfMoNbTaTi, HfMoNbTaZr, and HfMoTaTiZr.

Tensile and creep properties by small punch test

While most literature only reports compressive properties, tensile properties are of the most engineering importance. It was shown that small punch tests (SPT) provide tensile yield strength when the small-size membrane sample is punched through undergoing a stretching deformation74. Our SPT results show that the yield points for as-cast AL-A3, A6, and A8 are 150, 240, and 193 MPa at 1300 °C, respectively (Fig. 7). The SPT strength of ALA6 is considerably higher than that of TZM (Mo–0.5%Ti–0.1%Zr alloy, ASTM B386 TYPE 364) (~110 MPa) and H282 (~20 MPa) at 1300 °C. We also carried out SPT on equiatomic MoNbTaW and MoNbTaVW, they all fractured before any meaningful loading. The deformed regions of the as-cast samples in Fig. 7a assume half-dome shapes, indicating a significant degree of ductility. The fractured surfaces of the samples punch tested at 1300 °C (Fig. 7c–e) offer additional insights into the fracture behavior. The as-cast samples show river patterns, cleavage feathers, and shear ledges on the fractured surface, indicating a quasi-ductile transgranular fracture mechanism.

Fig. 7: Small punch test results at 1300 °C and its setup.
figure 7

Stress versus strain curves for AL-A3, A6, and A8 (a), and for AL-A6, TZM (Mo–0.5%Ti–0.1%Zr alloy), and Haynes-282 (b). SEM images of fracture surfaces after 1300 °C punch tests for ALA3 (c), ALA6 (d), and ALA8 (e); small-punch rig setup (f), die (g), and schematics (h). Samples were tested in an inert atmosphere and not coated.

AL-A3, A6, and A8 alloys are all highly creep resistant, with AL-A6 exceeding the properties of any known refractory alloys. Small punch creep (SPC) tests were conducted at 1000 and 1300 °C with 200 and 100 MPa for up to 100 h, respectively. The SPC creep deformation curves are shown in Fig. 8 for AL-A6, a typical fully age-hardened Ni-based superalloy Haynes 282, and a typical refractory alloy TZM, consolidated via powder metallurgy route. All the samples entered the steady-state creep stage without rupture at the 100th hour. At 1000 °C 200 MPa, the minimum creep rate of AL-A6 is 1.55 × 10–5 h−1, much lower than that of Haynes-282 (2.56 × 10–4 h−1) and TZM (2.83 × 10–4) tested under the same condition. Haynes-282 was not tested at 1300 °C due to its low melting point 1348 °C. AL-A6 maintained a low creep rate 2.99 × 10–5 h−1 even at 1300 °C 100 MPa. The creep resistance of AL-A6 is partially a result of its high melting points, as it is creep tested at 42–53% melting point where solid-solution strengthening still dominates. The large lattice strain and atomic-size mismatch from the multiple elements may also contribute to the increase in resistance to dislocation initiation and multiplication.

Fig. 8: Small punch creep tests of H282, TZM, and AL-A6 samples at 1000 and 1300 °C.
figure 8

Load is 7.88 kg (200 MPa) for 1000 °C and 3.94 kg (100 MPa) for 1300 °C. Thickness for the samples ranges from 0.42 to 0.43 mm. The minimum creep rates are 2.54 × 10–4, 2.83 × 10–4, 1.54 × 10–5, 3.00 × 10–5 h−1 for H282 at 1000 °C, TZM at 1000 °C, AL-A6 at 1000 °C, and AL-A6 at 1300 °C, respectively. The inflection at ~70 h for TZM reflects a partial crack in the sample.

AL-A6 exhibits a good combination of RT ductility and HT strength, but its oxidization resistance is poor, which is typical for Mo-based alloys. To address this concern when deploying the alloy for application, we developed a multilayer RMPEA-Si-B coating (including a Mo precoating and a Si–B pack cementation)75,76. The coating is chemically compatible with AL-A6 and provides significant improvement in oxidation resistance compared to the uncoated sample. At both isothermal and cyclic conditions at 1300 °C, the mass changes for the coated samples are 2.7 mg cm2 after 50 h in thermal condition and 1 mg cm−2 after 400 cycles in cyclic condition. The improved oxidation resistance by the coating is provided by its multilayered structure, where the top aluminoborosilicate acts as an oxygen barrier, a disilicide layer as source for silica formation, a boron-doped (metal)5Si3 T1 layer for oxidation resistance, and a (metal)5Si2B2 T2 and boride layer acts as a diffusion barrier and provides boron for T1 layer generation and for long term stability (see details in refs. 75,77).

In summary, we showcased accelerated (~1 year) discovery of RMPEAs with an improved combination of RT and HT mechanical properties by combining DFT-guided down-selection with rapid alloy synthesis and characterization. For rapid screening, we used four simple criteria, including the use of a universal metric for a priori down-selection of alloy compositions with higher strength, and identified Mo72.3W12.8Ta10.0Ti2.5Zr2.5 with well-balanced compressive, tensile, and creep properties up to 1300 °C. In combination with a borosilicate (self-healing) coating, this alloy is also resistant to oxidation at 1300 °C.

Methods

Calculation of mechanical properties

The reduced model for the zero-temperature shear yield-stress developed by Maresca and Curtin78 is represented by an analytical equation governed by the stiffness constants (Cij’s):

$${\tau }_{{\rm {y}}0}=0.04{\alpha }^{-\frac{1}{3}}\bar{\mu }{\left(\frac{1+\bar{\nu }}{1-\bar{\nu }}\right)}^{\frac{4}{3}}\left[\frac{{\sum }_{n}{c}_{n}\Delta {V}_{n}^{2}}{{\bar{b}}^{6}}\right]$$
(1)

where \(\alpha\) is the edge dislocation line tension parameter, \(\bar{\mu }\) is the shear moduli, \(\bar{\nu }\) is Poisson’s ratio, \({c}_{n}\) is the concentration of type n solute, \(\Delta {V}_{n}\) is the misfit volume of type n solute, and \(\bar{b}\) is the Burger’s vector. The \(\bar{\mu }\), \(\bar{\nu }\), and \(\bar{B}\) (bulk moduli) are related to Cij’s by

$$\bar{\mu }=\sqrt{\frac{1}{2}{C}_{44}\left({C}_{11}-{C}_{12}\right)}$$
(2)
$$\bar{B}=\left({C}_{11}+2{C}_{12}\right)$$
(3)
$$\bar{\nu }=\frac{3\bar{B}-2\bar{\mu }}{2(3\bar{B}+\bar{\mu })}$$
(4)

The DFT-calculated stiffness constants, elemental elastic parameters and alloy volumes were used to estimate the Cij’s. The misfit volume was estimated from the rule of mixture. The model was fitted with the rule-of-mixtures stiffness constants, and this fitting parameter is absorbed by the 0.04 value making it just a scaling factor.

At finite temperature and finite strain rate standard thermal activation theory then leads to the predicted tensile yield stress as: (with k Boltzmann’s constant, T (in Kelvin)

$${{\sigma }_{{\rm {y}}0}=3.0{\tau }}_{{\rm{y}}{0}}\left[1-{\left(\frac{{kT}}{\Delta {E}_{{\rm {b}}}}\mathrm{ln}\frac{\dot{{\epsilon }_{0}}}{\dot{\epsilon }}\right)}^{2/3}\right]$$
(5)

Here, 3.06 is the Taylor factor for isotropic bcc polycrystal strength controlled by edge dislocations, \(\dot{{\epsilon }_{0}}\) is a reference strain rate. The quantities \(\bar{\mu }\) and \(\bar{\nu }\) are the shear modulus and Poisson’s ratio of the alloy, and 0.0833 (bcc) is related to the edge dislocation line tension as \(\Gamma\) = \(\alpha \bar{\mu }{b}^{2}\). The experimental Poisson’s ratio was used in strength calculation.

Combinatorial synthesis and characterization

A set of 32 samples (1.5 g each) were synthesized with a generalized formula of (Mo1−αWα) × Tay(TiZr)zAlβCrγ using the combinatorial arc-melt system (SP-MAM32, MTI Corporation). High-purity elemental raw materials (at least 99.9% pure) were acquired from the Materials Preparation Center at Ames Laboratory. The starting materials were in the form of foil or wire for convenient handling. Synthesized samples were machined into a cylinder shape ~3 mm in height with parallel top and bottom surfaces using an electrical discharge machining (EDM). Elastic moduli were measured using the ultrasound pulse-echo technique. Modulus was calculated by from:

$${\rm{Poisson}}{{\mbox{'}}}{\rm{s}}\; {\rm{ratio}}\,v=\frac{1-2{({V}_{{\rm {T}}}/{V}_{{\rm {T}}})}^{2}}{2-2{({V}_{{\rm {T}}}/{V}_{{\rm {T}}})}^{2}}$$
(6)
$${\rm{Young}}{{\mbox{'}}}{\rm{s}}\; {\rm{modulus}}\,E=\frac{{V}_{{\rm {L}}}^{2}\rho (1+v)(1-2v)}{1-v}$$
(7)

where \({V}_{{\rm {T}}}\) is the shear (transverse) velocity, \({V}_{{\rm {L}}}\) is the longitudinal velocity, \(\rho\) is the density, and \(v\) is the Poisson’s ratio of the samples. The densities of the samples were measured by Archimedes’ method79. To compute for shear and longitudinal velocities (\({V}_{{\rm {T}}}\) and \({V}_{{\rm {L}}}\)), the respective round-trip transit time values were measured by a precision thickness gauge (38DL PLUS, Olympus) with shear wave (V156, Olympus) and longitudinal wave (M112, Olympus) contact transducer. The velocity was calculated by dividing sample thickness by one-half of the round-trip transit time. The compositions of samples were measured by energy-dispersive spectroscopy (EDS) on a scanning electron microscope (SEM) (Teneo, FEI Inc). X-ray diffraction patterns were acquired using a Bruker DaVinci D8 diffractometer with a Cu target equipped with an autosampler.

Standard-size sample preparation, testing, and characterization

Down-selected standard-sized samples (10–20 g) were synthesized by arc melting of elemental pure metal sheets (at least 99.9% pure) acquired from the Materials Preparation Center at Ames Laboratory. Samples were flipped five times during the process to ensure homogeneity. Cylinders (2 mm in diameter and 6mm in length) and disks (8 mm in diameter and 0.5 mm in thickness) were EDM machined from arc-melted buttons for RT compression tests and HT small punch tests, respectively. The RT compression tests were conducted on an Instron 5569 testing frame with a 0.2 mm min−1 loading rate following a modified ASTM E9-19 standard80 due to the small sample size. The small punch tests were conducted on a compact small punch rig (OTF-1500X-S-ST, MTI Corporation) co-designed by the authors and MTI corporation following European standards on small-punch testing of metallic materials81 with 0.5 mm/min loading rate in flowing Ar atmosphere. The punch die set is made of SiC (OTF-1500X-S-CR-SIC, MTI Corp.), the punch rod is made of Mo (OTF-1500X-S-CR-MO, MTI Corp.), and the half-ball that transfers the force to the sample is made from ZrO2 with a diameter of 2.5 mm. Yield point was determined from load–displacement curve via an empirical relationship:

$${{S}}_{{{\rm {y}}}}={\rm{\alpha }}\frac{{{P}}_{{\rm {y}}}}{{t}^{2}}$$
(8)

where Sy is the tensile yield, Py is the load at yield, t is the sample thickness, and α is a geometry-dependent correlation coefficient (0.479 in our case)74,81. For the small punch creep tests (SPC), the same mold and sample geometry was used, the load was applied by a dead load system when the target temperature has reached, and the creep stress was evaluated using Eq. (8). The SPC tests were conducted in a flowing Ar atmosphere for up to 100 h. The minimum creep rate was calculated using an empirical equation:

$${\dot{\varepsilon }}_{\min }=0.3922\,{\dot{v}}_{\min }^{1.191}\,\left[1/{\rm {h}}\right]$$
(9)

where \({\dot{\varepsilon }}_{\min }\) is the minimum creep rate, \({\dot{v}}_{\min }\) is the minimum creep deformation rate as determined from plotting the second derivative of the creep deformation curve.

The strain in the SPT and SPC plots were calculated assuming the stain accumulates in the radial direction of the sample via the following equation:

$$\varepsilon =\frac{\Delta l}{{l}_{0}+\Delta l}=\frac{{{\rm {sqrt}}}\left({0.5D}^{2}+{v}^{2}\right)-0.5D}{0.5D+{{\rm {sqrt}}}\left({0.5D}^{2}+{v}^{2}\right)-0.5D}\times 100 \%$$
(10)

where \(\varepsilon\), \(\Delta l\), and \({l}_{0}\) are the strain, change in gauge length, and initial gauge length, respectively. D is the diameter of the lower punch opening, and v is the displacement from the LVDT.

The HT compression test was conducted in a vacuum of 10−5 Torr or better with an initial strain rate of 10−3 s−1 using an Instron mechanical testing machine fitted with a vacuum furnace. Cylindrical specimens with a diameter of ~6.5 mm and a height of ~11mm were used for the test. The cross-sectional and fractured-surface imaging was done using a scanning electron microscope (SEM) (Teneo, FEI Inc.) equipped with an energy-dispersive X-ray spectroscopy (EDS) detector.

Density-functional theory (DFT)

The all-electron, DFT-based Korringa–Kohn–Rostoker (KKR) Greens’ function method combined with the coherent potential approximation (CPA) was used to calculate total energy, lattice parameters, universal strength metric, and moduli of arbitrary solid-solution alloys47. The KKR-CPA performs averaging overall all possible configurations simultaneously with DFT charge self-consistency, which properly includes alloy-induced Friedel impurity-charge screening and core-level shifts due to alloying. As such only small symmetry-equivalent unit cells are required, e.g., 1-atom (2-atom) cells for A1/A2 (A3). Hence, thousands of alloy compositions can be performed in an HTP manner. For defected cells (e.g., stacking faults) larger unit cells are required but the configurational averages are fully performed over the full cell. We used the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional for solids82. We employed a site-centered, spherical-harmonic basis with s, p, d, and f-orbital symmetries (i.e., Lmax = 3) in all calculations. The self-consistent charge density was obtained from Green’s function using a complex-energy contour integration and Gauss-Laguerre quadrature (with 24 energy points on a semi-circular mesh enclosing the bottom to the top of the valence states). An equally spaced Monkhorst–Pack k-space mesh of 24 × 24 × 24 for bcc/fcc was used for Brillouin zone integrations. Core electrons were treated fully relativistically (which includes spin–orbit coupling), while semi-core/valence electrons were treated scalar relativistically (neglecting spin–orbit coupling).