Ductility mechanisms in complex concentrated refractory alloys from atomistic fracture simulations

We begin by considering differences between elemental Nb, MoNbTaW, and $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ systems based on the Rice theory [10], which predicts a crossover condition, under pure mode I loading, for cleavage bond breaking (brittle behavior) and dislocation emission (ductile behavior) at an advancing crack tip. Specifically intrinsically brittle (ductile) behavior is predicted if the ratio ($\gamma_{\text{\text{USF}}}/\gamma_{\text{surface}}$) of USF energy ($\gamma_{\text{\text{USF}}}$) to surface energy ($\gamma_{\text{surface}}$) is larger (smaller) than a threshold value defined as:

where $\nu$ is the Poisson’s ratio, $\theta$ is the inclination angle between the fracture plane and the dislocation glide plane, and $\phi$ is the inclination angle between the slip direction with the normal to the crack tip.

For both $\{110\}$ and $\{112\}$ fracture planes, Table I presents a comparison of values for $\gamma_{\text{\text{USF}}}/{\gamma_{\text{surface}}}$, and the estimated threshold values $\left[\gamma_{\text{\text{USF}}}/{\gamma_{\text{surface}}}\right]_{\textrm{threshold}}$ calculated from the MLIP due to Yin et al. [22] for NbMoTaW and Nb, and the MLIP developed here for $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ (see Methods section). In calculating these values we derive $\nu$ from the calculated single-crystal elastic constants using Voigt-Reuss-Hill (VRH) averaging, and consider all possible $\{110\}$ and $\{112\}$ slip planes (see Table SIV), corresponding to those with a value of $\theta$ and $\phi$ that maximizes $\left[\gamma_{\text{\text{USF}}}/{\gamma_{\text{surface}}}\right]_{\textrm{threshold}}$. Further, calculated $\gamma_{\text{\text{USF}}}$ results represent average values that do not account for local variations due to compositional fluctuations in the RHEAs. Due to these approximations we do not expect the properties in Table I to give a rigorous threshold condition for brittle versus ductile behavior, but the relative values provide a guide to what is expected in the atomistic fracture simulations presented in the next sub-section. Specifically, the Rice criterion with the values in Table I predicts brittle cleavage bond breaking fracture behavior for NbMoTaW and Nb. For $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$, $\gamma_{\text{\text{USF}}}/{\gamma_{\text{surface}}}$ are close to threshold values, and ductile dislocation may be expected. Additionally, the numbers in Table I suggest that the effect of adding group VI elements (Mo and W) to Nb is to raise the tendency toward brittle fracture, while the addition of group IV elements (Ti and Hf) is to increase ductility, consistent with trends reported previously from first-principles calculations [11, 23, 24].

We consider next direct atomistic simulations of crack-tip propagation using the $K$-test geometry with a semi-infinite crack [15] (see Methods section). Simulation results are presented for four crack-tip orientations, corresponding to three crack line directions ($\langle 1\bar{1}0\rangle$, $\langle 1\bar{1}1\rangle$, and $\langle 1\bar{1}2\rangle$) with a $\{110\}$ fracture plane, and one crack line direction ($\langle 1\bar{1}0\rangle$) for a $\{112\}$ fracture plane.

Figure 1 contrasts simulation results for loading of $\{110\}\langle 1\bar{1}2\rangle$ oriented cracks for pure Nb ((a) and (b), and the RHEAs MoNbTaW ((b) and (c)) and $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ ((d) and (e)). For pure Nb, the crack is observed to propagate by bond breaking, indicating brittle fracture behavior, consistent with the Rice theory predictions from Table I. For NbMoTaW crack propagation for the first 50 $\AA$ proceeds through a mixed mode of bond breaking and dislocation emissions, but switches to bond breaking at later stages. Thus, despite the results in Table I, indicating more brittle behavior for NbMoTaW relative to Nb, plasticity is observed in the initial stages of crack propagation for the former. This behavior may represent an effect of fluctuations associated with compositional complexity, giving rise to local environments with lower barriers for dislocation nucleation than would be expected from the average $\gamma_{\text{\text{USF}}}$ value.

Turning to the results for $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ in Fig. 1 (d) and (e), the initial crack propagation is associated with a dislocation emission event, and further crack advances lead to more dislocation emissions, significant blunting of the crack tip and twin formation. These observations reflect a ductile fracture mode, correlating with the near threshold values for $\gamma_{\text{\text{USF}}}/{\gamma_{\text{surface}}}$ in Table I. For the remaining crack tip orientations we focus on comparisons of results for MoNbTaW and $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$.

Figures 2 and 3 show representative atomistic configurations for crack propagation in NbMoTaW and $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ for $\{110\}\langle 1\bar{1}0\rangle$ and $\{110\}\langle 1\bar{1}1\rangle$ oriented crack tips, respectively. Consistent with the Rice theory, the crack propagates via cleavage bond breaking and remains sharp in NbMoTaW for both of these crack orientations. For $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$, the results in Fig. 2 are similar to those in Fig. 1, with initial propagation accompanying the emission of dislocations (Fig. 2(c)) and eventual blunting of the crack tip. The behavior for $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ is somewhat more complex for $\{110\}\langle 1\bar{1}1\rangle$ oriented crack tips, with fewer dislocation emissions observed with a lower degree of crack blunting and cleavage bond breaking observed in the later stages.

We consider next the propagation of cracks normal to the $\{112\}$ surface, as shown in Fig. 4. For this $\{112\}\langle 1\bar{1}0\rangle$ crack orientation, MoNbTaW cracks no longer propagate in an atomically sharp manner, and display dislocation emission. Similar to the initial stages of crack propagation in this system for the $\{110\}\langle 1\bar{1}2\rangle$ crack orientation, this behavior for MoNbTaW is in contrast to the picture from the Rice theory, and may reflect the effects of compositional fluctuations arising from chemical complexity. This dislocation emission does not lead to a high degree of crack blunting, however, and the crack propagation switches to predominantly bond breaking in later stages. For $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ the behavior for $\{112\}$ fracture planes is similar to that observed for fracture on $\{110\}$ planes, with crack propagation accompanying dislocation emission, twin formation, and crack blunting.

To further analyze the evolution of crack behavior once it begins to advance, we analyze the crack tip opening displacement (CTOD) as a function of the crack extension $\Delta$a, reflecting a nanoscale version of the R-curve used in fracture mechanics. In this analysis the CTOD is derived as distance between the opposite surfaces of the crack, the shape of which is approximated by that of a hemisphere (see SI. Fig. S1-8 for details). Figure 5 shows the atomistic crack-resistance R-curve, plotting CTOD versus $\Delta$a. The CTOD in $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ is higher than or approximately the same as that of NbMoTaW for all $\Delta$a values, for all crack orientations considered. These results thus show good qualitative correspondence with the predictions of the Rice theory (Table I), suggesting more ductile behavior for the former RHEA. For all orientations considered, the crack tips in $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ involve a more irregular crack propagation path and a larger CTOD as the loading proceeds. Those results indicate that a higher degree of crack tip blunting and fracture resistance for $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$, also in qualitative agreement with the Rice theory.

A notable finding from the results presented in the previous sub-section is that the mechanisms for crack-tip propagation (plasticity versus bond breaking) and geometry (blunted versus sharp) can differ between initial and later stages of crack propagation, starting from an atomically sharp crack tip. We consider next the origins of these changes, facilitated by an analysis of the strain energy release ($G_{I}$) as a function of crack tip extension ($\Delta$a). A plot of $G_{I}$ versus $\Delta$a represents a nanoscale version of the $G_{R}$ curve considered in fracture mechanics, with higher slopes in this relationship corresponding to higher fracture resistance (i.e., more ductile behavior). To compute $G_{I}$ we make use of the following relationship with the stress intensity factor ($K_{I}$) under plane strain condition:

where $\nu$ is the Poisson’s ratio and $E$ is the Young’s modulus. As above, we determine values for the isotropic elastic moduli in Eq. 2 via the VRH averages derived from the calculated single-crystal elastic constants. The resulting $G_{R}$ curves are plotted in Fig. 6.

We consider first the results for pure Nb, illustrated for a $\{110\}\langle 1\bar{1}2\rangle$ crack orientation in Fig. 6(a). The constant slope reflects the same operating fracture mode, namely cleavage fracture, over the entire range of $\Delta$a considered. We note that these results for pure Nb were obtained using the MLIP due to Yin et al. [22], although a similar curve with a constant slope having a magnitude within 22 percent of that shown in Fig. 6 is obtained for Nb using the MLIP developed in the current work, see Supplemental Fig.S9. This result indicates that the differences between Nb and $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ in Fig. 6 (discussed below) are not an artifact of using two separate MLIP models in the fracture simulations, but rather reflect the role of alloy chemistry.

For $\{110\}\langle 1\bar{1}0\rangle$ and $\{110\}\langle 1\bar{1}1\rangle$ crack orientations, the results for MoNbTaW in Fig. 6 (b), (c) show behavior similar to pure Nb, with a constant slope reflecting cleavage fracture for all values of $\Delta$a plotted. By contrast, results for MoNbTaW in Fig. 6 (a), (d) for $\{110\}\langle 1\bar{1}2\rangle$ and $\{112\}\langle 1\bar{1}0\rangle$ orientations, as well as $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ for all orientations, show changes in slope reflecting different fracture modes for initial and later stages of crack propagation.

Considering the results for $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$, we find that although the crack tip propagates predominately by dislocation activity, the changes in slopes in Fig. 6 indicate different modes of plasticity. From an analysis of the atomic configurations before and after the change in slopes, we identify two dislocation-assisted crack growth modes: nucleation-limited extension and slip-limited extension. In the nucleation-limited mode, a large change in the loading stress intensity factor ($\Delta K_{I}$) is required to nucleate a dislocation at the crack tip, but once it is formed, the dislocation propagates into the bulk relatively easily. In the slip-limited mode, abundant dislocation sources are active at the crack tip, and the nucleated dislocations interact with already emitted dislocations, causing resistance to slip. It is noted that in this slip-limited mode, the emission of a dislocation into the bulk crystal is sometimes accompanied by the annihilation of others. The changes in the dominant mode of ductile fracture are captured by the sharp changes in the slope in the $G_{R}$ curves for $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$. With this insight, we detail below the changes in the fracture modes observed in the simulations.

For the $\{110\}\langle 1\bar{1}2\rangle$ (Fig. 6(a)) and $\{112\}\langle 1\bar{1}0\rangle$ (Fig. 6(d)) crack orientations, both RHEAs exhibit two distinct stages of crack extension. In the initial stage for NbMoTaW, for both crack orientations, the crack extends by a mixture of dislocation nucleation and slip limited crack extension, while in the later stage a bond-breaking cleavage mode prevails. For $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$, for both orientations, crack extension is initially dislocation-nucleation limited, and switches to slip-limited in later stages. For the nucleation limited regime, the average changes in the loading stress intensity factor are $\Delta K_{I}\approx 0.13\text{MPa}\sqrt{\text{m}}$ and $\Delta K_{I}\approx 0.12\text{MPa}\sqrt{\text{m}}$ to nucleate a dislocation for $\{110\}\langle 1\bar{1}2\rangle$ and $\{112\}\langle 1\bar{1}0\rangle$, respectively; subsequent dislocation emission occurs at intervals of $\Delta K_{I}\approx 0.05\text{MPa}\sqrt{\text{m}}$ and $\Delta K_{I}\approx 0.04\text{MPa}\sqrt{\text{m}}$, respectively. The later crack extension stage for $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ is dislocation-slip limited, requiring on average $\Delta K_{I}\approx 0.18\text{MPa}\sqrt{\text{m}}$ and $\Delta K_{I}\approx 0.26\text{MPa}\sqrt{\text{m}}$ for $\{110\}\langle 1\bar{1}2\rangle$ and $\{112\}\langle 1\bar{1}0\rangle$, respectively, to advance an already existing dislocation, and one or more dislocations are nucleated almost immediately after such an event; in the interval between dislocation emissions, the number, total length, and the shape of the dislocations change as a result of their interaction. In the case of the $\{112\}\langle 1\bar{1}0\rangle$ crack orientation, at a later stage a new grain is formed at the crack tip and the nearby grain boundary acts as a source for dislocations (see Fig. 4(f)).

Considering the $\{110\}\langle 1\bar{1}0\rangle$ (Fig. 6(b)), and $\{110\}\langle 1\bar{1}1\rangle$ (Fig. 6(c)) crack orientations, MoNbTaW displays a constant slope corresponding to a cleavage fracture mode. By contrast, $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ again displays two crack-extension stages. For the $\{110\}\langle 1\bar{1}0\rangle$ orientation, crack extension is initially dislocation-nucleation limited, and switches to slip-limited in later stages, similar to the orientations discussed above. For the nucleation limited regime, the average loading is $\Delta K_{I}\approx 0.11\text{MPa}\sqrt{\text{m}}$ to nucleate a dislocation. Subsequent dislocation emission occurs at intervals of $\Delta K_{I}\approx 0.13\text{MPa}\sqrt{\text{m}}$. The later crack extension stage for $\text{Nb}_{45}\text{Ta}_{25}\text{Ti}_{15}\text{Hf}_{15}$ is dislocation-slip limited: on average, requiring $\Delta K_{I}\approx 0.17\text{MPa}\sqrt{\text{m}}$ to advance an already existing dislocation. For the $\{110\}\langle 1\bar{1}1\rangle$ orientation, the initial stage (lower slope) corresponds to a slip-limited regime in which three dislocation emission events are observed. In between these emission events, multiple dislocations interact, causing resistance to slip, and a higher slope in Fig. 6(c). At later stages, the crack-growth mode switches to cleavage (bond-breaking), with a correspondingly lower slope in Fig. 6(c).

For the atomistic simulations we employed a previously published MLIP for NbMoTaW and Nb [22] based on the moment tensor potential (MTP) formalism [27]. For NbTaTiHf, we have developed a new MLIP based on the atomic cluster expansion (ACE) formalism [28, 29, 30, 31]. The ACE potential is trained on a large database of energies and forces derived from DFT calculations performed using the Vienna Ab-Initio Simulation Package (VASP) [32, 33]. The initial dataset contains energies and forces for the configurations (beside those corresponding to the ab-initio molecular dynamics simulations) used in the training of the MoNbTaW MLIP [34], as well as some configurations taken from the first-principles study of Borges et al. [35]. An initial ACE MLIP for NbTaTiHf was trained on these configurations and used to initiate an iterative refinement of the potential, employing an active learning process in which configurations with high values of the so-called extrapolation grade [31] were identified from finite-temperature molecular-dynamics (MD) simulations under periodic boundary conditions at elevated temperatures (1000K, 1600K, 2000K, and 2400K) and free surfaces at room temperature. These MD simulations and extrapolation-grade analysis were performed using the performant atomic cluster expansion (PACE) implementation [29] in the Large- scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software [36]. The configurations with high extrapolation grade were extracted for DFT calculations and added to the training set to refine the potential further. This iterative process continued until finite-temperature MD simulations produced no configurations with an extrapolation grade higher than 5. For the final version of the potential, we considered a total of 6493 structures, and 263,810 atoms. The MLIP was trained on 90% of the data, with the remaining used as a test set using the pacemaker software [29]. The results of the training are summarized in Fig. 7, showing parity plots for energy (a) and forces (b), for both training and test sets. Further details of the DFT calculations and MLIP training are provided in Table SI-III in the Supplemental Information (SI), where additional results are presented for properties derived from the potential, including elastic constants, generalized stacking faults, local lattice distortion, etc.

The MLIPs were used in fracture simulations employing the $K$-test geometry with a semi-infinite crack [15]. In these simulations, the response of a crack tip is analyzed for incrementally increasing values of an applied stress intensity factor $K_{I}$, which sets the boundary conditions far away from the crack tip. The simulation cells are approximately 55 $\AA$ thick along the crack line direction, 300 $\AA$ long in the direction of crack propagation, and 300 nm along the direction normal to the fracture plane. The simulation cells contain from 278,000 to 314,000 atoms. Results highlighted in the main text are found to be relatively insensitive to further increases in simulation size. As described in the text, four different orientations are considered for the fracture plane and crack line direction. These simulation geometries are chosen to align the crack surface along the bcc planes ($\{110\}$ and $\{112\}$). Periodic boundary conditions are used in the line direction of the crack tip with free surfaces employed normal to the other two directions. The crack is introduced by removing a layer of atoms to mimic the slightly blunted crack and to generate a traction free surface. The simulation cell is loaded via the anisotropic elastic displacement field solution of linear elastic fracture mechanics (LEFM) [37] under plane strain conditions, followed by the static relaxation of the atoms in the interior of the simulation cell; in these relaxations the atoms in the exterior (2$r_{c}$ from the boundary) are kept frozen at positions dictated by the LEFM solutions. The simulations were repeated multiple times, incrementing the loading by $\Delta K=0.01\text{MPa}\sqrt{\text{m}}$ in each subsequent iteration. All static relaxations were performed in LAMMPS using the “FIRE” algorithm [38] with energy tolerance of $1\times 10^{-9}$ and force tolerance of $1\times 10^{-3}$ eV/$\AA$.

This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, under Contract No. DE-AC02-05-CH11231 within the Damage Tolerance in Structural Materials (KC 13) program.

The study made use of resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under the same contract number, using NERSC Award No. BES-ERCAP0027535.

This research used the Lawrencium computational cluster resource provided by the IT Division at the Lawrence Berkeley National Laboratory (Supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231)

The Scientific Computing Group (also known as High-Performance Computing Services) under Science IT supports the mission of Lawrence Berkeley National Laboratory by providing technology and consulting support for science and technical programs, in the areas of data management, HPC cluster computing, and Cloud services.