Modeling of heterogeneous site energy distributions in precipitate nucleation

The simulation of heat changes resulting from phase transitions can help to interpret differential scanning calorimetry (DSC) measurements, e.g. of metallic alloy systems in which multiple reactions overlap during non-isothermal heat treatments. So far, simulated DSC curves mostly exhibit sharp reaction peaks as commonly just one mean energy value for a certain type of nucleation site is assumed. This work proposes an efficient model for treating heterogeneous nucleation site energy variations within the framework of classical nucleation theory (CNT). The site energies are assumed to vary according to a Rayleigh distribution and a scaling function. The effect on the nucleation behavior of precipitates is studied. A consideration of the distribution of heterogeneous site energies has the potential to significantly smoothen the numerical treatment of precipitation processes compared to the non-distributed case. The comparison to previously published simulations of DSC curves during the cooling of an AA6005 aluminum alloy demonstrates the advantages of this extension, especially for slow cooling rates.

In this regard, a recently published computational analysis modeled nucleation at secondphase particles to simulate precipitation during continuous cooling of an AA6005 aluminum alloy [20].The calculations therein are compared to differential scanning calorimetry (DSC) data from Milkereit et al [21,22] and good agreement for the peak positions and areas over a wide range of cooling rates is obtained.However, particularly at low cooling rates, the simulated exothermic reaction peaks exhibit a significantly narrower and sharper shape than the experimental data.The two precipitation reactions considered in this case both nucleate on previously existing Fe-containing precipitate particles (coarse primary AlFeMnSi precipitates at higher temperatures and Mn-rich dispersoids at lower temperatures [21]).These precipitates are present in a certain range of sizes, shapes, and probably orientation relationships with the matrix.It stands to reason that the differences between experimental data and simulations can be attributed to the fact that real material microstructures offer heterogeneous nucleation sites in a certain range of energies, while only one single energy is used in the simulations [20].
It is well known that, for example, the grain boundary energy can vary significantly depending on the crystal orientation and the segregation of solute atoms [23].Similarly, the dislocation line energy depends on the local configuration of a given dislocation segment, such as the dislocation character (edge and screw components) [24] and the dislocation spacing [25].
While the magnitude of energy values associated with a given type of defect is usually well explored in the literature, information on the distribution of defect energies is rather limited.A few studies [26][27][28] on grain boundary energies in recrystallized microstructures report an inverse relationship between the multiples of a random distribution and the relative grain boundary energy, meaning that grain boundaries with lower energy tend to be more frequent than those with higher energy.Based on considerations on energy minimization, this can be expected for any defect in a sufficiently equilibrated microstructure and, in a broader sense, suggests a right-skewed distribution of defect energies.
In contrast to the complexity of a real microstructure, computer simulations commonly use only one single energy value to describe heterogeneous nucleation sites.Solidification models are an exception in this regard, where pragmatic approaches to address this issue include using polynomial laws [29] or Gaussian distributions [30,31] to assign different activation undercooling values to a pre-determined number of nucleation sites.However, these formulations are not embedded into CNT and describe the nucleation of the solid phase from the liquid independent of time.Hence, an increment in undercooling is required to create new nuclei, and the nucleation rate becomes zero for constant undercooling.In this regard, one publication introduces a distribution of cavity angles into CNT to describe polymeric foaming [32].However, since it represents a purely geometrical contribution to the catalytic factor, this approach is not applicable in the present context of nucleation site energy distributions.
Inspired by the work of Zurob et al on recrystallization [33], a pragmatic treatment is introduced here to account for inhomogeneities within the framework of CNT.A variation of heterogeneous nucleation energies is considered based on a Rayleigh distribution and a scaling function.It is shown that considering nucleation site energy distributions improves the accuracy of DSC simulations, especially for slow cooling rates.

Setup
The starting point for the present work is the treatment of heterogeneous nucleation site energies introduced in [20], where the precipitation of stable and metastable phases has been analyzed during continuous cooling of an AA6005 aluminum alloy.Therein, the Gibbs energy change ∆G nuc upon formation of a nucleus with radius ρ, is expressed as with T being the temperature, X i the chemical compositions of the precipitating phase and the matrix, D f the driving force, γ eff the effective interfacial energy, and G het the heterogeneous nucleation energy.The latter can be separated into a geometrical part ξ geo and an energy τ , both of which are characteristic of the given nucleation site: The energy variable τ may correspond to, e.g. a grain boundary energy or a dislocation line energy.The scaling function m(x) changes the distribution characteristics and will be introduced and explained in a later paragraph.The variable x represents an abstract measure of the favorability of a given nucleation site.First, the Rayleigh distribution f(x) is introduced as with σ being the mode of the distribution.The Rayleigh distribution is mainly chosen for convenience because its cumulative distribution function F(x) is given by an exact analytical expression that can be easily inverted.The cumulative distribution function F(x) further represents the normalized number density of available nucleation sites (complement to the occupancy).Demanding that F(0) = 0, the expression for F(x) is given as where N 0 represents the total and N av the available number of nucleation sites.The variation in G het , as defined in equation ( 2), is realized with the scaling function m(x), which assigns every x-position in the distribution a particular value of G het .The form of m(x) is a priori undefined, but it should (i) provide a smooth transition from non-distributed to distributed energies and (ii) be unity at the mode of the Rayleigh distribution.The expression which is adopted in the present work is The first requirement is realized with the parameter d, which essentially controls the width of the distribution and adopts values between zero and one.For d equal to zero, every nucleation site is identical concerning the site energy G het .The parameter n mainly influences the skewness of the distribution (symmetrical or right-skewed).It is practical to define a cut-off x cut , for which a value is selected that ensures F(x cut ) > 0.999.The choice of σ, while irrelevant for the actual calculation of the nucleation rate, is important for the definition of the cut-off, and for σ = 2, a value of x cut = 7.5 is chosen.The functions describing the distribution of site energies for some choices of parameters are shown in figure 1.
The Rayleigh distribution with the chosen mode and the associated cut-off remains fixed in the evaluation process (figure 1(a)).Depending on the choice of the parameters d and n from equation ( 5), the scaling function m(x) is altered (figure 1(c)), which translates the ratio of available to total nucleation sites into site energy scaling values (figures 1(b) and (d) and equation ( 2)).A broad distribution (d = 1.0, n = 0.75) means that the density of nucleation sites with energies at the maximum and minimum of the scaling function is relatively high.In Heterogeneous nucleation naturally begins at the most favorable sites, that is, those with the largest value of G het .As nucleation commences, the sites with the highest energy are progressively occupied, and the remaining sites have an increasingly lower value of G het .This effect is accounted for by tracking the number of available nucleation sites and updating the energy accordingly.For this purpose, the actual value of x at any given point in time is evaluated from the occupancy using equation ( 4) with Thus, the value of m(x) for a given timestep is obtained and used to determine the current G het and the nucleation rate J.The reduction of available sites based on ultimately results in a lower scaling value for the subsequent timestep.The automatic timestep control is set up such that the nucleation rate does not decrease the value of x by more than a pre-defined fraction ε.The process is schematically shown in figure 2.

Parameter study
This section discusses the influence of the distribution width d on the evolution of the fundamental quantities that govern the precipitation process.The MatCalc [9] simulation for cooling of an AA6005 aluminum alloy at 0.1 K min −1 , as published in [20], is considered for this purpose (for details on chemical composition and parameters, see tables 1 and 2 in section 3).For simplicity, only precipitation of the stable β-Mg 2 Si phase is considered here.The default setup from Miesenberger et al [20] corresponds to a value of zero for parameter d, which is increased to 0.2 and 0.4 in the parameter study given in figure 3, meaning that the energy gain at the heterogeneous nucleation site is progressively increasing.The parameter n is held constant at a value of 1, meaning the relation between x and m(x) is linear.The nucleation process as described by CNT is always time-dependent and not all nucleation sites are activated at the exact same moment, even if their associated nucleation barrier is identical (d = 0).Consequently, the first nuclei have a certain timeframe for growth before full occupancy has been reached.Especially for slow cooling rates, growth may sufficiently affect the supersaturation to suppress further nucleation before all sites are occupied.This is why, during continuous cooling, nucleation may come to a halt temporarily until the temperature is low enough (i.e., the driving force is high enough) to activate the remaining nucleation sites in a second nucleation wave.Naturally, the comparably abrupt activation of identical nucleation sites for the non-distributed case makes multiple nucleation waves less likely to occur.However, the rate of 0.1 K min −1 chosen in the parameter study is sufficiently slow and leads to two distinct nucleation waves at 490 • C and 380 • C (figure 3(c)).In contrast, the distributed nucleation site energy buffers the initial nucleation wave, spreading it over a wider temperature range.Therefore, the reduction of the supersaturation is not as abrupt, and the temporary increment of the critical nucleation barrier and the critical radius around 490 • C is less pronounced (figures 3(a) and (b)).Overall, the distinct nucleation waves from the non-distributed case begin to merge with increasing distribution width (figure 3(c)), and the calculated DSC signal (figure 3(d)) is significantly smoother.The evolution of the scaling function m(x) during cooling is shown in figure 4.

Results and discussion
All simulations are performed with the thermokinetic software package MatCalc (version 6.04) [9] with the open thermodynamic database 'mc_al.tdb'(version 2.035) [34] and the open diffusion database 'mc_al.ddb'(version 2.004).The chemical composition used in the calculations is that of the simplified system given in table 1.The simulations are now carried out for both experimentally observed types of precipitates: β-Mg 2 Si and B'-Al 4 Mg 8 Si 7 [21].Both nucleate heterogeneously at second-phase particles.
The procedure for calculating DSC curves and the nucleation model for precipitation at second-phase particle interfaces based on the parameters κ i and κ s is reported in [20].The simulations therein are used as a benchmark to demonstrate the improvements offered by the present approach.The simulation parameters are summarized in table 2 and apply to all cooling rates.
The interface energy of the second-phase particle is assumed to be comparable to a grain boundary in Al, with an energy of 0.5 Jm −2 [28,36].The precise nucleation sites for β-Mg 2 Si and B'-Al 4 Mg 8 Si 7 are large intermetallic constituent particles and dispersoids, respectively.Hence, the number densities in table 2 are based on typical experimental values [37,38].The measured and simulated DSC curves are compared in figure 5.
During the cooling of AA6005, stable β-Mg 2 Si forms at high temperatures and low cooling rates, whereas faster cooling rates promote the nucleation of metastable B'-Al 4 Mg 8 Si 7 [21].Particularly for the high-temperature reaction peak (β-Mg 2 Si), the introduction of an energy distribution significantly improves the quality of the calculated DSC signal.The agreement with experimental data is excellent across all cooling rates.
The impact of the distribution model is amplified for decreasing cooling rates as the difference in nucleation onset temperatures between the distributed and non-distributed case translates into increasingly larger timespans.On the contrary, the differences between the two simulations for the low-temperature reaction (B'-Al 4 Mg 8 Si 7 ) are moderate because the metastable phase only becomes dominant at higher cooling rates.Overall the simulations are slightly less accurate for the low-temperature reaction compared with the high-temperature reaction.This could be connected to some challenges regarding, for instance, the thermodynamic description of metastable phases, which is more complex than that of a stoichiometric equilibrium phase such as β-Mg 2 Si.In the present work Cu is neglected and only B'-Al 4 Mg 8 Si 7 is considered because the experimental data [21] indicates this to be the dominant phase.It is well known that moderate changes in the Cu content will affect the precipitation sequence [39].Hence, certain amounts of other precipitates may also contribute to the DSC signal at lower temperatures, depending on the relative Cu content in the matrix, which is directly related to the cooling rate.This being said, the agreement is fairly reasonable given the simplifications used in the present work.
Regarding the simulation parameters used, the new setup with distributed nucleation site energies requires a small increment with respect to the precipitate interface energy inside the incoherent boundary (see κ i in table 1).This change moves the overall nucleation onset to a slightly lower temperature and compensates for the shift introduced by broadening the energy distribution.The remaining parameters are almost identical.Evaluation of equation ( 5) with the distribution parameters listed in table 2 shows that the first nucleation sites occupied for β-Mg 2 Si and B'-Al 4 Mg 8 Si 7 in the distributed case  [21,22] is compared to the MatCalc [9] simulations in the present work and those by Miesenberger et al [20].are roughly 2.5 and 1.8 times more favorable than the default values at the mode, respectively.

Summary
In the present work, we propose a pragmatic and efficient way to incorporate energy variations into the description of heterogeneous nucleation sites within the framework of CNT.The value for the heterogeneous nucleation site energy is varied by linking a Rayleigh distribution to a scaling function, which alters the nucleation barrier based on the current occupancy of nucleation sites.Nucleation of metastable and stable precipitates at Fe-containing constituents and dispersoids in an AA6005 aluminum alloy during cooling after solutionizing is selected as a use case.The simulations are compared to experimental DSC measurements.It is shown that the proposed extension, based on just two parameters, provides an efficient means for improving simulated DSC curves, especially at low cooling rates.

Figure 1 .
Figure 1.Overview of the functions used to describe a continuously distributed nucleation site energy with different parameters: (a) Rayleigh distribution f(x) with fixed mode σ and cut-off xcut, (b) Rayleigh distribution as a function of the scaling function m(x).(c) Scaling function m(x) (d) cumulative distribution F(x) as a function of m(x).

Figure 2 .
Figure 2.Link between the heterogeneous energy distribution and the nucleation rate J as calculated from CNT.In every timestep, the heterogeneous nucleation site energy is calculated based on the current value of m(x), which directly depends on the occupancy for the given type of nucleation site.

Figure 3 .
Figure 3. Influence of the distribution width d on the evolution of (a) the critical nucleation barrier G * /k B T (b) the critical nucleation radius r crit (c) the nucleation rate J and (d) the simulated DSC-signal during cooling of AA6005 at 0.1 K min −1 .For this comparison, only the heterogeneous nucleation of β-Mg 2 Si is considered.Calculations are performed with MatCalc [9].

Figure 4 .
Figure 4. Evolution of the scaling function m(x) during cooling of AA6005 at 0.1 K min −1 .For this comparison, only the precipitation of β-Mg2Si is considered.

Figure 5 .
Figure 5. DSC curves for cooling of an AA6005 aluminum alloy at various rates.The experimental data from Milkereit et al[21,22] is compared to the MatCalc[9] simulations in the present work and those by Miesenberger et al[20].

Table 1 .
Simplified chemical composition for aluminum alloy AA6005 used in the present work.All values in mass percent.

Table 2 .
Simulation parameters used in the precipitation kinetics calculations for cooling of AA6005.