- Email: [email protected]
Implementation of high interfacial energy anisotropy in phase field simulations
Scripta Materialia 54 (2006) 1919–1924 www.actamat-journals.com
Implementation of high interfacial energy anisotropy in phase field simulations Ning Ma, Qing Chen 1, Yunzhi Wang
*
Department of Materials Science and Engineering, The Ohio State University, 2041 College Road, Columbus, OH 43210, USA Received 10 May 2005; received in revised form 7 January 2006; accepted 3 February 2006 Available online 3 March 2006
Abstract We show changes of equilibrium shape and artifacts of enhanced particle coalescence and numerical pinning, introduced by boundary thickness variation associated with anisotropy in interfacial energy in phase field simulations. Methods to implement high anisotropy in interfacial energy without altering boundary thickness are presented. 2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Phase field; Anisotropy; Interfacial energy; Boundary width
1. Introduction In materials, it is more common than not, for the thermodynamic and kinetic properties of hetero-phase interfaces and grain boundaries to vary from one boundary to another, as well as from one inclination to another along the same boundary. Such anisotropies in boundary properties are important factors that determine the morphology of a multi-phase or a polycrystalline microstructure and the kinetics of its evolution [1,2]. Recently there have been increasing efforts to incorporate the boundary energy and mobility anisotropies in the widely used phase field method of microstructural evolution [3–13]. Caginalp and Fife [14] introduced a quadratic form for the gradient-energy coefficient as a function of interface inclination. McFadden and coworkers [15] employed a more general dependence of the gradient-energy coefficient on inclination. When the shape of a c-plot (the polar representation of the interfacial energy as a function of interface orientation) is non-convex, the time-dependent Ginzburg–Landau (TDGL) equa*
Corresponding author. E-mail address: [email protected] (Y. Wang). 1 Permanent address: Thermo-Calc Software, Stockholm Technology Park, SE-113 47, Stockholm, Sweden.
tion (also called the Allen–Cahn equation) becomes illposed. Eggleston and coworkers [16] have solved this problem by modifying the interfacial energy function without changing the equilibrium shape of a particle. Debierre and coworkers [17] simulated faceted interfaces by using an approximated c-plot with rounded cusps. In these models [14–17], only the gradient-energy coefficient is formulated as a function of the local inclination. One consequence of this is that the interface thickness increases linearly with the interfacial energy. This may not be a problem for solidification because the anisotropy of a solid–liquid interface is limited. For solid–solid interfaces, however, the energies of interfaces of different inclinations could differ by orders of magnitude [18]. In this case the boundary width will differ by orders of magnitude, leading to either artifacts of numerical pinning of the low energy interfaces or enhanced coalescence of the highenergy interfaces when uniform mesh is used [19,20]. In addition, curvatures that have a radius commensurable with the boundary width may not be accounted for accurately in the phase field method. Even though adaptive meshing may partially solve the problem, this will increase dramatically the amount of calculation and will not eliminate the unrealistically large variation in interface thickness. Therefore, proper implementation of interfacial
1359-6462/$ - see front matter 2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.scriptamat.2006.02.005
1920
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
energy anisotropy in the phase field method requires a selfconsistent adjustment of boundary width in accordance with the variation of the gradient-energy coefficient with boundary inclination. In this paper we present an extension of the previous models by allowing both the gradient-energy coefficient and the expansion coefficients in the Landau free energy polynomial to vary with boundary inclination. In principle, one may change the interfacial energy and boundary width independently. Without losing generality, we keep the interfacial width constant while varying the interfacial energy. Using simple examples we demonstrate that variation in interface thickness has a significant effect on the equilibrium particle shape and particle growth kinetics. Applications of the approach to grain growth in systems of anisotropic boundary energy with respect to both inclination and misorientation can be found in Refs. [9,12]. 2. Phase field model In the phase field model, an arbitrary microstructure is described by continuum fields. For simplicity, consider a system whose microstructure is characterized fully by a non-conserved field, g(r). The total free energy of the system is given by: Z j2 F ¼ Wf ðgÞ þ jrgj2 dr ð1Þ 2 v where f(g) is the local free energy, which is chosen in this study to have a simple double well form, i.e., f(g) = (1 g)2 g2, W is the height of the hump between the two wells, and j2 is the gradient-energy coefficient. The temporal evolution of g(r) follows the TDGL equation og dF ¼ L ð2Þ ot dg where L is the kinetic coefficient characterizing interface mobility. For a planar interface, the steady-state profile of g(r) is given by [15,16]: pffiffiffiffiffi 1 x W g ¼ tanh pffiffiffi þ1 ð3Þ 2 2 2j If we assume that within the interfacial region g(r) changes from 0.1 to 0.9, the interface thickness can be estimated as: pffiffiffi j k 2:2 2 pffiffiffiffiffi ð4Þ W This solution shows that the interface thickness is propor. The interfacial energy is given by [21]: tional to pjffiffiffi W Z 1 og r ¼ j2 dx ð5Þ ox 1 Substituting Eq. (3) into Eq. (5) and carrying out the integration, we have pffiffiffiffiffi Wj r ¼ pffiffiffi ð6Þ 3 2
Eq. (6) shows pffiffiffiffiffithat the interfacial energy depends on the product of W j. To describe anisotropy in interfacial energy or interface width, it is necessary to allow both the height of the double well hump, W, and the gradientenergy coefficient, j2, to vary as functions of the inclination angle given by: gy H ¼ arctan ð7Þ gx Then the total free energy variation with respect to g yields: h i h i oH oH o W ðHÞ f ðgÞ o W ðHÞ f ðgÞ ogx ogy dF ¼ W ðHÞfg ox oy dg j2 r2 g þ 2jjH rH rg
þ jjH j2 þ jjHH gx Hy gx Hy
ð8Þ
In the above equations and the ones that follow, the subscripts denote partial derivatives, and gy oH ¼ 2 ogx gx þ g2y oH g ¼ 2 x 2 ogy gx þ gy
ð9Þ
Mobility anisotropy can also be introduced in a similar way by introducing an inclination-dependent mobility coefficient, e.g., L = L(H) [9,15,22]. The above analysis is quite general and is not confined to any particular H-dependence of the interfacial energy and interface width. In the present study, we assume a two-fold symmetry for the interfacial energy as a function of H unless mentioned separately: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rðHÞ ¼ r0 sin2 ðHÞ þ a cos2 ðHÞ ð10Þ The invariance of the interface width with respect to the inclination then requires: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi W ðHÞ ¼ W 0 sin2 ðHÞ þ a cos2 ðHÞ ð11Þ and j2 ðHÞ ¼ j20
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sin2 ðHÞ þ a cos2 ðHÞ
ð12Þ
In the above equations, a, r0, W0, j0 are constants. The evolution of g can be computed by substituting Eqs. (8), (9), (11), (12) to Eq. (2). According to Wulff construction, the equilibrium shape is p anffiffiffi ellipse with the semi-major to semi-minor axis ratio of a. To examine the growth of precipitates during, for example, a massive transformation, one can introduce a volumetric driving force into the free energy, e.g., following Kobayashi [3], Eq. (1) becomes: Z j2 2 2 3 F ¼ Wf ðgÞ þ jrgj þ mð10 15g þ 6g Þg dr ð13Þ 2 v where m is a constant related to the magnitude of the driving force for precipitate growth.
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
In our simulations, the parameters in the above equations are chosen as the following (unless mentioned separately, all are in reduced units): j20 ¼ 4, W0 = 2, a = 25, L = 1 and m = 0.5. The time step and grid size are, respectively, Dt = 0.001 and d0 = 1. 3. Equilibrium shape of a precipitate We first validate our model by comparing the simulated equilibrium shape of an isolated particle with the one determined by the Wulff construction. The simulation was carried out for a two-dimensional (2D) system of 128 128d 20 with a periodic boundary condition. Initially, a circular particle with r = 20d0 was placed at the center of the computational cell. Following Ward [23], we used a volume conserved TDGL equation to compute the equilibrium shape. When the aspect ratio of the precipitate shape became invariant with time, equilibrium was assumed. The isograms of g profiles representing the equilibrium precipitate shape are shown in Fig. 1(a) by the solid lines. For comparison the equilibrium shape obtained from the Wulff construction is also plotted in the figure given by the dashed line. It can be seen that the result from the phase
1921
field simulation agrees well with that from the Wulff construction. Parallel simulations were carried out by adopting only anisotropy in the gradient-energy coefficient and the results for two different values of j2 are shown in Fig. 1(b) and (c). It can be readily seen that the interface thickness is non-uniform and the shape defined by the contour line at g = 0.5 deviates significantly from the equilibrium one given by the dashed line. In fact, when the interface thickness is comparable to the tip radius, it is no longer valid to describe the curvature of the tip using the phase field method. In principle one can reduce the ratio of interfacial width over tip radius by either refining the numerical mesh or sharpening the interface. This would either reduce the computational efficiency or introduce numerical pinning to interfacial locations of the smallest width. Any trial-and-error approach to balance these factors will become more and more impractical as the degree of anisotropy increases. 4. Kinetics of particle growth and coalescence In this example we consider the growth of two secondphase particles from a parent phase matrix. The system size
Fig. 1. Comparison of equilibrium particle shapes obtained by phase field simulations (thin solid lines) and by Wolff construction (thick dotted lines) with (a) constant and (b) and (c) varying interface thickness (j2 = 0.5 in (c)). Considering the two-fold shape of the particle, only a portion of the computational cell is shown.
1922
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
employed is 256 256d 20 . Initially, two spherical precipitates with r = 10d0 were introduced into the system. Two parallel simulations were carried out with constant and varying interfacial thickness, respectively. Snapshots of the temporal evolution of the two particles in the two cases are shown in Fig. 2. The corresponding growth kinetics of the particles is given in Fig. 3. It can be seen that at the same moment the interfaces of the two growing particles start to overlap with each other in the case of varying boundary width (Fig. 2(b)), while they are well separated in the case of constant boundary width (Fig. 2(a)). Therefore, particle coalescence occurs at an earlier moment in the case of varying interface width. This is also reflected in Fig. 3, P where the total area of the product phase (calculated by r gðrÞ) as a function of time shows an abrupt change when the two particles coalesce. Thus a large variation of interface thickness in phase field simulations could result in significant errors in the overall growth kinetics. Since in the simulation the driving force for the growth
Fig. 3. Growth and coalescence kinetics of the two particles shown in Fig. 2.
of the particles were provided by a difference in the bulk free energy between the precipitate phase and the matrix phase, the shapes of the precipitates are quite different from the equilibrium one shown in Fig. 1. 5. Particles of sharp corners and faceted shapes The equilibrium shape of a precipitate may exhibit sharp corners and facets if the c-plot of a system is non-convex or has cusps. Missing orientations and facets have been studied recently using the phase field method [16,17], where only anisotropy in the gradient-energy coefficient is considered. In the following examples we examine the variation in boundary width at sharp corners and its effect on equilibrium shape and kinetics of shape evolution. An interfacial energy anisotropy function with four-fold symmetry was employed: rðHÞ ¼ r0 ð1 þ 0:5 cos 4HÞ
Fig. 2. Contour plots of intermediate microstructures of two growing particles at reduce time s = 120 with (a) constant and (b) varying interface thickness.
ð14Þ
As shown in Fig. 4, the polar plot of Eq. (14) has round cusps at 45 orientation. Orientations smaller than 41.7 are missing within the range of 0 < H < p/4. For convexification of the polar plot of the inverse of the interfacial energy the approach proposed by Eggleston et al. [16] was used. Two parallel simulations were carried out, one with anisotropy in the gradient-energy term only and one with anisotropy in both the gradient-energy term and the freeenergy hump following Eqs. (8) and (9). Initially, a circular particle of r = 30d0 was placed in the center of a matrix of 128d0 · 128d0. The invariant particle shapes after long time relaxation in the two cases are shown in Fig. 5(a) and (b). The enlarged contours of corners in the two cases are shown in Fig. 5(c) and (d). While smoothly round corners are produced in the case of uniform boundary width, kinks are observed at the locations where the boundary width changes in the case of varying boundary width. These kinks pin the interface and prevent the particle evolving further after initial relaxation, leading to a shape significantly different from the one assumed by the Wulff construction (the dotted line). In the case of uniform interface width the
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
Fig. 4. Polar plot of normalized interfacial energy, r(H)/r0, as a function of interface inclination.
1923
equilibrium particle shape obtained is close to a square, i.e., the equilibrium shape following the Wulff construction, but with rounded corners. With the presence of the kinks, the plane normal is erroneously calculated by the centered finite differencing formula employed in the numerical calculations. As a consequence, the interfacial energy cannot be assigned correctly. This problem may be avoided either by aligning the inclination of the interfaces with the numerical grid orientation or by applying a one-sided finite difference scheme [16]. The former is employed in this study. The equilibrium particle shapes obtained for the two cases are given in Fig. 6(a) and (b). In both cases, the equilibrium shape obtained from the phase field simulation is close to the one given by the Wulff construction. It is obvious that simulations must be carried out with great caution in the case of varying boundary width because the simulation results become sensitive to the relative orientation of the interface inclination with respect to the numerical grids.
Fig. 5. Particle shape under (a) varying and (b) constant boundary width (j2 = 8 in both cases); (c) and (d) are the corresponding enlarged contours of corners. The dotted lines show the equilibrium shape from Wulff construction. The arrows indicate the orientation of the contour normal calculated using centered finite difference. In both cases the starting particle has a circular shape.
1924
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
be a straightforward way to alter the free energy hump [24,25]. In this case, the approach presented in Refs. [26– 28] that treats the interface as a mixture of the two equilibrium phases could be more convenient in implementing high interfacial energy anisotropy. A quantitative phase field study on side-plate development in Ti–6Al–4V using such an approach will be presented in a separate paper. Acknowledgement This work is supported by the US Air Force through the Metals Affordability Initiative under contract F33615-99-25215. References
Fig. 6. Particle shape under (a) varying and (b) constant boundary width (j2 = 8 in both cases) obtained when boundary inclinations are aligned with the numerical grid orientations.
6. Summary We have demonstrated through several simple examples that when simulating microstructural evolution in systems of high interfacial energy anisotropy using the phase field method, it is necessary to adjust the double well hump of the local free energy and the gradient-energy coefficient simultaneously to maintain a reasonable interface width. Otherwise the artificially large variation in boundary width will alter both the morphology of the microstructure and the kinetics of its evolution. For congruent transformations or grain growth where the local free energy is a function of non-conserved fields only, the method discussed seems to be able to deal effectively with high interfacial energy anisotropy. If the local free energy depends on both conserved and non-conserved fields, however, there may not
[1] Gottstein G, Shvindlerman LS. Grain boundary migration in metal. New York (NY): CRC Press; 1999. [2] Sutton AP, Balluffi RW. Interfaces in crystalline material. New York (NY): Oxford University Press; 1995. [3] Kobayashi R. Exp Math 1994;3:59. [4] Warren JA, Boettinger WJ. Acta Metall Mater 1995;43:689–703. [5] Wheeler AA, McFadden GB, Boettinger WJ. Proc Roy Soc London Ser A 1996;452:495–525. [6] Karma A, Rappel WJ. Phys Rev E 1996;53:R3017. [7] Khachaturyan AG. Philos Mag A 1996;74:3. [8] Braun RJ, Cahn JW, McFadden GB, Wheeler AA. Philos Trans R Soc London 1997;355:1787; Braun RJ, Cahn JW, McFadden GB, Wheeler AA. Acta Metall 1995; 43:689. [9] Kazaryan A, Wang Y, Dregia SA, Patton BR. Phys Rev B 2000;61: 14275. [10] Kazaryan A, Wang Y, Dregia SA, Patton BR. Phys Rev B 2001;63: 184102. [11] Kazaryan A, Patton BR, Dregia SA, Wang Y. Acta Mater 2002;50: 499. [12] Kazaryan A, Wang Y, Dregia SA, Patton BR. Acta Mater 2002;50: 2491. [13] Ma N, Kazaryan A, Dregia SA, Wang YZ. Acta Mater 2004;52:3869. [14] Caginalp G, Fife P. Phys Rev B 1986;34:4940. [15] McFadden GB, Wheeler AA, Braun RJ, Coriell SR, Sekerka RF. Phys Rev E 1993;48:2016. [16] Eggleston JJ, McFadden GB, Voorhees PW. Physica D 2001;150: 91–103. [17] Debierre JM, Karma A, Celestini F, Guerin R. Phys Rev E 2003;68: 041604. [18] Aaronson HI, Laird C. Ford Motor Co. Scientific Laboratory Report, 1967. [19] Vaithyanathan V, Wolverton C, Chen LQ. Acta Mater 2004;52: 2973–87. [20] Loginova I, Agren J, Amberg G. Acta Mater 2004;52:4055–63. [21] Allen SM, Cahn JW. Acta Metall 1979;27:1085–95. [22] Wang SL, Sekerka RF. Phys Rev B 1996;53:3760. [23] Ward MJ. SIAM J Appl Math 1996;56:1247–79. [24] Shen C, Chen Q, Wen YH, Simmons JP, Wang Y. Scripta Mater 2004;50:1023. [25] Shen C, Chen Q, Wen YH, Simmons JP, Wang Y. Scripta Mater 2004;50:1029. [26] Grafe U, Bottger B, Tiaden J, Fries SG. Scripta Mater 2000;42:1179. [27] Tiaden J, Nestler B, Diepers HJ, Steinbach I. Physica D 1998;115:73. [28] Kim SG, Kim WT, Suzuki T. Phys Rev D 1999;60:7186.
Implementation of high interfacial energy anisotropy in phase field simulations Ning Ma, Qing Chen 1, Yunzhi Wang
*
Department of Materials Science and Engineering, The Ohio State University, 2041 College Road, Columbus, OH 43210, USA Received 10 May 2005; received in revised form 7 January 2006; accepted 3 February 2006 Available online 3 March 2006
Abstract We show changes of equilibrium shape and artifacts of enhanced particle coalescence and numerical pinning, introduced by boundary thickness variation associated with anisotropy in interfacial energy in phase field simulations. Methods to implement high anisotropy in interfacial energy without altering boundary thickness are presented. 2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Phase field; Anisotropy; Interfacial energy; Boundary width
1. Introduction In materials, it is more common than not, for the thermodynamic and kinetic properties of hetero-phase interfaces and grain boundaries to vary from one boundary to another, as well as from one inclination to another along the same boundary. Such anisotropies in boundary properties are important factors that determine the morphology of a multi-phase or a polycrystalline microstructure and the kinetics of its evolution [1,2]. Recently there have been increasing efforts to incorporate the boundary energy and mobility anisotropies in the widely used phase field method of microstructural evolution [3–13]. Caginalp and Fife [14] introduced a quadratic form for the gradient-energy coefficient as a function of interface inclination. McFadden and coworkers [15] employed a more general dependence of the gradient-energy coefficient on inclination. When the shape of a c-plot (the polar representation of the interfacial energy as a function of interface orientation) is non-convex, the time-dependent Ginzburg–Landau (TDGL) equa*
Corresponding author. E-mail address: [email protected] (Y. Wang). 1 Permanent address: Thermo-Calc Software, Stockholm Technology Park, SE-113 47, Stockholm, Sweden.
tion (also called the Allen–Cahn equation) becomes illposed. Eggleston and coworkers [16] have solved this problem by modifying the interfacial energy function without changing the equilibrium shape of a particle. Debierre and coworkers [17] simulated faceted interfaces by using an approximated c-plot with rounded cusps. In these models [14–17], only the gradient-energy coefficient is formulated as a function of the local inclination. One consequence of this is that the interface thickness increases linearly with the interfacial energy. This may not be a problem for solidification because the anisotropy of a solid–liquid interface is limited. For solid–solid interfaces, however, the energies of interfaces of different inclinations could differ by orders of magnitude [18]. In this case the boundary width will differ by orders of magnitude, leading to either artifacts of numerical pinning of the low energy interfaces or enhanced coalescence of the highenergy interfaces when uniform mesh is used [19,20]. In addition, curvatures that have a radius commensurable with the boundary width may not be accounted for accurately in the phase field method. Even though adaptive meshing may partially solve the problem, this will increase dramatically the amount of calculation and will not eliminate the unrealistically large variation in interface thickness. Therefore, proper implementation of interfacial
1359-6462/$ - see front matter 2006 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.scriptamat.2006.02.005
1920
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
energy anisotropy in the phase field method requires a selfconsistent adjustment of boundary width in accordance with the variation of the gradient-energy coefficient with boundary inclination. In this paper we present an extension of the previous models by allowing both the gradient-energy coefficient and the expansion coefficients in the Landau free energy polynomial to vary with boundary inclination. In principle, one may change the interfacial energy and boundary width independently. Without losing generality, we keep the interfacial width constant while varying the interfacial energy. Using simple examples we demonstrate that variation in interface thickness has a significant effect on the equilibrium particle shape and particle growth kinetics. Applications of the approach to grain growth in systems of anisotropic boundary energy with respect to both inclination and misorientation can be found in Refs. [9,12]. 2. Phase field model In the phase field model, an arbitrary microstructure is described by continuum fields. For simplicity, consider a system whose microstructure is characterized fully by a non-conserved field, g(r). The total free energy of the system is given by: Z j2 F ¼ Wf ðgÞ þ jrgj2 dr ð1Þ 2 v where f(g) is the local free energy, which is chosen in this study to have a simple double well form, i.e., f(g) = (1 g)2 g2, W is the height of the hump between the two wells, and j2 is the gradient-energy coefficient. The temporal evolution of g(r) follows the TDGL equation og dF ¼ L ð2Þ ot dg where L is the kinetic coefficient characterizing interface mobility. For a planar interface, the steady-state profile of g(r) is given by [15,16]: pffiffiffiffiffi 1 x W g ¼ tanh pffiffiffi þ1 ð3Þ 2 2 2j If we assume that within the interfacial region g(r) changes from 0.1 to 0.9, the interface thickness can be estimated as: pffiffiffi j k 2:2 2 pffiffiffiffiffi ð4Þ W This solution shows that the interface thickness is propor. The interfacial energy is given by [21]: tional to pjffiffiffi W Z 1 og r ¼ j2 dx ð5Þ ox 1 Substituting Eq. (3) into Eq. (5) and carrying out the integration, we have pffiffiffiffiffi Wj r ¼ pffiffiffi ð6Þ 3 2
Eq. (6) shows pffiffiffiffiffithat the interfacial energy depends on the product of W j. To describe anisotropy in interfacial energy or interface width, it is necessary to allow both the height of the double well hump, W, and the gradientenergy coefficient, j2, to vary as functions of the inclination angle given by: gy H ¼ arctan ð7Þ gx Then the total free energy variation with respect to g yields: h i h i oH oH o W ðHÞ f ðgÞ o W ðHÞ f ðgÞ ogx ogy dF ¼ W ðHÞfg ox oy dg j2 r2 g þ 2jjH rH rg
þ jjH j2 þ jjHH gx Hy gx Hy
ð8Þ
In the above equations and the ones that follow, the subscripts denote partial derivatives, and gy oH ¼ 2 ogx gx þ g2y oH g ¼ 2 x 2 ogy gx þ gy
ð9Þ
Mobility anisotropy can also be introduced in a similar way by introducing an inclination-dependent mobility coefficient, e.g., L = L(H) [9,15,22]. The above analysis is quite general and is not confined to any particular H-dependence of the interfacial energy and interface width. In the present study, we assume a two-fold symmetry for the interfacial energy as a function of H unless mentioned separately: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rðHÞ ¼ r0 sin2 ðHÞ þ a cos2 ðHÞ ð10Þ The invariance of the interface width with respect to the inclination then requires: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi W ðHÞ ¼ W 0 sin2 ðHÞ þ a cos2 ðHÞ ð11Þ and j2 ðHÞ ¼ j20
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sin2 ðHÞ þ a cos2 ðHÞ
ð12Þ
In the above equations, a, r0, W0, j0 are constants. The evolution of g can be computed by substituting Eqs. (8), (9), (11), (12) to Eq. (2). According to Wulff construction, the equilibrium shape is p anffiffiffi ellipse with the semi-major to semi-minor axis ratio of a. To examine the growth of precipitates during, for example, a massive transformation, one can introduce a volumetric driving force into the free energy, e.g., following Kobayashi [3], Eq. (1) becomes: Z j2 2 2 3 F ¼ Wf ðgÞ þ jrgj þ mð10 15g þ 6g Þg dr ð13Þ 2 v where m is a constant related to the magnitude of the driving force for precipitate growth.
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
In our simulations, the parameters in the above equations are chosen as the following (unless mentioned separately, all are in reduced units): j20 ¼ 4, W0 = 2, a = 25, L = 1 and m = 0.5. The time step and grid size are, respectively, Dt = 0.001 and d0 = 1. 3. Equilibrium shape of a precipitate We first validate our model by comparing the simulated equilibrium shape of an isolated particle with the one determined by the Wulff construction. The simulation was carried out for a two-dimensional (2D) system of 128 128d 20 with a periodic boundary condition. Initially, a circular particle with r = 20d0 was placed at the center of the computational cell. Following Ward [23], we used a volume conserved TDGL equation to compute the equilibrium shape. When the aspect ratio of the precipitate shape became invariant with time, equilibrium was assumed. The isograms of g profiles representing the equilibrium precipitate shape are shown in Fig. 1(a) by the solid lines. For comparison the equilibrium shape obtained from the Wulff construction is also plotted in the figure given by the dashed line. It can be seen that the result from the phase
1921
field simulation agrees well with that from the Wulff construction. Parallel simulations were carried out by adopting only anisotropy in the gradient-energy coefficient and the results for two different values of j2 are shown in Fig. 1(b) and (c). It can be readily seen that the interface thickness is non-uniform and the shape defined by the contour line at g = 0.5 deviates significantly from the equilibrium one given by the dashed line. In fact, when the interface thickness is comparable to the tip radius, it is no longer valid to describe the curvature of the tip using the phase field method. In principle one can reduce the ratio of interfacial width over tip radius by either refining the numerical mesh or sharpening the interface. This would either reduce the computational efficiency or introduce numerical pinning to interfacial locations of the smallest width. Any trial-and-error approach to balance these factors will become more and more impractical as the degree of anisotropy increases. 4. Kinetics of particle growth and coalescence In this example we consider the growth of two secondphase particles from a parent phase matrix. The system size
Fig. 1. Comparison of equilibrium particle shapes obtained by phase field simulations (thin solid lines) and by Wolff construction (thick dotted lines) with (a) constant and (b) and (c) varying interface thickness (j2 = 0.5 in (c)). Considering the two-fold shape of the particle, only a portion of the computational cell is shown.
1922
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
employed is 256 256d 20 . Initially, two spherical precipitates with r = 10d0 were introduced into the system. Two parallel simulations were carried out with constant and varying interfacial thickness, respectively. Snapshots of the temporal evolution of the two particles in the two cases are shown in Fig. 2. The corresponding growth kinetics of the particles is given in Fig. 3. It can be seen that at the same moment the interfaces of the two growing particles start to overlap with each other in the case of varying boundary width (Fig. 2(b)), while they are well separated in the case of constant boundary width (Fig. 2(a)). Therefore, particle coalescence occurs at an earlier moment in the case of varying interface width. This is also reflected in Fig. 3, P where the total area of the product phase (calculated by r gðrÞ) as a function of time shows an abrupt change when the two particles coalesce. Thus a large variation of interface thickness in phase field simulations could result in significant errors in the overall growth kinetics. Since in the simulation the driving force for the growth
Fig. 3. Growth and coalescence kinetics of the two particles shown in Fig. 2.
of the particles were provided by a difference in the bulk free energy between the precipitate phase and the matrix phase, the shapes of the precipitates are quite different from the equilibrium one shown in Fig. 1. 5. Particles of sharp corners and faceted shapes The equilibrium shape of a precipitate may exhibit sharp corners and facets if the c-plot of a system is non-convex or has cusps. Missing orientations and facets have been studied recently using the phase field method [16,17], where only anisotropy in the gradient-energy coefficient is considered. In the following examples we examine the variation in boundary width at sharp corners and its effect on equilibrium shape and kinetics of shape evolution. An interfacial energy anisotropy function with four-fold symmetry was employed: rðHÞ ¼ r0 ð1 þ 0:5 cos 4HÞ
Fig. 2. Contour plots of intermediate microstructures of two growing particles at reduce time s = 120 with (a) constant and (b) varying interface thickness.
ð14Þ
As shown in Fig. 4, the polar plot of Eq. (14) has round cusps at 45 orientation. Orientations smaller than 41.7 are missing within the range of 0 < H < p/4. For convexification of the polar plot of the inverse of the interfacial energy the approach proposed by Eggleston et al. [16] was used. Two parallel simulations were carried out, one with anisotropy in the gradient-energy term only and one with anisotropy in both the gradient-energy term and the freeenergy hump following Eqs. (8) and (9). Initially, a circular particle of r = 30d0 was placed in the center of a matrix of 128d0 · 128d0. The invariant particle shapes after long time relaxation in the two cases are shown in Fig. 5(a) and (b). The enlarged contours of corners in the two cases are shown in Fig. 5(c) and (d). While smoothly round corners are produced in the case of uniform boundary width, kinks are observed at the locations where the boundary width changes in the case of varying boundary width. These kinks pin the interface and prevent the particle evolving further after initial relaxation, leading to a shape significantly different from the one assumed by the Wulff construction (the dotted line). In the case of uniform interface width the
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
Fig. 4. Polar plot of normalized interfacial energy, r(H)/r0, as a function of interface inclination.
1923
equilibrium particle shape obtained is close to a square, i.e., the equilibrium shape following the Wulff construction, but with rounded corners. With the presence of the kinks, the plane normal is erroneously calculated by the centered finite differencing formula employed in the numerical calculations. As a consequence, the interfacial energy cannot be assigned correctly. This problem may be avoided either by aligning the inclination of the interfaces with the numerical grid orientation or by applying a one-sided finite difference scheme [16]. The former is employed in this study. The equilibrium particle shapes obtained for the two cases are given in Fig. 6(a) and (b). In both cases, the equilibrium shape obtained from the phase field simulation is close to the one given by the Wulff construction. It is obvious that simulations must be carried out with great caution in the case of varying boundary width because the simulation results become sensitive to the relative orientation of the interface inclination with respect to the numerical grids.
Fig. 5. Particle shape under (a) varying and (b) constant boundary width (j2 = 8 in both cases); (c) and (d) are the corresponding enlarged contours of corners. The dotted lines show the equilibrium shape from Wulff construction. The arrows indicate the orientation of the contour normal calculated using centered finite difference. In both cases the starting particle has a circular shape.
1924
N. Ma et al. / Scripta Materialia 54 (2006) 1919–1924
be a straightforward way to alter the free energy hump [24,25]. In this case, the approach presented in Refs. [26– 28] that treats the interface as a mixture of the two equilibrium phases could be more convenient in implementing high interfacial energy anisotropy. A quantitative phase field study on side-plate development in Ti–6Al–4V using such an approach will be presented in a separate paper. Acknowledgement This work is supported by the US Air Force through the Metals Affordability Initiative under contract F33615-99-25215. References
Fig. 6. Particle shape under (a) varying and (b) constant boundary width (j2 = 8 in both cases) obtained when boundary inclinations are aligned with the numerical grid orientations.
6. Summary We have demonstrated through several simple examples that when simulating microstructural evolution in systems of high interfacial energy anisotropy using the phase field method, it is necessary to adjust the double well hump of the local free energy and the gradient-energy coefficient simultaneously to maintain a reasonable interface width. Otherwise the artificially large variation in boundary width will alter both the morphology of the microstructure and the kinetics of its evolution. For congruent transformations or grain growth where the local free energy is a function of non-conserved fields only, the method discussed seems to be able to deal effectively with high interfacial energy anisotropy. If the local free energy depends on both conserved and non-conserved fields, however, there may not
[1] Gottstein G, Shvindlerman LS. Grain boundary migration in metal. New York (NY): CRC Press; 1999. [2] Sutton AP, Balluffi RW. Interfaces in crystalline material. New York (NY): Oxford University Press; 1995. [3] Kobayashi R. Exp Math 1994;3:59. [4] Warren JA, Boettinger WJ. Acta Metall Mater 1995;43:689–703. [5] Wheeler AA, McFadden GB, Boettinger WJ. Proc Roy Soc London Ser A 1996;452:495–525. [6] Karma A, Rappel WJ. Phys Rev E 1996;53:R3017. [7] Khachaturyan AG. Philos Mag A 1996;74:3. [8] Braun RJ, Cahn JW, McFadden GB, Wheeler AA. Philos Trans R Soc London 1997;355:1787; Braun RJ, Cahn JW, McFadden GB, Wheeler AA. Acta Metall 1995; 43:689. [9] Kazaryan A, Wang Y, Dregia SA, Patton BR. Phys Rev B 2000;61: 14275. [10] Kazaryan A, Wang Y, Dregia SA, Patton BR. Phys Rev B 2001;63: 184102. [11] Kazaryan A, Patton BR, Dregia SA, Wang Y. Acta Mater 2002;50: 499. [12] Kazaryan A, Wang Y, Dregia SA, Patton BR. Acta Mater 2002;50: 2491. [13] Ma N, Kazaryan A, Dregia SA, Wang YZ. Acta Mater 2004;52:3869. [14] Caginalp G, Fife P. Phys Rev B 1986;34:4940. [15] McFadden GB, Wheeler AA, Braun RJ, Coriell SR, Sekerka RF. Phys Rev E 1993;48:2016. [16] Eggleston JJ, McFadden GB, Voorhees PW. Physica D 2001;150: 91–103. [17] Debierre JM, Karma A, Celestini F, Guerin R. Phys Rev E 2003;68: 041604. [18] Aaronson HI, Laird C. Ford Motor Co. Scientific Laboratory Report, 1967. [19] Vaithyanathan V, Wolverton C, Chen LQ. Acta Mater 2004;52: 2973–87. [20] Loginova I, Agren J, Amberg G. Acta Mater 2004;52:4055–63. [21] Allen SM, Cahn JW. Acta Metall 1979;27:1085–95. [22] Wang SL, Sekerka RF. Phys Rev B 1996;53:3760. [23] Ward MJ. SIAM J Appl Math 1996;56:1247–79. [24] Shen C, Chen Q, Wen YH, Simmons JP, Wang Y. Scripta Mater 2004;50:1023. [25] Shen C, Chen Q, Wen YH, Simmons JP, Wang Y. Scripta Mater 2004;50:1029. [26] Grafe U, Bottger B, Tiaden J, Fries SG. Scripta Mater 2000;42:1179. [27] Tiaden J, Nestler B, Diepers HJ, Steinbach I. Physica D 1998;115:73. [28] Kim SG, Kim WT, Suzuki T. Phys Rev D 1999;60:7186.