Frédéric De Geuser
email: frederic.de-geuser@simap.grenoble-inp.fr

Brian M Gable

Barry C Muddle

CALPHAD based kinetic Monte-Carlo simulation of clustering in binary Al-Cu alloy

Keywords: age-hardening, aluminium alloys, Monte-Carlo, simulation Kinetic Monte-Carlo simulations, precipitation, Guinier-Preston zones, aluminium alloys

published or not. The documents may come

CALPHAD based kinetic Monte-Carlo simulation of clustering in binary Al-Cu alloy

Introduction

Age hardening aluminium alloys rely on the formation of a supersaturated solid solution that will decompose into strengthening solute-rich precipitates. Although phase diagrams predict the emergence of the stable phase(s), kinetic effects often lead to the formation of a sequence of metastable phases that are more readily nucleated. These metastable precipitates are not equally desirable in terms of the mechanical properties they attribute to the alloy. The precipitation sequence in a given aluminium system can be greatly influenced by the heat treatment temperature, the introduction of defects through deformation and the use of microalloying elements.

Monte-Carlo and kinetic Monte-Carlo simulations have been used to study clustering and precipitation phenomena in such alloys. These methods have proven very successful from a fundamental perspective in validating the application of classical nucleation theory in substitutional alloys [1] and in improved understanding of how kinetic aspects such as solute mobility and solute-vacancy binding influence precipitation [2][3][4]. Kinetic Monte-Carlo simulations are also used extensively to study irradiation induced segregation and precipitation [5,6]. However, applications of the approach commonly assume a fixed lattice and pair-wise interactions, and are likely to be limited to nearest-neighbour interactions, and such factors can be limiting if the simulations are to be universally predictive. Although sometimes possible [7,8], realistic and reliable energetic calibration through experimental measurements can also be problematic because of the need for overly simplistic approximations. Energetic calibration via ab initio calculations have often proven more reliable [9][10][11], but are more complex to perform.

In this paper, we describe an atomistic model that determines the energy of the system through both two-and four-body interactions between atoms, and introduce a method for deriving the input parameters from a thermodynamic database. The model is applied to kinetic Monte-Carlo simulation of the clustering of Cu atoms in an Al-Cu alloy. Details of the benefits and the necessity for employing multi-body interactions are also addressed.

The atomistic model

The system is represented by a cubic box of L 3 unit cells of aluminium (i.e. N s = 4 L 3 sites). The lattice is considered to be rigid and each atom is situated at a node of the FCC lattice. Vacancies are considered as "pseudo-atoms" in the sense that the vacancy-solute or vacancy-Al interactions are expressed in a similar form to interatomic interactions. The Al-Cu binary system is thus treated as a ternary Al-Cu-vacancy system.

The energy of the system is assumed to be well described by an Ising model restricted to first nearest neighbour interactions. This is considered reasonable as the focus is the early stages of clustering, when crystal structure changes are not significant. In order to assess the influence of multi-body interactions, we compare two different descriptions of the energy of the system: one based on a pair-wise interaction model and one on a four-body interaction model.

Pair-wise interactions model

In the pair-wise interaction model, the total internal energy of the system in a configuration (1) is given by :

pair E 1 = 1 2 n,m i,j ij p i n p j m ( 1
)
where i, j represent the different species present in the alloy (including vacancies); n, m represent the 12 nearest neighbours sites and ij is an effective interaction parameter. The occupation function p i n is equal to 1 if an atom of the type i occupies the site n, or 0 otherwise.

Four-body interactions model

In a four-body interaction model, we can write similarly :

4-body

E 1 = 1 4 n,m,p,q i,j,k,l ijkl p i n p j m p k p p l q (2)
Instead of counting the number of i -j pairs to calculate the energy as in the pair-wise model, we count the number of i -j -k -l tetrahedrons. The two models, represented schematically in Fig. 1, are equivalent if an i -j -k -l tetrahedron can be expressed by summing the effective pair interactions involved in the tetrahedron. Taking into account multi-site interactions is equivalent to considering the dependence of the atomic interactions on the local concentration and leads to more realistic energetic properties [12]. It does not necessarily imply a symmetric phase diagram as is the case with the pair-wise interaction model.

Kinetics

In order to study the kinetics of a diffusional transition, it is important to relate the transition rate, from the initial state of the system to the product, to the physical mechanism of diffusion. Some authors have used transition rates based solely on the energy difference between initial and final configurations (e.g. [8]), essentially performing equilibrium Monte-Carlo simulations. This is potentially misleading because the evolution of the system towards an energy minimum does not have any physical meaning. Differences in elemental diffusivities, for example, may not be taken into account. In the present case, atomic diffusion is assumed to occur only by a vacancy mechanism, in which a vacancy exchanges its position with one of its 12 nearest neighbours.

During a vacancy jump (i.e. an atom jump into a vacant site), we do not consider explicitly the trajectory between the two equilibrium positions. As the jump time is short compared to the residence time at equilibrium sites, we can express the transition rate between a pre-jump configuration (1) and a post-jump configuration (2) by [13] :

W 1→2 = ν exp - E act kT ( 3
)
where ν is an attempt frequency and E act is the activation energy or the migration barrier. It corresponds to the energy difference between the saddle point (the maximum energy position of the migrating atom) and the initial position of the atom, as shown in Fig 2: .

E act = E sp -E 1 (4)
Following [14], we assume that the saddle point energy can be written as the average energy of the two configurations 1 and 2, plus a constant Q i that depends on the nature of the jumping atom, so that :

E act = Q i + E 1 + E 2 2 -E 1 = Q i + ∆E 2 (5)

Derivation of the parameters

Other groups, examining the early stages of phase separation in aluminium alloys, have derived their input thermodynamic parameters either from experimental solubility [7,8,15] or from ab initio simulations [9,16]. Here we use CALPHAD thermodynamic databases to derive the effective interaction parameters. There are two advantages to this approach: 1) it does not require any assumption as to the phase(s) that are likely to form, because it considers only the solid solution of solute on a FCC aluminium lattice, and 2) since these energy expressions for Al and Cu have been optimized to generate the Al-Cu binary phase diagram, the CALPHAD description already accounts implicitly for the short and long-range interactions.

The derivation is based on that used in [12] for determining the Cahn-Hilliard gradient energy coefficient κ. For the sake of simplicity, we will first consider a binary alloy A-B (e.g. A = Al and B = Cu), on a FCC lattice (with 12 nearest neighbours) and x A and x B the atomic fractions of A and B (where

x A + x B = 1).
The method is based on the expression of the internal energy1 of the system when all the atoms are randomly distributed. It is important to notice that this does not mean that an undersaturated solution is necessarily random, as interactions between the atoms will create a short-range correlation. What the method does imply, however, is that the effect of this short range correlation on the free energy is mainly carried by the entropy term. This term does not have to be known explicitly as it will be automatically adjusted by the stochastic process of the Monte-Carlo simulation. In terms of CALPHAD data, the free energy of the system can be expressed by the sum of the free energy contribution E A and E B of the pure elements A and B respectively, an entropic term corresponding to an ideal mixing, and the excess free energy of mixing (per atom) that can be written with Redlich-Kister polynomials [12]:

F = x B E B + x A E A + kT (x B ln(x B ) + x A ln(x A )) + x A x B F mix = x B E B + (1 -x B )E A + kT (x B ln(x B ) + (1 -x B ) ln(1 -x B )) + x B (1 -x B ) n i=0 L i AB (x A -x B ) i = x B E B + (1 -x B )E A + kT (x B ln(x B ) + (1 -x B ) ln(1 -x B )) + x B (1 -x B ) n i=0 L i AB (1 -2x B ) i ( 6
)
where n is the order of the polynomial used and L i AB is the ith order interaction parameter (n = 2 is the highest order found in the ThermoCalc software database for the Al-Cu system). The latter parameters can be expressed as

L i AB = E i AB + T ∂L i AB ∂T
, with E i AB being the internal energy contribution. We further assume that the short range order contribution is mainly carried by the entropic terms, and we write the internal energy per atom as:

E = x B E B + (1 -x B )E A + x B (1 -x B ) 2 i=0 E i AB (1 -2x B ) i = x B E B + (1 -x B )E A + x B (1 -x B ) [ (E 0 AB + E 1 AB + E 2 AB ) + x B (-2E 1 AB -4E 2 AB ) + x 2 B (4E 2 AB ) ] (7)

Four-body interactions

According to eq. 2, the internal energy per atom of the system can be calculated by counting the number of A-A-A-A, A-A-A-B, A-A-B-B, A-B-B-B and B-B-B-B tetrahedrons. Each atom is involved in 8 tetrahedrons and each tetrahedron contains 4 atoms, so that there are two tetrahedrons to be counted per atom. On a FCC lattice, a nearest neighbour tetrahedron is symmetric and A-A-A-B is equivalent to A-A-B-A or any other combination. The energy is then related to the probability of finding each tetrahedron: which can be written

4-body E = 2(x 4 A AAAA + 4x 3 A x B AAAB + 6x 2 A x 2 B AABB + 4x A x 3 B ABBB + x 4 B BBBB ) (8)
4-body E = x B 2 BBBB + (1 -x B ) 2 AAAA + x B (1 -x B ) [ 2(4 AAAB -3 AAAA -BBBB ) +x B 2(6 AABB -8 AAAB + 3 AAAA -BBBB ) +x 2 B 2(4 ABBB -6 AABB + 4 AAAB -AAAA -BBBB ) ] (9)
which,as for eq. 7,is in the form:

E = x B E B + (1 -x B )E A + x B (1 -x B )E mix (10)
where

x B E B + (1 -x B )E A = x B 2 BBBB + (1 -x B )
2 AAAA is the internal energy of an ideal solution and is equal to the energy of a fraction x B of atoms of B and a fraction (1 -x B ) of atoms of A physically separated. E mix is a measure of the nonideality of the solution: when E mix > 0, adding B atoms to the solid solution is not energetically favourable (i.e. increases the energy) and there will be a miscibility gap in the phase diagram. Both eq. 9 and eq. 7 are of the same order and it is now possible to derive the various interaction parameters by solving a system of equations corresponding to each term of the polynomials. The determination of all ijkl is then straightforward. Their values are given in appendix A. In particular, we have for the pure elements:

AAAA = 1 2 E A BBBB = 1 2 E B ( 11
)
It is of particular interest to define the B-B binding energy, which corresponds to the energy of a solute B-B pair. This represents the energy that one needs to separate two B atoms in solid solution. A balance of the two situations (2 isolated B atoms and 2 clustered B atoms) gives the expression of this energy:

4-body E binding B-B = 2 (2 AAAB -AAAA -AABB ) = 1 6 E 0 AB + 3E 1 AB + 5E 2 AB ( 12
)

Pair-wise interactions

In the pair-wise interaction model, there are 6 pairs per atom and counting the number of A-A, A-B and B-B pairs gives:

pair E = 6 x 2 A AA + 2x A x B AB + x 2 B BB (13)
which, by inserting x A = 1 -x B , can be written:

pair E = x B 6 BB + (1 -x B ) 6 AA + x B (1 -x B ) 6(2 AB -AA -BB ) (14)
Although both eq. 14 and eq. 7 follow the same pattern as eq. 10, namely an ideal term and an excess mixing term, a direct identification of the coefficients is not possible here. Whereas in eq 14 the excess energy term is constant, in the CALPHAD description E mix depends on the concentration x B . Equation 7 needs to be truncated in order to be used for calibrating the pair-wise interaction model. We will now describe two possible approaches for this truncation and thus approximation.

Approximation method I

A first approach involves keeping only the 0th order term of the excess energy of mixing from eq. 7. This gives :

AA = 1 6 E A BB = 1 6 E B ( 15
)
for the pure elements, and

AB = 1 12 (E A + E B + E 0 AB + E 1 AB + E 2 AB ) (16)
for the A-B cross terms. Similarly to eq. 12, we can calculate the B-B binding energy which reads :

pair E binding (I) B-B = 2 AB -AA -BB = 1 6 (E 0 AB + E 1 AB + E 2 AB ) (17)
The attraction of this approach is that the energy terms for the pure elements are untouched by the approximation and, as in the case of the for the four-body model, do not depend on terms originating from the second element. However, it would be more valid mathematically to expand the entire right-hand side of eq. 7 rather than only the term for the excess energy of mixing and truncate the full expression at the appropriate order. This is the basis of the second approximation method.

Approximation method II

For this approach, eq. 7 can be expanded to:

E = E A + x B (E B -E A + E AB + E 1 AB + E 2 AB ) + x 2 B (-E AB -3E 1 AB -5E 2 AB ) + x 3 B (2E 1 AB + 8E 2 AB ) + x 4 B (-4E 2 AB ) (18)
and eq. 14 to :

pair E = 6 AA + x B 6(2 AB -2 AA ) + x 2 B 6( AA + BB -2 AB ) (19)
Truncating eq. 18 to the quadratic term leads to : for the pure elements, and

AA = 1 6 E A BB = 1 6 (E A -E B + 2E 1 AB + 4E 2 AB ) (20)
AB = 1 12 (2E A -E B + E 0 AB + 5E 1 AB + 9E 2 AB ) (21)
for the A-B interaction term. Although this second approach is more rigorous mathematically, it is noteworthy that the effective interaction parameter BB for pure element B contains terms related to element A, which is in constrast to the outcome for the four-body model, and perhaps initially couter-intuitive.

However, adopting this approach, we can then write the B-B binding energy :

pair E binding (II) B-B = 1 6 (E 0 AB + 3E 1 AB + 5E 2 AB ) (22)
which is identical to the B-B binding energy for the four-body model, and which may thus facilitate a comparison between the two. The influence of the binding energy on the behaviour of the solution will be discussed in section 4.

An interesting feature of both approximation methods is that the terms that are neglected are Redlich-Kister polynomial terms of order >0. This means that for a Redlich-Kister polynomial of order 0 (i.e. a regular solution), the pair-wise model and the four-body model are equivalent for both approximations. If the polynomial is of higher order, the multi-body interactions will have an influence.

Interaction with the vacancy

The ij and ijkl involving a vacancy can be deduced from the formation energy of a vacancy in the solvent, E form vac and the binding energy between a solute and a vacancy, E binding i-vac . These energy terms can be written as functions of ij or ijkl by (i) expressing the energy needed to take an A atom from a position close to i to the surface for E form vac , and (ii) taking the energy difference between a system in which a vacancy is next nearest neighbour to a solute atom and a system where they are far apart, thus:

Pair E form vac = 12 Avac -6 AA 4-body E form vac = 8 AAAvac -6 AAAA (23)
and

Pair E binding i-vac = Avac + Ai -AA -ivac 4-body E binding i-vac = 2 ( AAAvac + AAAi -AAAA -AAivac ) (24)
It is further assumed that the effective four-body interactions involving a vacancy can be expressed as the sum of the pairwise effective interactions composing the four-body tetrahedron. This is equivalent to assuming that the four-body model is equivalent to the pairwise model for any given tetrahedron in which a vacancy is involved. For instance, we can write:

AAAvac = 1 2 (3 AA + 3 Avac ) (25)
The explicit expressions for all terms are listed in Appendix A. The kinetic parameters can be derived from diffusivity data. When ∆E = 0 (Fig. 2), the activation energy from eq. 5 is equal to the activation energy for tracer diffusion (or self diffusion) Q i of element i corrected by the formation energy of a vacancy in the neighbourhood of the atom i, so that:

Q i = Q i -(E form vac -E binding i-vac ) Pair Q i = Q i -(11 Avac + 5 ivac -5 AA -Ai ) 4-body Q i = Q i -(6 AAAvac + 2 AAivac -4 AAAA -2 AAAi ) (26)
with i being either A or B.

The attempt frequency ν can be related to the pre-exponential factor D i 0 of the diffusivity:

ν = D i 0 f i a 2 ( 27
)
where a is the lattice parameter and f i is the correlation factor related to the diffusion of species i. The latter takes into account the fact that successive atomic jumps are not independent (the reverse jump having a higher probability than jumps in any other directions). In dilute alloys, the correlation factor can be calculated using the 5-frequencies model [17]. This model requires that, in a dilute alloy, the vacancy is either neighbouring a solute atom or not. In the first case, it can jump to another equivalent position (with an exchange frequency W 0 ), or become a solute neighbour (W 4 ). In the second case, it can either jump to another solute neighbouring site (W 1 ), exchange with the solute itself (W 2 ), or leave the solute (W 3 ). The correlation factor can then be expressed by [17]:

f i = W 1 + 7 2 W 3 W 1 + W 2 + 7 2 W 3 (28)
Through eqs. 3, 5 and 27, we can write :

W 1 = D A 0 f A a 2 exp - Q A kT , ( 29a
)
W 2 = D i 0 f i a 2 exp - Q i kT , ( 29b
)
W 3 = D A 0 f A a 2 exp - Q A -E binding i-vac kT (29c)
so that

f i = 1 -f A D i 0 D A 0 exp -Qi-QA kT 1 + 7 2 exp E binding i-vac kT (30)
D Al 0 [20] Q Al [20] D Cu 0 [20]
Q Cu The simulation volume contains 64 3 lattice sites (128 3 for some simulations) with periodic boundary conditions (so that every site has 12 nearest neighbours). We use the rhombohedral primitive cell of the FCC aluminium lattice, which is easier to handle from an algorithmic point of view. As a consequence, the simulation volume is not a cube but a rhombohedron. However, for clarity, when showing snapshots of the simulation box, it is positioned within a cube defined by base vectors parallel to the < 001 > FCC directions, in order to more easily identify the orientation of the clusters. One single vacancy is kept in the box for all times. This apparent vacancy concentration has been taken into account by rescaling the simulated time accordingly [9]. The validity of this approach has been confirmed by comparing two simulations with different sample volumes (results not shown here). Depending on the desired simulation time and the amount of on-the-fly data processing, the simulations lasted for real times ranging from several minutes to a few days.

This model has been used to study the clustering of an Al-1.7at%Cu alloy during isothermal ageing at 100 • C. The starting point of each simulation is a random solid solution and we assume the vacancy concentration to be close to that which would be observed in pure aluminium at a reasonable solution heat treatment temperature of 530 • C. In other words, we assume that all the vacancies have been quenched in from the solution heat treatment temperature with no subsequent annihilation. This approximation is the main source of discrepancies between the real physical time and the simulated time. The amount of quenched-in vacancies and the rate at which they annihilate makes it difficult to consider the simulated time scale as absolute.

The numerical values for the Cu-Cu binding energies are

pair E binding (I) Cu-Cu ≈ -0.024eV 4-body E binding Cu-Cu = pair E binding (II) Cu-Cu ≈ 0.118eV (31)
It is noteworthy that, in the case of the pair-wise interaction model, the sign of the binding energy changes for the 2 approximation methods. The implication of these different Cu-Cu binding energies will be discussed in detail in a following section.

Pair-wise interactions (approximation method I)

In the approximation method I, the Cu-Cu binding energy is negative which corresponds to a situation where there is no miscibility gap. All Cu atom will tend to be surrounded by Al atoms, even more so than in an ideal solution. To characterise this, it is possible to count the average number of Cu atoms surrounding each Cu atom (in a first nearest neighbour position) n Cu 1 . In a perfectly random solid solution, this number will be equal to Z 1 C 0 , with Z 1 , the first coordination number (or number of first nearest neighbours, i.e. 12 for a FCC lattice) and C 0 , the overall Cu concentration of the alloy (C 0 =1.7%). At the beginning of the simulation,

n Cu 1 (t = 0) = Z 1 C 0 = 0.204.
As the alloy is dilute, we can make the assumption that Cu atoms exist in no bigger clusters than pairs, so that Cu can only be found as either a single Cu atom surrounded by 12 Al atoms, or a 2-atom cluster. If this is the case, the formation and dissolution of these pairs are thermally activated processes, and the equilibrium number of Cu-Cu pairs is determined by the energy needed to create a Cu-Cu pair (-E binding Cu-Cu ). The equilibrium value of n Cu 1 is given by:

n Cu 1 (t = ∞) = Z 1 C 0 exp E binding Cu-Cu kT (32)
In the case of the approximation method I, as the binding energy is negative, n Cu 1 (t = ∞) < Z 1 C 0 . The formation of a pair of Cu atoms costs energy, so that there will be fewer CuCu pairs than in a random solid solution. It is expected that n Cu 1 will decrease from Z 1 C 0 = 0.204 to its equilibrium value of

Z 1 C 0 exp E binding
Cu-Cu /kT = 0.097, which is what is shown by the simulation (Fig 3). Using pair-wise interactions in the approximation method I, the simulations demonstrate cluster dissolution rather than aggregation. The absence of a miscibility gap, arising from this approximation method, is of course, contrary to common theoretical and experimental evidence of phase equilibria for this alloy (e.g. [21]). This trend for dissolution indicates that the higher order CALPHAD terms are necessary for the phase energy description to replicate the binary Al-Cu phase diagram.

Pair-wise interactions (approximation method II) and four-body interactions

For the approximation method II and for the four-body model, the Cu-Cu binding energy is positive, indicating that there is a clustering tendency. In this case, the assumption that only isolated Cu atoms and Cu-Cu pairs existed is not valid, as clustering is very likely to occur. Figure 4 provides a sequence of 3D atom maps over the simulation volume calculated on the basis of both the pair-wise model in the approximation method II (above) and the four-body model (below) for three different simulation times up to 10 7 s (around 4 months).

For the sake of clarity, not all the atoms within the simulation volume are represented. We show only Cu atoms, and we represent non-isolated Cu atoms (i.e. those Cu atoms that have at least another Cu atom as a first nearest neighbour) by larger spheres. The viewing direction is slightly misoriented from the [100] Al direction. It is clear from the figure that both descriptions lead to the clustering of the Cu atoms. In the case of the pair-wise model (II), spherical Cu-rich clusters appear. In the four-body model, the Cu atom clusters adopt a disc-like shape along the 001 Al planes, and the disc thickness is one atomic plane. This monolayer aggregate of Cu atoms is the generally accepted structure of so-called Guinier-Preston (GP) zones [22,23].

Although the clustering tendency is the same (i.e. the binding energies are equal), the difference in behaviour between the two models is a direct consequence of the limitation of the pair-wise model in this context. In this model, the interactions between the atoms do not depend on their environment. The four-body model, on the other hand captures this dependence in the sense that a pair of Cu atoms has a different energy depending upon whether it is surrounded by other Cu atoms or by Al atoms.

This can be illustrated by considering the energy of a cluster comprising 3 nearest neighbour Cu atoms. The energy of such a cluster can be computed by comparing the situation where the 3 Cu atoms are isolated to the situation where they are clustered. In a pairwise model, the energy of such a cluster is simply the sum of the energy of three Cu-Cu pairs. In a four-body model, the three pairs influence each other and this cluster is demonstrably less favourable. The energy descriptions for these phenomena are quite complex and can not be simply defined in terms of the sum of pair-wise interactions.

Although no long-range interactions are explicitly defined in the model, the improvement in the simulations achieved by considering multi-body interactions is dramatic as the four-body model is, indeed, able to predict the formation of GP zones.

. Kinetics of GP-zone formation

We can now follow more closely the evolution of the clustering reaction for simulations based upon the four-body interaction, for it is possible to monitor the residual matrix Cu content, the number density and the average size of the GP-zones as a function of simulation time. In order to be able to compute these three parameters, we need to define a minimum cluster size below which a cluster is assumed to be part of the matrix. As the cluster size distribution does not show any discontinuity that we could exploit in order to determine the critical size of a GP-zone, we have chosen to use a value derived from the critical nucleus size for a Al-1.7at%Cu alloy, obtained by a modified embedded-atom method [24]. This previous work defined a critical radius of 0.4nm, which corresponds approximately to 10 Cu atoms. This critical size has been obtained for the initial Cu concentration and will then change as the supersaturation of the matrix decreases. It is not possible to account for this effect, so that simulations will eventually under-estimate the critical size of the GP-zones. However, this should not affect the results significantly, because as the matrix concentration decreases, the probability of a 'random' matrix cluster reaching the initial critical size decreases significantly.

Figure 5 describes the evolution of the solute content of the matrix, the number density and the average size of the GP zones as a function of time. The overall formation sequence for the GP zones can be conveniently, if somewhat arbitrarily, described in three stages. Although there are no abrupt transitions between these different stages, it is possible to define limits at key points in simulation time: i.e. when the matrix solute concentration first starts to decrease significantly (Fig. 5(a)), when the number density of clusters reaches a maximum and when the cluster size starts to increase significantly, as in Figs. 5(b) and 5(c), respectively.

During the first stage, only the number density changes significantly. Nuclei are created that all have approximately the same critical size. The overall volume fraction is not sufficient to significantly affect the matrix composition. This might be regarded as the nucleation stage.

When the nuclei start to grow, they increasingly consume the matrix Cu, so that the matrix solute composition decreases significantly. As new critical nuclei are created at the same time, the average size does not increase markedly. As the matrix Cu content decreases, so does the supersaturation and the nucleation rate, so that the rate of increase in the number density slows. This region constitutes a stage of combined nucleation + growth.

When the supersaturation of the matrix becomes very low, nucleation essentially ceases. Two phenomena occur that give rise to a steep increase in the average size of the GP-zones: the GP-zones tend to grow as driven by the remaining supersaturation, and the requirement to minimise interfacial energy sees the smaller GP-zones dissolve sacrificially to benefit the larger ones, resulting in a decreasing number density. This is a regime of combined growth + coarsening.

Temperature dependence

We have compared microstructures obtained for different ageing temperatures. In Fig. 6(a), the number density of GP zones is shown as a function of ageing time for four different temperatures ranging from 100 • C to 140 • C. It is clear from this figure that increasing the temperature not only accelerates the kinetics of the reaction, but also decreases the maximum number density reached by the system. This is due to an increase in the solute solubility and a consequent decrease in the An alternative representation would see the number density of zones plotted as a function of time multiplied by the ratio of the Cu diffusivities at a given temperature T and at 100 • C. Indeed, temperature can influence the reaction in two ways: by increasing the diffusivity and by changing the level of supersaturation and thus the corresponding driving force. Re-scaling the plot will suppress the effect of the increased diffusivity. This rescaled representation is shown in Fig. 6(b). It is observed that, although the maximum densities differ, the four curves overlap nearly perfectly. This confirms that the nucleation rate is largely transport controlled, and that the effect of the diffusivity is more important than that of the supersaturation. Indeed, the maximum number density reached by the alloy aged at 130 • C is only about 10% lower than that of the alloy aged at 100 • C.

Figure 7(a) shows the evolution of the average size (number of atoms) of the GP zones for different temperatures. Although the temperature does affect the growth rate of the particles, if the time is re-scaled to account for the diffusivity differences, the four curves are a near perfect match (Fig. 7(b)). From Fig. 6(b) and Fig. 7(b), we can conclude that the variation in supersaturation induced by changing the temperature from 100 • C to 120 • C does not affect the growth rate of the particles and only slightly modifies the number density. This is particularly noteworthy because, if we assume that mechanical properties are related to the number density of particles (of comparable size), it means that an increase in the ageing temperature of only 20 • C will accelerate ageing kinetics by an order of magnitude to deliver very similar results in terms of properties.

Discussion

Albeit simplistic in terms of interaction range, the present model is demonstrated to predict the formation of coherent, monolayer solute Cu zones that bear close similarity to Guinier-Preston zones. Although the form of such zones is usually attributed to elastic interactions, it is not necessary to explicityly introduce such longer-range interactions in the model. However, it should be noted that this does not mean that the formation of GP-zones can be attributed solely to chemical interaction between the atoms. While the model does consider the lattice as being fixed, so that the distance between atoms does not change locally (i.e. the atom positions are assumed to be the node of a perfect FCC lattice), the method of calibration of the model through semi-empirical CALPHAD thermodynamic potentials implicitly takes into account all forms of interaction between the atoms, including the elastic energy component. In any case, the formation of such elongated objects in an atomistic simulation with an interaction range limited to first nearest neighbour is fundamentally instructive.

The limitations of the model are severe enough to limit the quantitativity of the simulations. However, we believe that, through its simplicity, this model is an interesting prototype for the nucleation of GP-zones type objects.

The nucleation rate is determined primarily by the diffusivity, and by the local supersaturation, and it is very likely that it will not be severely affected by the limitation of the interaction range. The same is true for the growth regime. The coarsening regime, however, is very different. In the Lifschitz, Slyozov, Wagner (LSW) treatment [25,26], the supersaturation is negligible and the driving force for coarsening arises from the tendency of the system to minimise the interfacial energy. The local solute solubility in the matrix is higher around those precipitates with smaller radii of curvature, which tend to dissolve, and lower around larger precipitates, which tend to grow. If generalised for long range elastic interaction, this phenomenon might be even more effective, since the long-range chemical and strain environments of neighbouring precipitate may overlap. The lattice distortion around precipitates of different size is also likely to be very different. For larger precipitates of the same form as the GP zones (e.g. θ phase in the same system), it is also known that the presence of a precipitate can favour the nucleation of another in its own lattice distortion field (so-called autocatalytic nucleation [27]). By limiting the range of interactions to nearest neighbours, such effects, if they exist, cannot be taken into account within the current model. The special value of such simulations lies in the fact that, at these very early stages (i.e. at the nucleation stage), the scale of the features (typically 10 atoms) makes it very difficult, if not impossible, to observe experimentally.

Conclusion

We have demonstrated that the use of CALPHAD energy descriptions in kinetic Monte-Carlo simulations for the study of early stages of solute clustering in Al-Cu alloys is possible, and that a multi-body interactions model can lead to realistic simulations of observed behaviours. Pair-wise interactions are too simplistic to describe the clustering behaviour of the solute atoms, even at the early stages of clustering, as the example of Al-Cu has shown. Four-body interactions, limited to the first nearest neighbour distance, are sufficient to predict the formation of monoatomic layers of Cu atoms on the 001 planes of aluminium, which is the accepted form and structure of GP zones. Because of the range of the interactions, the validity of the description, however, does not extend to advanced stages of growth. Nevertheless, it gives an instructive prototype for the nucleation of this kind of coherent planar object.

An advantage of the model is that it can potentially be extended to any other solute of interest, provided the corresponding thermodynamic database and diffusivity values are reliably available. In particular, it can potentially be used as a design tool for the study of the influence of microalloying elements. For this purpose, however, it is limited to the early stages of clustering, when long range interactions can be neglected. It can, however, be used on a broader fundamental study of the role of microalloying elements and ageing temperature on the subsequent precipitation. These aspects are currently being explored.

F

3 Figure 1 .

 31 Figure 1. Schematic comparison of (a) the six different pairwise, and (b) the four-body atomic interactions.

4 Figure 2 .

 42 Figure 2. Schematic description of the energies involved in the calculation of the transition rate through eq (3).

4 .

 4 Influence of multi-body interactions

Figure 3 .

 3 Figure 3. Average number of Cu atom neighbours for each Cu atom, calculated according to the pair-wise model in the approximation method I.

13 Figure 4 .

 134 Figure 4. 3D images of the simulation box using both the pair-wise model (approx. II) (above) and the four-body model (below) after different ageing times at 100 • C. All Al atoms are omitted for the sake of clarity and clustered Cu atoms are represented with larger spheres. The black cube corresponds to the < 001 > directions of the Al FCC lattice. The empty corners of the cube are due to the rhombohedral shape of the simulation box and are no artefacts. Both insets on the 10 7 s volumes show a typical particle in a different orientation, emphasizing their spherical shape (pair-wise) or disk shape (4-body).

  supersaturation of the matrix for a given initial concentration.

Figure 5 .

 5 Figure 5. Reaction kinetics at 100 • Cshowing evolution of (a) the matrix solute concentration, (b) the GP zone number density, and (c) the average size of the GP zones.

3 )Figure 6 .Figure 7 .

 367 Figure 6. Evolution of the number density of GP zones as a function of time for different temperatures. (a) Normal time scale. (b) Re-scaled time scale accounting for the diffusivity difference.

Table 1

 1 includes the numerical values used in the present model and Appendix A provides the expressions for the actual interaction energies.

	une 28, 2010	13:12	Philosophical Magazine	kmc
		10
			F o r
			P
				e
				e r
				R
				e
				v i e
				w
				O n l
				y

Table 1 .

 1 Summary of the energy parameters (in eV) and diffusivities employed. See Appendix A for expressions for the interaction energies.

	E form vac [18] E binding Cu-vac [19]
	E 2 AlCu
	E 1 AlCu
	E 0 AlCu
	E Cu
	E Al

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

It is assumed that the enthalpy and the internal energy (and so the Gibbs free energy and the Helmholtz free energy) are equal, as the pressure term is usually negligible in metals.

Acknowlegdments

The authors gratefully acknowledge the support of the Australian Research Council (ARC) and Monash University. One of the authors (BCM) is recipient of an ARC Federation Fellowship. The authors also acknowledge gratefully the assistance of Dr C.R. Hutchinson with application of the CALPHAD database and of E.A. Lass with derivation of the interaction parameters.

Appendix A. Expression of the pair and 4-body energies used in this work

The general expressions of the pair and four-body energies used in this work are given here, in terms of CALPHAD data, vacancy formation energy in pure solvent (E form vac ) and solute-vacancy binding energy (E B-vac ) We assumed during the derivation of the expression that the pure elements free energy Redlich-Kister polynomials were of the order 0, the mixing free energy Redlich-Kister polynomials for solutes in the solvent were of the order 2 and the mixing free energy Redlich-Kister polynomials for solutes in one other were of the order 1. This was the case for our solvent (Al) and the considered solutes (Cu, Mg, Ag,...). Although we introduce a vacancy-vacancy binding energy (E vac-vac ), we have considered its value to be zero, which has no influence on the results, since there was only one vacancy at a time in our simulation box.

A.1. Pair wise model (approximation method I)

A.2. Pair wise model (approximation method II)