Phase field modeling of polycrystalline freezing T. Pusztai, G. Bortel, L. Gránásy

(1)

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.

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 emerg- ing 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:grana@szfki.hu (L. Gr´an´asy).

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 dis- ordered 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

doi:10.1016/j.msea.2005.09.057

(2)

field approaches relying on orientation fields represent a power- ful tool for modeling complex polycrystalline patterns.

2. Phase ﬁeld theory with orientation ﬁeld(s)

Our starting point is a free energy functional based equivalent of the standard binary phase field theory by Warren and Boet- tinger[18]. In this approach the free energy functional consists of the usual square gradient, double well, and driving force terms.

F =

d³r ε²_φT

2 |∇φ|²+f(φ, c, T)

, (1)

where the local phase state of matter (solid or liquid) is characterized by the phase fieldφ(a structural order parameter) and the solute concentrationc,ε_φ 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) φ²(1−φ)² andp(φ) =φ³(10 – 15φ+ 6φ²), respectively, while the free energy scale w(φ)=(1−c)w_A+cw_B [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:

φ˙ = −Mφ

δF

δφ +ζφ=Mφ

∇ ∂I

∂∇φ

−∂I

∂φ

+ζφ

˙

c= ∇M_c∇ δF

δc −ζ_j

= ∇ v_m

RTDc(1−c)∇ ∂I

∂c

−∇

∂I

∂∇c

−ζ_j , (2) whereIis the total free energy density (including the gradient terms),vmthe molar volume,Dthe diffusion coefficient in the liquid, andζjare the appropriate noise terms. The time scales for the two fields are determined by the mobility coefficients appear- ing in the coarse-grained equations of motion:MφandMc. These

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], wherek 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 disadvan- tages: It has divergences at the polesϑ= 0 and␲, and one has to use trigonometric functions that are time consuming in numeri- cal calculations. Therefore, we opt for the four symmetric Euler parameters,q0= cos(η/2),q1=c1sin(η/2),q2=c2sin(η/2), and q3=c3sin(η/2), a representation free of such difficulties. (Here ciare the components of the unit vectorcof the rotation axis.) These four parametersq= (q0,q1,q2,q3), often referred to as quaternions, satisfy the relationship

iq²_i =1, therefore, can be viewed as a point on the four-dimensional (4D) unit sphere [19]. (Here

istands for summation with respect toi= 0, 1, 2, and 3, a notation used throughout this paper.)

The angular difference δ between two orientations repre- sented by quaternionsq1andq2can be expressed as cos(δ) = 1/2 [Tr(R)−1], where the matrix of rotation Ris 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)⁻¹. 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, where²=(q2−q1)²=

iq_i², is the square of the Euclidian distance between the pointsq1and q2on the 4D unit sphere. Comparing this expression with the Taylor expansion of the function cos(δ), one finds that 2is 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 neighbor- ing 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 contributionforito 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)

(3)

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 ofS= 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 withH. 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]

f_ori=2HT[1−p(φ)]

i(∇q_i)² 1/2

, (4)

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) andq(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 chem- ical 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 man- ner 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.forirecovers this behavior by realizing a strong cou- pling between the orientation and phase fields.

Note thatforiconsists 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 relax- ational 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]

θ˙= −M_θδF

δθ +ζ_θ=M_θ

∇ ∂I

∂∇θ

−∂I

∂θ

+ζ_θ, (5) where the orientational mobilityMθ controls the rate at which regions reorient.

In 3D, we need to take into account the quaternion properties (

iq²_i =1) in the derivation of the equation of motion for the four orientational fieldsqi(r), which is done using the method of Lagrange multipliers,

∂q_i

∂t = −Mq

δF

δq_i +ζi=Mq

∇ ∂I

∂∇q_i

− ∂I

∂q_i

+ζi, (6) whereMqis 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]. Uti- lizing the relationship

iq_i(∂q_i/∂t)=0 that follows from the quaternion properties, and expressing the Lagrange multiplier in terms ofqiandqi, the equation of motion for the orientation (quaternion) fields can be expressed as[15]

∂qi

∂t =M_q





∇



HT[1−p(φ)] ∇qi

i(∇q_i)² _1/2





−q_i

kq_k∇



HT[1−p(φ)] ∇qk

i(∇q_i)² _1/2









+ζ_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 quater- nion properties of theqifields are retained. (ζL,iandζS,iare 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

(4)

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 dis- cretized as described by Karma and Rappel [20] and Plapp [21]. As the size of the 3D simulation box and the acces- sible 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 trans- port 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 mobilityM_φ=M_φ,0{1− 3ε_3D+4ε_3D[(∇φ)⁴_x+(∇φ)⁴_y+(∇φ)⁴_z]/|∇φ|⁴}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), wherecL= 0.466219,cS= 0.399112, and care the concentrations at the liquidus, solidus, and the initial homogeneous liquid mixture, respectively.

Since the physical thickness of the interface is in the nanome- ter 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 weret= 1.31 ns andx= 13.1 nm. The diffusion coefficient in the liquid was DL= 10⁻⁹m²/s. Unless stated otherwise, dimensionless mobil-

Mθ∝Drot Mθ/Mφ∝Drot/Dtr Drot/Dtr

enables the formation of new grains at the perimeter as detailed in Ref.[10].

5. Results and discussion

5.1. Polycrystalline solidiﬁcation 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 han- dled 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 solidiﬁcation in 3D

Dendritic solidification of four particles with random orientation grown with cubic crystal symmetries is shown inFig. 3(a) and (b). These morphologies are consistent with previous results for single-crystal dendritic growth[1,23].

(5)

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 thez-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 inFig. 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ánásy, T. Börzsönyi, T. Pusztai, Phys. Rev. Lett. 88 (2002) 206105.

(6)