Simulating Diffusion Induced Grain Boundary Migration in Binary Fe–Zn

: A recently developed phase-ﬁeld model for simulating diffusion-induced grain boundary migration (DIGM) is applied to binary Fe–Zn. The driving force for the boundary migration is assumed to come from the coherency strain energy mechanism suggested by Sulonen. The effect of the angle of the grain boundary with the surface on the velocity of the boundary migration is studied in detail. The simulation results compare favorably with experimental observations, such as the oscillatory motion of the grain boundary, velocity of the moving grain boundary during DIGM, and the maximum value of mole fraction of Zn at the surface after 20 h of heat treatment.


Introduction
The discontinuous precipitation (DP) reaction, where a grain boundary constitutes the growth front behind which alternating lamellae of the precipitating phase and a solute impoverished matrix phase grow, was for a long time discussed as a recrystallization reaction [1].As a consequence, the mechanism of grain boundary motion in DP reactions was not an issue.Somehow this view seems to have affected the early theories of DP [2,3] although it was by that time realized that grain boundary diffusivity was an important factor.Hillert [4] concluded that there would be no driving force for grain boundary motion in DP unless there were a deviation from the local equilibrium between the precipitating phase and the boundary.In a short note, Sulonen [5] suggested that a zone in front of the advancing boundary should form where the lattice misfit due to the inward diffusion of solutes in the bulk from the boundary leads a to coherency strain energy, acting as the driving force for grain boundary motion.Later it was also suggested that the boundaries moved due to a pulling effect exerted by the precipitates [6].Some recent studies have shown that the grain boundary movement in similar processes exhibits jerky motions [7] and corresponding simulations have been performed to explain DP reaction with reaction front tracking and eventual composition changes in the microstructure [8].There are other factors that could influence the grain boundary motion, one such factor could be the formation of grain boundary phases [9,10] which may eventually control grain growth.In this study, the formation of such phases is neglected for simplicity.
When further investigating the role of grain boundary diffusion of zinc during DP in Fe-Zn alloys, Hillert and Purdy [11] found that grain boundaries in iron samples, when exposed to an atmosphere of zinc and without forming a precipitate phase, moved.They called the phenomenon Chemically-Induced Grain Boundary Migration and suggested that the driving force for the boundary motion had a chemical origin.Nowadays, the phenomenon is often called DIGM (Diffusion Induced Grain Boundary Migration) [12][13][14][15][16][17][18] and is considered to be an important component of the DP reaction.
To model DIGM the phase-field method is a suitable choice.Cahn et al.
[14] suggested a phase-field model for DIGM, with which they could study different mechanisms for the driving force.They concluded that a driving force based on the coherency strain energy is ". . .consistent with the observed features of DIGM in nearly all cases" [14].The purpose of the present paper is to test this prediction on a real case, the binary Fe-Zn system, using a recently developed phase-field model for simulating DIGM by Mukherjee et al. [19].

Phase-Field Model for DIGM in Fe-Zn
The model was described in detail in [19] for a hypothetical binary system.For the convenience of the reader, the main parts are repeated here for the Fe-Zn case.

Governing Equations
The basic expression of the phase-field method is the total Gibbs energy of the system (G), which is written as a functional for a binary system, such as Fe-Zn: where g 0 is Gibbs energy per mole, κ and are the gradient energy coefficients for gradients in mole fraction (x Zn ) and phase field variable φ, respectively, and V m corresponds to the molar volume.Here, 0 ≤ x Zn ≤ 1 and 0 ≤ φ ≤ 1.
The equations governing the evolution of mole fraction and phase field variable are the Cahn-Hilliard equation [20], Equation ( 2), and the Allen-Cahn equation [21], Equation (3), respectively: ẋZn In Equation (2), L is a kinetic parameter related to the diffusional mobility and M φ in Equation ( 3) is also a kinetic parameter but related to the boundary mobility.The g 0 in Equation (1) will be divided into two parts, i.e., g 1 and g 2 , where g 1 is the weighted average of the molar Gibbs energy of the phases present locally while g 2 is a function of mole fraction and phase field variable.In addition, the gradient term κ(∇x Zn ) 2 is dropped out due to its insignificance in the model.However, the gradient in the phase field variable (φ) has been taken into account due to the nature of our simulation setup.We consider only a single phase ferrite, α, in our simulations, thus the term g 1 then becomes: The term g 2 contains information about additional interfacial energy contribution apart from the gradient term in φ in Equation (1), and the strain energy contribution which will serve as a driving force for DIGM.Therefore, g 2 can be divided into g 21 and g 22 .The term g 21 maintains the shape of the interface inspired by Finel et al. [22], see Equation (5), Here, λ refers to the energy barrier coefficient given by 2dσ which is similar to in Equation (1), σ is the interfacial energy in J/m 2 , and Θ is given by, where d is the grid spacing and w is controlling the width of the interface.The term g 22 contains the strain energy contribution due to the incorporation of Zn atoms in either of the two grains, as a result of diffusion of Zn atoms along the boundary.
The term was proposed as an interaction between concentration and phase field variable by Cahn et al. [14] and later modified by Mukherjee et al. [19] to model DIGM in a hypothetical system, which for the Fe-Zn system can be written as: Here, x 0 Zn is the mole fraction of Zn in the grain in front of the moving grain boundary which is not yet swept by it, x Zn is the mole fraction of the area of the grain swept by the grain boundary, the term φ(1 − φ) suggests that it is confined to the interface and the pre-factor K was given by Hillert [13] as, where Y refers to Young's modulus of α-Fe.This, however, is supposed to change due to the diffusion of Zn into the matrix phase.Therefore, Vegard's law was used to calculate Y for x Zn = 0.1 (corresponding to the maximum value of x Zn observed at the sample surface in Ref. [23]).Due to a very small difference in the Poisson's ratio ν for pure Fe and the alloy formed with x Zn = 0.1, the ratio was taken for α-Fe, i.e., 0.29 and the lattice misfit η is given by, Here, d(ln a) corresponds to the difference in the natural logarithm of the lattice parameter of pure Fe and the alloy.The lattice parameters for pure Fe and pure Zn were taken from [24,25] and Vegard's law was again used to calculate the lattice parameter of the alloy formed due to Zn diffusion with x Zn = 0.1.The interaction term G int is responsible for grain boundary migration, which is triggered when Zn atoms from the vapor phase diffuse, much faster than the bulk diffusion, along the grain boundary.

Simulation Setup
The simulations performed are compared with the experiment on the Fe-Zn system by Chongmo and Hillert [23], where pure Fe samples were exposed to a vapor of Zn from turnings of three different alloys: Fe-33wt%Zn, 18.8wt%Zn and 9wt%Zn.The alloy with 18.8wt%Zn and the temperature 873 K were chosen for our simulations.According to microprobe measurements in [23], the mole fraction Zn on the sample surface was close to the solubility limit for Zn, i.e., x Zn = 0.1.The initial microstructure was set up such that the grain boundary extends from the top of the surface to 10 µm below the surface, a schematic of such a setup was shown in our previous work [19].
The enhanced diffusivity in the boundary was controlled by a factor F which simply was multiplied by the bulk diffusivity.One of the most striking results by Hillert and Purdy [11] was that the diffusivity of moving grain boundaries was up to four orders of magnitude larger than for stationary boundaries.This effect we cannot account for in the simulations, and the same high diffusivity was used for the whole simulation time.The factor F, however, will be affected by the thickness of the grain boundary.As per Kaur et al. [26], the boundary diffusivity is inversely proportional to its thickness.Since the boundary thickness for metals or metallic alloys is generally between 0.5 nm to 1 nm the factor F is changed from 10 4 to 200 corresponding to the thickness of the boundary used in our simulations, which is of the order of grid spacing, i.e., 0.1 µm.The thickness of the grain boundary is controlled by w in Equation (6).In our simulations, the value of w is taken, such that d/w = 2, where, d is the grid spacing.
The kinetic parameter M φ in Equation ( 3) is related to the boundary mobility M by where ξ = δ (dφ/dz) 2 dz, and it can be shown that ξ = 0.235/δ [27].The boundary mobility M for α-α grain boundaries in Fe has been suggested to be [28]: k M is a fudge factor that was set to around 200; it turned out that the value obtained from Equation (11) was too low to yield a reasonable agreement with experimental results.
In our simulations, σ, the interfacial energy, is used as a fitting parameter.All parameter values used are collected in Table 1.The thermodynamic and kinetic data were taken from Thermo-Calc Software TCFE12 Steels/Fe-alloys database [29] and MOBFE7 Steels/Fe-Alloys Mobility database [30], respectively.The mobility values for Fe and Zn at temperature T = 873 K are also listed in Table 1.

Numerical Details
The simulations were performed using YAPFI [31], a software for phase-field simulations based on the Wheeler-Boettinger-McFadden phase field model [32].In YAPFI, the governing equations are solved in a fully implicit manner, and the implementation is based on a finite volume approach that allows for multiple phase-field variables and multiple components.

Results and Discussion
The simulations were performed on a setup based on Chongmo and Hillert's experiments.The initial stage at t = 0 h shows that the grain boundary is at an angle with the surface of higher Zn activity, Figure 1.The results show DIGM after 3 h, where the grain boundary was observed to migrate faster at the junction of boundary and surface, Figure 2. The grain boundary on the other side of the plot in Figure 2, however, moves only due to its curvature as no activity of Zn is present here, therefore, no DIGM occurs.In Figure 2, the velocity of the grain boundary is higher close to the surface which is why it is more bent near the surface.This happens due to the driving force from the interaction term in Equation ( 7), which acts only when an asymmetric concentration gradient across the grain boundary is present.The gradient is created by the diffusion of Zn atoms along the grain boundary.Concurrently, the energy is generated corresponding to the interaction term and it results in a change in the overall energy of one of the grains.The grain with lower energy grows over the other resulting in a grain boundary migration.In the end, it leaves behind an enriched Zn surface, as shown in Figure 3 (blue curve).Due to grain boundary migration, a steep but continuously increasing profile of X Zn is left behind the growing front, Figure 3. On a closer look at the profile, it was noted that the effect of the grain boundary curvature was only active at the beginning, and the effect of G int was relatively low during that time.This is evident from Figure 3 where the red curve corresponding to G int = 0 and the blue curve with a non-zero G int overlap with each other till 3 µm.After that, the X Zn achieves a sufficiently high value which generates a significant concentration gradient across the interface.As a result, the interaction term becomes active and aids in pushing the grain boundary to 5.6 µm (blue curve).As per the experiments, the surface concentration of Zn should reach around 11 wt% after 16.5 h (see Figure 17 in Ref. [23]) which corresponds to X Zn = 0.095.Compared to the simulation results, X Zn reaches 0.077 after 16 h, Figure 4, which is in reasonable agreement with the experimental measurement.The extent of DIGM was tested with three types of grain boundaries, depending on their contact angle at the surface of high Zn activity, Figure 5. Here, the Zn activity is active on the right-hand side of the plot, same as in Figure 1.It was observed that the migrating boundary showed a dependence on the angle of contact (here, the angle of contact is calculated anti-clockwise from the surface of high Zn activity to the grain boundary).The position of the interface in the course of their migration is shown in Figure 6.The high-angled boundary, corresponding to Figure 5a, moved slower than a lowangled boundary Figure 5c.This is clear from Figure 6, where the position of the interface with a lower angle (SL 3 ) moves faster than SL 2 when it crosses over each other at t = 2000 s.The comparison could also be done with SL 1 when the starting point of SL 2 and SL 3 are set to the same starting point as that of SL 1 .A relative comparison showed that SL 3 covers the longest distance, followed by Sl 2 and then SL 1 , Figure 7.According to Chongmo and Hillert [23] the angle with which the grain boundary meets the surface should be almost 90 • due to initial annealing at 850 • C, however, a straight grain boundary meeting at a right angle at the surface would hardly show DIGM in our simulations.Therefore, three different angles were taken accordingly and the angle close to a right angle (SL 1 ) moved to 6.8 µm, SL 2 to 4.9 µm, and SL 3 to 5.3 µm.When the velocity of the interface was calculated at the surface, SL 3 seemed to have a higher velocity in the beginning, i.e., v = 4 nm/s and then quickly reduced to 1 nm/s.Whereas SL 2 had an initial velocity of v = 2 nm/s and SL 1 had v = 1 nm/s.On further analysis of the velocity of SL 3 , see Figure 8, it was noted that it achieved a high value in the beginning because of the combined effort from the two forces, i.e., one due to G int and the other due to boundary curvature, which was aiding in grain boundary migration, however, it was observed to drop quickly to 1 nm/s and followed a stop and go motion because of a constant struggle between the two forces.In the end, the velocity was zero as the forces ended up balancing each other.On comparison with the velocity calculated for such interfaces in the experiments, it was reported as v = 2.8 nm/s at 873 K, see Table 1 in Ref. [23], which is well under the range obtained in our simulations.[23] captured a clear image of the oscillations of the grain boundaries at a surface.In the simulations, this was observed when a grain boundary beneath the surface moved forward and backward, i.e., it reverted back after moving forward for a period of time.In Figure 9, the traces (where φ = 0.5) of the grain boundary are shown where, due to the slant nature of the boundary, the direction of motion was fixed, however, due to the faster movement of the interface the angle of the boundary at the surface changed quickly, thereby, leading to a force that would stop the boundary from moving and straighten out instead.As can be seen, the grain boundary has reverted back after t = 3 h and then stopped moving after 4 h.Although the two forces, one due to G int and the other due to the curvature, would like to ultimately balance each other, the curvature of the boundary changes quickly after 3 h and the bent boundary, corresponding to Figure 2, has a force to straighten out in the opposite direction.Therefore, the two forces are again aligned but in a different direction.However, the alignment changes just after the fourth hour, and the two forces are in constant struggle with one another, thereby leading to zero velocity.Such an observation was made by Chongmo and Hillert [23], where the grain boundary migration seemed to be balanced by the growth of Zn-rich grains.The growth of grains corresponds to the curvature effects and the chemical driving force which aids in grain boundary migration corresponds to the elastic strain energy (G int ).A demonstration of such an oscillation effect could be explained by Figure 10 where the position and the velocity of the interface are plotted with time.An unusual kink was observed between the third and fourth hour and, thereafter, a zero velocity is maintained till the twentieth hour.The unusual kink in Figure 10a corresponds to the change in position of the interface from forward to backward.According to Ref. [23] the grain boundary in a Zn-rich region "should eventually come in contact with the surface. . ." and "establish a right angle", i.e., it should stop moving after a few oscillations.In the results from our simulations, this can be validated by plotting the velocity of the interface.In Figure 10b, it is seen that the grain boundary velocity continuously decreases with time and ends up achieving a negative value after t = 3 h until it immediately reverts back to zero at the fourth hour.Later, the grain boundary was observed to be stationary between 4 to 20 h, and the rest of the boundary straightened out during this time.Another comparable feature with the experiments was the non-uniform composition profile along the surface where the grain boundary has swept away more than once.It was observed that the grain boundary reverted back when the curvature effect took precedence over the strain energy effect, thereby leading to an unusual mole fraction profile of Zn at the surface as shown in Figure 11.The mole fraction profile of Zn after 20 h of simulation shows a ditch between two peaks.This corresponds to the position where the boundary was stuck after oscillation.At this position, there is an equal driving force due to G int on either side of the boundary as the grain boundary is almost straight and there is no excess driving force due to curvature.

Conclusions
The phase field model proposed by Mukherjee et al. [19] was used successfully to explain most of the observable features of DIGM in the Fe-Zn system.The following conclusions can be made:

•
The assumption that the driving force for DIGM comes from the coherency strain energy (G int ), that is generated due to the concentration gradient across the interface, seems to be valid for Fe-Zn.This gradient was achieved by swift grain boundary diffusion that was controlled by the factor F, see Ref. was also observed that the grain boundary with a higher angle of contact between the boundary and the surface (close to 90 • ) moved slower than a low-angled grain boundary.When evaluated, the velocity of the low-angled grain boundary was very high at the beginning of migration and then gradually dropped to follow a stop-and-go motion till it reached zero velocity.

Figure 1 .
Figure 1.Plot of φ showing the initial setup with a slant grain boundary meeting the surface of high Zn activity.

Figure 2 .
Figure 2. Plot of variable φ shows bent grain boundary near the surface due to DIGM.

Figure 3 .
Figure 3. X Zn along the surface referring to Figure 2 at x = 10 µm, with and without G int after 3 h.

Figure 4 .
Figure 4. X Zn profile at the surface left behind by the moving grain boundary.

Figure 5 .
Figure 5. Plot of the variable φ at t = 0 h for three slant boundaries (a) SL 1 , (b) SL 2 and (c) SL 3 , based on different angle of contact at the surface, i.e., right side of the plots.

Figure 6 .
Figure 6.Position of the interface after 2 h for SL 1 , SL 2 , and SL 3 .

Figure 7 .
Figure 7. Relative position of the interface after 2 h for SL 1 , SL 2 , and SL 3 when the starting point is adjusted for SL 2 and SL 3 so that it is same for all three cases, i.e., 4 µm.A clear difference is observed when SL 3 moved further than Sl 2 and SL 1 .

Figure 8 .
Figure 8. Velocity of the interface for SL 3 setup where the velocity is observed to drop with time.Another feature of DIGM from the experiments was the oscillatory motion of grain boundaries, see Figure6in Ref.[23].The grain boundary was observed to go through an oscillatory motion until it stopped moving in the end.The micrograph corresponding to Figure7in Ref.[23] captured a clear image of the oscillations of the grain boundaries at a surface.In the simulations, this was observed when a grain boundary beneath the surface moved forward and backward, i.e., it reverted back after moving forward for a period of time.In Figure9, the traces (where φ = 0.5) of the grain boundary are shown where, due to the slant nature of the boundary, the direction of motion was fixed, however, due to the faster movement of the interface the angle of the boundary at the surface changed quickly, thereby, leading to a force that would stop the boundary from moving and straighten out instead.As can be seen, the grain boundary has reverted back after t = 3 h and then stopped moving after 4 h.Although the two forces, one due to G int and the other due to the curvature, would like to ultimately balance each other, the curvature of the boundary changes quickly after 3 h and the bent boundary, corresponding to Figure2, has a force to straighten out in the opposite direction.Therefore, the two forces are again aligned but in a different direction.However, the alignment changes just after the fourth hour, and the two forces are in constant struggle with one another, thereby leading to zero velocity.Such an observation was made byChongmo and Hillert [23], where the grain boundary migration seemed to be balanced by the growth of Zn-rich grains.The growth of grains corresponds to the curvature effects and the chemical driving force which aids in grain boundary migration corresponds to the elastic strain energy (G int ).A demonstration of such an oscillation effect could be explained by Figure10where the position and the velocity of the interface are plotted with time.An unusual kink was observed between the third and fourth hour and, thereafter, a zero velocity is maintained till the twentieth hour.The unusual kink in Figure10acorresponds to the change in position of the interface from forward to backward.According to Ref.[23]  the grain boundary in a Zn-rich region "should eventually come in contact with the surface. . ." and "establish a right angle", i.e., it should stop moving after a few oscillations.In the results from our simulations, this can be validated by plotting

Figure 10 .
Figure 10.(a) Position of the interface and (b) Velocity of the interface.These plots show the unusual kink (encircled in the figure) after 3 h which refers to the oscillation of grain boundary for SL 2 .

Figure 11 .
Figure 11.X Zn profile at the surface after 20 h.Therefore, the model for DIGM suggested byMukherjee et al. [19]  also explains most of the observed features of DIGM from experiments but there is a lack of full quantitative analysis, such as the Zn profile from the surface to the interior of the grains.A plausible explanation for this could be the inability to reproduce the nucleation of new grains at the surface as observed byChongmo and Hillert [23].Nevertheless, the current model is successful in reproducing most of the features of DIGM in a real system, such as Fe-Zn.
[19,33].• The observed features from the experiments by Chongmo and Hillert [23], such as a distorted grain boundary at the junction of the boundary and the surface, and the oscillatory motion of the grain boundary at the surface agreed qualitatively well with the simulation results.The reason for the oscillatory motion is thought to be due to a competitive relationship between two forces, i.e., one due to the interaction term (G int ) and the other due to the curvature of the boundary.During the course of migration, these two forces either aligned, opposed, or balanced each other multiple times.• The influence of the contact angle of a grain boundary with the surface on the velocity of the interface, and the mole fraction of Zn at the surface agreed quantitatively with the previous experiments [23].The velocity of the interface from Chongmo and Hillert [23] was under the range of the values obtained from the simulation.It