Phase-field modelling of transformation pathways and microstructural evolution in multi-principal element alloys

Multi-principal element alloys (MPEAs), generally known as high entropy alloys (HEAs), have gained much attention in the past two decades owing to their potentially better functional, electrochemical, and mechanical properties.^1–3 Since single phase MPEAs are not suitable for high temperature structural applications, a significant amount of recent efforts has been devoted to the development of multiphase MPEAs.^4–9 In this study, we focus on a two-phase MPEA, AlMo_0.5NbTa_0.5TiZr, that has some exceptional physical and mechanical properties superior to Ni-based superalloys such as low density and high high-temperature yield strength.¹⁰ The microstructure of the alloy was initially characterized to have two β (bcc) phases.¹¹ However, Jensen et al.¹⁰ found that the matrix (topologically continuous) phase is actually an ordered β′ (B2) phase and the β precipitate phase forms a self-organized, periodic array of cuboidal particles aligned along the ⟨001⟩ directions. At the solutionizing temperature (1400 °C), only the $β$ (bcc) phase is present in the microstructure.¹² This microstructure is uncommon because most precipitation-hardened alloys have disordered solution phases as matrix phases and ordered intermetallic phases as precipitate phases. For example, in Ni-base and Co-base superalloys, the disordered (⁠$γ$⁠, fcc) phase is the matrix phase while the ordered phase (⁠$γ′$⁠, L1₂) is the precipitate phase. The reason why we expect the matrix phase to be disordered is that at the solutionizing temperature, a single disordered phase is present and the ordered phase forms upon cooling. It is also interesting to note that the unique microstructural features observed in AlMo_0.5NbTa_0.5TiZr, including (i) narrow particle size distribution and (ii) highly regular and aligned cuboidal particles arrays, resemble well those of the $γ$ + $γ′$ microstructure in Ni-base and Co-base superalloys. The formation of the $γ$ + $γ′$ microstructures has been explained and modeled via the conventional nucleation and growth mechanism, i.e., the $γ′$ phase precipitates out of a supersaturated parent $γ$ phase solid solution. It should be noted that under specific heat treatment conditions (such as rapid-cooling), it is possible to obtain $γ+γ′$ two-phase microstructures through a spinodal-mediated phase separation pathway.^13,14 However, what lead to the “inverted superalloy-like” β′ + β microstructure observed in AlMo_0.5NbTa_0.5TiZr is still unclear. For a similar alloy system (Al_0.5NbTa_0.8Ti_1.5V_0.2Zr), two possible PTPs mediated through spinodal decomposition¹⁵ have been proposed, which could lead to the $β′+β$ microstructure. In order to tailor and optimize the microstructure for specific engineering applications, it becomes critical to understand the PTP leading to this β′ + β microstructure. In general, materials with an ordered phase as the continuous matrix phase may not have the well-balanced properties that the superalloys can offer and, thus, there is a need to invert the β′ (matrix) + β (precipitate) microstructure into a β (matrix) + β′ (precipitate) microstructure by adjusting the alloy composition and heat treatment schedule. Experimental efforts are on-going in this area to tailor the β′ + β microstructures in MPEAs.^9,13–18

In this work, using a prototype model system and phase-field simulations, we demonstrate a possible PTP leading to the β′ + β microstructure and study effects of the equilibrium volume fractions of the β′ and β phases and their elastic properties on inverting the β′ + β microstructure into a β + β′ one. The proposed PTP consists of (i) congruent ordering of β into β′ upon cooling, (ii) spinodal decomposition of β′ into solute-lean and solute-rich B2 phases, and (iii) congruent disordering of the solute-lean β′ phase into β, leading to a β′ + β microstructure.

Our prototype binary system has two phases: (i) a disordered phase that represents the bcc (⁠$β$⁠) phase and (ii) an ordered phase that represents the B2 (⁠$β′$⁠) phase. The free energy curves of the B2 (⁠$fB2$⁠) and bcc (⁠$fbcc$⁠) phases as function of solute concentration $(c)$ of the prototype model binary system at a temperature within the BCC+B2 field of the phase diagram are shown in Fig. 1. Both the free energies (⁠$fB2, fbcc$⁠) and solution concentration (⁠$c$⁠) are presented in non-dimensional forms. The non-dimensional solute concentration (⁠$c$⁠) can be related to real solute concentration (⁠$x$⁠) using $c=(x−xbcc)$/$(xB2−xbcc)$ where $xbcc$ and $xB2$ are the equilibrium solute concentrations of the bcc and B2 phases, respectively. The B2 phase has a miscibility gap and the spinodal region within which, $∂2fB2/∂c2<$ 0, is in between the two cross symbols in Fig. 1(a). The dashed line in Fig. 1(a) represents the absolute instability of the bcc phase with respect to spontaneous ordering, i.e., a disorder–order transition without an energy barrier. For both alloy 1 and alloy 2 indicated in Fig. 1(a), spinodal decomposition in the B2 phase will take place following congruent ordering of the starting bcc phase. This kind of spinodal decomposition is called conditional spinodal^19,20 and is known to occur in the Al-Fe system.²¹ The equilibrium microstructures in both alloys are mixtures of B2 and bcc phases with different volume fractions of these phases. We chose the solute atoms of our prototype system to represent $β′$ forming elements, such as Al and Zr. Hence, the equilibrium B2 phase is solute-rich (⁠$c=1)$ and the equilibrium bcc phase is solute-lean (⁠$c=0$⁠).

The phase-field model consists of two order parameters: (i) concentration $c(r,t)$ and (ii) long-range order (lro) parameter $η(r,t)$ that distinguishes the disordered bcc phase (⁠$η=0$⁠) from the ordered B2 phase (⁠$η=1$⁠). The total free energy of the system $F[cr,t,η(r,t)]$ is formulated as a functional of these order parameters,

where $f(c,η)$ is the local chemical free energy density, $κc$ is the gradient energy coefficient for concentration, $κη$ is the gradient energy coefficient for the lro parameter, and $Eel[cr,t,η(r,t)]$ is the elastic strain energy functional of the system. The local free density $f(c,η)$ is approximated by

where $hη=η2(3−2η)$ is an interpolation function, $gη=η21−η2$ is a double-well function and $ω$ is the barrier height between the ordered and disordered phases. A contour plot of the free energy surface, $f(c,η)$⁠, is shown in Fig. 1(b) and its projection onto the f-c plane is shown in Fig. 1(a). The color bar in Fig. 1(b) represents the magnitude of the local free energy density, $f(c,η)$⁠.

For simplicity, we assume that the lattice misfit is only a functional of concentration following the Vegard law.^22–24 The elastic energy functional is given by

where $Cijkl$ is the elastic modulus field, $ϵklT$ is the transformation strain field, and $ϵij$ is the total strain field. Both the $Cijkl$ and $ϵklT$ are functions of concentration $cr,t,$

where $ϵoo=1adadc, a$ is the lattice parameter, $Cijklbcc$ is the modulus of the disordered phase, and $CijklB2$ is the modulus of the ordered phase. The total strain field $ϵr$ is determined by the mechanical equilibrium equation: $∇.σ=0$⁠. When the modulus of the system is homogeneous (⁠$Cijklbcc=CijklB2$⁠), the elastic energy $Eel$ can be directly written as an explicit functional of the order parameters as derived by Khachaturyan et al.^25,26 using the Green's function solution in the reciprocal space. This closed form expression of $Eel$ greatly increases the numerical efficiency of the model. When the modulus of the system is inhomogeneous, an iterative procedure^27,28 is used to compute the strain field.

The time-evolution of the concentration field, $c(r,t)$⁠, and the lro order parameter field, $η(r,t)$⁠, follows the Cahn–Hilliard equation²⁹ and Allen–Cahn equation³⁰ (also known as time-dependent Ginsburg–Landau equation), respectively,

where $M$ is the chemical mobility, L is the kinetic coefficient characterizing the order-disorder transition, and $ζc$ and $ζη$ are the Langevin-noise terms for concentration and lro parameter, respectively. The Langevin-noise terms are added to the governing equations in order to incorporate thermal fluctuations in our model.^31,32 The governing equations [Eqs. (6) and (7)] are solved using the semi-implicit spectral method.³³ When the disordered or ordered phase is metastable, the congruent disorder ↔ order transition occurs by nucleation, which is simulated by the Langevin noise term $ζη$ in the Allen–Cahn [Eq. (7)] equation.^34,35 For simplicity, the Langevin noise for concentration is chosen to be zero, $ζc=0$⁠. All the numerical values of the model parameters used in the simulation are listed in Table I. All simulations are performed in two dimensions with a system size of $512lo×512lo$ where $lo$ is the grid size. The length scale in the phase field model is determined by the ratio of the interfacial energy to the elastic energy. Assuming an interfacial energy of $γ=50 mJ$/$m2$ for the coherent interface between the B2 and bcc phases and a typical shear modulus (⁠$C44$⁠) of 182 GPa for the bcc phase, the grid size $lo$ is determined to be $4.2 nm$⁠.

The simulation results obtained for alloy 1 are shown in Fig. 2, along with the PTP marked on the free-energy contour. In this case, we assumed homogeneous modulus (⁠$Cijklbcc=CijklB2$⁠) and a lattice misfit of $0.2%$ between the B2 and bcc phases. The evolution of the system in Fig. 2 is represented through three columns where the first column shows the concentration field $cr,t$⁠, the second column shows the schematic PTP on the free energy surface where the current state of the system is marked by the white open squares, and the third column shows the order parameter field $η(r,t)$⁠. The figure shows eight snapshots of the microstructural evolution process marked as (1–8) below the PTP in the second column. The initial condition of the simulation [snapshot (1)] is a homogeneous bcc phase corresponding to the state of the system at the homogenization temperature before cooling.

The snapshot (2) shows the congruently ordered (B2) state formed during isothermal aging at a lower temperature, where $ηr,t=1$ everywhere in the system without any concentration modulations. Since the alloy composition is located in the absolutely unstable region [shown as blue dashed line in Fig. 1(a)] for alloy 1, the congruent ordering is a barrierless transformation and the transformation can proceed without the assistance of the Langevin noise term [i.e., $ζηr,t=0$]. For alloy 2, the alloy composition is located in the metastable region [shown as blue solid line in Fig. 1(a)] and then the Langevin noise term $ζη(r,t)$⁠, is required to initiate the congruent ordering reaction. However, in both cases, the microstructure will reach the congruently ordered (B2) single-phase state first before spinodal decomposition. Note that the B2 phase has two anti-phase domains and can form anti-phase domain boundaries (APBs) during the congruent ordering. We assume that the anti-phase domain size is much larger than our simulated system size and ignore the presence of APBs. This assumption is justified from the experiments where the size of the anti-phase domains is much larger than the size of the disordered precipitates.¹⁰ 

The snapshots (3) and (4) capture the spinodal decomposition process in the B2 phase. The early stage of the spinodal decomposition shown in snapshot (3) resembles a basket-weave microstructure with two orthogonal concentration waves along the elastically soft directions (i.e., ⟨001⟩) superimposed on top of each other. After the concentration modulation develops completely, we observe ordered solute-lean domains embedded in an ordered solute-rich matrix, aligned along the ⟨001⟩ elastically soft directions as per the $Cijkl$ chosen in our model. The solute-lean ordered region is metastable, which can be seen from the free energy surface [Fig. 1(b)] and undergo congruent disordering. The snapshots (5)–(8) show the congruent disordering process. In the solute-lean B2 region, the bcc particles are formed by nucleation and growth (assisted by thermal fluctuations in the lro parameter that is captured by the Langevin noise term in the simulations), as indicated by the blue patches (snapshot (5), the right column) in the $η(r,t)$ field. Finally, most of the solute-lean regions become disordered domains embedded in a continuous, solute-rich ordered B2 phase matrix, as shown in snapshot (8).

Thus, our proposed PTP leads to a microstructure where discrete bcc precipitates embedded in a continuous B2 matrix, which agrees with that observed in the experiment.¹⁰ The comparison between the microstructure of the annealed sample and the phase-field simulated microstructure is shown in Fig. B1 (supplementary material). We have also simulated the microstructural evolution in alloy 2 with parameter given in Table I. The microstructural evolution in Alloy 2 leads to the inverted microstructure $β$(matrix) + $β′$(precipitate) as shown in Fig. 3.

From the above simulation results (Figs. 2 and 3), we realize that whether the ordered or disordered phase forms the continuous matrix phase is determined by the spinodal decomposition step of the PTP [marked in Fig. 1(b)]. At the end of the congruent ordering step, the system is a homogeneous single B2 phase without any concentration modulations. In the congruent disordering step [snapshots (5)–(8) in both Figs. 2 and 3], only the solute-lean B2 regions can transform to the bcc phase and the morphology of the solute-lean and solute rich regions do not change. Hence, the morphologies (i.e., discrete, continuous, or bi-continuous) of the ordered and disordered phases are predetermined by the spinodal decomposition process.

A clear understanding of the phase transformation pathway (PTP) is essential to design the alloy composition and heat-treatment, which will help us to tailor the microstructure for specific engineering applications. Below, we show an example of tailoring the microstructure under our proposed PTP by modifying the alloy composition and elastic modulus.

The as-cooled sample from the experiments has discrete precipitates consisting of cuboidal particles and small square shaped platelets, well-aligned in the ⟨001⟩ direction, which are in contrast to the elongated plate-like (rod in 2D) precipitates observed in our simulation (Fig. 2). A more quantitative definition of the discreteness can be found in the supplementary material (Section A). This could be due to the modulus mismatch between the two phases. We further studied the effect of modulus mismatch and alloy composition (i.e., volume fractions of the two phases) on microstructural evolution during spinodal decomposition. It was shown in the literature that an elastically more compliant (i.e., soft) phase tends to wrap around an elastically stiffer (i.e., hard) phase in order to minimize the elastic energy during phase separation.^36,37 In a complex multicomponent system, such as our MPEA, we do not have reliable measurements of the elastic modulus $Cijkl$ of the individual phases, namely, the solute-lean B2 (which will congruently disorder into bcc phase in the final microstructure) and solute rich B2 phases, to identify the elastically soft phase. Hence, we consider two extreme cases: (i) case A: soft solute-rich phase $Relas=CijklB2$/$Cijklbcc$=0.6; (ii) case B: hard solute-rich phase $Relas=1.4$⁠. In order to obtain the desired $Relas$⁠, the modulus of the solute-lean B2 phase was kept constant (values given in Table I), while changing the modulus of the solute-rich B2 phase. We simulated the spinodal decomposition process in these two cases for both alloy 1 and alloy 2. Note that the equilibrium volume fraction of the solute-rich phase in alloy 1 (∼60%) is higher than alloy 2 (∼40%). The simulated microstructures are shown in Fig. 4. Even though these microstructures are similar to those obtained by Li et al.,³⁸ the modulus mismatch considered in the current study (⁠$Relas=1.4 or 0.6$⁠) is much larger than that used in Li et al.²⁸ (⁠$Relas≈1.03 or 0.97$⁠). What we find is that such a large modulus mismatch is necessary to achieve significantly increased discreteness of the precipitate phase as compared to those obtained in the homogeneous modulus counterpart. We also performed a systematic study to quantify discreteness of the microstructure as a function of the inhomogeneity ratio $Relas$ (see the supplementary material, Fig. A1).

In case A, we observe a highly discrete solute-lean precipitates in alloy 1, as shown in Fig. 4(a-3). As the volume fraction of the solute-lean phase is low in this alloy, it becomes the precipitate phase (i.e., topologically discontinuous phase) and, because the matrix phase is elastically softer, the modulus mismatch further increased the discreteness of the precipitates. In alloy 2, where the volume fraction of the solute-lean phase is high, it becomes the matrix at the early stage of the decomposition process [Fig. 4(b-2)]. As the decomposition proceeds, the solute-rich phase becomes more topologically connected because it is elastically soft, leading to a nearly bi-continuous microstructure [Fig. 4(b-3)].

In case B, we observe a nearly bi-continuous microstructure in alloy 1 [Fig. 4(c-3)] and discrete solute-lean particles in a solute-rich matrix in alloy 2 [Fig. 4(d-3)]. The formation mechanisms of these microstructures can be explained with similar reasoning as that mentioned in case A. The microstructure observed in alloy 2 is interesting as it will give the desired inverted microstructure (β matrix + β′ precipitate) after congruent disordering of the solute-lean phase.

We have noted that the current simulations do not reproduce exactly the microstructures observed in the experiment (see, e.g., Fig. 4 of Ref. 10), i.e., the presence of both cube-shaped precipitates and plate-shaped satellite precipitates in between. Research is ongoing to identify the underlying mechanisms for the formation of such an interesting microstructure.

Recently, Li et al.³⁸ developed a phase field model to simulated the B2/bcc microstructures observed in HEAs. Even though the simulated microstructure resembles the microstructure observed in the experiments, they considered spinodal decomposition into a B2 and a bcc phase, which is not an isostructural transformation. Thus, the actual PTP in these alloys could be similar to what we have considered in the current study, i.e., congruent ordering of bcc into B2, followed by spinodal decomposition of the B2 phase into two B2 phases having different compositions, and one of them congruently disorders into bcc,²¹ or follow other non-conventional PTPs such as pseudospinodal^39,40 and mechanochemical spinodal⁴¹ that occur between phases having different crystal structures.

In summary, we have proposed and analyzed a phase transformation pathway (PTP) for phase separation in AlMo_0.5NbTa_0.5TiZr MPEA. Phase field simulations of microstructural evolutions along this PTP for different alloy compositions and different elastic modulus mismatch between the two co-existing phases were performed. They show that both the $β′$(matrix) + $β$ (precipitate) microstructure, which is observed in experiments and its inverted counterpart $β$ (matrix) + $β′$ (precipitate) can be obtained from the proposed PTP. They also show that the modulus mismatch between the phases can be used to increase the discreteness of the precipitate phase.

See the supplementary material for the effect of modulus inhomogeneity on the microstructure and the experimental microstructure.

The authors would like to acknowledge the financial support from Air Force Office of Scientific Research (AFOSR) under Grant No. FA9550-20-1-0015. Computational resources for this work is provided by Ohio Supercomputer Center⁴² under Project No. PAS0971.