- Email: [email protected]
Unexpected selection of growing dendrites by very-large-scale phase-field simulation
Journal of Crystal Growth 382 (2013) 21–25
Contents lists available at ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro
Unexpected selection of growing dendrites by very-large-scale phase-field simulation Tomohiro Takaki a,n, Takashi Shimokawabe b, Munekazu Ohno c, Akinori Yamanaka d, Takayuki Aoki b a
Mechanical and System Engineering, Kyoto Institute of Technology, Kyoto 606-8585, Japan Global Scientific Information and Computing Center, Tokyo Institute of Technology, Tokyo 152-8550, Japan c Division of Materials Science and Engineering, Hokkaido University, Sapporo, Hokkaido 060-8628, Japan d Mechanical Systems Engineering, Tokyo University of Agriculture and Technology, Koganei, Tokyo 184-8588, Japan b
art ic l e i nf o
a b s t r a c t
Article history: Received 17 June 2013 Received in revised form 17 July 2013 Accepted 24 July 2013 Communicated by M. Plapp Available online 6 August 2013
Dendrites are typical growth morphologies in alloy solidification. However, the ways in which dendrites grow and form solidification microstructures remain poorly understood. Here we show unexpected findings that reveal the survival of unfavorably oriented dendrites and highly complicated dendrite– dendrite interactions in three-dimensional space during the directional solidification of a binary alloy. These results are observed for the first time through very-large-scale phase-field computations performed by a graphics processing unit (GPU) supercomputer and a high-performance algorithm developed for parallel computing. & 2013 Elsevier B.V. All rights reserved.
Keywords: A1. Dendrites A1. Directional solidification A1. Computer simulation
1. Introduction Dendrites, which are prime examples of spontaneously generated patterns found in nature, are formed during the solidification of a number of materials such as non-metallic and metallic systems [1,2]. The prediction and control of dendritic microstructures are among the most important issues in the field of solidification science and engineering [3]. This is because the morphologies and sizes of inhomogeneous entities such as grains, segregation patterns, and secondary precipitates in solidification microstructures have a profound influence on several material properties. Advances in experimental technologies such as in-situ observation techniques of non-metallic [4–6] and metallic systems [7–10] have improved our understanding of dendritic growth. However, research to date has been centered on thin samples, which are essentially two-dimensional (2D) in terms of dendrite– dendrite interactions, small systems where only a small number of dendrites grow, or both. Yet the collective behavior of many dendrites in three-dimensional (3D) space is of primary importance in controlling the solidification microstructures of material, because this determines the size and morphology of dendritic microstructures.
n
Corresponding author. Tel.: +81 75 724 7317; fax: +81 75 724 7300. E-mail address: [email protected] (T. Takaki).
0022-0248/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jcrysgro.2013.07.028
Directional solidification is a useful way of evaluating interactions between dendrites. It is generally accepted that a favorably oriented dendrite, for which the inclination angle between the preferred growth direction and heat flow direction is small, can stop the growth of an unfavorably oriented dendrite, for which the inclination angle between the preferred growth direction and heat flow direction is larger than that for the favorably oriented dendrite [11–13]. Hereafter, this is referred to as “expected” selection behavior. Recently, an unexpected phenomenon, in which an unfavorably oriented dendrite can stop the growth of a favorably oriented dendrite, has been observed in both experiments and simulations [14–18]. Although this new unexpected phenomenon is very interesting, it was observed under simple bicrystal conditions. Because actual solidification microstructures are much more complicated, we need to investigate the interactions among multiple dendrites. The phase-field method has emerged as a powerful tool to simulate the complex dynamics of dendritic growth [19–28]. Its ability to describe the formation of complex dendritic patterns is well acknowledged. However, phase-field computations of dendritic growth have thus far been limited to 2D phenomena or 3D growth of a single dendrite [22,23,29,30], because of the large computational cost stemming from the diffuse interface model, where the order parameter that describes the different phases changes smoothly in a thin interface region [31,32]. Although some phase-field simulations of multiple grain growth in 3D have
22
T. Takaki et al. / Journal of Crystal Growth 382 (2013) 21–25
been reported for equiaxed grains [33,34], dendrites [35], and columnar cells [36,37] in directional solidification, these are limited to relatively small spaces and have simple morphologies. Therefore, as far as we know, there are no phase-field simulations for multiple and complicated dendrites in sufficient 3D space to evaluate realistic microstructures. Hence, the key bottleneck here is the development of large-scale computation technologies. We have recently enabled the first-ever petascale phase-field computation in single precision for the directional dendritic solidification of a binary alloy [38]. This has been realized by using 4000 graphics processing units (GPU) along with 16 000 central processing unit (CPU) cores of the TSUBAME 2.0 supercomputer at the Tokyo Institute of Technology. To achieve high performance, we have written our simulation code from scratch in the GPU-specific language CUDA and introduced an overlapping technique that avoids performance degradation due to inter-GPU communication. In this study, we show a very-large-scale 3D phase-field simulation of competitive dendritic growth in the directional solidification of an Al–Si binary alloy by using the GPU supercomputer TSUBAME2.0 and a custom-developed high-performance algorithm. We focus on phenomena captured only in a very large system during a very long time period. Specifically, our simulation were carried out in a system with the dimensions of 3.072 mm 3.078 mm 3.072 mm (4096 4104 4096 meshes) for a total time period of more than 100 s (four million computational steps)—the largest simulation of dendrite growth that have ever been reported—by using 768 GPUs with 768 CPUs in TSUBAME2.0. This very-large-scale simulation in both spatial and temporal domains is indispensable for elucidating the nature of dendritic competitive growth. In fact, we demonstrate an unexpected finding that is obtained only by computations on this scale.
2. Phase-field model The existence of a liquid or solid phase is expressed by a phasefield variable ϕ, which is defined as 0 in a liquid or 1 in a solid and changes smoothly at the interface region. By using ϕ, we define the following Ginzburg–Landau type free energy functional F [20]: Z 2 a F¼ ð∇ϕÞ2 þ WqðϕÞ þ pðϕÞf S þ ð1pðϕÞÞf L dV ð1Þ V 2 where a is the gradient energy coefficient, W is the height of the double-well potential, and fS and fL are the free energy densities in solid and liquid, respectively. The functions pðϕÞ and qðϕÞ are the energy density function and double-well function, respectively; here we use pðϕÞ ¼ ϕ3 ð1015ϕ þ 6ϕ2 Þ and qðϕÞ ¼ ϕ2 ð1ϕÞ2 [19]. The concentration in a binary alloy is defined as c ¼ ϕcS þ ð1ϕÞcL , where cS and cL are the concentrations for the solid and liquid, respectively [24]. These concentrations are made mutually dependent according to cS ¼ kcL with the partition coefficient k. Time evolutions for ϕ and c are described by the Allen–Cahn equation [32] and Cahn–Hilliard equation [31], respectively, as follows: 9 9 8 8 2 > > > > > > > > = = < < 6 2 ∂ϕ ∂ ∂a ∂ ∂a 2 2 þ ð ∇ϕ Þ ð∇ϕÞ ¼ Mϕ 6 a a ∇ a ∇ϕ þ 4 > > ∂ϕ ∂ϕ ∂t ∂x> ∂y> > > > > ; ; : ∂ : ∂ ∂x ∂y 9 8 3 > > > > = ∂< ∂a dpðϕÞ dqðϕÞ7 7; þ a ð∇ϕÞ2 SΔT W ð2Þ > ∂ϕ ∂z> dϕ dϕ 5 > > ; : ∂ ∂z ∂c ð1kÞc ¼ ∇ D ∇c þ ∇ϕ : ∂t 1ϕ þ kϕ
ð3Þ
Here, M ϕ is the phase-field mobility; S is the entropy of fusion; and ΔT is the undercooling expressed by ΔT ¼ T m ðcÞT, where T is the temperature and Tm(c) is the temperature on the liquidus line for concentration c. The diffusion coefficient is expressed by D ¼ DS þ ðDL DS Þð1ϕÞ=ð1ϕ þ kϕÞ where DS and DL are the diffusion coefficients in the solid and liquid, respectively. T changes with time t according to T ¼ GzRt þ T 0 where G is the temperature gradient, R is the cooling rate, z is the spatial coordinate in the heat flow or G direction, and T0 is the temperature at the bottom surface (z¼ 0). To introduce anisotropy of the kinetic coefficient and interface energy, we use 8 4 4 4 9 > ∂ϕ ∂ϕ ∂ϕ > > > > > þ þ < ~ ~ 4ζ ∂ x ∂ y ∂z~ = ~ ð4Þ η ∇ϕ ¼ ð13ζ Þ 1 þ ~ 4 > > 13ζ j∇ϕj > > > > ; : ~ ¼ M ϕ ηð∇ϕÞ ~ ~ ¼ aηð∇ϕÞ, ~ for M ϕ and a to obtain M ϕ ð∇ϕÞ and að∇ϕÞ respectively [29]. Here, M ϕ and a are the standard values of M ϕ ~ y; ~ z~ Þ is a material and a, respectively. The coordinate system ðx; coordinate system in which each axis corresponds to the 〈100〉 direction of a cubic lattice. The following coordinate transforma~ y; ~ z~ Þ: tion is used between the coordinate systems of ðx; y; zÞ and ðx; 8 ∂ϕ 9 2 ∂x ∂y ∂z 38 ∂ϕ 9 > > > ∂x > > = 6 ∂x~ ∂x~ ∂x~ 7> = < ∂x~ > < > ∂y ∂ϕ ∂ϕ ∂x ∂z 7 6 ¼ 4 ∂y~ ∂y~ ∂y~ 5 ∂y ∂y~ > > > > > ∂y ; ; : ∂ϕ > : ∂ϕ > ∂x ∂z > ∂z ∂z~ ∂z~ ∂z~ ∂z~ 2 3 1 0 0 6 cos θx sin θx 7 ¼40 5 0 sin θx cos θx 38 ∂ϕ 9 2 32 > cos θy 0 sin θy cos θz sin θz 0 > > = < ∂x > 7 ∂ϕ 6 76 0 1 0 ð5Þ 4 54 sin θz cos θz 0 5 ∂y > > > > ; sin θy 0 cos θy 0 0 1 : ∂ϕ ∂z where θx , θy , and θz are the rotation angles around the x, y and z axes, respectively. The preferred growth directions, which are different for every nucleus, are set to the lattice points with ϕ Z105 under the initial conditions. During growth, if the phase-field variable ϕ of a lattice point without the preferred growth direction, which is the nearest lattice point to a lattice point with a preferred growth direction, becomes larger than 10 5, the lattice point has the same preferred growth direction as that of the nearest lattice point. In this way, although two neighboring dendrites with different preferred growth directions sometimes impinge on each other at the low-temperature side with respect to the dendrite tip positions of the primary arm, this does not affect competitive dendrite growth, on which we focus in this study. a, W, and M ϕ in Eq. (2) can correspondingly be related to the interface thickness interface pffiffiffiffiffiffiffiffiffiffiffiffiffi δ, interface energy γ, and pffiffiffiffiffiffiffiffi mobility M by a ¼ 3δγ=b, W ¼ 6γb=δ, and M ϕ ¼ μ 2W =ð6aSÞ, respectively, where μ is the kinetic coefficient and b is a constant related to the interface thickness. Because, in this study, we aimed at very-large-scale phase-field simulation, we used a conventional model without an antitrapping current term [22,23,27,28], to reduce the computational costs associated with the antitrapping current. Therefore, the following results are not quantitative. However, we believe that the results give us important information on competitive growth between multiple dendrites and open up new roads for the investigation of microstructural growth mechanisms by large-scale phase-field simulations. 3. Numerical conditions Directional solidification of the Al–Si binary alloy is simulated in a computational domain of 3.072 mm 3.078 mm 3.072 mm
T. Takaki et al. / Journal of Crystal Growth 382 (2013) 21–25
divided into 4096 4104 4096 finite difference meshes with a uniform grid spacing of Δx ¼ Δy ¼ Δz ¼ 0:75 μm. The initial liquid concentration is set to an atomic fraction of 0.05 and the corresponding liquidus temperature is calculated as 903.47 K from T m ðcÞ ¼ mcþ T m , where m ¼ –600 K/at.frac. is the liquidus slope and Tm ¼933.47 K is the melting temperature of pure Al. Other material properties and computational conditions used in the present simulation are as follows: temperature gradient in z-direction G¼10 K/mm; cooling rate R¼ 0.3 K/s; temperature at z¼0 surface T0 ¼901.47 K; interface thickness δ ¼ 6Δx ¼ 4:5 μm; entropy of fusion S¼ 8.0 105 J/Km3; interface energy γ ¼ 0:16 J=m2 ; interface kinetic coefficient μ ¼ 0:32 103 m=Ks; diffusivity of solid DS ¼1 10 12 m2/s; diffusivity of liquid DL ¼3 10 9 m2/s; partition coefficient k¼0.13; and strength of anisotropy ζ ¼ 0:03. As for the boundary conditions, the periodic boundary condition was set to the edges of x and y, and the Neumann condition was employed at the edges of z for the phase-field and concentration. The time increment Δt was set to Δt ¼ 2:68 105 s. Four million computation steps, which correspond to 107 s in real time, were then performed. The 128 spherical solid seeds of radii 3Δx with random preferred growth directions are placed on the bottom surface (z ¼0). Here, the inclination angle θ for the heat flow direction is expressed by θy because all θx are set to zero just for convenience. We confirm that special textures were not formed by setting θx ¼ 0. Higher-order branching is formed by 10% fluctuation of the concentration. The present simulation typically takes 2 h 59 min to perform a 40 000 time step simulation, including data output. The total simulation runtime for obtaining numerical results is about 12 days. The performance we have achieved in this simulation is approximately 187 TFlops. Note that this performance is lower than we expected from the results of a previous simulation, i.e., 2PFlops [38]. Runtime configurations, such as the decomposed domain size, and the GPU kernel configuration used in this simulation are not sufficiently optimized to achieve the desired performance, and network resources are shared with other TSUBAME users, which probably result in degradation of the simulation performance.
23
4. Numerical results Fig. 1 shows the time slices of competitive dendritic growth. In the simulation, 128 solid seeds initially nucleated at the bottom wall grow into star-like particles (Fig. 1a), and then cover the bottom surface by lateral growth (Fig. 1b). After covering the entire bottom surface, the dendrites begin growing in a longitudinal direction that corresponds to the heat flow direction with simultaneous growth of branched dendrites, viz. the secondary arms (Fig. 1c). Competitive growth of dendrites occurs, as is reflected in the selection of growing dendrites (Fig. 1d). The number of seeds gradually decreases as the solidification front proceeds upward. When the front is located at 1.2 mm from the bottom (Fig. 1e), almost half of the original 128 seeds stop growing. The number of seeds further decreases and only 46
Fig. 2. The tip positions of 128 dendrites in the z-direction at 96.5 s. The dendrites are categorized into three groups A, B and C: (A) those that ceased growing in a region where the number of stopped dendrites increases rapidly with the tip distance along the z-coordinate, (B) those that ceased growing in a region where the number of stopped dendrites increases almost proportionally with tip distance, and (C) those that survived to 96.5 s.
Fig. 1. Time slices of competitive dendritic growth. The color indicates the inclination angle θ of the preferred growth direction 〈100〉 relative to the x-axis. (a) 128 spherulitic seeds change in morphology to star-like dendrites (1.1 s). (b) Dendrites preferentially grow in the lateral direction or wet the bottom wall (5.4 s). (c) Dendritic arms start to grow in heat flow direction (10.7 s). (d) Competitive growth occurs where the number of solid seeds decreases approximately in proportion to the square of the distance from the bottom (26.8 s). (e) Competitive growth continuously occurs where the number of seeds decreases approximately in proportion to the distance from the bottom (53.6 s). (f) Final morphology of the present simulation (107.1 s). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)
24
T. Takaki et al. / Journal of Crystal Growth 382 (2013) 21–25
ultimately survive after 107.1 s (Fig. 1f). The 3D simulations thus reproduce the full range of directional solidification phenomena from heterogeneous nucleation on the mold wall to competitive growth between multiple dendrites. Fig. 2 shows the dendrite tip position (growth front in heat flow direction) at 96.5 s for all 128 dendrites, each of which is specified by the seed number on the left. The dendrites are categorized into three groups: A, B and C. The dendrites that stop growing before 96.5 s are categorized separately from those still growing after 96.5 s. Of the former, dendrites categorized in A that do not grow beyond about 1.2 mm from the bottom rapidly decrease in number as the distance from the bottom increases. As shown in Table 1, in this category are 65 dendrites, half (32 dendrites) of which have a preferential growth direction largely inclined from the heat flow direction (A4) (i.e., inclination angle θ ¼ 301–451). This is predictable from a generally accepted dendritic selection mechanism in
Table 1 The number of dendrites categorized by dendrite tip position and inclination angle θ. Categories A, B, and C correspond to those in Fig. 2. θ ð1Þ
A
0–10 10–20 20–30 30–45 Total
(A1) (A2) (A3) (A4) 65
B 9 11 13 32
(B1) (B2) (B3) (B4) 17
C 2 5 6 4
(C1) (C2) (C3) (C4) 46
Total 23 15 6 2
34 31 25 38 128
which only dendrites with small inclination angles grow to survive [11–13]. At the same time, more than 10% of the dendrites in this first category A1 (9 dendrites) actually have a small inclination angle (θ ¼ 01–101). Among these, 3 dendrites are prevented from growing by neighboring dendrites with a smaller inclination angle; however, the remaining 6 are stopped by those with a larger inclination angle. Table 1 reveals that at the very beginning of solidification, the growth of more than 15% of all dendrites with a small inclination angle (θ ¼ 01–101) is actually blocked by neighboring dendrites with a larger angle. A second category B of dendrites that stop growing comprise 17 dendrites whose tip distances are below about 2.7 mm, where the solidification process should be close to the steady-state regime. Of these, 15 have an inclination angle of 101–451 and are stopped from growing by dendrites with a smaller inclination angle. The 2 dendrites in B1 with a small inclination angle (θ ¼ 01–101) are stopped by dendrites with a larger inclination angle. The growth of these dendrites is blocked through the interaction of more than three dendrites. This phenomenon is beyond the explanatory scope of the recently proposed 2D selection mechanism, which accounts for how unfavorably oriented dendrites can overgrow favorably oriented dendrites in both experimental and simulation studies [14–18]. Of the dendrites in the third category C (i.e., those still growing after 96.5 s), half have an inclination angle of 01–101 (C1). The number of dendrites decreases with increasing inclination angle. This is a reasonable result assuming a conventionally accepted dendrite selection mechanism; yet even in this category, 2 dendrites in C4 with a large inclination angle of 301–451 survive.
Fig. 3. Dendritic interactions around (a) dendrite A1 and (b) dendrite B1 categorized as C4 in Table 1. The inclination angle (in degrees) for each dendrite is indicated in parentheses.
T. Takaki et al. / Journal of Crystal Growth 382 (2013) 21–25
Why and how the 2 dendrites (A1, θ ¼ 35:11; B1, θ ¼ 42:51) continuously grow are shown in Fig. 3, which presents side and top views of these dendrites along with their surrounding dendrites. As shown in Fig. 3a, the leftward growth of A1 is stopped by A2 and A3 in the upper and lower parts, respectively; meanwhile, A1 grows rightward, successively stopping A4 (θ ¼ 12:91), A5 (θ ¼ 34:11), and A6 (θ ¼ 13:81), which have a smaller inclination angle. This selection process is also seen in the growth of A2 (θ ¼ 27:71), which stops the growth of a dendrite with a smaller inclination angle (A3, θ ¼ 7:31). The favorably oriented dendrites that grow adjacent to A1 (i.e., A7–A15) are unable to stop the growth of A1 because they grow almost parallel to or away from it. In Fig. 3b, B1 is stopped on both sides by favorable dendrites B3, B4, and B5. The surrounding dendrites B5, B6, B8, and B9 grow almost parallel to B1 and do not prevent its growth. As shown above, these 2 dendrites continue to grow in a similar way with a plate-like structure. In addition, B2 is observed to penetrate B1; A2 and A16 also penetrate one another. Thus, dendrite–dendrite interaction is a highly complicated phenomenon. Our simulations demonstrate that contrary to commonly held assumptions, the dendrite selection process cannot be discussed solely on the basis of single parameter, viz., the inclination angle of the preferential growth orientation with respect to the direction of heat flow. The detailed 3D morphology of the dendrites is a non-trivial factor that strongly affects dendrite–dendrite interactions. 5. Conclusions Our results demonstrate a series of directional solidification steps from nucleation to competitive growth in 3D for the first time. This very-large-scale simulation was achieved by multi-GPU computation with a high-performance algorithm best suited to phase-field computation. Our large-scale simulations have shown that a relatively large amount of unfavorably oriented dendrites can survive during unidirectional solidification. This finding has profound significance in controlling as-solidified microstructures, especially in terms of their crystallographic features. We believe that the present work constitutes an important step in understanding and developing solidification science, based on verylarge-scale phase-field simulations. Acknowledgments We are grateful for the opportunity to use the TSUBAME2.0 at the Tokyo Institute of Technology. We thank Dr. A. Nukada, Tokyo Institute of Technology, for fruitful comments on supercomputing. This study was financially supported in part by JSPS KAKENHI Grant Number 25289006 and the ISIJ Research Promotion Grant from the Iron and Steel Institute of Japan.
25
Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version of http://dx.doi.org/10.1016/j.jcrysgro.2013.07.028.
References [1] J.S. Langer, Reviews of Modern Physics 52 (1980) 1. [2] S.A. David, T. DebRoy, Science 257 (1992) 497. [3] M. Asta, C. Beckermann, A. Karma, W. Kurz, R. Napolitano, M. Plapp, G. Purdy, M. Rappaz, R. Trivedi, Acta Materialia 57 (2009) 941. [4] R. Trivedi, K. Somboonsuk, Acta Metallurgica 33 (1985) 1061. [5] H. Esaka, W. Kurz, Journal of Crystal Growth 72 (1985) 578. [6] J.S. Langer, Science 243 (1989) 1150. [7] R.H. Mathiesen, L. Arnberg, R. Mo, T. Weitkamp, A. Snigirev, Physical Review Letters 83 (1999) 5062. [8] R.H. Mathiesen, L. Arnberg, K. Ramsøskar, T. Weitkamp, C. Rau, A. Snigirev, Metallurgical and Materials Transactions B 33 (2002) 613. [9] H. Yasuda, I. Ohnaka, K. Kawasaki, A. Sugiyama, T. Ohmichi, J. Iwane, K. Umetani, Journal of Crystal Growth 262 (2004) 645. [10] H. Yasuda, T. Nagira, M. Yoshiya, N. Nakatsuka, A. Sugiyama, K. Uesugi, K. Umetani, ISIJ International 51 (2011) 402. [11] D. Walton, B. Chalmers, Transactions of the Metallurgical Society of the American Institute of Mining, Metallurgical and Petroleum Engineers 215 (1959) 447. [12] M. Rappaz, C.-A. Gandin, Acta Metallurgica et Materialia 41 (1993) 345. [13] C.A. Gandin, M. Rappaz, Acta Metallurgica et Materialia 42 (1994) 2233. [14] N. D'Souza, M.G. Ardakani, A. Wagner, B.A. Shollock, M. McLean, Journal of Material Science 37 (2002) 481. [15] A. Wagner, B.A. Shollock, M. McLean, Material Science Engineering A 374 (2004) 270. [16] Y.Z. Zhou, A. Volek, N.R. Green, Acta Materialia 56 (2008) 2631. [17] X.B. Meng, Q. Lu, X.L. Zhang, J.G. Li, Z.Q. Chen, Y.H. Wang, Y.Z. Zhou, T. Jin, X.F. Sun, Z.Q. Hu, Acta Materialia 60 (2012) 3965. [18] J. Li, Z. Wang, Y. Wang, J. Wang, Acta Materialia 60 (2012) 1478. [19] A.A. Wheeler, W.J. Boettinger, G.B. McFadden, Physical Review A 45 (1992) 7424. [20] R. Kobayashi, Physica D 63 (1993) 410. [21] A. Karma, W.-J. Rappel, Physical Review E 53 (1996) R3017. [22] A. Karma, W.-J. Rappel, Physical Review E 57 (1998) 4323. [23] A. Karma, Y.H. Lee, M. Plapp, Physical Review E 61 (2000) 3996. [24] S.G. Kim, W.T. Kim, T. Suzuki, Physical Review E 60 (1999) 7186. [25] M. Ode, S.G. Kim, T. Suzuki, ISIJ International 41 (2001) 1076. [26] W.J. Boettinger, J.A. Warren, C. Beckermann, A. Karma, Annual Review of Materials Science 32 (2002) 163. [27] M. Ohno, K. Matsuura, Physical Review E 79 (2009) 031603. [28] M. Ohno, Physical Review E 86 (2012) 051603. [29] W.L. George, J.A. Warren, Journal of Computational Physics 177 (2002) 264. [30] A. Yamanaka, T. Aoki, S. Ogawa, T. Takaki, Journal of Crystal Growth 318 (2011) 40. [31] J.W. Cahn, J.E. Hilliard, Journal of Chemical Physics 28 (1958) 258. [32] S.M. Allen, J.W. Cahn, Acta Metallurgica 27 (1979) 1085. [33] I. Steinbach, C. Beckermann, B. Kauerauf, Q. Li, J. Guo, Acta Materialia 47 (1999) 971. [34] T. Pusztai, G. Bortel, L. Gránásy, Europhysics Letters 71 (2005) 131. [35] I. Steinbach, Acta Materialia 56 (2008) 4965. [36] S. Gurevich, A. Karma, M. Plapp, R. Trivedi, Physical Review E 81 (2010) 011603. [37] Y.L. Tsai, C.C. Chen, C.W. Lan, International Journal of Heat Mass Transfer 53 (2010) 2272. [38] T. Shimokawabe, T. Aoki, T. Takaki, A. Yamanaka, A. Nukada, T. Endo, N. Maruyama, S. Matsuoka, in: Proceedings of SC11, 2011, p. 1.
Contents lists available at ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro
Unexpected selection of growing dendrites by very-large-scale phase-field simulation Tomohiro Takaki a,n, Takashi Shimokawabe b, Munekazu Ohno c, Akinori Yamanaka d, Takayuki Aoki b a
Mechanical and System Engineering, Kyoto Institute of Technology, Kyoto 606-8585, Japan Global Scientific Information and Computing Center, Tokyo Institute of Technology, Tokyo 152-8550, Japan c Division of Materials Science and Engineering, Hokkaido University, Sapporo, Hokkaido 060-8628, Japan d Mechanical Systems Engineering, Tokyo University of Agriculture and Technology, Koganei, Tokyo 184-8588, Japan b
art ic l e i nf o
a b s t r a c t
Article history: Received 17 June 2013 Received in revised form 17 July 2013 Accepted 24 July 2013 Communicated by M. Plapp Available online 6 August 2013
Dendrites are typical growth morphologies in alloy solidification. However, the ways in which dendrites grow and form solidification microstructures remain poorly understood. Here we show unexpected findings that reveal the survival of unfavorably oriented dendrites and highly complicated dendrite– dendrite interactions in three-dimensional space during the directional solidification of a binary alloy. These results are observed for the first time through very-large-scale phase-field computations performed by a graphics processing unit (GPU) supercomputer and a high-performance algorithm developed for parallel computing. & 2013 Elsevier B.V. All rights reserved.
Keywords: A1. Dendrites A1. Directional solidification A1. Computer simulation
1. Introduction Dendrites, which are prime examples of spontaneously generated patterns found in nature, are formed during the solidification of a number of materials such as non-metallic and metallic systems [1,2]. The prediction and control of dendritic microstructures are among the most important issues in the field of solidification science and engineering [3]. This is because the morphologies and sizes of inhomogeneous entities such as grains, segregation patterns, and secondary precipitates in solidification microstructures have a profound influence on several material properties. Advances in experimental technologies such as in-situ observation techniques of non-metallic [4–6] and metallic systems [7–10] have improved our understanding of dendritic growth. However, research to date has been centered on thin samples, which are essentially two-dimensional (2D) in terms of dendrite– dendrite interactions, small systems where only a small number of dendrites grow, or both. Yet the collective behavior of many dendrites in three-dimensional (3D) space is of primary importance in controlling the solidification microstructures of material, because this determines the size and morphology of dendritic microstructures.
n
Corresponding author. Tel.: +81 75 724 7317; fax: +81 75 724 7300. E-mail address: [email protected] (T. Takaki).
0022-0248/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jcrysgro.2013.07.028
Directional solidification is a useful way of evaluating interactions between dendrites. It is generally accepted that a favorably oriented dendrite, for which the inclination angle between the preferred growth direction and heat flow direction is small, can stop the growth of an unfavorably oriented dendrite, for which the inclination angle between the preferred growth direction and heat flow direction is larger than that for the favorably oriented dendrite [11–13]. Hereafter, this is referred to as “expected” selection behavior. Recently, an unexpected phenomenon, in which an unfavorably oriented dendrite can stop the growth of a favorably oriented dendrite, has been observed in both experiments and simulations [14–18]. Although this new unexpected phenomenon is very interesting, it was observed under simple bicrystal conditions. Because actual solidification microstructures are much more complicated, we need to investigate the interactions among multiple dendrites. The phase-field method has emerged as a powerful tool to simulate the complex dynamics of dendritic growth [19–28]. Its ability to describe the formation of complex dendritic patterns is well acknowledged. However, phase-field computations of dendritic growth have thus far been limited to 2D phenomena or 3D growth of a single dendrite [22,23,29,30], because of the large computational cost stemming from the diffuse interface model, where the order parameter that describes the different phases changes smoothly in a thin interface region [31,32]. Although some phase-field simulations of multiple grain growth in 3D have
22
T. Takaki et al. / Journal of Crystal Growth 382 (2013) 21–25
been reported for equiaxed grains [33,34], dendrites [35], and columnar cells [36,37] in directional solidification, these are limited to relatively small spaces and have simple morphologies. Therefore, as far as we know, there are no phase-field simulations for multiple and complicated dendrites in sufficient 3D space to evaluate realistic microstructures. Hence, the key bottleneck here is the development of large-scale computation technologies. We have recently enabled the first-ever petascale phase-field computation in single precision for the directional dendritic solidification of a binary alloy [38]. This has been realized by using 4000 graphics processing units (GPU) along with 16 000 central processing unit (CPU) cores of the TSUBAME 2.0 supercomputer at the Tokyo Institute of Technology. To achieve high performance, we have written our simulation code from scratch in the GPU-specific language CUDA and introduced an overlapping technique that avoids performance degradation due to inter-GPU communication. In this study, we show a very-large-scale 3D phase-field simulation of competitive dendritic growth in the directional solidification of an Al–Si binary alloy by using the GPU supercomputer TSUBAME2.0 and a custom-developed high-performance algorithm. We focus on phenomena captured only in a very large system during a very long time period. Specifically, our simulation were carried out in a system with the dimensions of 3.072 mm 3.078 mm 3.072 mm (4096 4104 4096 meshes) for a total time period of more than 100 s (four million computational steps)—the largest simulation of dendrite growth that have ever been reported—by using 768 GPUs with 768 CPUs in TSUBAME2.0. This very-large-scale simulation in both spatial and temporal domains is indispensable for elucidating the nature of dendritic competitive growth. In fact, we demonstrate an unexpected finding that is obtained only by computations on this scale.
2. Phase-field model The existence of a liquid or solid phase is expressed by a phasefield variable ϕ, which is defined as 0 in a liquid or 1 in a solid and changes smoothly at the interface region. By using ϕ, we define the following Ginzburg–Landau type free energy functional F [20]: Z 2 a F¼ ð∇ϕÞ2 þ WqðϕÞ þ pðϕÞf S þ ð1pðϕÞÞf L dV ð1Þ V 2 where a is the gradient energy coefficient, W is the height of the double-well potential, and fS and fL are the free energy densities in solid and liquid, respectively. The functions pðϕÞ and qðϕÞ are the energy density function and double-well function, respectively; here we use pðϕÞ ¼ ϕ3 ð1015ϕ þ 6ϕ2 Þ and qðϕÞ ¼ ϕ2 ð1ϕÞ2 [19]. The concentration in a binary alloy is defined as c ¼ ϕcS þ ð1ϕÞcL , where cS and cL are the concentrations for the solid and liquid, respectively [24]. These concentrations are made mutually dependent according to cS ¼ kcL with the partition coefficient k. Time evolutions for ϕ and c are described by the Allen–Cahn equation [32] and Cahn–Hilliard equation [31], respectively, as follows: 9 9 8 8 2 > > > > > > > > = = < < 6 2 ∂ϕ ∂ ∂a ∂ ∂a 2 2 þ ð ∇ϕ Þ ð∇ϕÞ ¼ Mϕ 6 a a ∇ a ∇ϕ þ 4 > > ∂ϕ ∂ϕ ∂t ∂x> ∂y> > > > > ; ; : ∂ : ∂ ∂x ∂y 9 8 3 > > > > = ∂< ∂a dpðϕÞ dqðϕÞ7 7; þ a ð∇ϕÞ2 SΔT W ð2Þ > ∂ϕ ∂z> dϕ dϕ 5 > > ; : ∂ ∂z ∂c ð1kÞc ¼ ∇ D ∇c þ ∇ϕ : ∂t 1ϕ þ kϕ
ð3Þ
Here, M ϕ is the phase-field mobility; S is the entropy of fusion; and ΔT is the undercooling expressed by ΔT ¼ T m ðcÞT, where T is the temperature and Tm(c) is the temperature on the liquidus line for concentration c. The diffusion coefficient is expressed by D ¼ DS þ ðDL DS Þð1ϕÞ=ð1ϕ þ kϕÞ where DS and DL are the diffusion coefficients in the solid and liquid, respectively. T changes with time t according to T ¼ GzRt þ T 0 where G is the temperature gradient, R is the cooling rate, z is the spatial coordinate in the heat flow or G direction, and T0 is the temperature at the bottom surface (z¼ 0). To introduce anisotropy of the kinetic coefficient and interface energy, we use 8 4 4 4 9 > ∂ϕ ∂ϕ ∂ϕ > > > > > þ þ < ~ ~ 4ζ ∂ x ∂ y ∂z~ = ~ ð4Þ η ∇ϕ ¼ ð13ζ Þ 1 þ ~ 4 > > 13ζ j∇ϕj > > > > ; : ~ ¼ M ϕ ηð∇ϕÞ ~ ~ ¼ aηð∇ϕÞ, ~ for M ϕ and a to obtain M ϕ ð∇ϕÞ and að∇ϕÞ respectively [29]. Here, M ϕ and a are the standard values of M ϕ ~ y; ~ z~ Þ is a material and a, respectively. The coordinate system ðx; coordinate system in which each axis corresponds to the 〈100〉 direction of a cubic lattice. The following coordinate transforma~ y; ~ z~ Þ: tion is used between the coordinate systems of ðx; y; zÞ and ðx; 8 ∂ϕ 9 2 ∂x ∂y ∂z 38 ∂ϕ 9 > > > ∂x > > = 6 ∂x~ ∂x~ ∂x~ 7> = < ∂x~ > < > ∂y ∂ϕ ∂ϕ ∂x ∂z 7 6 ¼ 4 ∂y~ ∂y~ ∂y~ 5 ∂y ∂y~ > > > > > ∂y ; ; : ∂ϕ > : ∂ϕ > ∂x ∂z > ∂z ∂z~ ∂z~ ∂z~ ∂z~ 2 3 1 0 0 6 cos θx sin θx 7 ¼40 5 0 sin θx cos θx 38 ∂ϕ 9 2 32 > cos θy 0 sin θy cos θz sin θz 0 > > = < ∂x > 7 ∂ϕ 6 76 0 1 0 ð5Þ 4 54 sin θz cos θz 0 5 ∂y > > > > ; sin θy 0 cos θy 0 0 1 : ∂ϕ ∂z where θx , θy , and θz are the rotation angles around the x, y and z axes, respectively. The preferred growth directions, which are different for every nucleus, are set to the lattice points with ϕ Z105 under the initial conditions. During growth, if the phase-field variable ϕ of a lattice point without the preferred growth direction, which is the nearest lattice point to a lattice point with a preferred growth direction, becomes larger than 10 5, the lattice point has the same preferred growth direction as that of the nearest lattice point. In this way, although two neighboring dendrites with different preferred growth directions sometimes impinge on each other at the low-temperature side with respect to the dendrite tip positions of the primary arm, this does not affect competitive dendrite growth, on which we focus in this study. a, W, and M ϕ in Eq. (2) can correspondingly be related to the interface thickness interface pffiffiffiffiffiffiffiffiffiffiffiffiffi δ, interface energy γ, and pffiffiffiffiffiffiffiffi mobility M by a ¼ 3δγ=b, W ¼ 6γb=δ, and M ϕ ¼ μ 2W =ð6aSÞ, respectively, where μ is the kinetic coefficient and b is a constant related to the interface thickness. Because, in this study, we aimed at very-large-scale phase-field simulation, we used a conventional model without an antitrapping current term [22,23,27,28], to reduce the computational costs associated with the antitrapping current. Therefore, the following results are not quantitative. However, we believe that the results give us important information on competitive growth between multiple dendrites and open up new roads for the investigation of microstructural growth mechanisms by large-scale phase-field simulations. 3. Numerical conditions Directional solidification of the Al–Si binary alloy is simulated in a computational domain of 3.072 mm 3.078 mm 3.072 mm
T. Takaki et al. / Journal of Crystal Growth 382 (2013) 21–25
divided into 4096 4104 4096 finite difference meshes with a uniform grid spacing of Δx ¼ Δy ¼ Δz ¼ 0:75 μm. The initial liquid concentration is set to an atomic fraction of 0.05 and the corresponding liquidus temperature is calculated as 903.47 K from T m ðcÞ ¼ mcþ T m , where m ¼ –600 K/at.frac. is the liquidus slope and Tm ¼933.47 K is the melting temperature of pure Al. Other material properties and computational conditions used in the present simulation are as follows: temperature gradient in z-direction G¼10 K/mm; cooling rate R¼ 0.3 K/s; temperature at z¼0 surface T0 ¼901.47 K; interface thickness δ ¼ 6Δx ¼ 4:5 μm; entropy of fusion S¼ 8.0 105 J/Km3; interface energy γ ¼ 0:16 J=m2 ; interface kinetic coefficient μ ¼ 0:32 103 m=Ks; diffusivity of solid DS ¼1 10 12 m2/s; diffusivity of liquid DL ¼3 10 9 m2/s; partition coefficient k¼0.13; and strength of anisotropy ζ ¼ 0:03. As for the boundary conditions, the periodic boundary condition was set to the edges of x and y, and the Neumann condition was employed at the edges of z for the phase-field and concentration. The time increment Δt was set to Δt ¼ 2:68 105 s. Four million computation steps, which correspond to 107 s in real time, were then performed. The 128 spherical solid seeds of radii 3Δx with random preferred growth directions are placed on the bottom surface (z ¼0). Here, the inclination angle θ for the heat flow direction is expressed by θy because all θx are set to zero just for convenience. We confirm that special textures were not formed by setting θx ¼ 0. Higher-order branching is formed by 10% fluctuation of the concentration. The present simulation typically takes 2 h 59 min to perform a 40 000 time step simulation, including data output. The total simulation runtime for obtaining numerical results is about 12 days. The performance we have achieved in this simulation is approximately 187 TFlops. Note that this performance is lower than we expected from the results of a previous simulation, i.e., 2PFlops [38]. Runtime configurations, such as the decomposed domain size, and the GPU kernel configuration used in this simulation are not sufficiently optimized to achieve the desired performance, and network resources are shared with other TSUBAME users, which probably result in degradation of the simulation performance.
23
4. Numerical results Fig. 1 shows the time slices of competitive dendritic growth. In the simulation, 128 solid seeds initially nucleated at the bottom wall grow into star-like particles (Fig. 1a), and then cover the bottom surface by lateral growth (Fig. 1b). After covering the entire bottom surface, the dendrites begin growing in a longitudinal direction that corresponds to the heat flow direction with simultaneous growth of branched dendrites, viz. the secondary arms (Fig. 1c). Competitive growth of dendrites occurs, as is reflected in the selection of growing dendrites (Fig. 1d). The number of seeds gradually decreases as the solidification front proceeds upward. When the front is located at 1.2 mm from the bottom (Fig. 1e), almost half of the original 128 seeds stop growing. The number of seeds further decreases and only 46
Fig. 2. The tip positions of 128 dendrites in the z-direction at 96.5 s. The dendrites are categorized into three groups A, B and C: (A) those that ceased growing in a region where the number of stopped dendrites increases rapidly with the tip distance along the z-coordinate, (B) those that ceased growing in a region where the number of stopped dendrites increases almost proportionally with tip distance, and (C) those that survived to 96.5 s.
Fig. 1. Time slices of competitive dendritic growth. The color indicates the inclination angle θ of the preferred growth direction 〈100〉 relative to the x-axis. (a) 128 spherulitic seeds change in morphology to star-like dendrites (1.1 s). (b) Dendrites preferentially grow in the lateral direction or wet the bottom wall (5.4 s). (c) Dendritic arms start to grow in heat flow direction (10.7 s). (d) Competitive growth occurs where the number of solid seeds decreases approximately in proportion to the square of the distance from the bottom (26.8 s). (e) Competitive growth continuously occurs where the number of seeds decreases approximately in proportion to the distance from the bottom (53.6 s). (f) Final morphology of the present simulation (107.1 s). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)
24
T. Takaki et al. / Journal of Crystal Growth 382 (2013) 21–25
ultimately survive after 107.1 s (Fig. 1f). The 3D simulations thus reproduce the full range of directional solidification phenomena from heterogeneous nucleation on the mold wall to competitive growth between multiple dendrites. Fig. 2 shows the dendrite tip position (growth front in heat flow direction) at 96.5 s for all 128 dendrites, each of which is specified by the seed number on the left. The dendrites are categorized into three groups: A, B and C. The dendrites that stop growing before 96.5 s are categorized separately from those still growing after 96.5 s. Of the former, dendrites categorized in A that do not grow beyond about 1.2 mm from the bottom rapidly decrease in number as the distance from the bottom increases. As shown in Table 1, in this category are 65 dendrites, half (32 dendrites) of which have a preferential growth direction largely inclined from the heat flow direction (A4) (i.e., inclination angle θ ¼ 301–451). This is predictable from a generally accepted dendritic selection mechanism in
Table 1 The number of dendrites categorized by dendrite tip position and inclination angle θ. Categories A, B, and C correspond to those in Fig. 2. θ ð1Þ
A
0–10 10–20 20–30 30–45 Total
(A1) (A2) (A3) (A4) 65
B 9 11 13 32
(B1) (B2) (B3) (B4) 17
C 2 5 6 4
(C1) (C2) (C3) (C4) 46
Total 23 15 6 2
34 31 25 38 128
which only dendrites with small inclination angles grow to survive [11–13]. At the same time, more than 10% of the dendrites in this first category A1 (9 dendrites) actually have a small inclination angle (θ ¼ 01–101). Among these, 3 dendrites are prevented from growing by neighboring dendrites with a smaller inclination angle; however, the remaining 6 are stopped by those with a larger inclination angle. Table 1 reveals that at the very beginning of solidification, the growth of more than 15% of all dendrites with a small inclination angle (θ ¼ 01–101) is actually blocked by neighboring dendrites with a larger angle. A second category B of dendrites that stop growing comprise 17 dendrites whose tip distances are below about 2.7 mm, where the solidification process should be close to the steady-state regime. Of these, 15 have an inclination angle of 101–451 and are stopped from growing by dendrites with a smaller inclination angle. The 2 dendrites in B1 with a small inclination angle (θ ¼ 01–101) are stopped by dendrites with a larger inclination angle. The growth of these dendrites is blocked through the interaction of more than three dendrites. This phenomenon is beyond the explanatory scope of the recently proposed 2D selection mechanism, which accounts for how unfavorably oriented dendrites can overgrow favorably oriented dendrites in both experimental and simulation studies [14–18]. Of the dendrites in the third category C (i.e., those still growing after 96.5 s), half have an inclination angle of 01–101 (C1). The number of dendrites decreases with increasing inclination angle. This is a reasonable result assuming a conventionally accepted dendrite selection mechanism; yet even in this category, 2 dendrites in C4 with a large inclination angle of 301–451 survive.
Fig. 3. Dendritic interactions around (a) dendrite A1 and (b) dendrite B1 categorized as C4 in Table 1. The inclination angle (in degrees) for each dendrite is indicated in parentheses.
T. Takaki et al. / Journal of Crystal Growth 382 (2013) 21–25
Why and how the 2 dendrites (A1, θ ¼ 35:11; B1, θ ¼ 42:51) continuously grow are shown in Fig. 3, which presents side and top views of these dendrites along with their surrounding dendrites. As shown in Fig. 3a, the leftward growth of A1 is stopped by A2 and A3 in the upper and lower parts, respectively; meanwhile, A1 grows rightward, successively stopping A4 (θ ¼ 12:91), A5 (θ ¼ 34:11), and A6 (θ ¼ 13:81), which have a smaller inclination angle. This selection process is also seen in the growth of A2 (θ ¼ 27:71), which stops the growth of a dendrite with a smaller inclination angle (A3, θ ¼ 7:31). The favorably oriented dendrites that grow adjacent to A1 (i.e., A7–A15) are unable to stop the growth of A1 because they grow almost parallel to or away from it. In Fig. 3b, B1 is stopped on both sides by favorable dendrites B3, B4, and B5. The surrounding dendrites B5, B6, B8, and B9 grow almost parallel to B1 and do not prevent its growth. As shown above, these 2 dendrites continue to grow in a similar way with a plate-like structure. In addition, B2 is observed to penetrate B1; A2 and A16 also penetrate one another. Thus, dendrite–dendrite interaction is a highly complicated phenomenon. Our simulations demonstrate that contrary to commonly held assumptions, the dendrite selection process cannot be discussed solely on the basis of single parameter, viz., the inclination angle of the preferential growth orientation with respect to the direction of heat flow. The detailed 3D morphology of the dendrites is a non-trivial factor that strongly affects dendrite–dendrite interactions. 5. Conclusions Our results demonstrate a series of directional solidification steps from nucleation to competitive growth in 3D for the first time. This very-large-scale simulation was achieved by multi-GPU computation with a high-performance algorithm best suited to phase-field computation. Our large-scale simulations have shown that a relatively large amount of unfavorably oriented dendrites can survive during unidirectional solidification. This finding has profound significance in controlling as-solidified microstructures, especially in terms of their crystallographic features. We believe that the present work constitutes an important step in understanding and developing solidification science, based on verylarge-scale phase-field simulations. Acknowledgments We are grateful for the opportunity to use the TSUBAME2.0 at the Tokyo Institute of Technology. We thank Dr. A. Nukada, Tokyo Institute of Technology, for fruitful comments on supercomputing. This study was financially supported in part by JSPS KAKENHI Grant Number 25289006 and the ISIJ Research Promotion Grant from the Iron and Steel Institute of Japan.
25
Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version of http://dx.doi.org/10.1016/j.jcrysgro.2013.07.028.
References [1] J.S. Langer, Reviews of Modern Physics 52 (1980) 1. [2] S.A. David, T. DebRoy, Science 257 (1992) 497. [3] M. Asta, C. Beckermann, A. Karma, W. Kurz, R. Napolitano, M. Plapp, G. Purdy, M. Rappaz, R. Trivedi, Acta Materialia 57 (2009) 941. [4] R. Trivedi, K. Somboonsuk, Acta Metallurgica 33 (1985) 1061. [5] H. Esaka, W. Kurz, Journal of Crystal Growth 72 (1985) 578. [6] J.S. Langer, Science 243 (1989) 1150. [7] R.H. Mathiesen, L. Arnberg, R. Mo, T. Weitkamp, A. Snigirev, Physical Review Letters 83 (1999) 5062. [8] R.H. Mathiesen, L. Arnberg, K. Ramsøskar, T. Weitkamp, C. Rau, A. Snigirev, Metallurgical and Materials Transactions B 33 (2002) 613. [9] H. Yasuda, I. Ohnaka, K. Kawasaki, A. Sugiyama, T. Ohmichi, J. Iwane, K. Umetani, Journal of Crystal Growth 262 (2004) 645. [10] H. Yasuda, T. Nagira, M. Yoshiya, N. Nakatsuka, A. Sugiyama, K. Uesugi, K. Umetani, ISIJ International 51 (2011) 402. [11] D. Walton, B. Chalmers, Transactions of the Metallurgical Society of the American Institute of Mining, Metallurgical and Petroleum Engineers 215 (1959) 447. [12] M. Rappaz, C.-A. Gandin, Acta Metallurgica et Materialia 41 (1993) 345. [13] C.A. Gandin, M. Rappaz, Acta Metallurgica et Materialia 42 (1994) 2233. [14] N. D'Souza, M.G. Ardakani, A. Wagner, B.A. Shollock, M. McLean, Journal of Material Science 37 (2002) 481. [15] A. Wagner, B.A. Shollock, M. McLean, Material Science Engineering A 374 (2004) 270. [16] Y.Z. Zhou, A. Volek, N.R. Green, Acta Materialia 56 (2008) 2631. [17] X.B. Meng, Q. Lu, X.L. Zhang, J.G. Li, Z.Q. Chen, Y.H. Wang, Y.Z. Zhou, T. Jin, X.F. Sun, Z.Q. Hu, Acta Materialia 60 (2012) 3965. [18] J. Li, Z. Wang, Y. Wang, J. Wang, Acta Materialia 60 (2012) 1478. [19] A.A. Wheeler, W.J. Boettinger, G.B. McFadden, Physical Review A 45 (1992) 7424. [20] R. Kobayashi, Physica D 63 (1993) 410. [21] A. Karma, W.-J. Rappel, Physical Review E 53 (1996) R3017. [22] A. Karma, W.-J. Rappel, Physical Review E 57 (1998) 4323. [23] A. Karma, Y.H. Lee, M. Plapp, Physical Review E 61 (2000) 3996. [24] S.G. Kim, W.T. Kim, T. Suzuki, Physical Review E 60 (1999) 7186. [25] M. Ode, S.G. Kim, T. Suzuki, ISIJ International 41 (2001) 1076. [26] W.J. Boettinger, J.A. Warren, C. Beckermann, A. Karma, Annual Review of Materials Science 32 (2002) 163. [27] M. Ohno, K. Matsuura, Physical Review E 79 (2009) 031603. [28] M. Ohno, Physical Review E 86 (2012) 051603. [29] W.L. George, J.A. Warren, Journal of Computational Physics 177 (2002) 264. [30] A. Yamanaka, T. Aoki, S. Ogawa, T. Takaki, Journal of Crystal Growth 318 (2011) 40. [31] J.W. Cahn, J.E. Hilliard, Journal of Chemical Physics 28 (1958) 258. [32] S.M. Allen, J.W. Cahn, Acta Metallurgica 27 (1979) 1085. [33] I. Steinbach, C. Beckermann, B. Kauerauf, Q. Li, J. Guo, Acta Materialia 47 (1999) 971. [34] T. Pusztai, G. Bortel, L. Gránásy, Europhysics Letters 71 (2005) 131. [35] I. Steinbach, Acta Materialia 56 (2008) 4965. [36] S. Gurevich, A. Karma, M. Plapp, R. Trivedi, Physical Review E 81 (2010) 011603. [37] Y.L. Tsai, C.C. Chen, C.W. Lan, International Journal of Heat Mass Transfer 53 (2010) 2272. [38] T. Shimokawabe, T. Aoki, T. Takaki, A. Yamanaka, A. Nukada, T. Endo, N. Maruyama, S. Matsuoka, in: Proceedings of SC11, 2011, p. 1.