Numerical model of solid phase transformations governed by nucleation and growth. Microstructure development during isothermal crystallization

A simple numerical model which calculates the kinetics of crystallization involving randomly distributed nucleation and isotropic growth is presented. The model can be applied to different thermal histories and no restrictions are imposed on the time and the temperature dependencies of the nucleation and growth rates. We also develop an algorithm which evaluates the corresponding emerging grain size distribution. The algorithm is easy to implement and particularly flexible making it possible to simulate several experimental conditions. Its simplicity and minimal computer requirements allow high accuracy for two- and three-dimensional growth simulations. The algorithm is applied to explore the grain morphology development during isothermal treatments for several nucleation regimes. In particular, thermal nucleation, pre-existing nuclei and the combination of both nucleation mechanisms are analyzed. For the first two cases, the universal grain size distribution is obtained. The high accuracy of the model is stated from its comparison to analytical predictions. Finally, the validity of the Kolmogorov-Johnson-Mehl-Avrami model is verified for all the cases studied.


I. INTRODUCTION
Transformed volume and grain morphology development due to solid-phase crystallization depend on two kinetic parameters: growth and nucleation rates. These parameters can be obtained from the morphological evolution observed by microscopy. 1 In addition, thermoanalytical techniques provide a simple and rapid way to measure the crystallization kinetics. 2 The aim of these kinds of analyses is to predict the crystallization behavior in order to define thermal treatments suitable to achieve a particular microstructure.
(1) gives: Although some authors 10 have cast doubts on the correctness of the KJMA theory, the relationship between ) (t α and ) (t ext α of Eq. (4) is exact. 11 Recent numerical simulations 8,9,[12][13][14][15] have confirmed it for several particular cases (a noteworthy analysis is given in ref. 12). The KJMA theory also holds in case of anisotropic growth provided that the grains have a convex shape and are aligned in parallel. 16 Moreover, the KJMA theory provides a good approximation when the anisotropy is moderate 13 or for soft-impingement and nonrandom nucleation. 17 Unfortunately, as far as we know, analytically exact solutions for the transformed fraction, ) (t α , are restricted to three particular situations under isothermal conditions: time-independent growth and nucleation rates, time-independent growth rate and nucleation rate proportional to a power of time, 18 and preexisting nuclei (site saturation). Recently, a quasi-exact solution of the KJMA model has been obtained under continuous heating conditions. 19 In contrast, many transformations are governed by time-dependent nucleation and growth rates and non-isothermal heat treatments are a common practice. Thus, numerical calculations are needed to simulate these general cases. The main difficulty in numerically solving Eq. (2) is the dependence on the time history through τ . A common solution is the analytical development and numerical integration of Eq. (2) for a particular set of conditions. Conversely, the number of general numerical solutions is quite sparse. Yinnon et al. 20 developed a simple method under the assumption of linear cooling or heating rate. Besides, Krüger 21 followed a different approach for non-isothermal transformations. The validity of the latter numerical solution is limited to some particular cases, as will be commented on in this paper.
The particular kinetic conditions of a phase transformation have an essential effect on the emerging grain morphology and, therefore, the material properties. Thus, several computer simulations have been developed to predict the resulting microstructure. 22 These simulations can be classified into two main groups: those based on a Monte Carlo method 23,24 and those based on cellular automata. 25,26 A common drawback of both approaches is the accumulative error at each evaluation step related to the spatial resolution resulting from space discretization. In general, the spatial resolution should be chosen high enough to reduce this error and the simulated volume should be high enough to reduce the statistical error related to the finite number of nuclei. 13 The problem is that a high spatial resolution limits the space extent. The problems related to the finite extent can be diminished by using periodic boundary conditions and by performing several runs to minimize the statistical error. conditions. In addition, the model is applied to the analysis of the effect of partial crystallization prior to isothermal treatments.
Afterwards, we introduce an algorithm that calculates the microstructure and grain-size distribution. The algorithm computes the microstructure from the previous numerical calculation of the nuclei extended volume. In contrast with the previous calculation, the actual transformed fraction is not obtained from Eq. (4). Indeed, with this algorithm, the microstructure and transformed fraction are calculated directly from nucleation, growth and impingement of the individual grains. Thus, its applicability is not restricted by the KJMA assumptions. Since our approach reaches high accuracy in short computation times with minimal computer requirements, the microstructure for 3D growth can be easily calculated. Indeed, we will compute the evolution of the transformed fraction and the grain-size distribution under isothermal conditions and for several nucleation mechanisms. The accuracy of our approach is tested against some analytical results. In particular, the prediction of the mean grain size is excellent and indicates that the grain-size distributions obtained in this paper are very accurate.
Finally, the correctness of the KJMA theory is verified in all the cases.

II. NUMERICAL CALCULATION OF THE KJMA KINETIC EQUATIONS
Our numerical approach follows Kolmogorov's 7 development; for a finite volume V of the parent phase, the extended volume is: The previous summation covers all the grains. Since the volume of the parent phase is finite, the number of grains N(t) is also finite. Thus, the actual transformed fraction can be derived from Eq. (4) provided that the volume of any arbitrary grain is much smaller than V. The latter condition is fulfilled if the number of grains is large enough. (According to appendix C the number of grains is equal to the ratio between V and the mean grain volume) The numerical calculation consists in the creation of an array which, for each single nucleus, i, stores its radius r i . To avoid an accumulative rounding error related to the calculation of N(t j ), at each time step, j, first we calculate the total number of grains, N(t j ) , and then the number of nuclei created, ΔN(t j ): (6) then, for all the previous created grains, i, their radius, r i (t j ), are updated: is the radius growth in the time interval Note that r is the radius of the extended volume of a nucleus, i.e., the radius supposing that the grain grows free.
For the grains created during the time interval an average radius is assumed: where 0 r is the critical germ size. In our calculations we will assume that 0 r is negligible.
A more accurate calculation is obtained if the actual radius is calculated for each new nucleus: where i τ is the creation time for the nucleus i. Besides, it is not necessary for the time interval to be constant; on the contrary, a more efficient simulation is obtained when the time interval is chosen such that a constant growth, r Δ , is imposed.
Finally, the extended fraction can be obtained according to Eq. (5): The simplicity of our approach relies on the fact that we solve the time history dependence by storing the radius for each nucleus. Moreover, no assumptions have been made on the time dependence of both I and G, so the numerical solution is general.

III. ACCURACY OF THE NUMERICAL CALCULATION
To check the accuracy of our numerical calculation, test runs were done for 3D growth during isothermal and continuous heating, for which there exist an exact 7 and a quasi-exact solution, 19 respectively. The calculations have been done for the particular G and I values of amorphous silicon: where T is the temperature in Kelvin and k B is the Boltzmann constant. In Table I we summarize the parameters found in the literature. 1 When nucleation and growth rates follow an Arrhenius dependence on temperature, the exact solution for the isothermal regime is: is an average activation energy and 0 t is the initial time. On the other hand, the quasi-exact solution for the continuous heating case is: 19 For the whole range of α and when the total number of nuclei is of the order of 10 5 , the deviations from the exact solutions are smaller than 10 -5 for the isothermal case, while for the non-isothermal case the discrepancies are smaller than 0.01. Actually, the larger discrepancy for the non-isothermal case is due to the analytical solution and not to the present calculation. In fact, simulations performed when the quasi-exact solution becomes exact (taking the same activation energy for growth and nucleation 19 ) resulted in deviations smaller that 10 -6 . In addition, similar minor discrepancies were also obtained for the case where all nuclei appear at 0 t t = -the site saturation case. 2 Often a thermal treatment consists of an initial constant heating period followed by an isothermal step. A widespread approximation to this problem consists in introducing a virtual initial time ' 0 t : where ' 0 t is the time necessary for the isothermal regime to reach its initial transformed fraction, i.e. the transformed fraction after the constant heating regime: where 0 , ex α is the corresponding initial extended transformed fraction. This analysis is the basis of Krüger's numerical approach. 21 Nevertheless, this approach is generally wrong since the state of the system depends on the thermal history, 20,28 i.e. a given value of α will correspond to a different state and consequently it will evolve at a different rate. From Indeed, the number of nuclei previously formed has a minor effect on the transformed fraction but has a pronounced effect on the subsequent crystallization evolution.
It is worth mentioning that our numerical approach keeps all the information related to the system such as temperature, transformed fraction, number of nuclei and size of the extended grains. Consequently, it can be applied to any arbitrary thermal history.

IV. ALGORITHM FOR EVALUATION OF THE MORPHOLOGY
In this section, we present an algorithm which evaluates the grain morphology grain. If the cell already belongs to another grain, then we are dealing with a "phantom nucleus". The concept of phantom nuclei was introduced by Avrami. The phantom nuclei do not appear in the lattice and have no effect on the calculation of the grain morphology, so they are discarded. However, as pointed out by Sessa et al., 12 they must be included in the calculation of the total extended transformed fraction performed in Section II.
The algorithm used for the grain morphology evolution will be explained with the help of Fig. 2 where black cells correspond to the nuclei, the circumferences indicate the size of the extended grains, the grey regions are the cells associated to a particular grain and white cells represent the untransformed volume. The extended radius of the circumferences, r i , is calculated from Eqs. (7) and (11). When the extended radius increment approaches x Δ all untransformed cells are checked to verify whether they have been incorporated to a neighboring grain or not. To save computing time, only the grains that have at least one transformed cell "in the vicinity" of the center of the untransformed cell are checked. The vicinity analyzed is determined by x Δ . If there is a grain, i, such that the distance from the nucleus (center of the black cell) to the center of the untransformed cell is smaller or equal to r i , then the number i is assigned to the cell.
When the cell is in the range of more than one grain, then the grain which reaches the cell first is assigned to the cell. Otherwise the cell remains untransformed. For example, the dashed circumferences represent the extended volume grown during the next step.
The cell a remains untransformed since there are no grains in the vicinity. Cells b and c only have a grain in the vicinity, which is the dark grey. Cell b will turn dark grey since its distance to the nucleus is smaller than the circumference radius, and conversely, cell c will remain untransformed. Finally two grains are located in the vicinity of cell d and the distance to both nuclei is smaller than the circumference radius, so cell d belongs to the grain that reaches it first.
Once the grain morphology has been built, the grain size distribution can be calculated directly from the radius of a grain. For 3D grain growth: where i v is the actual grain volume. And the transformed fraction is: Besides, nucleation and growth calculations are separate from the grain morphology evaluation. Indeed, x Δ is usually several orders of magnitude larger than r Δ and one single microstructure evaluation involves a large number of time steps.
Therefore, our algorithm can easily deal with time dependent nucleation and growth rates. Actually, the computer time required to evaluate nucleation and growth (Section II) is practically negligible when compared to the microstructure calculation. Hence, handling complex time dependences for nucleation and growth rates does not represent a significant amount of the computing time.

V. ACCURACY OF THE ALGORITHM
In the Monte Carlo and cellular automata calculation methods, the principal source of error is the space discretization. The growth at each step is associated to one cell. To minimize this error, the linear growth is adjusted to be an integer fraction of x Δ . The error can be avoided also by using particular growth modes. 29 However, in other cases it is unavoidable that, at each step, some cells are incorrectly assigned. The result is an accumulative error that progressively reduces the accuracy. Consequently, x Δ must be as small as possible. V should be high enough to reduce the boundary effects and to allow a high population of nuclei to minimize statistical errors. Taking into account the computer memory limitations, the latter restriction is especially dramatic in 3D simulations where the number of cells is 3 / x V Δ . Usually, satisfactory accuracy is reached by performing the calculation several times. 13 Indeed, averaging for n calculations is equivalent to involve n times more nuclei, thus the statistical error is reduced.
In our algorithm, since the grain growth is driven by the evolution of its extended volume, the accumulative error associated to x Δ is suppressed. A grain boundary is not allowed to grow unless its distance to the nuclei is less than or equal to  Fig. 1 in ref. 14).
The number of cells in our simulation is larger than that of Crespo et al.
However, they have performed the average for 32 simulations while we have done only one calculation, i.e. their calculation is equivalent to a single calculation with the same discretization but a number of cells of 32 times greater (537 million cells). Despite the fact that their total number of cells is slightly larger, our calculation is significantly more accurate. For instance, when comparing the evolution of the transformed fraction (Fig. 3); the largest deviation from the theoretical value in our case is 3 10 -4 , while Crespo et al. reported an accuracy of 10 -3 . For the final grain radius distribution, Fig.   4.c, our simulation exhibits an error clearly smaller than the one obtained by Crespo et al. (Fig. 2.f in ref. 14). Moreover, from Fig. 4 one can observe that the accuracy of the grain radius distribution improves as the transformation proceeds. The reason is that, as the transformation proceeds, the number of nuclei increases and consequently the error diminishes. In opposition, in the case of Crespo et al. (Fig. 2 in ref. 14), the accuracy diminishes as the transformation proceeds. The reason is the accumulative error related to the space discretization in the Monte Carlo method. Thus, one can conclude that the accumulative error due to the space discretization has been eliminated. Concerning the computer requirements, the program was run on a standard personal computer and lasted 36 minutes and the amount of memory required was 2 Gbytes.
Although our algorithm is able to work with larger values of x Δ , its accuracy is still limited by the space discretization. Indeed, the smaller the grain is, the greater the error introduced by x Δ . To minimize this effect, we have introduced the following condition: if the grain radius of Eq. (18) is larger than the extended radius r i , then we take, i i r r = , otherwise the initial value of Eq. (18) is taken. This condition significantly reduces the error introduced by x Δ at the first stages of nuclei growth but has no effect when the grain impingement takes place. When the grains impinge, the precise calculation of r i allows us to accurately establish which grain reaches the center of a particular cell first. However, the center of a cell belonging to a grain does not necessarily apply to the rest of the cell. Thus the error in the calculation of the grain radius is about 50% of x Δ . Therefore, x Δ must be, at least, one order of magnitude smaller than > < r . On the other hand, the larger x Δ is, the larger the number of nuclei and the smaller the statistical error is. Actually, the analysis of the accuracy of the algorithm with respect to x Δ confirmed that the optimum value of x Δ is approximately

ISOTHERMAL CONDITIONS
In this Section we will analyze the different microstructures that emerge depending on the nucleation mechanism under isothermal conditions. We will focus our attention on the case of amorphous silicon. In this case, nucleation is continuous. In addition, the nucleation mechanism can be modified by introducing pre-existing nuclei, e.g. by ion implantation prior to isothermal annealing 30 or by pre-annealing the sample. 31 Furthermore, the simultaneous nucleation by both mechanisms has also been observed in the crystallization of metallic glasses. 32 Thus, we will analyze the site saturation nucleation case (crystallization of preexisting nuclei) alone and mixed with the continuous nucleation. As will be stated, the results obtained here can be extrapolated to any system featuring these nucleation mechanisms. In particular, the results obtained for continuous nucleation and pre-existing nuclei are universal.

Continuous nucleation
Under isothermal conditions, both I and G are constant. Hence, the system has an exact dimensional scaling, Eq. (20) (see Appendix B), so the behavior is universal. 33 In Figs . The corresponding grain radius distribution for 2D growth is shown in Fig. 8.
Concerning the grain radius dispersion, defined by the standard deviation of the distribution:

Pre-existing nuclei
In this case, with the natural time and space scaling, 27 ( ) the dimensionless solution is universal, as well. The parameter n 0 is the pre-existing nuclei density. Fig. 7 shows the corresponding final grain radius distribution. Its average grain size is obviously ' . From the resulting microstructure we have calculated its standard deviation: Therefore the size distribution is significantly narrower than in the preceding case because all nuclei appear simultaneously.
This distribution fits a Gaussian distribution with remarkable accuracy (the square correlation coefficient is 0.9998): The final grain radius distribution for 2D growth is plotted in Fig. 8. Here again we fit the calculated distribution to a Gaussian distribution, though the fit is not as fine As outlined in appendix C, there is not a universal solution in this case. Instead, the grain size distribution depends on the relative contribution of both nucleation mechanisms. Indeed, in appendix C it is shown that the final mean grain size depends only on two parameters: n n n / ' 0 0 ≡ (where n is the density of grains that would result without preexisting nuclei) and the space scaling factor λ . In Fig. 9  continuous nucleation, it is proportional to a factor that depends on both nucleation and growth rates. Thus, once we know their temperature dependence, we can easily control the final grain size by choosing the appropriate temperature. However, since the solution is universal, we cannot control the grain size distribution (shape and width). On the other hand, the distribution is significantly narrower when only pre-existing nuclei are present. Here again the mean grain size can be easily controlled with temperature, but the distribution is also universal. Thus, the distribution can be modified by mixing both nucleation mechanisms. However, in this case, the resulting distribution is very similar to a superposition of the distributions corresponding to both nucleation mechanisms acting independently.
Although the results have not been detailed for all the simulations reported, the agreement between the transformed fraction calculated numerically or analytically from Eq. (4) and from the constructed microstructure is very good (the error is smaller than 0.003). These results, in addition to those of Sessa et al., 12 Crespo et al. 14 and Pusztai et al., 13 represent a direct confirmation of the Avrami theory (in particular of Eqs. (1) and (4)) for a number of particular cases.

VII. CONCLUSIONS
To sum up, we have introduced a new numerical method for solving the kinetic equations of KJMA theory. Our method, shares the main Avrami assumptions: random nucleation and isotropic growth. The numerical model takes into account all the parameters that define the system state and no restrictions are imposed either to the thermal history or to the growth and nucleation rate dependencies on temperature. A comparison between the numerical results and the analytical solutions speaks for the high accuracy of the method (relative error smaller than 10 -6 ).
In addition, we have introduced an algorithm to calculate the grain morphology from the extended microstructure. Since the extended microstructure is calculated according the previous numerical model, this algorithm keeps its flexibility and can be easily adapted to any conditions. In contrast, the calculation of the transformed fraction is not based on KJMA theory. Thus, its applicability is not restricted by the assumptions of KJMA theory. This fact means that our algorithm can be used to test the correctness and range of applicability of KJMA's theory. Compared to existing numerical methods, our algorithm features high accuracy in relatively short computational times and low computer memory requirements.
The algorithm has been used to obtain the universal grain size distribution for constant growth rate and constant continuous nucleation as well as for pre-existing nuclei. As far as we know, both results are the most accurate ever published. In addition, the case of continuous nucleation combined with pre-existing nuclei has been analyzed. In this case, a universal solution does not exist. However, the grain size distribution depends on the relative contribution of both mechanisms (apart from a size scaling factor).
In all cases, the agreement between numerical results and theoretical predictions is excellent. Table II summarizes the calculated average grain size and the standard deviation of the grain radius distribution.

This work has been supported by the Spanish Programa Nacional de Materiales
under contract number MAT2006-11144.

SUBSEQUENT TO CONSTANT A HEATING STEP.
During the isothermal regime, the extended transformed fraction has two contributions, one from the nuclei created during the heating step and a second one from the nuclei created during the isothermal step: where t 0 is the time where the isothermal regime starts and ) , where β is the heating rate and Let's assume that two systems, characterized by their particular values of the growth and nucleation rates (G 1 , I 1 ,G 2 , and I 2 ), differ only by a scale factor at given times, t 1 and t 2 , respectively (see Fig. 11). Their microstructure will maintain this scale relation provided that at Thus, the behavior of both systems is the same when scaled by τ and λ .
Finally, note that under the previous scaling, the dimensionless growth and nucleation rates are equal to 1. Indeed, the fact that dimensionless parameters do not depend on the particular values of I and G proves that the dimensionless system is universal, i.e. the evolution of any particular system can be obtained from the dimensionless system simply by multiplying the dimensionless time and space by τ and λ respectively.

MATERIAL.
For the sake of simplicity, the calculations are done for 3D growth only. The mean grain size is defined as: where N is the total number of grains. Likewise, when the transformation is over, the total volume is the sum of the volume of the individual grains, i v : and taking into account the definition of the grain radius (Eq. (18) In other words, nuclei become real grains only when they appear in the untransformed volume. Then, for isothermal and continuous nucleation, Eq. (C.7) reduces to: ( )