Large-scale phase-field simulation of three-dimensional isotropic grain growth in polycrystalline thin films

The MPF model proposed by Steinbach and Pezzolla [32] is employed for grain growth simulations. This model describes a polycrystalline system of N grains through N phase-field variables ϕ_i (i = 1, 2, ..., N), which take a value of 1 in the ith grain, 0 in the other grains, and 0 < ϕ_i < 1 at the grain boundaries. The sum of the phase fields at any spatial point in the system must be conserved as $\displaystyle {\sum }_{i=1}^{N}{\phi }_{i}=1.$ The migration of grain boundaries is reproduced by calculating the time evolution of ϕ_i for each spatial point under the constraints of free energy minimization and phase-field conservation. When considering pure curvature-driven grain growth, the specific form of the time evolution equation is given as:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\phi }_{i}}{\partial t}=-\displaystyle \frac{2}{n}\displaystyle \sum _{j=1}^{N}{s}_{i}{s}_{j}{M}_{ij}^{\phi }\displaystyle \sum _{k=1}^{N}\left\{\left({W}_{ik}-{W}_{jk}\right){\phi }_{k}+\displaystyle \frac{1}{2}\left({a}_{ik}^{2}-{a}_{jk}^{2}\right){{\rm{\nabla }}}^{2}{\phi }_{k}\right\},\end{eqnarray} \tag{ 1 }$$

where n is the number of non-zero phase fields at a given point; s_i is a step function that takes a value of 1 if ϕ_i > 0 and vanishes otherwise; and ${M}_{ij}^{\phi },$ W_ij, and a_ij denote the phase-field mobility, barrier height, and gradient coefficient of the boundary between the ith and jth grains, respectively. The parameters ${M}_{ij}^{\phi },$ a_ij, and W_ij are related to the width (δ ), energy (σ_ij), and mobility (M_ij) of the grain boundary through the following equations:

$$\begin{eqnarray}&&{M}_{ij}^{\phi }=\displaystyle \frac{{\pi }^{2}}{8\,\delta }{M}_{ij},{a}_{ij}=\displaystyle \frac{2}{\pi }\sqrt{2\delta {\sigma }_{ij}},{W}_{ij}=\displaystyle \frac{4{\sigma }_{ij}}{\delta }.\end{eqnarray} \tag{ 2 }$$

Note that the grain boundary width δ is just an artificial quantity in the current simulations, and must be sufficiently large to resolve the boundary regions. Here, we set δ to six times the grid spacing, which is reported to be a good compromise between computational efficiency and negligible grid anisotropy [33].

To examine the effects of film thickness on grain growth, multiple simulations were performed while changing the thickness of the computational domain. Figure 1 shows an example of the polycrystalline systems employed for simulating grain growth. Here, grains are distinguished according to their colors. The domain with periodic boundaries on x- and y-directions was divided into 8192 × 8192 × n_z points by regular cubic grids of size Δx, with n_z being varied as n_z = 48, 64, 80, 96; the system depicted in figure 1 corresponds to n_z = 64. For simply modeling ideal surfaces without thermal grooving, zero-flux boundaries were used in the thickness (z-) direction, which is equivalent to the conditions where surface energies are isotropic and much larger than grain boundary energies (i.e. grain boundaries migrate so as to form local contact angles of 90° with the surfaces). The zero-flux boundaries were realized by applying the Neumann condition with zero normal gradients to the phase-field variables on the boundaries, namely, as ∂ϕ_i/∂z = 0. In this case, the actual domain thickness, h, equals (n_z − 1)Δx. Similar treatment of free surfaces has been also employed in the Monte-Carlo Potts simulations in [15]. Here, it is worth noting that there are other possible choices of boundary condition for modeling surface effects. For instance, a Dirichlet boundary condition that fixes phase-field values can be used to represent fully pinned grain boundaries on surfaces or substrates, which is expected to result in the retardation of grain growth. The effects of boundary conditions on thin-film grain growth will be investigated in detail elsewhere. The time increment (Δt), grid spacing (Δx), grain boundary energy (σ), and mobility (M) for the simulations were set as follows, being non-dimensionalized with typical scales of 1 s, 10⁻⁶ m, 1 J m⁻², and 10⁻¹² m⁴ J⁻¹ s⁻¹, respectively: Δt = 0.075, Δx = 1, σ = 1, and M = 1.

Zoom In Zoom Out Reset image size

For each simulation, the initial polycrystalline structures with a fully random, equiaxed shape were created as follows: first, a polycrystalline structure with 1800 000 grains was prepared in a relatively thick domain (n_z = 144) by growing randomly distributed nuclei under a constant driving force and full periodic boundary conditions. Then, a partial region with a desired thickness was cut out from the thick domain. The number of grains contained in the obtained structures was approximately 800 000–1400 000 in proportion to the domain thickness.

All simulations were carried out on the GPU-rich supercomputer TSUBAME3.0 at the Tokyo Institute of Technology by utilizing our own CUDA C code [55] that was specifically developed for parallel GPU computation. The number of GPUs (NVIDIA Tesla P100) used for each simulation was 128, constant for all domain sizes. Technical detail of the parallel GPU computation is provided in our previous papers [55, 56]. To ensure effective MPF simulations, the APT algorithm [33–35] proposed by Kim et al [33] was employed, storing only non-zero phase-field variables. The maximum number of stored variables at each grid point was set to seven.

This section investigates the kinetic and morphological characteristics of grain growth while focusing on their transitions from 3D to 2D growth modes. First, to quantitatively evaluate the 3D–2D transition of microstructures, we calculated the temporal variations in the number fraction of 2D-like columnar grains, f_colum, for each film thickness h. Here, grains that appear both on the top and bottom surfaces of the domain are defined as columnar ones. The results are shown in figure 3, from which we can see that all of the curves appear to be sigmoidal shaped. However, with increasing h, the rate of increase in f_colum significantly decreases. If we compare the results for h = 47Δx and h = 95Δx, the time required for the development of a fully columnar structure (i.e. f_colum ≈ 1) is almost four times larger in the latter case than in the former case.

Zoom In Zoom Out Reset image size

Next, we examined grain growth kinetics. Figure 4 shows temporal variations in squared average grain size, 〈r〉², for different film thickness. Here, grain sizes on the top and bottom surfaces (r_sur) and those in the middle section (r_mid) of the domain are measured as their in-plane area (A)-equivalent radius (A/π)^1/2. The dotted lines indicate the times at which the columnar grain fraction f_colum shown in figure 3 reaches 0.5 and 0.99. For comparison, we also plotted the results for the cross sections of a 3D bulk system and for a full 2D system, which were obtained from the previous large-scale MPF simulations of ideal grain growth [56]. For steady-state ideal grain growth, it is well known that the kinetics of the system follows the parabolic growth law [57–60]:

$$\begin{eqnarray}&&{\left\langle r(t)\right\rangle }^{2}-{\left\langle r({t}_{0})\right\rangle }^{2}=kM\sigma \left(t-{t}_{0}\right),\end{eqnarray} \tag{ 3 }$$

where t₀ is the initial time and k is a constant. In all panels of figure 4, both 〈r_mid〉² and 〈r_sur〉² are observed to follow the parabolic law at relatively early and late stages of grain growth. However, for the early stages, while the temporal changes in 〈r_mid〉² coincide with that of the cross sections of a 3D bulk, the 〈r_sur〉² values increase visibly faster compared to them. Subsequently, the increases in 〈r_mid〉² are significantly accelerated around f_colum = 0.5, after which the grain growth rates in the surfaces and internal regions gradually decrease and converge to the identical one in the vicinity of f_colum = 1. Then, the slopes of size-time curves of thin-film grain growth approach that of the full 2D growth. The origin of the decrease in grain growth rates during the formation of columnar microstructures can be understood as follows: the zero-flux boundary condition used in the simulations guarantees that the grain boundaries migrate so as to form local contact angles of 90° with the surfaces. Therefore, once grain boundaries reach both the top and bottom surfaces, their through-thickness curvatures gradually decrease to zero due to interfacial tension. Finally, 3D grains become equivalent to 2D polygonal ones. Because the driving force (interfacial curvature) for 2D grain growth is generally smaller than that for 3D growth [56, 58], such change in the geometrical characteristics of grains results in the decrease in the slopes of grain size-time curves.

Zoom In Zoom Out Reset image size

The above results confirm that the 3D–2D transition of growth kinetics takes places around f_colum = 0.5 and gets accomplished near f_colum = 1 for all conditions of film thickness. We observed a similar transition behavior for the morphological characteristics of grain growth, namely, grain size distribution. As examples, figure 5 shows the temporal variations in normalized grain size distributions on the surfaces (r_sur/〈r_sur〉) and in the middle section (r_mid/〈r_mid〉), as calculated for film thickness of h = 47Δx and h = 79Δx. Here, the results are given for the times at which the columnar grain fraction f_colum is approximately 0.3, 0.5, 0.8, and 0.99. For comparison, the steady-state size distributions [56] for the cross section of a 3D bulk system and for a full 2D system are also depicted. As can be seen in the panels of figure 5, interestingly, grain size distributions on the surfaces exhibit a symmetrical shape similar to that of a full 2D system, even in the 3D growth regimes (f_colum ≈ 0.3, 0.5). On the other hand, at early stages (f_colum ≈ 0.3), size distributions in the middle section show common features with the bulk cross sections in terms of their right-skewed shapes and relatively long tails toward the large grain sizes. After f_colum ≈ 0.5, the distributions in the middle section gradually converge to those of the surfaces, and become almost identical around f_colum ≈ 1. It was confirmed from the results in figure 5 that, when the values of f_colum are almost the same, grain size distributions for films with different thickness take very similar shapes.

Zoom In Zoom Out Reset image size

The results presented in the previous section indicate that the thin-film grain growth behavior is uniquely determined by the columnar grain fraction f_colum, irrespective of film thickness. Herein, we discuss a way of comprehensively describing the grain growth in thin films with different thickness.

First, let us consider the theory of columnar structure formation reported by Walton et al [18]. According to their theory, if films have fully random initial structures with a similar average grain size, their columnar grain fraction f_colum during grain growth can be expressed by a function using only the ratio of time t to squared film thickness h², namely, as f_colum = f_colum (t/h²). That is, if we use a scaled time τ = t/h² instead of t, f_colum for films with different thickness is expected to take identical values at a given time τ. Note that, although the theory of Walton et al [18] is constructed using a 2D model of a through-thickness plane of thin film (corresponding to the side views in figure 2), the outline of the derivation may also hold for 3D by replacing 'length' with 'area' in their arguments. To test the above scaling law, we calculated and plotted the temporal variations of f_colum for different film thickness as a function of τ = t/h², as shown in figure 6. In the figure, the f_colum versus τ curves for different thickness are almost perfectly overlapped. It is evident, therefore, that f_colum of 3D films follows the scaling law of Walton et al.

Zoom In Zoom Out Reset image size

Next, we focus on the scaling of grain growth kinetics. Provided the mode of grain growth is dependent solely on the f_colum value, growth kinetics (variation in grain size) is expected to also follow some scaling rule, as with the case of f_colum; here, we test this idea. Dividing both sides of the parabolic law (equation (3)) by h² roughly estimates that, if time t is scaled by τ = t/h², squared average grain size 〈r〉² also scales as 1/h²:

$$\begin{eqnarray}&&\displaystyle \frac{{\left\langle r(\tau )\right\rangle }^{2}}{{h}^{2}}-\displaystyle \frac{{\left\langle r({\tau }_{0})\right\rangle }^{2}}{{h}^{2}}=kM\sigma \left(\tau -{\tau }_{0}\right),\end{eqnarray} \tag{ 4 }$$

the validity of this prediction can be observed in figures 7(a) and (b), where the increases in the square of the scaled average grain sizes, 〈r〉²/h², are given as functions of τ, respectively for the surfaces and middle section of the films with different thickness. As can be seen in these figures, the scaled grain sizes exhibit a very good agreement with each other, similarly to the case of f_colum.

Zoom In Zoom Out Reset image size

The results summarized above enable us to conclude that grain growth behaviors in thin films with random initial structures are identical under uniform magnification of 1/h² and, thus, can be comprehensively described irrespective of film thickness.