- Email: [email protected]
Phase field modeling of polycrystalline freezing
Materials Science and Engineering A 413–414 (2005) 412–417
Phase field modeling of polycrystalline freezing T. Pusztai, G. Bortel, L. Gr´an´asy ∗ Research Institute for Solid State Physics and Optics, H-1525 Budapest, POB 49, Hungary Received 5 July 2005
Abstract The formation of two and three-dimensional polycrystalline structures are addressed within the framework of the phase field theory. While in two dimensions a single orientation angle suffices to describe crystallographic orientation in the laboratory frame, in three dimensions, we use the four symmetric Euler parameters to define crystallographic orientation. Illustrative simulations are performed for various polycrystalline structures including simultaneous growth of randomly oriented dendritic particles, the formation of spherulites and crystal sheaves. © 2005 Elsevier B.V. All rights reserved. Keywords: Solidification; Polycrystalline matter; Phase field theory; Dendrites; Spherulites
1. Introduction Many of the crystalline materials in use are polycrystalline, i.e., they are composed of a large number of small crystallites, whose size, shape, and composition distributions determine their physical properties and failure characteristics [1]. These materials normally form by freezing an undercooled liquid, and their polycrystalline features are determined by the conditions of solidification and the subsequent heat treatments. The emerging polycrystalline structures can formally be divided into three main categories: (i) Grain structures formed by impinging single crystals. (ii) Polycrystalline growth forms that grow via the formation of new crystal grains at the perimeter of the structure. (iii) Impinging polycrystalline growth forms. In such processes nucleation of new grains plays an essential role. Nucleation of grains may happen homogeneously (via internal fluctuations of the liquid) or heterogeneously (with the assistance of foreign particles, container walls, etc.). We restrict the present study to the homogeneous formation of polycrystals. Modeling of such structures is an essential challenge that computational materials science has to face. In the past decade, the phase field theory became a method of choice when mod∗
Corresponding author. Tel.: +36 1 392 2222x3155; fax: +36 1 392 2219. E-mail address: [email protected] (L. Gr´an´asy).
0921-5093/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2005.09.057
eling complex solidification patterns [1]. While various aspects of polycrystalline solidification have been addressed previously within the phase field theory [2,3], description of systems with anisotropic growth needs further attention. (For a review on phase field approach to polycrystalline solidification, see Ref. [4].) The main difficulty here is that for describing an anisotropic growth of individual particles, one needs to consider the crystallographic orientation of all grains. While in the multi-phase field models [5] handling of this problem and the mimicking of nucleation via inserting supercritical crystallites into the simulation window “by hand” are straightforward, the possibility for a Langevin-type physical modeling of crystal nucleation is less obvious. Building on previous work by Kobayashi et al. [6], we recently developed a phase field theory that successfully addressed the noise-induced nucleation of crystallites with different crystallographic orientation in two-dimensional (2D) and quasi-2D systems [7]. Our approach describes the formation of such complex polycrystalline growth patterns as dendrites disordered by foreign particles [4,8], fractal like aggregates [9], spherulites [4,9,10], and “quadrites” [4,9] observed in polymeric thin films. (Similar models have been used to describe grain boundary dynamics, including grain boundary wetting and grain rotation [6,11,12].) Besides their importance in technical alloys, polymers, and minerals, polycrystalline structures are also interesting from a biological viewpoint (semicrystalline amyloid spherulites have been associated with the Alzheimer and Kreutzfeld-Jacob diseases, and the type II diabetes [13,14]). Until recently, simulation of complex polycrystalline growth
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417
patterns has been restricted to two dimensions. While in 2D a single field is sufficient to describe the local crystallographic orientation, relative to which the anisotropies of the interfacial and dynamic properties are counted, in 3D at least three parameter fields (e.g., the Euler angles) are needed to formulate the theory. We have recently developed a 3D model of polycrystalline solidification (nucleation and growth) [15]. At essentially the same time, a similar approach has been independently developed by Kobayashi and Warren for describing grain coarsening in 3D polycrystals [16,17]. In this paper, we address polycrystalline solidification in two and three dimensions. We are going to demonstrate that phase field approaches relying on orientation fields represent a powerful tool for modeling complex polycrystalline patterns. 2. Phase field theory with orientation field(s) Our starting point is a free energy functional based equivalent of the standard binary phase field theory by Warren and Boettinger [18]. In this approach the free energy functional consists of the usual square gradient, double well, and driving force terms. ε2φ T 3 2 |∇φ| + f (φ, c, T ) , F= d r (1) 2 where the local phase state of matter (solid or liquid) is characterized by the phase field φ (a structural order parameter) and the solute concentration c, εφ is a constant, and T is the temperature. The gradient term for the phase field leads to a diffuse crystal–liquid interface, a feature observed both in experiment and computer simulations [1]. The local free energy density is assumed to have the form f (φ, c, T ) = w(φ)Tg(φ) + [1 − p(φ)] fS (c) + p(φ) fL (c), where the “double well” and “interpolation” functions have the forms g(φ) = (1/4) φ2 (1 − φ)2 and p(φ) = φ3 (10 – 15φ + 6φ2 ), respectively, while the free energy scale w(φ) = (1 − c)wA + cwB [18]. The respective free energy surface has two minima (φ = 0 and 1, corresponding to the crystalline and liquid phases), whose relative depth is the driving force for crystallization and is a function of both temperature and composition as specified by the free energy densities in the bulk solid and liquid, fS,L (c, T), taken here from the ideal solution model. Time evolution is described by the following equations of motion: ∂I δF ∂I φ˙ = −Mφ + ζ φ = Mφ ∇ − + ζφ δφ ∂∇φ ∂φ δF c˙ = ∇Mc ∇ − ζj δc vm ∂I ∂I =∇ , (2) Dc(1 − c)∇ −∇ − ζj RT ∂c ∂∇c where I is the total free energy density (including the gradient terms), vm the molar volume, D the diffusion coefficient in the liquid, and ζ j are the appropriate noise terms. The time scales for the two fields are determined by the mobility coefficients appearing in the coarse-grained equations of motion: Mφ and Mc . These
413
coarse-grained mobilities can be taken from experiments and/or evaluated from atomistic simulations [1]. For example, the mobility Mc , is directly proportional to the classic inter-diffusion coefficient for a binary mixture, the mobility Mφ dictates the rate of crystallization. Solving these equations numerically in three dimensions with an anisotropic interface free energy γ = γ,0 {1 − 3ε3D + 4ε3D [(∇φ)4x + (∇φ)4y + (∇φ)4z ]/|∇φ|4 }, one obtains single crystal growth forms as exemplified in Fig. 1. The local crystallographic orientation is considered as the relative orientation of a local coordinate system (fixed to the crystal lattice) with respect to a reference system (fixed to the laboratory frame). In 2D, a single orientation field, θ ∈ [0, 1], is used. It is a normalized orientation angle, which is related to the angle Ξ relative to one of the axes of the (Cartesian) laboratory frame as Ξ = (2/k)(θ – 1/2) ∈ [−/k, /k], where k is the symmetry index. (For four-fold symmetry, k = 4.) In 3D, the relative orientation with respect to the laboratory system is uniquely defined by a single rotation of angle η around a specific axis, and can be expressed in terms of the three Euler angles. However, this representation has disadvantages: It has divergences at the poles ϑ = 0 and , and one has to use trigonometric functions that are time consuming in numerical calculations. Therefore, we opt for the four symmetric Euler parameters, q0 = cos(η/2), q1 = c1 sin(η/2), q2 = c2 sin(η/2), and q3 = c3 sin(η/2), a representation free of such difficulties. (Here ci are the components of the unit vector c of the rotation axis.) These four parameters q = (q0 , q1 , q2
, q3 ), often referred to as quaternions, satisfy the relationship i qi2 = 1, therefore, can be viewed as
a point on the four-dimensional (4D) unit sphere [19]. (Here i stands for summation with respect to i = 0, 1, 2, and 3, a notation used throughout this paper.) The angular difference δ between two orientations represented by quaternions q1 and q2 can be expressed as cos(δ) = 1/2 [Tr(R) − 1], where the matrix of rotation R is related to the individual rotation matrices R(q1 ) and R(q2 ), that rotate the reference system into the corresponding local orientations, as R = R(q1 )·R(q2 )−1 . After lengthy but straightforward algebraic manipulations one finds that the angular difference can be expressed in terms of the differences of quaternion coordinates:
cos(δ) = 1 – 2 2 + 4 /2, where 2 = (q2 − q1 )2 = i qi2 , is the square of the Euclidian distance between the points q1 and q2 on the 4D unit sphere. Comparing this expression with the Taylor expansion of the function cos(δ), one finds that 2 is indeed an excellent approximation of δ. Relying on this approximation, we express the orientational difference as δ ≈ 2 . The free energy of small-angle grain boundaries increases approximately linearly with the misorientation of the neighboring crystals, saturating at about twice of the free energy of the solid–liquid interface. Our goal is to reproduce this behavior of the small angle grain boundaries. To penalize spatial changes in the crystal orientation, in particular the presence of grain boundaries, we introduce an orientational contribution fori to the integrand in Eq. (1), which is invariant to rotations of the whole system. In 2D, it is assumed to have the form [7] fori = HT [1 − p(φ)]|∇θ|,
(3)
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417
414
Fig. 1. 3D single-crystal growth forms predicted by the phase field theory (without orientation fields) with the properties of Ni–Cu, for isotropic (left) and anisotropic (right, ε3D = 0.03) interfacial free energies, at a supersaturation of S = 0.7 (1/8 of the object has been calculated on a 800 × 800 × 800 grid, which has been then mirrored to obtain the full particles). The φ = 0.5 surface is shown.
where the grain boundary energy scales with H. This specific choice ensures a narrow grain boundary [4,11,12], and successfully describes polycrystalline solidification and grain boundary dynamics [6–12]. In 3D, we postulate [15]
1/2 fori = 2HT [1 − p(φ)] (∇qi )2 , (4) i
a definition that boils down to the 2D model, provided that the orientational transition across the grain boundary has a fixed rotation axis as in 2D. To model crystal nucleation in the liquid, the orientation fields θ(r) and q(r) are extended to the liquid, where they are made to fluctuate in time and space [7,15]. Assigning local crystal orientation to liquid regions, even a fluctuating one, may seem artificial at first sight. However, due to geometrical and/or chemical constraints, a short-range order exists even in simple liquids, which is often similar to the one in the solid. Rotating the crystalline first-neighbor shell so that it aligns optimally with the local liquid structure, one may assign a local orientation to every atom in the liquid. The orientation obtained in this manner fluctuates in time and space. The correlation of the atomic positions/angles shows how good this fit is. (In the model, the fluctuating orientation fields and the phase field play these roles.) Approaching the solid from the liquid, the orientation becomes more definite (the amplitude of the orientational fluctuations decreases) and matches to that of the solid, while the correlation between the local liquid structure and the crystal structure improves. fori recovers this behavior by realizing a strong coupling between the orientation and phase fields. Note that fori consists of the factor (1 – p(φ)) to avoid double counting the orientational contribution in the liquid, which is per definitionem incorporated into the free energy of the bulk liquid. With appropriate choice of the model parameters, an ordered liquid layer surrounds the crystal as seen in atomistic simulations [1]. Time evolution of the orientation fields is governed by relaxational dynamics both in 2D and 3D, and noise terms are added to model the thermal fluctuations. In 2D, the equation of motion
for the (non-conserved) orientation field reads as [7] ∂I ˙θ = −Mθ δF + ζθ = Mθ ∇ ∂I − + ζθ , δθ ∂∇θ ∂θ
(5)
where the orientational mobility Mθ controls the rate at which regions reorient.
In 3D, we need to take into account the quaternion properties ( i qi2 = 1) in the derivation of the equation of motion for the four orientational fields qi (r), which is done using the method of Lagrange multipliers, ∂I ∂I ∂qi δF − + ζi , = −Mq + ζ i = Mq ∇ (6) ∂t δqi ∂∇qi ∂qi where Mq is the common mobility coefficient for the symmetric quaternion fields, I is here the total free energy density (including the gradient terms and the term with the Lagrange multiplier), while ζ i are the appropriate noise terms [15]. Utilizing the relationship i qi (∂qi /∂t) = 0 that follows from the quaternion properties, and expressing the Lagrange multiplier in terms of qi and qi , the equation of motion for the orientation (quaternion) fields can be expressed as [15] ∂qi ∇qi = Mq ∇ HT [1 − p(φ)] 1/2
∂t 2 i (∇qi ) ∇qk + ζi . −qi qk ∇ HT [1 − p(φ)] 1/2
k 2 (∇q ) i i (7) Gaussian white noises of amplitude ζ i = ζ S,i + (ζ L,i − ζ S,i ) p(φ) have been added to the orientation fields so that the quaternion properties of the qi fields are retained. (ζ L,i and ζ S,i are the amplitudes in the liquid and solid.) This formulation of the model is valid for triclinic lattice without symmetries (space group P1). In the case of other crystals, the crystal symmetries yield equivalent orientations that do
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417
not form grain boundaries. These symmetries can be taken into account, when discretizing the differential operators used in the equations of motions for the quaternion fields. Calculating the angular difference between a central cell and its neighbors, all equivalent orientations of the neighbor have to be considered, the respective angular differences δ be calculated (using matrices of rotation R = R·S·R−1 , where S is a symmetry operator), of which the smallest δ value shall be used in calculating the differential operator. 3. Numerical solution The equations of motion have been solved numerically using an explicit finite difference scheme. Periodic boundary conditions were applied. The time and spatial steps were chosen to ensure stability of our solutions. The noise has been discretized as described by Karma and Rappel [20] and Plapp [21]. As the size of the 3D simulation box and the accessible time window are seriously restricted by the computational power needed, to observe nucleation, we enhanced the amplitude of the noise relative to the value from the fluctuation–dissipation theorem. A parallel code has been developed and run on two PC clusters, consisting of 60 and 100 nodes, respectively. 4. Physical properties In the present computations, we use the physical properties of the Ni–Cu system [12]. This choice is not particularly restrictive, as it is formally equivalent to a pure material, where thermal diffusion replaces solute diffusion as the dominant transport mechanism [12]. Note that our model is no way restricted to metals. We fix the temperature to be 1574 K, as in previous studies. Dendritic growth in our 3D simulations originates from the anisotropy of the phase field mobility Mφ = Mφ,0 {1 − 3ε3D + 4ε3D [(∇φ)4x + (∇φ)4y + (∇φ)4z ]/|∇φ|4 } that reflects the orientation dependence of the molecular attachment kinetics [1]. In the simulations for needle crystals, a rotational ellipsoidal symmetry of the phase field mobility has been assumed. In the 2D simulations Mφ = Mφ,0 {1 + ε2D cos[k(ϑ − 2πθ/k)]} Our calculations were performed at various supersaturations S = (cL − c)/(cL − cS ), where cL = 0.466219, cS = 0.399112, and c are the concentrations at the liquidus, solidus, and the initial homogeneous liquid mixture, respectively. Since the physical thickness of the interface is in the nanometer range and the typical solidification structures are far larger (m to mm), a full simulation of polycrystalline solidification from nucleation to particle impingement cannot be performed even with the fastest of the present supercomputers. Since we seek here a qualitative understanding, following previous work [6–10,12,15], the interface thickness has been increased (by a factor of about 20, d = 20.6 nm), while the interface free energy of the pure constituents has been reduced (by a factor of 6). This allows us to follow the life of crystallites from birth to impingement on each other. The time and spatial steps were t = 1.31 ns and x = 13.1 nm. The diffusion coefficient in the liquid was DL = 10−9 m2 /s. Unless stated otherwise, dimensionless mobil-
415
ities of Mφ,0L = 3.55 × 10−1 m3 /Js and Mq,L = 8.17 m3 /Js, and Mq,S = 0 were applied, while DS = 0 was taken in the solid. The time scale of orientational ordering is determined by the mobilities Mθ and Mq that can be related to the rotational diffusion coefficient Drot . In highly undercooled liquids, dynamic heterogeneities exist at the nanometer scale, leading to a reduction of the ratio Drot /Dtr of the rotational and translational diffusion coefficients. Recent experiment has shown that the rate of crystallization in highly supercooled liquids is proportional to Dtr [22]. In our model, the growth velocity scales linearly with Mφ , so consistency requires Mφ ∝ Dtr . Since we also expect that Mθ ∝ Drot , we arrive at Mθ /Mφ ∝ Drot /Dtr . Reduction of Drot /Dtr enables the formation of new grains at the perimeter as detailed in Ref. [10]. 5. Results and discussion 5.1. Polycrystalline solidification in 2D We performed illustrative simulations for anisotropic growth of polycrystalline patterns crystals in 2D. Noise-induced nucleation and growth in a binary alloy leads to poly-dendritic structures shown in Fig. 2(a), demonstrating that the formation of category (i) polycrystals (impinging single crystals) can be handled within the present theoretical framework. Reducing the orientational mobility makes it difficult for newly forming crystal regions to reorient with the parent crystal to lower its free energy at the growth front that is advancing with a velocity scaling with the translational diffusion coefficient. Thus epitaxy cannot keep pace with solidification, i.e., the orientational order that freezes in is incomplete, leading to a continuous formation of new grains at the perimeter as discussed recently in detail by us [9,10]. Under such conditions category (ii) polycrystalline structures form, illustrated by the polycrystalline spherulite displayed in Fig. 2(b). In many systems the formation of crystal fibers is preferred (polymers, biopolymers, etc.), that can be mimicked in our 2D approach by assuming a large two-fold (k = 2) anisotropy of the kinetic coefficient. Starting the solidification of such a system from a single crystal under reduced orientational mobility, first fibrils form and then secondary fibrils nucleate at the growth front to form crystal ‘sheaves’. The diverging ends of these sheaves subsequently fan out with time to form “eyes” (uncrystallized holes on the two sides of the initial fiber), and finally a roughly spherical growth form emerges. Fig. 2(c) shows the impingement of such polycrystalline growth forms that yields a category (iii) structure. This progression of crystallization is nearly universal in polymeric materials [22]. 5.2. Polycrystalline solidification in 3D Dendritic solidification of four particles with random orientation grown with cubic crystal symmetries is shown in Fig. 3(a) and (b). These morphologies are consistent with previous results for single-crystal dendritic growth [1,23].
416
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417
Fig. 2. Polycrystalline freezing in 2D. (a) Impinging dendritic single crystals. (Part of a simulation on a 7000 × 7000 grid, S = 0.8, and with an interfacial free energy anisotropy of ε2D = 0.15, k = 4. Such large anisotropies are not unusual in polymeric systems.) Composition map is shown (cL : blue, cL : yellow). (b) Polycrystalline spherulite formed by trapping of orientational disorder (4000 × 4000 grid, S = 0.95, and with a kinetic anisotropy of ε2D = 0.5, k = 4). Orientation field is shown. Different colors correspond to different crystallographic orientations. The liquid is denoted by black. (c) Impinging crystal sheaves formed by branching of needle crystals (2000 × 2000 grid, S = 0.9, interfacial free energy anisotropy of ε2D = 0.5, k = 2). Orientation field is shown. Coloring is as for panel (b).
Fig. 3. Polycrystalline freezing in 3D. (a and b) Snapshots showing the growth of four randomly oriented dendrites assuming cubic crystal symmetry (400 × 400 × 400 grid, S = 0.75, and with a kinetic anisotropy of ε3D = 0.3). Note the effect of periodic boundary conditions: the branches that grow out of the simulation on one side of the simulation box, enter on the opposite side. (c) Polycrystalline spherulite formed by trapping of orientational disorder calculated triclinic crystal symmetry (300 × 300 × 300 grid, S = 0.85, and with a kinetic anisotropy of ε3D = 0.3). The growth front is shaded according to the angular difference between the z-axis of the local and laboratory frames (see color bar). (d) A crystal sheaf formed by branching of a needle crystal simulated with triclinic crystal symmetry (250 × 250 × 500 grid, S = 0.85, and with an ellipsoidal symmetry of the phase field mobility). The φ = 0.5 surface is shown.
As in 2D, reducing the orientational mobility, a uniform crystallographic orientation cannot be established at the growth front, and new crystal grains form induced by orientational defects frozen into the solid [4,9,10]. At large driving forces space filling polycrystalline forms called spherulites appear [4,9,10]. A 3D polycrystalline spherulite obtained with triclinic crystal structure is displayed in Fig. 3(c). In the case of a large ellipsoidal anisotropy of the phase field mobility, needle crystals grow at small supersaturations. With increasing supersaturation random branching occurs, and more space filling patterns such as crystal sheaves [Fig. 3(d)] and polycrystalline spherulites form. 6. Summary We demonstrated that the phase field models with appropriate orientation fields, we proposed recently, are capable for describing various polycrystalline solidification processes including impinging dendritic particles, branching needle crystals, and spherulitic polycrystals both in 2D and 3D. These methods open up the way for phase field modeling of complex polycrystalline morphologies in three dimensions.
Acknowledgments The authors thank M. Tegze for the enlightening discussions. This work has been supported by contracts OTKA-T-037323 and ESA PECS No. 98005, and by the EU FP6 Integrated Project IMPRESS under Contract No. NMP3-CT-2004-500635, and forms part of the ESA MAP Projects No. AO-99-101 and AO-99-114. T.P. and G.B. acknowledge support by the Bolyai J´anos Scholarship of the Hungarian Academy of Sciences.
References [1] J.J. Hoyt, M. Asta, A. Karma, Mater. Sci. Eng. Rep. R 41 (2003) 121. [2] K.R. Elder, F. Drolet, J.M. Kosterlitz, M. Grant, Phys. Rev. Lett. 72 (1994) 677. [3] B. Morin, K.R. Elder, M. Sutton, M. Grant, Phys. Rev. Lett. 75 (1995) 2156. [4] L. Gr´an´asy, T. Pusztai, J.A. Warren, J. Phys. Condens. Matter 16 (2004) R1205. [5] I. Steinbach, F. Pezzola, B. Nestler, M. Seesselberg, R. Prieler, G.J. Schmitz, J.L.L. Rezende, Physica D 94 (1996) 135. [6] R. Kobayashi, J.A. Warren, W.C. Carter, Physica D 119 (1998) 415. [7] L. Gr´an´asy, T. B¨orzs¨onyi, T. Pusztai, Phys. Rev. Lett. 88 (2002) 206105.
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417 [8] L. Gr´an´asy, T. Pusztai, J.A. Warren, T. B¨orzs¨onyi, J.F. Douglas, V. Ferreiro, Nat. Mater. 2 (2003) 92. [9] L. Gr´an´asy, T. Pusztai, G. Tegze, J.A. Warren, J.F. Douglas, Phys. Rev. E 72 (art. no. 011605) (2005). [10] L. Gr´an´asy, T. Pusztai, T. B¨orzs¨onyi, J.A. Warren, J.F. Douglas, Nat. Mater. 3 (2004) 635. [11] A.E. Lobkovsky, J.A. Warren, Phys. Rev. E 63 (2001) 051605. [12] J.A. Warren, R. Kobayashi, A.E. Lobkovsky, W.C. Carter, Acta Mater. 51 (2003) 6035. [13] L.-W. Jin, K.A. Claborn, M. Kurimoto, M.A. Geday, I. Maezawa, F. Sohraby, M. Estrada, W. Kaminsky, B. Kahr, Proc. Natl. Acad. Sci. 100 (2003) 15297.
417
[14] M.R.H. Krebs, E.H.C. Krebs, S.S. Roger, A.M. Donald, Biophys. J. 88 (2005) 2013. [15] T. Pusztai, G. Bortel, L. Gr´an´asy, Europhys. Lett. 71 (2005) 131. [16] R. Kobayashi, J.A. Warren, TMS Lett. 2 (2005) 1. [17] R. Kobayashi, J.A. Warren, Physica A 356 (2005) 127. [18] J.A. Warren, W.J. Boettinger, Acta Metall. Mater. 43 (1995) 689. [19] G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers, McGraw-Hill, 1970, Chapter 14.10. [20] A. Karma, W.-J. Rappel, Phys. Rev. E 60 (1999) 3614. [21] M. Plapp, personal communication. [22] J.H. Magill, J. Mater. Sci. 36 (2001) 3143. [23] J. Bragard, A. Karma, Y.H. Lee, M. Plapp, Interf. Sci. 10 (2002) 121.
Phase field modeling of polycrystalline freezing T. Pusztai, G. Bortel, L. Gr´an´asy ∗ Research Institute for Solid State Physics and Optics, H-1525 Budapest, POB 49, Hungary Received 5 July 2005
Abstract The formation of two and three-dimensional polycrystalline structures are addressed within the framework of the phase field theory. While in two dimensions a single orientation angle suffices to describe crystallographic orientation in the laboratory frame, in three dimensions, we use the four symmetric Euler parameters to define crystallographic orientation. Illustrative simulations are performed for various polycrystalline structures including simultaneous growth of randomly oriented dendritic particles, the formation of spherulites and crystal sheaves. © 2005 Elsevier B.V. All rights reserved. Keywords: Solidification; Polycrystalline matter; Phase field theory; Dendrites; Spherulites
1. Introduction Many of the crystalline materials in use are polycrystalline, i.e., they are composed of a large number of small crystallites, whose size, shape, and composition distributions determine their physical properties and failure characteristics [1]. These materials normally form by freezing an undercooled liquid, and their polycrystalline features are determined by the conditions of solidification and the subsequent heat treatments. The emerging polycrystalline structures can formally be divided into three main categories: (i) Grain structures formed by impinging single crystals. (ii) Polycrystalline growth forms that grow via the formation of new crystal grains at the perimeter of the structure. (iii) Impinging polycrystalline growth forms. In such processes nucleation of new grains plays an essential role. Nucleation of grains may happen homogeneously (via internal fluctuations of the liquid) or heterogeneously (with the assistance of foreign particles, container walls, etc.). We restrict the present study to the homogeneous formation of polycrystals. Modeling of such structures is an essential challenge that computational materials science has to face. In the past decade, the phase field theory became a method of choice when mod∗
Corresponding author. Tel.: +36 1 392 2222x3155; fax: +36 1 392 2219. E-mail address: [email protected] (L. Gr´an´asy).
0921-5093/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2005.09.057
eling complex solidification patterns [1]. While various aspects of polycrystalline solidification have been addressed previously within the phase field theory [2,3], description of systems with anisotropic growth needs further attention. (For a review on phase field approach to polycrystalline solidification, see Ref. [4].) The main difficulty here is that for describing an anisotropic growth of individual particles, one needs to consider the crystallographic orientation of all grains. While in the multi-phase field models [5] handling of this problem and the mimicking of nucleation via inserting supercritical crystallites into the simulation window “by hand” are straightforward, the possibility for a Langevin-type physical modeling of crystal nucleation is less obvious. Building on previous work by Kobayashi et al. [6], we recently developed a phase field theory that successfully addressed the noise-induced nucleation of crystallites with different crystallographic orientation in two-dimensional (2D) and quasi-2D systems [7]. Our approach describes the formation of such complex polycrystalline growth patterns as dendrites disordered by foreign particles [4,8], fractal like aggregates [9], spherulites [4,9,10], and “quadrites” [4,9] observed in polymeric thin films. (Similar models have been used to describe grain boundary dynamics, including grain boundary wetting and grain rotation [6,11,12].) Besides their importance in technical alloys, polymers, and minerals, polycrystalline structures are also interesting from a biological viewpoint (semicrystalline amyloid spherulites have been associated with the Alzheimer and Kreutzfeld-Jacob diseases, and the type II diabetes [13,14]). Until recently, simulation of complex polycrystalline growth
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417
patterns has been restricted to two dimensions. While in 2D a single field is sufficient to describe the local crystallographic orientation, relative to which the anisotropies of the interfacial and dynamic properties are counted, in 3D at least three parameter fields (e.g., the Euler angles) are needed to formulate the theory. We have recently developed a 3D model of polycrystalline solidification (nucleation and growth) [15]. At essentially the same time, a similar approach has been independently developed by Kobayashi and Warren for describing grain coarsening in 3D polycrystals [16,17]. In this paper, we address polycrystalline solidification in two and three dimensions. We are going to demonstrate that phase field approaches relying on orientation fields represent a powerful tool for modeling complex polycrystalline patterns. 2. Phase field theory with orientation field(s) Our starting point is a free energy functional based equivalent of the standard binary phase field theory by Warren and Boettinger [18]. In this approach the free energy functional consists of the usual square gradient, double well, and driving force terms. ε2φ T 3 2 |∇φ| + f (φ, c, T ) , F= d r (1) 2 where the local phase state of matter (solid or liquid) is characterized by the phase field φ (a structural order parameter) and the solute concentration c, εφ is a constant, and T is the temperature. The gradient term for the phase field leads to a diffuse crystal–liquid interface, a feature observed both in experiment and computer simulations [1]. The local free energy density is assumed to have the form f (φ, c, T ) = w(φ)Tg(φ) + [1 − p(φ)] fS (c) + p(φ) fL (c), where the “double well” and “interpolation” functions have the forms g(φ) = (1/4) φ2 (1 − φ)2 and p(φ) = φ3 (10 – 15φ + 6φ2 ), respectively, while the free energy scale w(φ) = (1 − c)wA + cwB [18]. The respective free energy surface has two minima (φ = 0 and 1, corresponding to the crystalline and liquid phases), whose relative depth is the driving force for crystallization and is a function of both temperature and composition as specified by the free energy densities in the bulk solid and liquid, fS,L (c, T), taken here from the ideal solution model. Time evolution is described by the following equations of motion: ∂I δF ∂I φ˙ = −Mφ + ζ φ = Mφ ∇ − + ζφ δφ ∂∇φ ∂φ δF c˙ = ∇Mc ∇ − ζj δc vm ∂I ∂I =∇ , (2) Dc(1 − c)∇ −∇ − ζj RT ∂c ∂∇c where I is the total free energy density (including the gradient terms), vm the molar volume, D the diffusion coefficient in the liquid, and ζ j are the appropriate noise terms. The time scales for the two fields are determined by the mobility coefficients appearing in the coarse-grained equations of motion: Mφ and Mc . These
413
coarse-grained mobilities can be taken from experiments and/or evaluated from atomistic simulations [1]. For example, the mobility Mc , is directly proportional to the classic inter-diffusion coefficient for a binary mixture, the mobility Mφ dictates the rate of crystallization. Solving these equations numerically in three dimensions with an anisotropic interface free energy γ = γ,0 {1 − 3ε3D + 4ε3D [(∇φ)4x + (∇φ)4y + (∇φ)4z ]/|∇φ|4 }, one obtains single crystal growth forms as exemplified in Fig. 1. The local crystallographic orientation is considered as the relative orientation of a local coordinate system (fixed to the crystal lattice) with respect to a reference system (fixed to the laboratory frame). In 2D, a single orientation field, θ ∈ [0, 1], is used. It is a normalized orientation angle, which is related to the angle Ξ relative to one of the axes of the (Cartesian) laboratory frame as Ξ = (2/k)(θ – 1/2) ∈ [−/k, /k], where k is the symmetry index. (For four-fold symmetry, k = 4.) In 3D, the relative orientation with respect to the laboratory system is uniquely defined by a single rotation of angle η around a specific axis, and can be expressed in terms of the three Euler angles. However, this representation has disadvantages: It has divergences at the poles ϑ = 0 and , and one has to use trigonometric functions that are time consuming in numerical calculations. Therefore, we opt for the four symmetric Euler parameters, q0 = cos(η/2), q1 = c1 sin(η/2), q2 = c2 sin(η/2), and q3 = c3 sin(η/2), a representation free of such difficulties. (Here ci are the components of the unit vector c of the rotation axis.) These four parameters q = (q0 , q1 , q2
, q3 ), often referred to as quaternions, satisfy the relationship i qi2 = 1, therefore, can be viewed as
a point on the four-dimensional (4D) unit sphere [19]. (Here i stands for summation with respect to i = 0, 1, 2, and 3, a notation used throughout this paper.) The angular difference δ between two orientations represented by quaternions q1 and q2 can be expressed as cos(δ) = 1/2 [Tr(R) − 1], where the matrix of rotation R is related to the individual rotation matrices R(q1 ) and R(q2 ), that rotate the reference system into the corresponding local orientations, as R = R(q1 )·R(q2 )−1 . After lengthy but straightforward algebraic manipulations one finds that the angular difference can be expressed in terms of the differences of quaternion coordinates:
cos(δ) = 1 – 2 2 + 4 /2, where 2 = (q2 − q1 )2 = i qi2 , is the square of the Euclidian distance between the points q1 and q2 on the 4D unit sphere. Comparing this expression with the Taylor expansion of the function cos(δ), one finds that 2 is indeed an excellent approximation of δ. Relying on this approximation, we express the orientational difference as δ ≈ 2 . The free energy of small-angle grain boundaries increases approximately linearly with the misorientation of the neighboring crystals, saturating at about twice of the free energy of the solid–liquid interface. Our goal is to reproduce this behavior of the small angle grain boundaries. To penalize spatial changes in the crystal orientation, in particular the presence of grain boundaries, we introduce an orientational contribution fori to the integrand in Eq. (1), which is invariant to rotations of the whole system. In 2D, it is assumed to have the form [7] fori = HT [1 − p(φ)]|∇θ|,
(3)
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417
414
Fig. 1. 3D single-crystal growth forms predicted by the phase field theory (without orientation fields) with the properties of Ni–Cu, for isotropic (left) and anisotropic (right, ε3D = 0.03) interfacial free energies, at a supersaturation of S = 0.7 (1/8 of the object has been calculated on a 800 × 800 × 800 grid, which has been then mirrored to obtain the full particles). The φ = 0.5 surface is shown.
where the grain boundary energy scales with H. This specific choice ensures a narrow grain boundary [4,11,12], and successfully describes polycrystalline solidification and grain boundary dynamics [6–12]. In 3D, we postulate [15]
1/2 fori = 2HT [1 − p(φ)] (∇qi )2 , (4) i
a definition that boils down to the 2D model, provided that the orientational transition across the grain boundary has a fixed rotation axis as in 2D. To model crystal nucleation in the liquid, the orientation fields θ(r) and q(r) are extended to the liquid, where they are made to fluctuate in time and space [7,15]. Assigning local crystal orientation to liquid regions, even a fluctuating one, may seem artificial at first sight. However, due to geometrical and/or chemical constraints, a short-range order exists even in simple liquids, which is often similar to the one in the solid. Rotating the crystalline first-neighbor shell so that it aligns optimally with the local liquid structure, one may assign a local orientation to every atom in the liquid. The orientation obtained in this manner fluctuates in time and space. The correlation of the atomic positions/angles shows how good this fit is. (In the model, the fluctuating orientation fields and the phase field play these roles.) Approaching the solid from the liquid, the orientation becomes more definite (the amplitude of the orientational fluctuations decreases) and matches to that of the solid, while the correlation between the local liquid structure and the crystal structure improves. fori recovers this behavior by realizing a strong coupling between the orientation and phase fields. Note that fori consists of the factor (1 – p(φ)) to avoid double counting the orientational contribution in the liquid, which is per definitionem incorporated into the free energy of the bulk liquid. With appropriate choice of the model parameters, an ordered liquid layer surrounds the crystal as seen in atomistic simulations [1]. Time evolution of the orientation fields is governed by relaxational dynamics both in 2D and 3D, and noise terms are added to model the thermal fluctuations. In 2D, the equation of motion
for the (non-conserved) orientation field reads as [7] ∂I ˙θ = −Mθ δF + ζθ = Mθ ∇ ∂I − + ζθ , δθ ∂∇θ ∂θ
(5)
where the orientational mobility Mθ controls the rate at which regions reorient.
In 3D, we need to take into account the quaternion properties ( i qi2 = 1) in the derivation of the equation of motion for the four orientational fields qi (r), which is done using the method of Lagrange multipliers, ∂I ∂I ∂qi δF − + ζi , = −Mq + ζ i = Mq ∇ (6) ∂t δqi ∂∇qi ∂qi where Mq is the common mobility coefficient for the symmetric quaternion fields, I is here the total free energy density (including the gradient terms and the term with the Lagrange multiplier), while ζ i are the appropriate noise terms [15]. Utilizing the relationship i qi (∂qi /∂t) = 0 that follows from the quaternion properties, and expressing the Lagrange multiplier in terms of qi and qi , the equation of motion for the orientation (quaternion) fields can be expressed as [15] ∂qi ∇qi = Mq ∇ HT [1 − p(φ)] 1/2
∂t 2 i (∇qi ) ∇qk + ζi . −qi qk ∇ HT [1 − p(φ)] 1/2
k 2 (∇q ) i i (7) Gaussian white noises of amplitude ζ i = ζ S,i + (ζ L,i − ζ S,i ) p(φ) have been added to the orientation fields so that the quaternion properties of the qi fields are retained. (ζ L,i and ζ S,i are the amplitudes in the liquid and solid.) This formulation of the model is valid for triclinic lattice without symmetries (space group P1). In the case of other crystals, the crystal symmetries yield equivalent orientations that do
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417
not form grain boundaries. These symmetries can be taken into account, when discretizing the differential operators used in the equations of motions for the quaternion fields. Calculating the angular difference between a central cell and its neighbors, all equivalent orientations of the neighbor have to be considered, the respective angular differences δ be calculated (using matrices of rotation R = R·S·R−1 , where S is a symmetry operator), of which the smallest δ value shall be used in calculating the differential operator. 3. Numerical solution The equations of motion have been solved numerically using an explicit finite difference scheme. Periodic boundary conditions were applied. The time and spatial steps were chosen to ensure stability of our solutions. The noise has been discretized as described by Karma and Rappel [20] and Plapp [21]. As the size of the 3D simulation box and the accessible time window are seriously restricted by the computational power needed, to observe nucleation, we enhanced the amplitude of the noise relative to the value from the fluctuation–dissipation theorem. A parallel code has been developed and run on two PC clusters, consisting of 60 and 100 nodes, respectively. 4. Physical properties In the present computations, we use the physical properties of the Ni–Cu system [12]. This choice is not particularly restrictive, as it is formally equivalent to a pure material, where thermal diffusion replaces solute diffusion as the dominant transport mechanism [12]. Note that our model is no way restricted to metals. We fix the temperature to be 1574 K, as in previous studies. Dendritic growth in our 3D simulations originates from the anisotropy of the phase field mobility Mφ = Mφ,0 {1 − 3ε3D + 4ε3D [(∇φ)4x + (∇φ)4y + (∇φ)4z ]/|∇φ|4 } that reflects the orientation dependence of the molecular attachment kinetics [1]. In the simulations for needle crystals, a rotational ellipsoidal symmetry of the phase field mobility has been assumed. In the 2D simulations Mφ = Mφ,0 {1 + ε2D cos[k(ϑ − 2πθ/k)]} Our calculations were performed at various supersaturations S = (cL − c)/(cL − cS ), where cL = 0.466219, cS = 0.399112, and c are the concentrations at the liquidus, solidus, and the initial homogeneous liquid mixture, respectively. Since the physical thickness of the interface is in the nanometer range and the typical solidification structures are far larger (m to mm), a full simulation of polycrystalline solidification from nucleation to particle impingement cannot be performed even with the fastest of the present supercomputers. Since we seek here a qualitative understanding, following previous work [6–10,12,15], the interface thickness has been increased (by a factor of about 20, d = 20.6 nm), while the interface free energy of the pure constituents has been reduced (by a factor of 6). This allows us to follow the life of crystallites from birth to impingement on each other. The time and spatial steps were t = 1.31 ns and x = 13.1 nm. The diffusion coefficient in the liquid was DL = 10−9 m2 /s. Unless stated otherwise, dimensionless mobil-
415
ities of Mφ,0L = 3.55 × 10−1 m3 /Js and Mq,L = 8.17 m3 /Js, and Mq,S = 0 were applied, while DS = 0 was taken in the solid. The time scale of orientational ordering is determined by the mobilities Mθ and Mq that can be related to the rotational diffusion coefficient Drot . In highly undercooled liquids, dynamic heterogeneities exist at the nanometer scale, leading to a reduction of the ratio Drot /Dtr of the rotational and translational diffusion coefficients. Recent experiment has shown that the rate of crystallization in highly supercooled liquids is proportional to Dtr [22]. In our model, the growth velocity scales linearly with Mφ , so consistency requires Mφ ∝ Dtr . Since we also expect that Mθ ∝ Drot , we arrive at Mθ /Mφ ∝ Drot /Dtr . Reduction of Drot /Dtr enables the formation of new grains at the perimeter as detailed in Ref. [10]. 5. Results and discussion 5.1. Polycrystalline solidification in 2D We performed illustrative simulations for anisotropic growth of polycrystalline patterns crystals in 2D. Noise-induced nucleation and growth in a binary alloy leads to poly-dendritic structures shown in Fig. 2(a), demonstrating that the formation of category (i) polycrystals (impinging single crystals) can be handled within the present theoretical framework. Reducing the orientational mobility makes it difficult for newly forming crystal regions to reorient with the parent crystal to lower its free energy at the growth front that is advancing with a velocity scaling with the translational diffusion coefficient. Thus epitaxy cannot keep pace with solidification, i.e., the orientational order that freezes in is incomplete, leading to a continuous formation of new grains at the perimeter as discussed recently in detail by us [9,10]. Under such conditions category (ii) polycrystalline structures form, illustrated by the polycrystalline spherulite displayed in Fig. 2(b). In many systems the formation of crystal fibers is preferred (polymers, biopolymers, etc.), that can be mimicked in our 2D approach by assuming a large two-fold (k = 2) anisotropy of the kinetic coefficient. Starting the solidification of such a system from a single crystal under reduced orientational mobility, first fibrils form and then secondary fibrils nucleate at the growth front to form crystal ‘sheaves’. The diverging ends of these sheaves subsequently fan out with time to form “eyes” (uncrystallized holes on the two sides of the initial fiber), and finally a roughly spherical growth form emerges. Fig. 2(c) shows the impingement of such polycrystalline growth forms that yields a category (iii) structure. This progression of crystallization is nearly universal in polymeric materials [22]. 5.2. Polycrystalline solidification in 3D Dendritic solidification of four particles with random orientation grown with cubic crystal symmetries is shown in Fig. 3(a) and (b). These morphologies are consistent with previous results for single-crystal dendritic growth [1,23].
416
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417
Fig. 2. Polycrystalline freezing in 2D. (a) Impinging dendritic single crystals. (Part of a simulation on a 7000 × 7000 grid, S = 0.8, and with an interfacial free energy anisotropy of ε2D = 0.15, k = 4. Such large anisotropies are not unusual in polymeric systems.) Composition map is shown (cL : blue, cL : yellow). (b) Polycrystalline spherulite formed by trapping of orientational disorder (4000 × 4000 grid, S = 0.95, and with a kinetic anisotropy of ε2D = 0.5, k = 4). Orientation field is shown. Different colors correspond to different crystallographic orientations. The liquid is denoted by black. (c) Impinging crystal sheaves formed by branching of needle crystals (2000 × 2000 grid, S = 0.9, interfacial free energy anisotropy of ε2D = 0.5, k = 2). Orientation field is shown. Coloring is as for panel (b).
Fig. 3. Polycrystalline freezing in 3D. (a and b) Snapshots showing the growth of four randomly oriented dendrites assuming cubic crystal symmetry (400 × 400 × 400 grid, S = 0.75, and with a kinetic anisotropy of ε3D = 0.3). Note the effect of periodic boundary conditions: the branches that grow out of the simulation on one side of the simulation box, enter on the opposite side. (c) Polycrystalline spherulite formed by trapping of orientational disorder calculated triclinic crystal symmetry (300 × 300 × 300 grid, S = 0.85, and with a kinetic anisotropy of ε3D = 0.3). The growth front is shaded according to the angular difference between the z-axis of the local and laboratory frames (see color bar). (d) A crystal sheaf formed by branching of a needle crystal simulated with triclinic crystal symmetry (250 × 250 × 500 grid, S = 0.85, and with an ellipsoidal symmetry of the phase field mobility). The φ = 0.5 surface is shown.
As in 2D, reducing the orientational mobility, a uniform crystallographic orientation cannot be established at the growth front, and new crystal grains form induced by orientational defects frozen into the solid [4,9,10]. At large driving forces space filling polycrystalline forms called spherulites appear [4,9,10]. A 3D polycrystalline spherulite obtained with triclinic crystal structure is displayed in Fig. 3(c). In the case of a large ellipsoidal anisotropy of the phase field mobility, needle crystals grow at small supersaturations. With increasing supersaturation random branching occurs, and more space filling patterns such as crystal sheaves [Fig. 3(d)] and polycrystalline spherulites form. 6. Summary We demonstrated that the phase field models with appropriate orientation fields, we proposed recently, are capable for describing various polycrystalline solidification processes including impinging dendritic particles, branching needle crystals, and spherulitic polycrystals both in 2D and 3D. These methods open up the way for phase field modeling of complex polycrystalline morphologies in three dimensions.
Acknowledgments The authors thank M. Tegze for the enlightening discussions. This work has been supported by contracts OTKA-T-037323 and ESA PECS No. 98005, and by the EU FP6 Integrated Project IMPRESS under Contract No. NMP3-CT-2004-500635, and forms part of the ESA MAP Projects No. AO-99-101 and AO-99-114. T.P. and G.B. acknowledge support by the Bolyai J´anos Scholarship of the Hungarian Academy of Sciences.
References [1] J.J. Hoyt, M. Asta, A. Karma, Mater. Sci. Eng. Rep. R 41 (2003) 121. [2] K.R. Elder, F. Drolet, J.M. Kosterlitz, M. Grant, Phys. Rev. Lett. 72 (1994) 677. [3] B. Morin, K.R. Elder, M. Sutton, M. Grant, Phys. Rev. Lett. 75 (1995) 2156. [4] L. Gr´an´asy, T. Pusztai, J.A. Warren, J. Phys. Condens. Matter 16 (2004) R1205. [5] I. Steinbach, F. Pezzola, B. Nestler, M. Seesselberg, R. Prieler, G.J. Schmitz, J.L.L. Rezende, Physica D 94 (1996) 135. [6] R. Kobayashi, J.A. Warren, W.C. Carter, Physica D 119 (1998) 415. [7] L. Gr´an´asy, T. B¨orzs¨onyi, T. Pusztai, Phys. Rev. Lett. 88 (2002) 206105.
T. Pusztai et al. / Materials Science and Engineering A 413–414 (2005) 412–417 [8] L. Gr´an´asy, T. Pusztai, J.A. Warren, T. B¨orzs¨onyi, J.F. Douglas, V. Ferreiro, Nat. Mater. 2 (2003) 92. [9] L. Gr´an´asy, T. Pusztai, G. Tegze, J.A. Warren, J.F. Douglas, Phys. Rev. E 72 (art. no. 011605) (2005). [10] L. Gr´an´asy, T. Pusztai, T. B¨orzs¨onyi, J.A. Warren, J.F. Douglas, Nat. Mater. 3 (2004) 635. [11] A.E. Lobkovsky, J.A. Warren, Phys. Rev. E 63 (2001) 051605. [12] J.A. Warren, R. Kobayashi, A.E. Lobkovsky, W.C. Carter, Acta Mater. 51 (2003) 6035. [13] L.-W. Jin, K.A. Claborn, M. Kurimoto, M.A. Geday, I. Maezawa, F. Sohraby, M. Estrada, W. Kaminsky, B. Kahr, Proc. Natl. Acad. Sci. 100 (2003) 15297.
417
[14] M.R.H. Krebs, E.H.C. Krebs, S.S. Roger, A.M. Donald, Biophys. J. 88 (2005) 2013. [15] T. Pusztai, G. Bortel, L. Gr´an´asy, Europhys. Lett. 71 (2005) 131. [16] R. Kobayashi, J.A. Warren, TMS Lett. 2 (2005) 1. [17] R. Kobayashi, J.A. Warren, Physica A 356 (2005) 127. [18] J.A. Warren, W.J. Boettinger, Acta Metall. Mater. 43 (1995) 689. [19] G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers, McGraw-Hill, 1970, Chapter 14.10. [20] A. Karma, W.-J. Rappel, Phys. Rev. E 60 (1999) 3614. [21] M. Plapp, personal communication. [22] J.H. Magill, J. Mater. Sci. 36 (2001) 3143. [23] J. Bragard, A. Karma, Y.H. Lee, M. Plapp, Interf. Sci. 10 (2002) 121.