Atomistic and continuums modeling of cluster migration and coagulation in precipitation reactions.

The influence of vacancy preference towards one of the constituents in a binary system on the formation of precipitates was investigated by atomistic and continuums modeling techniques. In case of vacancy preference towards the solute atoms, we find that the mobility of individual clusters as well as entire atom clusters is significantly altered compared to the case of vacancy preference towards the solvent atoms. The increased cluster mobility leads to pronounced cluster collisions, providing a precipitate growth and coarsening mechanism competitive to that of pure solute evaporation and adsorption considered in conventional diffusional growth and Ostwald ripening. A modification of a numerical Kampmann-Wagner type continuum model for precipitate growth is proposed, which incorporates the influence of both mechanisms. The prognoses of the modified model are validated against the growth laws obtained with lattice Monte Carlo simulations and a growth simulation considering solely the coalescence mechanism.

Introduction

The vacancy/atom exchange mechanism is the dominant physical mechanism for diffusive transport in crystalline solids. This process is thermally activated, with a jump probability given by exp(−E/k), E being the migration energy barrier, k being the Boltzmann constant and T being temperature. E is often represented as a superposition of pair-wise bond energies between atoms, the magnitude of which depends on the sort of atoms occupying the connected sites. In the present work, in addition to the atoms occupying regular crystal lattice sites, vacancies are treated as a separate atomic species, with the energies εVa for the bonds between the lattice atoms X and the vacancy Va.

If the vacancy/atom interactions εVa have different values for different atoms, vacancy preference toward particular sorts of alloying elements occurs, i.e. the vacancy will spend more or less mean residence time in the neighborhood of specific atoms. Apart from the influence on the long-range diffusion kinetics of the single atoms (monomers), such a preference can also increase or reduce the mobility of entire atom clusters. Although cluster migration is observed in systems without vacancy preference at sufficiently low temperatures [1] also, its occurrence is significantly pronounced when the solute atom/vacancy complex is energetically more stable than the solvent/vacancy complex [2], [3], [4]. The random movement of clusters occasionally leads to cluster collisions, resulting in cluster coagulation and, as such, constitutes another precipitate growth and coarsening mechanism, in addition to the absorption and desorption of monomers. Using classical lattice Monte Carlo simulation [5], the combined absorption/desorption and cluster collision kinetics is investigated in the first part of this work.

In a complementary modeling approach, a numerical Kampmann–Wagner-type model (NKW) for simulation of multi-component multi-particle precipitation kinetics, as implemented in the software MatCalc [6], [7], is applied for simulation of these processes on the continuum scale. With this approach, computationally efficient studies can be performed on the interaction of a large ensemble of precipitates taking into account nucleation, growth and coarsening of all particles involved in the process. In the original implementation of the corresponding models in MatCalc, absorption and emission of monomers at the precipitate surface are considered as the only mechanisms controlling the kinetic evolution of the system.

Analytically, coarsening of precipitates governed by the evaporation and absorption of monomers can be described by the Ostwald-ripening model, during which smaller precipitates dissolve by preferential emission (evaporation) of monomers and larger particles keep growing by preferential absorption of these. The Lifshitz–Slyozow–Wagner (LSW) theory predicts an evolution of the particle number density with time as t−1 and the mean article radius as t1/3 [8], [9]. If cluster coagulation is also considered as a possible mechanism for precipitate evolution, a mathematical treatment of the combined processes becomes rather involved [10]. Simple closed analytical solutions are lacking.

In the present paper, the role of cluster coagulation in precipitation reactions is examined, first, by atomistic Monte Carlo simulation of nucleation, growth and coarsening on the example of bcc-Cu precipitates in a bcc Fe matrix. The differences in the evolution of precipitate number density and mean radius are studied for various vacancy/solute interaction energies εVa and discussed on basis of the resulting cluster evolution characteristics. Subsequently, complementary simulations are performed with a simple numerical model where clusters of different size and mobility are treated as single molecules, sticking to each other after collision, and evolving the system purely on basis of the coagulation mechanism. This simulation allows for a clear separation of the effect of coagulation from the combined precipitation reaction.

Finally, based on these observations, a semi-empirical continuum model is suggested, which is implemented into the precipitation kinetics framework of MatCalc. It is shown that this simple model is fully capable of reproducing the system behavior observed in the atomistic simulations, for both limiting situations of either negligible influence of coalescence, or cases where particle collisions represent the dominant growth and coarsening mechanism.

Modeling solid-state precipitate coagulation – state of the art

The numerical Kampmann–Wagner [11] model (NKW), which is widely used for precipitation kinetics simulations, allows for an efficient and accurate description of precipitation processes. In the NKW model, the vast number of precipitates in a representative volume element is grouped into classes of precipitates with same size and chemical composition. Each size class is determined by the number of identical precipitates, which are formed in a nucleation step based on the nucleation rate J and time step Δt [6], [7]. The further evolution of this size class is determined by diffusional growth. Eventually, the precipitates of a certain size class dissolve again, which causes the particular size class to be removed from the system. The classical NKW model does not involve the evolution of precipitates (or precipitate size classes) caused by coagulation.

Cluster coagulation has been treated several times in literature. The work by Smoluchowski [10] shows that the evolution of the precipitate number density N of a size class i is given asβ(i, j) is the collision rate of particles with sizes i and j. In some special cases, this set of equations can be solved analytically, depending on the functional form of β(i, j), which in turn depends on the particle collisions being of either Stokes [12], turbulent [13], Brownian [14] or gradient type [15].

In the case of particle migration in a solid solution, coagulation is governed by Brownian-type motion. For this, Smoluchowski derived an expression for β(i, j), which is a function of the cluster diffusion coefficients D and the precipitate radii r as

In the liquid phase, where the cluster diffusion coefficient is inversely proportional to its radius (Einstein–Stokes equation), β(i, j) can be reasonably well approximated with values independent of the sizes i and j. Under these conditions, analytical solution of the equation set is possible [10].

In contrast, the values of β(i, j) describing collisions in solid-state systems are strongly size-dependent. Binder and Stauffer [16] suggest that the cluster migration rate scales with n−4/3, n being the number of atoms in the cluster. Binder [17] remarks that this dependency might evolve towards n−1 with continued cluster growth. Fratzl and Penrose [1] confirmed the latter function in Monte Carlo simulations with the vacancy exchange mechanism. Consequently, β(i, j) is found to scale with cluster radius r in the range between r−4 and r−3 (as r ∼ n1/3).

In the absence of closed analytical expressions for the evolution of the number densities of particles (or particle size classes) subjected to simultaneous monomer attachment and detachment and particle coalescence, this problem is investigated numerically. For this task, three different approaches are used, involving atomistic (Monte Carlo) simulations and continuums simulations based on the mean-field model of precipitate nucleation, growth and coarsening within the MatCalc precipitation kinetics framework [18]. The corresponding approaches are introduced in the following section.

Simulation methods

Lattice Monte Carlo (LMC)

The atomistic simulations of the precipitation kinetics are performed with the Lattice Monte Carlo (LMC) method as implemented recently into MatCalc [18]. For the present investigation, the Fe–3 at.% Cu system at 500 °C is considered, which has already been analyzed in previous studies [18], [19]. A simulation box with an edge length of 50 body centered cubic cells (2.5 × 105 lattice sites) and periodic boundary conditions are used. Atom migration proceeds with the vacancy exchange mechanism.

The atom/vacancy exchange is performed according to the Metropolis algorithm, i.e. the acceptance probability P of a single exchange event is given as P = exp(−ΔE/k), with k – Boltzmann constant, T – temperature (for ΔE < 0, P = 1). The total energy of the system, E, is represented by the sum of all pair-wise bond energies extending to the first (ε(1)) and the second (ε(2)) coordination sphere according ton represents the number of corresponding bonds. εFeFe and εCuCu are set to zero, thus creating the energy reference line. εCuFe was set such as to reproduce the solubilities of bcc-Fe and bcc-Cu phases at 500 °C as given from the thermodynamic assessments of Turchanin et al. [20]. The corresponding values are 0.12 at.% Cu in Fe-bcc and 0.26 at.% Fe in Cu-bcc. The asymmetry of the solubilities was reproduced by introducing a linear dependence of the εFeCu parameter on the local chemical environment of the neighbor shell around the atomic bond with

The value of is assessed with 4.61 kJ/mol, and it represents the CuFe bond energy in an Fe-rich chemical environment. is 3.15 kJ/mol, representing the same energy in Cu-rich environment. The actual value of bond energy according to Eq. (4) is evaluated after determining the site fractions x of Fe and Cu atoms in first and second coordination shell. Various combinations of εFeVa and εCuVa are examined in the simulations. The energies are set as ε(2) = ε(1) · (r2/r1)−6, with r being the interatomic distance to the atom in the i-th coordination sphere (for a bcc lattice (r2/r1)−6 is 27/64).

In the simulations, the real-time increment corresponding to a single exchange event is evaluated aswith f being the vacancy diffusion correlation factor (0.727 for bcc lattice [21]), a being the nearest neighbor distance in bcc-Fe (2.485 Å), D being the diffusion coefficient of the jumping atom A in bcc-Fe (1.5 × 10−22 m2 s−1 for Cu [22], 7.4 × 10−23 m2 s−1 for Fe [23]) and xVa,MC being the vacancy site fraction in the simulation box (xVa,MC = 4 × 10−6).

Cluster coagulation simulation

In this simulation approach, the solute atoms/clusters are treated as randomly moving gas particles, which can only grow by coalescence after collision with another particle. In each simulation step, the particles move in random direction and the probability for movement is evaluated as described below. The results of this calculation (the time evolution of the cluster sizes and numbers) are taken for comparison with the LMC simulations. The simulation of the coagulation process is performed on a system built in the same manner as in the LMC model, defined as cube with an edge length of 14.35 nm (50 bcc cells with cell parameter 2.87 Å). The system composition and temperature was set equal to the ones of the LMC simulation (3 at.% Cu, 500 °C).

Initially, all Cu atoms are placed randomly on sites corresponding to a bcc lattice. From that moment on, the atoms are treated as clusters, whereby atoms adjacent to each other with the distance lower than 2.485 Å (the initial nearest neighbor distance) are considered as belonging to the same cluster. Furthermore, the clusters are treated as spherical particles concentrated in their center of mass, and a radius r = (3πnV/4)1/3, with n being the number of atoms in the cluster and V being the Cu atom volume (V = 1.18 × 10−29 m3). In one time interval, the centers of mass of the clusters are assumed to be free to move in any direction (randomly generated vector components) by the distance of d = 0.25 nm and the probability P = n−4/3, where n is the number of atoms in the cluster (cluster size). This kinetics is representative for the model suggested by Binder and Staufer [16]. If clusters collide, i.e. the distance of their centers of mass becomes smaller than the sum of their radii, they are considered as one cluster afterwards. The center of gravity of the new (bigger) cluster is associated with the center of gravity of the cluster pair at the moment of collision. According to random walk theory, the time increment in these simulations is given bywith DCu = 1.5 × 10−22 m2 s−1 taken from experiment [22].

Continuum model

In the continuum simulation, the precipitation kinetics framework of the software MatCalc [6], [7] is utilized. The nucleation kinetics of precipitates is evaluated from classical nucleation theory [24]. Accordingly, the rate J at which new precipitates are formed in the system isN is the number of available nucleation sites, Z is the Zeldovich factor (lowering the stability of the critical nucleus due to the thermal vibrations), β is the atomic attachment rate, τ is the induction period, t is the time and G* is the energy barrier for formation of a critical nucleus. The latter is given aswith F being the driving force for precipitation and γ being the interfacial energy. The interfacial energy γ is obtained from the generalized broken bond model [25], [26] being dependent on the composition of both matrix and precipitate phases. A mean-field model based on the thermodynamic extremum principle [27] is used for the further evolution of precipitate size and chemical composition. The numerical time integration is performed according to the numerical Kampmann–Wagner approach [11], with an extension for consideration of cluster collisions, which is described below. Details on the models implemented in the MatCalc software are given in the corresponding references above.

Results and discussion

LMC simulation

In analogy to previous studies [2], [3] on the effect of the vacancy-solute binding energy on cluster evolution, the system kinetics is investigated in dependence of an asymmetry parameter a* defined asa* is a quantity describing the affinity of the vacancy to either the matrix-forming Fe atoms or the precipitate-forming Cu atoms/clusters. With negative values of a*, the vacancy prefers to be surrounded by as many Fe atoms as possible to minimize the total system energy by maximizing the number of FeVa bonds. Positive values of a* will maximize the number of CuVa bonds.

The corresponding system behavior can be quantified by inspecting the ‘local chemical environment’ LCE of the vacancy during its migration through the crystal on a statistical basis. Fig. 1 shows an analysis of the LCE for different values of a* starting with a random solution (full circles), representing the starting conditions in the simulations, and after Cu-clusters have formed with a mean radius of 1 nm. The LCE, defined here as the chemical composition of the first coordination sphere around the vacancy, has been measured 1000 times after every 1000 vacancy jumps. The figure emphasizes that, with increasing cluster sizes, the mean LCE becomes increasingly pronounced either Fe-rich or Cu-rich.

Fig. 1

Mean local chemical environment (LCE) around the vacancy for various values of the asymmetry parameter a* and different times during the precipitation process.

Fig. 2 shows the influence of a* on the entire precipitation process. The characteristics of the system are visualized by displaying the precipitate cluster sizes during the simulation in the form of “stair-fountain” diagrams (compare Ref. [3]). In this analysis, Cu atoms adjacent to each other as nearest neighbors are considered as a precipitate cluster if the cluster size exceeds nine Cu atoms. In both diagrams (Fig. 2), a continuous drop in cluster size indicates the occurrence of Ostwald ripening, where a small precipitate dissolves and the mean precipitate radius of all other particles increases continuously. For a* = 2, the sudden appearance of some considerably larger cluster sizes in the later steps is a clear manifestation of cluster collisions, since this sudden increase of an individual cluster size cannot be reasoned on basis of the monomer absorption/desorption process. The right image in Fig. 2 clearly shows that both mechanisms, monomer absorption/desorption as well as particle coalescence, can occur simultaneously, the dominant mechanism being controlled by the vacancy solute interaction energies εFeVa and εCuVa.

Fig. 2

Stair-fountain diagram of the cluster sizes observed during the Monte Carlo simulation for various values of the asymmetry parameter a*.

The effect of various a* values on the precipitation kinetics is investigated next by observing the number density and mean radius of precipitates throughout the precipitation process. From the MC simulations, we observe that negative values of a* inhibit the precipitation process compared to a* = 0 such that the corresponding evolution curves are shifted to longer simulation times. The general shape of the curves remains almost unaffected (Fig. 3a and b), though.

Fig. 3

Evolution of the particle number density and mean radius of the precipitates in the LMC simulation depending on the value of the asymmetry parameter a*.

Increasing the value of a* from 0 to +2 shifts the onset of the curves to the left, indicating an acceleration of the nucleation process (Fig. 3c and d). Increasing the a* value further, however, leads to a flattening of the curves, while the precipitation onset remains almost unchanged (Fig. 3e and f). This observation indicates that positive values of a* do not significantly affect the first stages of precipitation, however, they act as strong decelerators of the precipitation process in the later stages when larger clusters of Cu have already formed.

These features can be reasoned on basis of the vacancy preference toward either the Fe (a* < 0) or Cu (a* > 0) atoms. The precipitation process proceeds fastest for a* being approximately +1, which is interpreted as the condition for the most effective transport of Cu atoms through the matrix. Decreasing the value of a* leads to an increase of the vacancy presence towards the Fe atoms, delaying the migration of Cu atoms through the Fe-matrix and, thus, the entire precipitation process. Increasing the value of a* leads to vacancy trapping on the Cu atoms, which again retards their movement due to repeated Cu atom/vacancy exchange in the same spatial configuration.

The different effects of positive and negative a* values on the shape of the evolution curves are an indication of different operating mechanism. The horizontal shift of the curve without shape modification (in the logarithmic time scale) observed for a* < 0 suggests that the trapping of vacancies at the solvent Fe atoms slows down the precipitation process due to a reduced transport kinetics of Cu atoms. Since vacancies avoid the vicinity to Cu atoms, the probability of Cu atom/vacancy exchanges becomes increasingly lower with higher values of a*. Therefore, the effective transport of Cu atoms through the bulk volume is retarded and the integral precipitation process is slowed down. This process does not affect the shape of the precipitation curves since the dominant precipitation-controlling mechanisms remain the same.

The flattening of the curve observed for a* > 0 is attributed to the affinity and trapping of vacancies at Cu atoms. In the initial stages, the precipitation process is only slightly affected by this attractive behavior, however, in the later stages, the vacancies are more or less trapped inside the Cu clusters and they can, thus, not contribute to the long-range transport of Cu outside the clusters. This is necessary, however, to support the absorption/desorption processes that commonly drive the coarsening process. The precipitation process slows down and the number density and mean radius curves become flatter.

Simultaneously, the trapping of vacancies at Cu atoms increases the effective cluster mobility, since the vacancies stay in the vicinity of the clusters and cluster interfaces. The pronounced transport of Cu-atoms along the cluster surface facilitates the Brownian-type movement of the entire cluster, particle collisions are becoming more likely and coalescence of precipitates becomes more prominent. This effect is clearly visible in the stair fountain diagram in Fig. 2 and has already been discussed.

Movies of the precipitation process for different values of a∗ are provided as complementary material in the electronic version of this publication.

Cluster coagulation model

In order to study the characteristics of cluster coalescence isolated from the absorption/desorption mechanism, a coagulation simulation was performed as described above. The analysis of the evolution curves allows for identification of the system kinetics after some initial induction period (Fig. 4). From the results of this simulation, the number density is found to evolve in time within t−1/2 − t−3/5, while the mean particle radius evolves within t1/6 − t1/5. The mean radius dependence of t1/5 was predicted by Fratzl and Penrose [1], but these authors used the dependence D = D1n−1 in their analysis. Nevertheless, they showed also that the D dependence on the cluster size applied in this simulation would result in the mean radius growth law with t1/6. The time evolution of the number densities assures volume conservation of the solute atoms. The trends identified in this coagulation calculation are observed also in the Monte Carlo simulations with the a* ⩾ +2 condition (evolution curves from LMC simulation with a* = +2 are shown in Fig. 4 for comparison), confirming thus the migration of the entire clusters for those cases. The qualitative difference of the particle number trends (Fig. 4a) observed during the initial stage of the cluster coagulation simulation (particle number is high and constant) and the LMC simulation (particle number increases) is a result of the different particle definitions used in the both methods (cluster coagulation – all particles; LMC – precipitates larger than nine adjacent atoms only).

Fig. 4

Evolution of the particle number density and mean radius of the precipitates in the cluster coagulation simulation. The results of Monte Carlo simulation for a* = +2 are also shown.

The MatCalc continuum model used in this study describes the integral precipitation process by the evolution of discrete size classes rather than dealing with individual clusters. In the original formulation [11], particle growth or shrinkage is governed by the diffusional transport of monomers through the matrix. In the modified extended numerical algorithm, cluster coagulation can be introduced into the evolution equations by incorporating the collision processes into the evolution of the number density N of size class i.

In a simple, semi-empirical approximation, the rate of precipitate coalescence in size class i is proposed to be proportional to the total cluster density N, the mean cluster cross-sectional area represented by the square of the mean particle radius , the diffusion coefficient D of the clusters of size class i and a dimensionless proportionality constant k with

The cluster diffusion coefficient D for a given precipitate size class scales with the diffusion coefficient of the monomers as (r/r1)−4, with r being the radius of size class i, in agreement with the Binder–Stauffer model [16]. In the numerical integration of the evolution equations, the rate is evaluated in each time interval and the number of precipitates N of size class i is reduced accordingly. For the sake of computational efficiency, the amount of solute atoms corresponding to the volume of the collided particles is released homogeneously into the matrix and is available for attachment to the remaining particles. By this procedure, the mean radius and number density evolution of the entire precipitate population is reasonably well simulated although the individual collision processes are not treated explicitly.

The influence of the proportionality coefficient k introduced in Eq. (10) on the evolution of the precipitate number density and mean radius is investigated in Fig. 5. For k = 0, the curves represent the results of the traditional NKW model. The number density of precipitates increases in the beginning of the reaction as determined by the transient nucleation rate J given in Eq. (7). The mean radius, at this time, is determined by the critical nucleation radius given by classical nucleation theory.

Fig. 5

Evolution of the particle number density and mean radius of the precipitates in the continuum precipitation kinetics simulation depending on the value of k parameter in Eq. (10). Various growth law lines are also shown (see text).

Without precipitate coalescence, the number density reaches a plateau after supersaturation has decreased due to solute atom depletion of the matrix. Simultaneously, the mean radius increases due to precipitate growth and reaches a plateau for the same reason. The plateaus for both quantities finally evolve into the well-known LSW coarsening kinetics according to t−1 and t1/3.

If coalescence is accounted for according to Eq. (10) and the procedure described above, we observe a gradual decrease of the number density of precipitates after reaching a peak value. For higher values of k, both curves for the number density and the mean radius approach the time dependence observed earlier in the pure particle coalescence simulation. The corresponding lines are shown in the diagram. The graphs also show that the peak density of precipitates reduces, while the mean radius evolution is significantly enhanced. As expected, in the limit of long simulation times, the evolution of number density and mean radius becomes controlled again by the LSW kinetics. Moreover, a transition between the cluster collision and evaporation associated regimes can be observed, which is delayed with increasing value of the coefficient k. This is in agreement with the analysis of Fratzl and Penrose [1], who postulated such a transition to occur at some critical particle mean radius.

Fig. 6, finally, compares the results of the atomistic LMC simulations for a* = +2 with the extended continuums simulations with and without precipitate coalescence. The graphs show that, using a value of k = 0.5 for the present simulation conditions, the evolution of number density and mean radius computed by the atomistic simulation and the continuum model are in excellent qualitative agreement. The shapes of both curves can reasonably well be reproduced over the entire precipitation process including the nucleation, growth and coarsening stages.

Fig. 6

Comparison of the LMC simulation for a* = +2 with the classic (k = 0) and modified (k = 0.5) NKW model results, together with the limiting growth law lines.

Summary

In the present work, the influence of precipitate coalescence on the integral precipitation process is investigated using three different computational approaches (lattice Monte Carlo, coagulation simulation, modified NKW approach). In atomistic lattice Monte Carlo simulations, the effect of the vacancy/solute atom bond energy on the mobility of clusters as well as the overall precipitation kinetics is investigated. In agreement with previous work, we find that the vacancy/solute atom bond energy controls the vacancy preference toward the matrix or precipitate phase and, thus, alters the observed precipitation kinetics. Vacancy preference towards the solvent atoms generally retards precipitate growth by slowing down the monomer diffusion kinetics in the matrix. If the vacancy is trapped on the solute atoms, the early stage precipitation kinetics temporarily shifts from a dependency on monomer evaporation and absorption to a kinetics governed by precipitate cluster collision and coalescence, until, in the asymptotic limit, the classical LSW kinetics is observed again.

A numerical analysis of the cluster collision kinetics delivers a t1/6 − t1/5 growth law for the mean radius and t−1/2 − t−3/5 dependency of the particle number density for purely coalescence-controlled reactions.

In order to incorporate these phenomena in NKW-type continuum-scale modeling, a semi-empirical evolution equation for the number density of precipitate size classes is suggested. This model is implemented into the MatCalc precipitation kinetics framework and it shows excellent qualitative agreement with the lattice Monte Carlo results. Both, the evolution of the number density as well as the mean radius of the precipitates, are well reproduced.