Hindawi Publishing Corporation The Scientific World Journal Volume 2013, Article ID 564272, 8 pages http://dx.doi.org/10.1155/2013/564272

Research Article Phase-Field Simulations at the Atomic Scale in Comparison to Molecular Dynamics Marco Berghoff, Michael Selzer, and Britta Nestler Institute of Applied Materials, Karlsruhe Institute of Technology, 76133 Karlsruhe, Germany Correspondence should be addressed to Marco Berghoff; [email protected] Received 23 August 2013; Accepted 2 October 2013 Academic Editors: O. Balitskii and N. Sekido Copyright © 2013 Marco Berghoff et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Early solidification is investigated using two different simulation techniques: the molecular dynamics (MD) and the phase-field (PF) methods. While the first describes the evolution of a system on the basis of motion equations of particles, the second grounds on the evolution of continuous local order parameter field. The aim of this study is to probe the ability of the mesoscopic phase-field method to make predictions of growth velocity at the nanoscopic length scale. For this purpose the isothermal growth of a spherical crystalline cluster embedded in a melt is considered. The system in study is Ni modeled with the embedded atom method (EAM). The bulk and interfacial properties required in the PF method are obtained from MD simulations. Also the initial configuration obtained from MD data is used in the PF as input. Results for the evolution of the cluster volume at high and moderate undercooling are presented.

1. Introduction The phase-field is a powerful method to describe solidification phenomena [1, 2] on the mesoscopic length scale. It has been used to model homogeneous and heterogeneous nucleation, microstructure formation in solids, and motion of grain boundaries. The phase-field models includes formulations for pure substances [3], for multicomponent systems [4], and for polycrystalline structure [5] and solidification in eutectic [6, 7], peritectic [7], and monotectic [8] systems. Other simulation techniques for dendrite growth are cellular automata [9] and hybrid methods such as the multiscale diffusion Monte Carlo (DMC) [10, 11]. The phase-field method requires previous knowledge of the material properties of the system in study. The input includes bulk properties such as density, heat capacity and latent heat, and others such as interfacial and kinetic growth coefficients, being the latter properties which are hardly accessible in experiments. Here molecular simulations play a fundamental role, since they provide a link between an interaction potential and all the required properties. Kinetic coefficients, interfacial free energies, and their dependence on the crystal orientation, of pure metals and alloys, can be

directly obtained from simulations of inhomogeneous liquidcrystal systems [12–15]. There are previous phase-field studies in which the material parameters have been obtained from molecular simulations to model dendritic growth in pure Ni [16, 17], CO2 hydrates [18], and binary alloys. In this work we compare MD and PF simulation results for the isothermal growth velocity of a nanoscopic spherical crystalline cluster of Ni embedded in the melt. For this purpose we use precise thermodynamic and kinetic data obtained from molecular dynamics and use equivalent initial configuration in the phase-field simulations.

2. Properties of Ni (EAM) The properties of Ni obtained from molecular dynamic simulations are given in this section. The embedded atom potential of Foiles [19] is used. The thermal dependence of bulk properties is determined from independent simulations of equilibrated homogeneous systems at pressure 𝑝 = 0. Liquid phase is prepared by melting of an initial crystalline phase, by setting the temperature of the thermostat well above the melting temperature. Crystalline phase is obtained by

2

The Scientific World Journal

relaxation of the perfect crystalline phase. The samples are prepared in NpT simulation runs and the production in NVE ensemble. Densities of the bulk phases are fitted to the function 𝜌𝛼 (𝑇) = 𝑎 + 𝑏 ⋅ 𝑇 + 𝑐 ⋅ 𝑇2 , for the fcc solid phase (𝛼 = s) 𝑎 = 8901.6 kg/m3 , 𝑏 = −0.20379 kg/m3 K, and 𝑐 = −0.0000614202 kg/m3 K2 while for the liquid phase (𝛼 = ℓ) 𝑎 = 8992.26 kg/m3 , 𝑏 = −0.667037 kg/m3 K, and 𝑐 = 0.0000331612 kg/m3 K2 . The parameters were obtained in the temperature intervals 1000 K to 3000 K for the liquid and 300 K to 1900 K for the solid. In both cases parts of the metastable phases are included in the fit procedure. The latent heat is obtained from the difference between ̃ 0 (𝑇) = the enthalpy of the liquid and the solid and fitted to 𝐿 2 𝑎 + 𝑏 ⋅ 𝑇 + 𝑐 ⋅ 𝑇 with 𝑎 = −15980.2 J/kg, 𝑏 = 324.745 J/kg K, and 𝑐 = −0.0810984 J/kg K2 . These values are determined in the range from 1000 K to 1900 K. The latent heat expressed in units of energy per volume of the solid phase is given by ̃ 0 (𝑇) ⋅ 𝜌s (𝑇). 𝐿s (𝑇) = 𝐿 The specific heat capacity at constant volume is given by 𝑐̃V𝛼 (𝑇) = 𝑎 + 𝑏 ⋅ 𝑇 where 𝑎 = 419.452 J/kg K and 𝑏 = 0.020388 J/kg K2 for the solid and 𝑎 = 563.024 J/kg K and 𝑏 = −0.06952 J/kg K2 for the liquid. The thermal diffusivity for the model is 2.1⋅10−7 m2 /s [20] while the experimental values range from 120 ⋅ 10−7 m2 /s [21] to the recent value 8.7⋅10−7 m2 /s [22] (obtained in microgravity conditions). The value obtained from simulations is at least 4 times lower than that in experiments, and the difference is due to the electronic contribution to thermal transport in metals which are not accounted in classical molecular dynamics. The melting temperature of the model is precisely estimated from simulations of inhomogeneous solid-liquid system at the point where the velocity of growth is zero. For the model 𝑇𝑚 = 1748 K. The kinetic coefficient is estimated from linear relationship between planar growth velocity and undercooling close to coexistence. Growth velocities are obtained from independent simulations for the crystal orientations [100], [110], and [111]. The kinetic growth coefficient depends also on the orientation of the crystal. Its orientational dependence is represented here by the expansion 𝑘 (̂ 𝑛) = 1 − 3𝜖𝑘 + 4𝜖𝑘 𝑄 + 𝛿𝑘 (𝑃 + 30𝑆) 𝑘0

(1)

with 𝑄 = 𝑛14 + 𝑛24 + 𝑛34 , 𝑃 = 𝑛16 + 𝑛26 + 𝑛36 , and 𝑆 = 𝑛12 𝑛22 𝑛32 . 𝑛̂ = (𝑛1 , 𝑛2 , 𝑛3 ) is the unit vector normal to the interface. The parameter 𝑘0 is a value of the magnitude of the kinetic coefficient; 𝜖𝑘 and 𝛿𝑘 are the strength of the anisotropy of the kinetic coefficient. The values for the model EAM F85 are 𝑘0 = 0.319205 m/s K, 𝜖𝑘 = −0.196511, and 𝛿𝑘 = 0.230331, estimated from the kinetic coefficients for different crystal orientations; 𝑘100 = 0.33 m/s K, 𝑘110 = 0.23 m/s K, and 𝑘111 = 0.12 m/s K. The interfacial stiffness and interfacial free energy were obtained from the analysis of capillary waves spectrum of large solid-liquid interfaces at coexistence (i.e., at 𝑇 = 𝑇𝑚 for

𝑃 = 0). The orientational dependence of the interfacial free energy 𝛾 is described by a cubic harmonic expansion 𝛾 (̂ 𝑛) 3 17 = 1 + 𝜖1 (𝑄 − ) + 𝜖2 (3𝑄 + 66𝑆 − ) 𝛾0 5 7 94 33 + 𝜖3 (5𝑄 − 16𝑆 − 𝑄 + ) 13 13

(2)

2

with 𝛾0 being an orientational averaged value of the interfacial free energy, while the coefficients 𝜖𝑖 describe the strength of the anisotropy. Interfacial stiffness is defined as 𝛾̃𝛼𝛽 = 𝛾 + 𝜕2 𝛾/𝜕̂ 𝑛𝛼 𝜕̂ 𝑛𝛽 , where 𝑛𝛼 and 𝑛𝛽 indicate unit vectors tangent to the interface plane [23]. The parameters of (2) are obtained from the system of equations for the stiffnesses 𝛾̃𝛼𝛽 and their numerical values are obtained from the analysis of the capillary waves spectrum [14]. Here we consider the crystal orientations [100], [110][110] , [110][001] , and [111]; the interface energy for the [110] orientation depends on the parallel direction denoted by the subscripts. In respect to the four different interface energies we used a cubic harmonic expansion with four fitting parameters. Accurate values of these variables and parameters for the model EAM F85 are 𝛾0 = 0.302 J/m2 , 𝜖1 = 0.10191, 𝜖2 = −0.00134, and 𝜖3 = 0.00876.

3. Phase-Field Model for Pure Material Isothermal crystal growth for pure material is modeled by the variables of the internal energy density 𝑒 and of two phasefields 0 ≤ 𝜙s ≤ 1 (solid) and 0 ≤ 𝜙ℓ ≤ 1 (liquid). The phase-field variables fulfill the constraint 𝜙s + 𝜙ℓ = 1, so a single phase-field variable 𝜙 = 𝜙s is sufficient to describe the evolution of the phase boundaries in the system. Note that 𝜙ℓ = 1 − 𝜙. The variable 𝜙(𝑥,⃗ 𝑡) denotes the local fraction at 𝑥⃗ of the considered phase at time 𝑡. The phase-field model is based on an entropy functional to ensure consistency with classical irreversible thermodynamics 1 S (𝑒, 𝜙) = ∫ 𝑠 (𝑒, 𝜙) − (𝜖𝑎 (𝜙, ∇𝜙) + 𝑤 (𝜙)) 𝑑𝑥. 𝜖 Ω

(3)

The bulk entropy density 𝑠 depends on the internal energy density 𝑒 and the phase-field variable 𝜙(𝑥,⃗ 𝑡). The functions 𝑎(𝜙, ∇𝜙) and 𝑤(𝜙) reflect the thermodynamics of the interfaces and 𝜖 = 4 is a small length scale parameter related to the thickness of the diffuse interface. In detail the gradient entropy with anisotropy is adopted from (2) and modeled by the factor (𝐴(̂ 𝑛))2 depending on the orientation of the interface. Consider the following: 𝑎 (𝜙, ∇𝜙) =

𝛾0 󵄨 󵄨2 󵄨 󵄨2 𝑛))2 󵄨󵄨󵄨∇𝜙󵄨󵄨󵄨 = 𝛾 (̂ 𝑛) 󵄨󵄨󵄨∇𝜙󵄨󵄨󵄨 , (𝐴 (̂ 𝑇

(4)

where 𝑛̂ = −∇𝜙/|∇𝜙| is the normalized gradient vector. The function 𝜔(𝜙) = (16/𝜋2 )(𝛾0 /𝑇)𝜙(1−𝜙) is the obstacle potential. It is set that 𝜔(𝜙) = ∞ for 𝜙 ∉ [0, 1]. Note that 𝛾0 /𝑇 is the interface entropy density.

The Scientific World Journal

3

From (3) one derives the equations for the nonconserved phase-field variable 𝜙 by taking the functional derivatives in the form 𝜏𝜖

𝜕𝜙 𝛿S = , 𝜕𝑡 𝛿𝜙

(5)

where 𝜏 = 𝜏(𝜙, ∇𝜙) is an anisotropic relaxation parameter 𝑛) with (1), dependant on temperature, 𝜏(̂ 𝑛) = 𝐿(𝑇)/𝑇𝑇𝑚 𝑘(̂ and it is related to the kinetic coefficient 𝑘0 . The evolution of the phase-field variables is described by 𝜏𝜖

𝜕𝑎 (𝜙, ∇𝜙) 𝜕𝑎 (𝜙, ∇𝜙) 𝜕𝜙 1 𝜕𝑤 (𝜙) = 𝜖 (∇ ⋅ − )− 𝜕𝑡 𝜕∇𝜙 𝜕𝜙 𝜖 𝜕𝜙 +

𝜕𝑠 (𝑒, 𝜙) . 𝜕𝜙 (6)

The thermodynamic relation 𝑒 = 𝑓+𝑇𝑠 gives us 𝜕𝑠(𝑒, 𝜙)/𝜕𝜙 = −(1/𝑇)(𝜕𝑓(𝑇, 𝜙)/𝜕𝜙), detailed in [24]. Here is the bulk free energy density 𝑓 (𝑇, 𝜙) = ∑ 𝐿𝛼 (𝑇) 𝛼∈{s,ℓ} 𝑇

+ ∑ (∫ 𝛼∈{s,ℓ}

𝑇𝑚

by increasing the temperature well above the melting point while the atoms in the solid region are constrained to remain at fixed positions. The temperature is linearly increased while an isotropic barostat is applied. In the next step, the stable melt at high temperature is cooled down to the temperature of the crystal. The melt reaches a metastable (undercooled) state. At the end of this run, a crystalline cluster embedded in the liquid phase at the same temperature is obtained. In the last step, the relaxation of the atoms at the interface is made by an NpT simulation of some picoseconds in which all particles are allowed to move. 4.2. Local Order Parameter. At microscopic level a configuration is described by positions and velocities of particles at a given time. The crystalline and liquid regions in this configuration can be identified by a suitable definition of the local order parameter, a function of the coordinates of a particle and its neighbors which adopts different values in the crystalline and in the disordered liquid phases. Here we use the bond local order parameter 𝑞6 𝑞6 [25, 26] 𝑍

𝑞6 𝑞6 (𝑖) :=

𝑇 − 𝑇𝑚 ℎ (𝜙𝛼 ) 𝑇𝑚 𝑠V𝛼

𝑇 ̃ ̃ 𝑑𝑇−𝑇 ̃ ̃ 𝑑𝑇 ) ℎ (𝜙𝛼 ) , (𝑇) ∫ 𝑠V𝛼 (𝑇) ̃ 𝑇𝑚 𝑇

4. Initial Configuration and Calibration 4.1. Preparation of the Sample. Growth simulations begin from a configuration made of a stabilized crystalline spherical cluster embedded in the melt. The preparation of the initial configuration involves a sequence of NpT simulation runs. First the crystalline solid phase is equilibrated at 𝑝 = 0 starting from atoms arranged in a perfect fcc structure in a cubic box with an initial density close to the experimental value at room temperature (𝜌 = 8.9 g/cm3 ). The temperature in this step is chosen equal to the one used later in the simulation of growth. In the second step two regions are defined in the simulation box: the melting zone and the atoms in the crystalline phase (which is chosen here as a sphere of radius ˚ The region outside the crystalline zone is melted 𝑟 = 25 A).

(8)

The inner product in this sum is given by 6

q6 (𝑖) ⋅ q6 (𝑗) = ∑ 𝑞̃6𝑚 (𝑖) 𝑞̃6𝑚 (𝑗)

(7) where ℎ(𝜙) is a monotonous function on [0, 1] with ℎ(0) = 0, ℎ(1) = 1 and ℎ󸀠 (0) = ℎ󸀠 (1) = 0. We chose ℎ(𝜙) = 𝜙3 (6𝜙2 − 15𝜙 + 10). 𝐿𝛼 (𝑇) is the latent heat and 𝑠V𝛼 (𝑇) = 𝜌𝛼 (𝑇)̃𝑐V𝛼 (𝑇) is the volumetric heat capacity both depending on the temperature 𝑇. 𝑇𝑚 is the melting temperature. For the isothermal case we can define two constants 𝑓(𝑇, 𝜙) =: ∑𝛼∈{s,ℓ} 𝐶𝑇𝛼 ℎ(𝜙𝛼 ) and consider here that the free energy contribution to the equation of motion (6) is a constant term 𝜕𝑓/𝜕𝜙 = (𝐶𝑇s − 𝐶𝑇ℓ )ℎ󸀠 (𝜙); note that ℎ󸀠 (𝜙s ) = −ℎ󸀠 (𝜙ℓ ). Equation (6) is solved by a time dependent forward euler scheme with a second order spatial discretisation on a regular Cartesian grid.

1 𝑖 ∑ q (𝑖) ⋅ q6 (𝑗) . 𝑍𝑖 𝑗=1 6

∗

(9)

𝑚=−6

with −1/2

6 󵄨2 󵄨 𝑞̃6𝑚 (𝑖) = 𝑄6𝑚 (𝑖) ( ∑ 󵄨󵄨󵄨󵄨𝑄6𝑚 (𝑖)󵄨󵄨󵄨󵄨 )

,

(10)

𝑚=−6

𝑍

𝑖 where 𝑄6𝑚 (𝑖) = (1/𝑍𝑖 ) ∑𝑗=1 𝑌6𝑚 (𝜃(𝑟𝑖𝑗⃗ ), 𝜙(𝑟𝑖𝑗⃗ )), 𝑍𝑖 is the number of neighboring atoms within a cut-off radius of ˚ and 𝑌6𝑚 are 6th order spherical harmonic functions. 3.36 A, In order to obtain an equivalent initial configuration for the PF method, the local order parameter values given at atom positions are mapped on a regular grid with cell-width of ˚ The order parameter in a given grid point is given by 1 A. the average of all particles contained in cell u = (𝑖, 𝑗, 𝑘) and its first two neighbor cells k = (𝑚, 𝑛, 𝑙) with |u − k| < 2. The order parameter values obtained from MD simulation were scaled so that the volume of the crystal cluster is equal in the coarse-grained and in the original initial configurations. The distribution of the order parameter in the bulk phases is Gaussian centered at 𝜇 with standard deviation 𝜎. The definition of order parameter in the phase-field method is as follows: 𝜙 = 0 if 𝑞6 𝑞6 < 𝑎 := 𝜇ℓ − 2𝜎ℓ , and 𝜙 = 1 if 𝑞6 𝑞6 > 𝑏 := 𝜇s + 2𝜎s and 𝜙 = (𝑞6 𝑞6 − 𝑎)/(𝑏 − 𝑎) otherwise.

4.3. Planar Growth. An analytic formulation of the phasefield method in the sharp interface limit assumes a temperature gradient at the interface. The formulation reads Δ𝑇 = 𝛽V + 𝜅Γ;

(11)

4

The Scientific World Journal

Table 1: Kinetic coefficient of Ni in m/sK from molecular-dynamics simulations and phase-field simulation. Interface [100] [110] [111]

From (1) 0.33 0.23 0.12

PF 0.32 0.21 0.11

Sun et al. [12] 0.358 ± 0.022 0.255 ± 0.016 0.241 ± 0.04

here Δ𝑇 denotes the undercooling in the interface, V is the velocity, and 𝛽 is a kinetic coefficient. The curvature 𝜅 of a two dimensional interface is the sum of the two principal curvatures 1/𝑟1 , 1/𝑟2 ; so 𝜅 = 2/𝑟 if 𝑟1 = 𝑟2 = 𝑟. The GibbsThomson-coefficient Γ = 𝜎𝑇/𝐿 with 𝜎 surface tension and 𝐿 being the latent heat. In equilibrium, that is, for V = 0, (11) reads as Δ𝑇 = 𝜅Γ = Γ/𝑟 (in 2D). It follows that the critical radius 𝑟𝑐 = Γ/Δ𝑇. A small spherical cluster with radius 𝑟𝑐 is in equilibrium. For ˚ A cluster example, at 𝑇 = 1550 K the critical radius 𝑟𝑐 = 9.7 A. smaller than 𝑟𝑐 melts. The phase-field simulations corroborate ˚ and this result: for 𝑇 = 1550 K the nucleus grows for 𝑟 = 10 A ˚ melts for 𝑟 = 9.5 A. For a planar interface (11) is simplified to Δ𝑇 = 𝛽V, taking account of the absence of curvature. The kinetic coefficient adapted from MD is 𝑘 = 1/𝛽. So the phase-field relaxation parameter is 𝜏 (̂ 𝑛) =

𝜏 𝐿 (𝑇) = 0 . 𝑇𝑇𝑚 𝑘 (̂ 𝑛) 𝑘 (̂ 𝑛)

(12)

Simulations of a planar front in comparison to the analytic solution and molecular dynamics simulations are shown in Figure 1. Additionally we list the kinetic coefficient in Table 1. Both methods converge at low undercooling to the solution given by (11). At higher undercooling velocities obtained from PF simulations still exhibit a linear dependence with undercooling Δ𝑇; this is given in the range from 1400 K to 𝑇𝑚 ; the fit of the kinetic coefficient is made in this range and agrees with the expected kinetics from (1). The MD simulations indicate that this linearity is not more valid; the velocity is lower than that expected from extrapolation of the linear relation V𝐼 = 𝑘Δ𝑇. This effect can also be recognized in the results of Hoyt et al. [15]. Many effects influence the growth velocity of MD at higher undercooling, with the increase of order in the liquid bulk being the most important which leads to the formation of small crystalline clusters in the bulk liquid. The phase-field shows this nonlinearity only for very high undercoolings Δ𝑇 > 500 K.

5. Spherical Cluster Growth Growth simulations starting from the samples prepared as described in Section 4.1 are performed using both MD and PF methods. At microscopic level NpT simulations are performed at 𝑝 = 0 using the Andersen thermostat and barostat. The number of particles is 𝑁 = 256000 which ˚ The corresponds to a cubic box of length of about 145 A. ˚ initial radius of the spherical crystalline cluster is 25 A which corresponds to about 5.5 ⋅ 103 particles (see Figure 2(a));

that is much larger than the size of the critical cluster at 𝑇 = 1500 K which is made of about 103 particles [27]. For 𝑇 ≥ 1600 K the initial size of the cluster is lower than the critical value and it melts as expected. PF simulations are performed imposing the isothermal condition starting with the same configuration as in MD but converted as described in Section 4.2. Results for the evolution of the cluster volume are shown in Figure 3; snapshots from the MD and PF simulations are shown in Figure 2. At moderate undercooling (𝑇 > 1550 K) the MD and PF show comparable growth rates. For higher undercooling the growth in the phase-field model is notably faster than in the MD method. The discrepancy becomes more pronounced at lower temperature and higher simulation time. Assuming that the cluster conserves its spherical shape during the first stage of growth we define the radius as 𝑟 = (3𝑉/4𝜋)1/3 and the radial velocity as V = 𝑑𝑟/𝑑𝑡. Figure 4 shows the evolution of radial velocity in the MD simulation; three regimes can be identified; initially the velocity tends to decrease, then exhibits a linear increase, and finally decreases again. The first decrease is due to the preparation method of the sample; the system relaxes while the interface is formed. Growth and relaxation cannot be decoupled; both occur simultaneously in this initial regime. In the second regime the growth velocity increases linearly in time. The decrease of velocity in the last regime is due to the finite size of the sample; when the cluster is large enough it interacts with its periodic images through the periodic boundaries of the simulation box and loses its spherical shape. In the PF simulation no extra relaxation regime is observed; the dynamic is determined by the initial definition after the interface is relaxated. Moreover, in the PF the cluster is not affected by finite size effects as in the MD simulation since no periodic boundaries are required. In comparison the PF exhibits an apparent linear increase from the beginning. This means that the starting radius of the nucleus is not the same for both simulation methods, so a comparison for the same radius instead of the time is more adequate. Figure 5 shows the radial velocity versus cluster radius at different temperatures. Results obtained from PF simulations are well described by Gibbs-Thomson equation (11). The analytical velocity is V = 𝑘avg (Δ𝑇 − 2Γ/𝑟). So we use 𝑉𝑇 (𝑟) = 𝑎𝑇 + 𝑏𝑇 /𝑟 to fit the data with the parameters 𝑎𝑇 and 𝑏𝑇 for each temperature 𝑇. We applied the fit in the region starting ˚ where the PF has relaxated the interface. with 𝑟 = 30 A, These fits are also applied for the MD results in the second regime of growth, that is, for crystalline clusters of size of ˚ ≾ 𝑟 ≾ 60 A ˚ (the points in Figure 5). According about 40 A to the Gibbs-Thomson equation the radial velocity reaches an asymptotic value (V𝑇,∞ = V𝑇,𝑟 → ∞ ), an ideal situation where the cluster is large enough and conserves its shape. This limit cannot be directly obtained from simulations first because the simulations are performed in a finite volume and also because the crystalline cluster would adopt a nonspherical shape due to the anisotropy of the interfacial energy and growth coefficient. We estimate this hypothetical limit as a basis of comparison between the PF and MD simulations. For this purpose we extrapolate the fit and obtain V𝑇,∞ = 𝑎𝑇 ; one value is obtained at each temperature. The asymptotic

5

250

250

200

200 Growth velocity (m/s)

Growth velocity (m/s)

The Scientific World Journal

150

100

50

0 1000

150

100

50

1100

1200

1300 1400 1500 Temperature (K)

PF Fitted PF Sun et al.

1600

0 1000

1700

1100

1200

PF Fitted PF Sun et al.

Hoyt et al. k100 ΔT MD

1600

1700

1600

1700

Hoyt et al. k110 ΔT (b) [110]

250

250

200

200 Growth velocity (m/s)

Growth velocity (m/s)

(a) [100]

150

100

150

100

50

50

0 1000

1300 1400 1500 Temperature (K)

1100

1200

1300 1400 1500 Temperature (K)

PF Fitted PF Sun et al.

1600

1700

Hoyt et al. k111 ΔT (c) [111]

0 1000

1100

PF [100] PF [110] PF [111]

1200

1300 1400 1500 Temperature (K) Fitted [100] Fitted [110] Fitted [111] (d) Phase-field

Figure 1: Growth velocity of a planar front at different temperatures for the interface orientations (a) [100], (b) [110], and (c) [111] in comparison to molecular dynamic results. The values of the kinetic coefficient can be found in Table 1. (d) shows the phase-field results for all orientations.

radial velocities obtained from the PF and MD simulations are shown in Figure 6. The results of both methods lie between the linear relations V100 = 𝑘100 Δ𝑇 and V111 = 𝑘111 Δ𝑇. For the MD the results show a similar behavior as observed in the planar growth case; the deviation from the linear relation at high undercooling observed in the PF simulations is due to the correction introduced in the kinetic coefficient. As shown in Section 4.3 the crystal grows in different crystal orientations with a different velocity. In the PF the kinetic anisotropy is responsible for this. The planar front simulations from Figure 1 show that the kinetics are very

well reproducible from the PF simulations. If we assume that we have a perfect sphere the mean kinetic can be calculated as a surface area integral of the unit sphere over the kinetic 𝑛)𝑑𝐴/ ∫𝑆2 1𝑑𝐴 ≈ 0.175. So coefficient; we get 𝑘avg = ∫𝑆2 𝑘(̂ V(Δ𝑇) = 𝑘avg Δ𝑇 is the excepted growth velocity for a sphere; it is shown as the solid line in Figure 6. At this early stage of growing the crystal is very similar to a sphere. The preferred growing in [100] directions and the surface anisotropy formed the crystal to its equilibrium shape; for Ni this is similar to an octahedron; in the limit this shape has only [111] directions. So the area of the [100]

6

The Scientific World Journal

(a)

(b)

(e)

(c)

(f)

(d)

(g)

(h)

Figure 2: Snapshots of the simulations at 1550 K. MD simulations at (a) initial time 80 ps, (c) 100 ps, and (d) 220 ps. The evolution of the crystalline cluster from PF simulations with the same volume is shown in the sequence from (e) to (h).

12 40

 = dr/dt (m/s)

10

3

˚ ) V (105 A

8 6 4

30

20

10 2 0

0 0

50

1400 K: MD 1450 K: MD 1500 K: MD 1525 K: MD 1550 K: MD 1575 K: MD

100 t (ps)

150

200

PF PF PF PF PF PF

Figure 3: Evolution of the crystalline cluster volume at different temperatures. MD results are the average over 8 independent simulation runs.

orientations shrinks and the area of the [111] orientations increases; it follows that the growth rate is between 𝑘avg and 𝑘111 , specially below 𝑘avg .

0

50

MD 1400 K MD 1450 K MD 1500 K

100 t (ps)

150

200

MD 1525 K MD 1550 K MD 1575 K

Figure 4: Evolution of the radial velocity at different temperatures obtained from MD simulations.

We use the sphericity 𝑠 as the measure for the sphericality of the crystal as follow. For a given volume 𝑉 and a given surface area 𝑆 of the crystal we calculate the radius 𝑟3 = 3𝑉/4𝜋 from the volume assuming a sphere and calculate the surface area 𝑆𝑉 = 4𝜋𝑟2 for a sphere. Then 𝑠 := 𝑆𝑉/𝑆 [28]. The volume of the crystal is defined as 𝑉 := ∫Ω 𝜙𝑑𝑉 and the surface area as 𝑆 := ∫Ω ‖∇𝜙‖𝑑𝑉; this surface

7

60

60

50

50

40

40

 = dr/dt (m/s)

 = dr/dt (m/s)

The Scientific World Journal

30 20 10

30 20 10

0 30

35

40

45

50

55

60

0 30

35

40

˚ r (A) MD 1400 K MD 1450 K MD 1500 K

45

50

55

60

˚ r (A)

MD 1525 K MD 1550 K MD 1575 K

PF 1525 K PF 1550 K PF 1575 K

PF 1400 K PF 1450 K PF 1500 K

(a)

(b)

Figure 5: Evolution of the radial velocity at different temperatures. Molecular dynamic (a) and phase-field (b). In lines the fits V𝑇 (𝑟) := 𝑎𝑇 + 𝑏𝑇 /𝑟.

80

1

70 0.99

50 Sphericity s

 = dr/dt (m/s)

60

40 30

0.98

0.97

20 0.96

10 0 1350

1400

1450

PF MD k100 ΔT

1500

1550 T (K)

1600

1650

1700

1750

k110 ΔT k111 ΔT kavg ΔT

Figure 6: Asymptotic radial velocity versus temperature. Values of V∞ are obtained from the asymptotic limit of the fits shown in Figure 5.

area measurement method assumes a smooth nondeformed interface; this is not given for PF interfaces if there are driving forces like curvatures and undercoolings. We use a factor ˚ to correct this. The estimated for a sphere of radius 25 A surface area is multiplied by this factor 1.089. Figure 7 shows the sphericity over the radius; the crystal becomes non spherical for longer growing and faster for higher undercoolings. This means that the velocity of growing

0.95 25

30

35

40

45

50

55

60

65

˚ Radius r (A) 1400 K 1450 K

1500 K 1550 K

Figure 7: Sphericity of the crystal for different undercoolings. The contour shows the shape of the [100] cutting plane through the center of the crystal for 𝑇 = 1400 K and 1550 K.

is below 𝑘avg especially for higher undercoolings. The nonspherical shape is shown in Figure 2(h).

6. Conclusion Many have shown the linking between MD and PF [16, 17]. So it is a common way to obtain MD parameters for PF simulations. The PF simulations on the natural PF scale, for example, mesoscopic scale, produce correct results compared

8 with experiments [29, 30]. We have confirmed that a PF with MD parameters is correct even on the atomistic scale. This is not surprising, since we are very close to the sharp interface limits. The mesoscopic view at the atomistic scale provides the same growth rates as an atomic simulation method, although the atoms are not resolved in the PF. The fluctuations of the interface on the atomistic length scale will be investigated in a consecutive work. However, you will not simulate too high undercooled melts, where the mesoscopic view of the PF will not consider all the nonlinearity from the MD interface. We have used the sphericity to measure the shape and explain the discrepancy to the analytic linear velocity. While it is known that PF works correctly on mesoscopic scale and with this study it is shown that PF also works on the atomistic scale, one is able to do multiscale simulations starting on the atomistic scale with one simulation method.

Acknowledgments The authors gratefully acknowledge the financial support by the German Research Foundation (DFG) in the Priority Program SPP 1296. They like to thank J. Horbach and R. Rozas from Heinrich Heine University of D¨usseldorf for providing the MD data.

References [1] L.-Q. Chen, “Phase-field models for microstructure evolution,” Annual Review of Materials Science, vol. 32, no. 1, pp. 113–140, 2002. [2] W. J. Boettinger, J. A. Warren, C. Beckermann, and A. Karma, “Phase-field simulation of solidification,” Annual Review of Materials Science, vol. 32, no. 1, pp. 163–194, 2002. [3] A. Karma and W.-J. Rappel, “Quantitative phase-field modeling of dendritic growth in two and three dimensions,” Physical Review E, vol. 57, no. 4, pp. 4323–4349, 1998. [4] B. Nestler, H. Garcke, and B. Stinner, “Multicomponent alloy solidification: phase-field modeling and simulations,” Physical Review E, vol. 71, no. 4, Article ID 041609, 2005. [5] L. Gr´an´asy, T. Pusztai, and J. A. Warren, “Modelling polycrystalline solidification using phase field theory,” Journal of Physics: Condensed Matter, vol. 16, no. 41, pp. R1205–R1235, 2004. [6] A. Karma, “Phase-field model of eutectic growth,” Physical Review E, vol. 49, no. 3, pp. 2245–2250, 1994. [7] B. Nestler and A. A. Wheeler, “A multi-phase-field model of eutectic and peritectic alloys: numerical simulation of growth structures,” Physica D, vol. 138, no. 1-2, pp. 114–133, 2000. [8] B. Nestler, A. A. Wheeler, L. Ratke, and C. St¨ocker, “Phase-field model for solidification of a monotectic alloy with convection,” Physica D, vol. 141, no. 1-2, pp. 133–154, 2000. [9] M. F. Zhu and C. P. Hong, “A modified cellular automaton model for the simulation of dendritic growth in solidification of alloys,” ISIJ International, vol. 41, no. 5, pp. 436–445, 2001. [10] P. Zhu and R. W. Smith, “Dynamic simulation of crystal growth by Monte Carlo method-I. Model description and kinetics,” Acta Metallurgica et Materialia, vol. 40, no. 4, pp. 683–692, 1992. [11] M. Plapp and A. Karma, “Multiscale finite-difference-diffusion Monte-Carlo method for simulating dendritic solidification,” Journal of Computational Physics, vol. 165, no. 2, pp. 592–619, 2000.

The Scientific World Journal [12] D. Y. Sun, M. Asta, and J. J. Hoyt, “Kinetic coefficient of Ni solid-liquid interfaces from molecular-dynamics simulations,” Physical Review B, vol. 69, no. 2, Article ID 024108, 2004. [13] A. Kerrache, J. Horbach, and K. Binder, “Molecular-dynamics computer simulation of crystal growth and melting in Al50 Ni50 ,” EPL, vol. 81, no. 5, Article ID 58001, 2008. [14] R. E. Rozas and J. Horbach, “Capillary wave analysis of rough solid-liquid interfaces in nickel,” EPL, vol. 93, no. 2, Article ID 26006, 2011. [15] J. J. Hoyt, M. Asta, and A. Karma, “Atomistic and continuum modeling of dendritic solidification,” Materials Science and Engineering R, vol. 41, no. 6, pp. 121–163, 2003. [16] J. Bragard, A. Karma, Y. H. Lee, and M. Plapp, “Linking phasefield and atomistic simulations to model dendritic solidification in highly undercooled melts,” Interface Science, vol. 10, no. 2-3, pp. 121–136, 2002. [17] B. Nestler, M. Selzer, and D. Danilov, “Phase-field simulations of nuclei and early stage solidification microstructures,” Journal of Physics: Condensed Matter, vol. 21, no. 46, Article ID 464107, 2009. [18] G. Tegze, T. Pusztai, G. T´oth et al., “Multiscale approach to CO2 hydrate formation in aqueous solution: phase field theory and molecular dynamics. Nucleation and growth,” The Journal of Chemical Physics, vol. 124, no. 23, Article ID 234710, 2006. [19] S. M. Foiles, “Application of the embedded-atom method to liquid transition metals,” Physical Review B, vol. 32, no. 6, pp. 3409–3415, 1985. [20] R. E. Rozas and J. Horbach, “A comparison of EAM models of Ni close to the melting temperature,” Submitted. [21] V. Y. Zinov’yev, V. F. Polev, S. G. Taluts, G. P. Zinov’yeva, and S. A. Il’inykh, “Diffusivity and thermal conductivity of 3D-transition metals in solid and liquid states,” Physics of Metals and Metallography, vol. 61, no. 6, pp. 85–92, 1986. [22] K. Nagata, H. Fukuyama, K. Taguchi, H. Ishii, and M. Hayashi, “Thermal conductivity of molten Al, Si and Ni measured under microgravity,” High Temperature Materials and Processes, vol. 22, no. 5-6, pp. 267–273, 2003. [23] J. J. Hoyt, M. Asta, T. Haxhimali et al., “Crystal-melt interfaces and solidification morphologies in metals and alloys,” MRS Bulletin, vol. 29, no. 12, pp. 935–939, 2004. [24] H. Garcke, B. Nestler, and B. Stinner, “A diffuse interface model for alloys with multiple components and phases,” SIAM Journal on Applied Mathematics, vol. 64, no. 3, pp. 775–799, 2004. [25] P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, “Bondorientational order in liquids and glasses,” Physical Review B, vol. 28, no. 2, pp. 784–805, 1983. [26] P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, “Numerical evidence for bcc ordering at the surface of a critical fcc nucleus,” Physical Review Letters, vol. 75, no. 14, pp. 2714–2717, 1995. [27] J. Bokeloh, R. E. Rozas, J. Horbach, and G. Wilde, “Nucleation barriers for the liquid-to-crystal transition in Ni: experiment and simulation,” Physical Review Letters, vol. 107, no. 14, Article ID 145701, 2011. [28] H. Wadell, “Volume, shape, and roundness of quartz particles,” The Journal of Geology, vol. 43, no. 3, pp. 250–280, 1935. [29] B. Echebarria, R. Folch, A. Karma, and M. Plapp, “Quantitative phase-field model of alloy solidification,” Physical Review E, vol. 70, no. 6, Article ID 061604, 2004. [30] S. Gurevich, A. Karma, M. Plapp, and R. Trivedi, “Phasefield study of three-dimensional steady-state growth shapes in directional solidification,” Physical Review E, vol. 81, no. 1, Article ID 011603, 2010.