The total Gibbs functional for the class of considered systems has the form: G = 1 V m ∫ V G c h x μ , c d + G ∇ ∇ x μ , ∇ c d + G e l d r , (1) where the term G _ch relates to the chemical Gibbs energy depending on the concentrations of the alloy components and point defects; G _∇ is a gradient energy part and the last term G _el corresponds to an elastic contribution; V _m is the molar volume; integration is taken over the whole volume V.

The free energy density G _ch can be decomposed into two terms, describing atomic and defect subsystems as follows: G _ch = G _Fe−Cr−Al + G _d . The molar Gibbs energy of the atomic subsystem G _Fe−Cr−Al is defined as: G F e − C r − A l = G F e − C r − A l ref + G F e − C r − A l i d + G F e − C r − A l e x , where G F e − C r − A l ref is the reference Gibbs energy, G F e − C r − A l i d denotes the ideal Gibbs energy contribution due to random mixing of atoms, G ^ex is the excess term defining a deviation from ideality. The reference Gibbs energy is written in the standard form G F e − C r − A l ref = ∑ μ G μ 0 x μ . The corresponding expressions for the reference energies G μ 0 can be obtained with the help of phase diagrams by using CALPHAD method (Dinsdale, 1991). The ideal part G ^id relates to entropic contribution G ^id = RT∑ _μ x _μ  ln x _μ , where R is the universal gas constant, T is the temperature. The excess term G ^ex is defined by temperature dependent interaction coefficients L _μν = L _μν (T) in the following form: G ^ex = ∑ _μ≠ν x _μ x _ν L _μ,ν . The Gibbs energy G _d for a point defect subsystem d = {i, v} can be written as follows: G d = ∑ d G d f + G d i d + G d int . Here G d f is the defect formation energy. It is defined through the formation energies G d f , μ of defects in pure materials Fe, Cr and Al and mean concentrations of atoms in the alloy as G d f = ∑ μ G d f , μ x ̄ μ , where x ̄ μ relates to nominal concentration of μ. The corresponding entropic contribution is G d i d = R T c d ⁡ ln c d . Interactions between defects with host atoms and alloying elements ( G i int and G v int ) can be described by the component G d int = c d ∑ μ x μ G d − μ int , where for the defect-atom interaction energies we use definitions (Was, 2016): G d − μ int = ( G coh μ + G d f , μ ) / Z . Here G coh μ denotes the corresponding cohesive energy, Z relates to the coordination number. Point defects formation energies for pure materials ( G d f , μ ) and the corresponding cohesive energies ( G coh μ ) can be taken from the literature or can be calculated by ab initio calculations or molecular dynamics simulations.

By combining all above expressions, the total Gibbs energy for the ternary alloy Fe − Cr − Al with point defects takes the form: G c h x F e , x C r , x A l , c v , c i = G F e 0 x F e + G C r 0 x C r + G A l 0 x A l + G v f c v + + G i f c i + R T x F e ⁡ ln x F e + x C r ⁡ ln x C r + x A l ⁡ ln x A l + c v ⁡ ln c v + c i ⁡ ln c i + x F e x C r L F e , C r + x F e x A l L F e , A l + x C r x A l L C r , A l + c v x F e G v − F e int + x C r G v − C r int + x A l G v − A l int + c i x F e G i − F e int + x C r G i − C r int + x A l G i − A l int . (2)

The gradient term G _∇ has the form: G ∇ = ∑ μ ( κ μ / 2 ) ( ∇ x μ ) 2 + ∑ d ( κ d / 2 ) ( ∇ c d ) 2 . The energy constants {κ _μ } are defined according to the Hilliard approach Ref. (Hilliard, 1970). Next we assume κ _μ = κ, where κ = L _Fe,Cr a ²/6, a = ∑ _μ a _μ x _μ , a _μ is the lattice constant for pure materials. In order to avoid negative gradient energy and ensure the disorder of the system, for the energy constants κ _i and κ _v we use the formalism discussed in Refs. (Liang et al., 2018; Yan et al., 2020) and define them as κ _d = κ _v = κ _i .

The α − α′-interface can be considered as a semi-coherent due to small difference in lattice parameters of two main elements Fe and Cr, as usual (Ribis and Lozano-Perez, 2012). To make our description more general we consider the elastic contribution G _el in the Gibbs energy has the standard form G e l = V m 2 C ijkl ( e i j − e i j 0 ) ( e k l − e k l 0 ) / 2 , here C _ijkl denotes the elastic constants tensor, e i j = ∂ x i u j + ∂ x j u i / 2 is the elastic strain tensor, where u = (u _x , u _y , u _z ) is the displacement vector. By using the Vegard law one defines the lattice parameters misfit e μ 0 = ( 1 / a ) ( d a / d x μ ) caused by introducing alloying element μ. Hence, the eigenstrain tensor e i j 0 acquires the general form e i j 0 ( r , t ) = δ i j ∑ μ e μ 0 δ x μ ( r , t ) , here δ _ij is the Kronecker delta-symbol, δ x μ ( r , t ) = ( x μ ( r , t ) − x ̄ μ ) . By taking into account that α′-phase is hard phase, one can consider only Chromium contribution to elastic effects. In such a case eigenstrain and corresponding stress tensors acquire the form e i j 0 ( r , t ) = δ i j e C r 0 δ x C r ( r , t ) , σ i j 0 = C ijkl e k l 0 with e C r 0 = 0.00614 . The elastic energy can be decomposed into four parts: G e l 0 = ( V m / 2 ) σ i j 0 e i j 0 related to eigenstrain components (proportional to squared concentration of alloying elements); G e l * ( { e s } s = 1 6 ) related to elastic components defined through elementary strains e _s ; G e l c is responsible for coupling between deformation elastic deformations e ₁ and alloying element concentration; G e l a corresponds to cubic symmetry.

As was shown previously (Cahn, 1962; Onuki, 2002; Khachaturyan, 2013) in BCC metals/alloys the minimum free energy is along the soft direction ⟨100⟩. The contribution to elastic energy coming from the lattice mismatch takes the form (Cahn, 1961): G e l 0 = Y [ 100 ] V m ( e C r 0 δ x C r ) 2 , Y [ 100 ] = C 11 + C 12 − 2 ( C 12 ) 2 / C 11 , where Y _[100] is the averaged stiffness parameter for cubic crystals related to the free energy minimum along ⟨100⟩ direction.

By using standard relations between stiffness components and elastic moduli one can write a contribution coming from elastic continuum as follows (Onuki, 2002): G e l * = 1 2 K e 1 2 + μ 2 Φ e 2 , e 3 + μ 3 Ψ e 4 , e 5 , e 6 , (3) where the stretching (Φ) and shear (Ψ) contributions are defined by Φ = 1 4 π 2 3 − cos 2 π e 2 − − cos 2 π e 2 + − cos 4 π 6 e 3 , Ψ = 1 4 π 2 3 − cos 2 π e 4 − cos 2 π e 5 − cos 2 π e 6 . (4)

Here elementary strain components are: e ₁ = e _xx + e _yy + e _zz , e 2 = e x x − e y y , e 3 = ( 2 e z z − e x x − e y y ) / 3 e 4 = e x y , e ₅ = e _yz , e ₆ = e _zx , e 2 ± = e 2 / 2 ± e 3 / 6 . For small strains we get Φ ≈ μ 2 ( e 2 2 + e 3 2 ) / 2 , and Ψ ≈ μ 3 ( e 4 2 + e 5 2 + e 6 2 ) / 2 , as usual.

The component describing coupling between deformation and concentration has the form (Onuki, 2002): G e l c = α V m δ x C r ∇ ⋅ u , where α = − e C r 0 3 K . Elastic misfit results in anisotropy with a parameter ξ _a = (C ₁₁ − C ₁₂ − 2C ₄₄)/C ₄₄ defining anisotropy contribution Onuki (2002). G e l a = V m 2 τ a ∑ i ≠ j | ∇ i 2 ∇ j 2 w | , τ a = e C r 0 3 K 2 ξ a C 44 C 11 2 , (5) where ∇ _i ∇ _j w = n _i n _j δx _Cr  (r, t), n _i denotes component of the interface normal n.

In our calculations we will follow approach developed by Onuki [see Ref. (Onuki, 2002) and citations therein] and results of molecular dynamics (Dai et al., 2022) and assume that all elastic constants depend on concentration of alloying elements: C i j = ∑ μ C i j μ x μ . A condition of mechanical equilibrium implies ∂ x j C ijkl e k l − δ k l e C r 0 δ x C r = 0 .

By combining all above expressions, the elastic contribution G _el to the total Gibbs energy for the ternary alloy Fe − Cr − Al becomes the form: G e l = Y 100 V m 2 e C r 0 δ x C r 2 + V m τ a 2 ∑ i ≠ j | ∇ i ∇ j w | − V m 3 K e C r 0 δ x C r ∇ ⋅ u + Φ e 2 , e 3 ; x C r + Ψ e 4 , e 5 , e 6 ; x C r . (6)

In our approach we use concentration dependent stiffness constants K = K ₀ + K ₁ x _Cr , μ ₂ = μ ₂₀ + μ ₂₁ x _Cr , μ ₃ = μ ₃₀ + μ ₃₁ x _Cr , ξ a = ξ 0 a + ξ 01 a x C r , where bulk, shear moduli and anisotropy components are collected in Table 1. For stiffness parameter Y _[100], next, we use elastic constants defined through the nominal compositions, i.e., Y [ 100 ] = Y [ 100 ] x ̄ C r . By using the Hook law for elastic stress tensor σ i j = C ijkl e k l = C ijkl ( e k l − e k l 0 ) one finds for diagonal components σ _ii : σ _xx = σ _yy = σ _zz = (K ₀ + K ₁ x)e ₁, for off–diagonal ones we get σ _xy = (μ ₃₀ + μ ₃₁ x _Cr ) sin (2πe ₄)/2π, σ _yz = (μ ₃₀ + μ ₃₁ x _Cr ) sin (2πe ₅)/2π, σ _xy = (μ ₃₀ + μ ₃₁ x _Cr ) sin (2πe ₆)/2π.

By taking into account the mass conservation law and defining the concentration of the Iron as the host element as: x _Fe = 1 − x _Cr − x _Al dynamics of the complete system will be governed by the following equations of the Cahn-Hilliard type (Cahn, 1962; Hilliard, 1970): ∂ t x C r = ∇ ⋅ M CrCr ∇ δ G δ x C r + M CrAl ∇ δ G δ x A l , ∂ t x A l = ∇ ⋅ M AlAl ∇ δ G δ x A l + M CrAl ∇ δ G δ x C r , ∂ t c v = ∇ ⋅ L v ∇ δ G δ c v , ∂ t c i = ∇ ⋅ L i ∇ δ G δ c i (7)

Mobility coefficients M _Cr,Cr , M _Al,Al and M _Cr,Al are as follows (Huang et al., 1995; Wu et al., 2001; Wu et al., 2004; Lee et al., 2020): M C r , C r = x C r 1 − x C r 2 M C r + x C r x A l M A l + x C r x F e M F e , M A l , A l = x A l 1 − x A l 2 M A l + x A l x F e M F e + x A l x C r M C r , M C r , A l = x C r x A l x F e M F e − 1 − x C r M C r − 1 − x A l M A l , (8) where M _Fe , M _Al and M _Cr are the corresponding mobilities for pure elements, which are defined in the standard manner: M _μ = D _μ /RT, D _μ is the corresponding atomic diffusivity. The mobility for the ensembles of vacancies and interstitials we choose in the standard manner: L _d = D _d c _d /RT, defined through the corresponding point defect diffusivity D _d . Here we use a simplified model for point defects mobilities.

In the Cahn-Hilliard-Cook equation Eq. 7 with a variable mobility Eq. 8, it is not convenient to accurately discretize the gradient operator and divergence operator by using the finite-difference method. However, the Fourier spectral method developed in Refs. (Chen and Shen, 1998; Zhu et al., 1999; Kim, 2007; Canuto et al., 2012; Biner, 2017) allows us to solve numerically Eq. 7 without much more difficulty than the case with a constant mobility. In the framework of the Fourier spectral method the system of coupled non-linear differential equations takes the form: ∂ t x C r k , t = i k ⋅ M C r , C r r i k δ G δ x C r k r + M C r , A l r i k δ G δ x A l k r k , ∂ t x A l k , t = i k ⋅ M A l , A l r i k δ G δ x A l k r + M C r , A l r i k δ G δ x C r k r k , ∂ t c v k , t = i k ⋅ L v r i k δ G δ c v k r k , ∂ t c i k , t = i k ⋅ L i r i k δ G δ c i k r k . (9)

Here, x _Cr,Al  (k, t), c _v,i  (k, t), δ G / δ x C r , A l k and δ G / δ c v , i k represent the Fourier transforms of x _Cr,Al (r, t), c _v,i (r, t), δ G / δ x C r , A l and δ G / δ c v , i in r-space, respectively. k is the vector in the Fourier space with the magnitude k = k x 2 + k y 2 ; { } _k is the Fourier transform and [ ] _r means the inverse Fourier transform to real space.

In further modeling, we use thermodynamic and elastic parameters characterizing ternary systems collected in Table 1.

In order to define ranges of concentrations of alloying additions Cr and Al and temperature, when phase decomposition is possible we exploit the linear stability analysis of a homogeneous state to inhomogeneous perturbations. According to obtained kinetic equations one can find phase diagram defining critical values of concentration of alloying elements at fixed temperatures shown in Figure 1. It follows that phase separation related to instability of the solid solution to inhomogeneous perturbations is possible at fixed range of Chromium and Aluminum concentrations. With the temperature growth the instability domain bounded by spinodal curves shrinks. Outside of the corresponding spinodal the system is stable to homogeneous perturbations. In this domains the solid solution is realized. It is clear from this diagram that in the case of high values of the annealing temperature (see, for example, the curve at T = 850 K), when precipitation is impossible in the binary Fe-Cr systems, the addition of Aluminum with a concentration > 2 % induces such processes. At high temperatures (see curve at T = 900 K) with the temperature increase phase decomposition related to formation of α′-precipitates enriched by Cr is possible only at elevated concentration of Al and with content of Cr around 50 at%.

To study process of phase decomposition at annealing numerically we choose values of nominal concentration of two alloying elements and fix the annealing temperature according to stability diagram. We take solid solution as initial microstructure of the corresponding alloy with homogeneously distributed concentration fields of all components including point defects. As initial configurations we use: ⟨ x μ ( 0 ) ⟩ = x ̄ μ , ⟨ c d ( 0 ) ⟩ = c d e q ; ⟨ ( x μ ( 0 ) − x ̄ μ ) 2 ⟩ = 1 0 − 3 x ̄ μ , ⟨ ( c d ( 0 ) − c d e q ) 2 ⟩ = 1 0 − 3 c d e q , where c d e q relate to equilibrium point defect concentration. All simulations were done on square lattice L × L of the linear size L = NΔx with N (N = 128) nodes in each direction and effective lattice constant for simulations Δx = 1 by exploiting Fourier spectral method (Chen and Shen, 1998; Zhu et al., 1999; Canuto et al., 2012; Biner, 2017). We take Δt′ = 10⁻³ as dimensionless time step of integration. Our simulations were done in dimensionless time t′ = tℓ ²/D _Al , for spatial coordinate we use r′ = r/ℓ, where ℓ is of the lattice parameter a. Boundary conditions were periodic. We discuss dynamics of α′-phase formation and growth and analyze effect of alloying elements and annealing temperature onto main statistical properties of precipitates.

According to provided stability analysis in order to study α′-phase precipitation next we consider alloy Fe − 30%Cr − 5%Al as a model system with annealing temperature T = 710 K. The main indicator of phase separation is the growth of the dispersion of the concentration field, determined in the standard way: ⟨ ( δ x μ ) 2 ⟩ = ⟨ x μ 2 ⟩ − ⟨ x μ ⟩ 2 . This quantity plays the role of an effective order parameter during phase decomposition. Indeed, if ⟨ ( δ x μ ) 2 ⟩ grows in time, it means that enriched and depleted domains of μ-component are formed and can be distinguished by concentration values, therefore phase decomposition occurs (García-Ojalvo and Sancho, 2012).

Typical scenario of phase decomposition is shown in Figure 2. Snapshots of the concentration field of Cr and Al are shown in Figure 2A at time instants related to markers indicated at protocols of effective order parameters in Figure 2B. Form Figure 2A it follows that after a certain incubation period (t < t _c ) domains enriched and depleted in the concentration of impurities emerge from the initial homogeneous distribution (solid solution). Thus, after 22 h of annealing, fine Chromium precipitates can be observed. At the same time, Aluminum is homogeneously distributed in the matrix (outside the precipitates). The number of precipitates and their size increase with time. During this nucleation stage, all formed precipitates are characterized by almost the same size. With further annealing one passes to the coarsening stage, where precipitates continue to grow and interact with each other: small precipitates dissolve, while large ones take material from the matrix and continue to grow following Ostwald ripening scenario. Detailed analysis of α′-phase precipitation can be seen from protocols of the order parameters shown in Figure 2B. One finds that after a certain incubation period dispersion of both x _Cr and x _Al grow rapidly, meaning a passage to the nucleation stage. At the coarsening stage the growth rate of dispersions ⟨ ( δ x C r ) 2 ⟩ and ⟨ ( δ x A l ) 2 ⟩ decreases, they continue to increase slowly and the system goes to a quasi-stationary regime. Here the total amount of precipitates decreases with time, and the dispersion in their sizes increases.

Results of statistical analysis of precipitation are shown in Figure 3. Protocols for the Chromium concentration in the α′-phase and the volume fraction of precipitates are shown in Figure 3A. Both values grow rapidly after the incubation period t _c at the nucleation stage whilst their slow growth is observed at coarsening stage. At the physical time of ∼ 230 hours almost 70% of Chromium is concentrated in the α′-phase, whilst the volume fraction of precipitates is about 20%. Obtained data relate well to experimental observations and results of numerical modeling Yang et al. (2020), Chen et al. (2019), Lee et al. (2020). Figure 3B illustrates dynamics of the mean precipitate radius and their number density recalculated from obtained data. It can be seen that the number of precipitates and their mean radius rapidly increase with time at the nucleation stage. At the coarsening stage the mean precipitate size (radius) ⟨R _p ⟩ increases according to the power-law ⟨R _p ⟩∝ t ⁿ with the exponent n close to the classical Lifshits-Slyozov-Wagner (LSW) value 1/3; the total amount of precipitates at this stage decreases exponentially (see power-law and exponential decay fitting shown by dash and dash-dot lines, with fitting parameter τ). Data obtained for the mean precipitate size (around 1–1.2 nm at annealing time 300 h) and precipitate number density (1.5 × 10¹⁹ cm ⁻³) for the studied alloy relate well to experimental and numerical data for class of system related to Fe-Cr-Al alloys (Messoloras et al., 1984; Chen et al., 2019; Lee et al., 2020; Yang et al., 2020). Scaling dynamics of mean precipitate size in Fe-Cr alloys was observed experimentally [see, for example, Refs. (Brenner et al., 1982; Messoloras et al., 1984; Bley, 1992; Hyde et al., 1995; Soriano-Vargas et al., 2010)], where the scaling exponent n depends on the nominal composition and annealing temperature.

Size distribution α′-phase precipitates f (R _p /⟨R _p ⟩) is shown in Figure 3C at different time instants by markers. It can be seen that data obtained from the numerical experiment at different times relate to each other meaning an universality of the size distribution of precipitate sizes at the coarsening stage. The solid curve shows the classical universal LSW distribution (Lifshitz and Slyozov, 1961; Wagner, 1961). The obtained numerical data are in good agreement with the LSW distribution. Hence, LSW theory can be used to predicate particle size distribution in Fe-Cr-Al systems at thermal treatment.

Let us consider a behavior of defect subsystem during precipitation. From Figure 4A it follows that the dispersion of the vacancy concentration decreases with time, while the dispersion of the concentration of interstitials slowly increase after the decreasing. Spatial rearrangement of point defects is well seen from snapshots in Figure 4A (see insets). Here vacancies are quasi-homogeneously distributed in a bulk, while interstitials are mostly localized in small precipitates (R _p < ⟨R _p ⟩) and near the phase interface in large precipitates (R _p > ⟨R _p ⟩). Detailed information about defect distribution around precipitates can be found in profiles of alloying elements and defects shown in Figure 4B. Here one finds the following: the interstitial concentration is higher in precipitates compared to its value in the matrix phase; they are mostly localized at regions of precipitates with elevated curvature (see the snapshot of c _i (r) in Figure 4A); vacancy concentration is locally increased around interfaces and has lower values in precipitates comparing to their values in matrix that relates to thermodynamics; α′-phase precipitates contain extremely small amount of Al atoms comparing to matrix phase.

By studying phase decomposition of alloys with different content of alloying elements and different temperatures one finds that an increase in one of the component can control time scale for precipitation kinetics. By studying α′-phase precipitation it was found that the incubation time t _c for precipitation decreases with content of Cr and Al as is shown in Figures 5A, B. Here we denote critical value for concentration of alloying elements as x C r , A l c (from phase diagram in Figure 1). It means that an increase in alloying element content accelerates driven forces for phase decomposition decreasing the incubation time t _c needed for emergence of precipitates. From theory of phase decomposition it is known that the quantity t _c depends on the annealing temperature (Khachaturyan, 2013). The corresponding dependence is shown in Figure 5C for the Fe-30%Cr-5%Al system (as an example), where one gets an infinite time needed to initiate phase decomposition in the vicinity of the critical temperature related to spinodal depicted in the phase diagram in Figure 1.

Effect of alloying elements and temperature influence onto dynamics of mean precipitate size, volume fraction of precipitates and precipitate size distribution is shown in Figure 6. In Figures 6A, B we plot the corresponding protocols in physical time scale according to the time related to phase decomposition starting from t _c . From obtained results it follows that an increase Al content results in slight increase in precipitate size and volume fraction (cf. curves 1 and 2), whereas Cr addition increases the corresponding values of above quantities crucially (cf. curves 1 and 3). These results are related well to previous studies [see for example, Ref. (Lee et al., 2020)]. With the temperature decrease from 710 K down to 690 K (cf. curves 1 and 4) one gets low values of the mean precipitate size and low values for volume fraction of α’ precipitates that is constituent with results discussed in Ref. (Chen et al., 2019). Distributions over precipitate size in alloys with different content of Al and Cr and two different temperatures are shown in Figure 6C. One finds that obtained data are well predicted by classical LSW distribution. Therefore, it follows that additional introduction of Cr or Al does not change universality of precipitate size distribution at coarsening stage at different temperatures. By comparing dynamics of the precipitate size and precipitate size distributions in alloys with different content of doping one concludes that both dynamics of the system and statistical properties of growing precipitates remain universal. The difference in measured quantities (mean precipitate size, precipitate number density, precipitate volume fraction) is related mostly to scaling factors determined by nominal composition of the studied alloys.