The Scientific World Journal

Volume 2013 (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

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.

Abstract.

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 liquid-crystal 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], CO₂ 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 . 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 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 , for the fcc solid phase ()  kg/m³,  kg/m³ K, and  kg/m³ K² while for the liquid phase ()  kg/m³,  kg/m³ K, and  kg/m³ K². 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 the enthalpy of the liquid and the solid and fitted to with  J/kg,  J/kg K, and  J/kg K². 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 .

The specific heat capacity at constant volume is given by where  J/kg K and  J/kg K² for the solid and  J/kg K and  J/kg K² for the liquid.

The thermal diffusivity for the model is  m²/s [20] while the experimental values range from  m²/s [21] to the recent value  m²/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  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

with , , and . is the unit vector normal to the interface. The parameter 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  m/s K, , and , estimated from the kinetic coefficients for different crystal orientations;  m/s K,  m/s K, and  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 ). The orientational dependence of the interfacial free energy is described by a cubic harmonic expansion

with being an orientational averaged value of the interfacial free energy, while the coefficients describe the strength of the anisotropy. Interfacial stiffness is defined as , 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], , , 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  J/m², , , and .

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 phase-fields (solid) and (liquid). The phase-field variables fulfill the constraint , so a single phase-field variable is sufficient to describe the evolution of the phase boundaries in the system. Note that . 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

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 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 depending on the orientation of the interface. Consider the following:

where is the normalized gradient vector.

The function is the obstacle potential. It is set that for . Note that is the interface entropy density.

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

where is an anisotropic relaxation parameter dependant on temperature, with (1), and it is related to the kinetic coefficient .

The evolution of the phase-field variables is described by

The thermodynamic relation gives us , detailed in [24]. Here is the bulk free energy density

where is a monotonous function on with , and . We chose . is the latent heat and is the volumetric heat capacity both depending on the temperature . is the melting temperature.

For the isothermal case we can define two constants and consider here that the free energy contribution to the equation of motion (6) is a constant term ; note that .

Equation (6) is solved by a time dependent forward euler scheme with a second order spatial discretisation on a regular Cartesian grid.

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 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 ( g/cm³). 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 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 [25, 26]

The inner product in this sum is given by

with

where , is the number of neighboring atoms within a cut-off radius of 3.36 Å, and are 6th order spherical harmonic functions.

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 1 Å. The order parameter in a given grid point is given by the average of all particles contained in cell and its first two neighbor cells with .

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: if , and if and otherwise.

4.3. Planar Growth

An analytic formulation of the phase-field method in the sharp interface limit assumes a temperature gradient at the interface. The formulation reads

here denotes the undercooling in the interface, is the velocity, and is a kinetic coefficient. The curvature of a two dimensional interface is the sum of the two principal curvatures ,; so if . The Gibbs-Thomson-coefficient with surface tension and being the latent heat.

In equilibrium, that is, for , (11) reads as (in 2D). It follows that the critical radius . A small spherical cluster with radius is in equilibrium. For example, at  K the critical radius  Å. A cluster smaller than melts. The phase-field simulations corroborate this result: for  K the nucleus grows for  Å and melts for  Å.

For a planar interface (11) is simplified to , taking account of the absence of curvature. The kinetic coefficient adapted from MD is . So the phase-field relaxation parameter is

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).

Table 1: Kinetic coefficient of Ni in m/sK from molecular-dynamics simulations and phase-field simulation.

Interface
From (1)
PF
Sun et al. [12]

100
0.33
0.32

110
0.23
0.21

111
0.12
0.11

(a) [100]

(b) [110]

(c) [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.

The MD simulations indicate that this linearity is not more valid; the velocity is lower than that expected from extrapolation of the linear relation . 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  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 using the Andersen thermostat and barostat. The number of particles is which corresponds to a cubic box of length of about 145 Å. The initial radius of the spherical crystalline cluster is 25 Å which corresponds to about particles (see Figure 2(a)); that is much larger than the size of the critical cluster at  K which is made of about 10³ particles [27]. For  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 ( 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 and the radial velocity as . 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 . So we use to fit the data with the parameters and for each temperature . We applied the fit in the region starting with  Å, where the PF has relaxated the interface. These fits are also applied for the MD results in the second regime of growth, that is, for crystalline clusters of size of about (the points in Figure 5). According to the Gibbs-Thomson equation the radial velocity reaches an asymptotic value (), 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 ; one value is obtained at each temperature. The asymptotic radial velocities obtained from the PF and MD simulations are shown in Figure 6. The results of both methods lie between the linear relations and . 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 coefficient; we get . So is the excepted growth velocity for a sphere; it is shown as the solid line in Figure 6.

(a)

(b)

(c)

(d)

(e)

(f)

(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).

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

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

(a)

(b)

Figure 5: Evolution of the radial velocity at different temperatures. Molecular dynamic (a) and phase-field (b). In lines the fits .

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

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] orientations shrinks and the area of the [111] orientations increases; it follows that the growth rate is between and , specially below .

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 from the volume assuming a sphere and calculate the surface area for a sphere. Then [28].

The volume of the crystal is defined as and the surface area as ; this surface 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 estimated for a sphere of radius 25 Å to correct this. The 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 is below especially for higher undercoolings. The nonspherical shape is shown in Figure 2(h).

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  K and  K.

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 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üsseldorf for providing the MD data.

References

L.-Q. Chen, “Phase-field models for microstructure evolution,” Annual Review of Materials Science, vol. 32, no. 1, pp. 113–140, 2002.

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.

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.

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.

L. Gránásy, 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.

A. Karma, “Phase-field model of eutectic growth,” Physical Review E, vol. 49, no. 3, pp. 2245–2250, 1994.

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.

B. Nestler, A. A. Wheeler, L. Ratke, and C. Stöcker, “Phase-field model for solidification of a monotectic alloy with convection,” Physica D, vol. 141, no. 1-2, pp. 133–154, 2000.

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.

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.

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.

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.

A. Kerrache, J. Horbach, and K. Binder, “Molecular-dynamics computer simulation of crystal growth and melting in Al₅₀Ni₅₀,” EPL, vol. 81, no. 5, Article ID 58001, 2008.

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.

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.

J. Bragard, A. Karma, Y. H. Lee, and M. Plapp, “Linking phase-field and atomistic simulations to model dendritic solidification in highly undercooled melts,” Interface Science, vol. 10, no. 2-3, pp. 121–136, 2002.

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.

G. Tegze, T. Pusztai, G. Tóth et al., “Multiscale approach to CO₂ 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.

S. M. Foiles, “Application of the embedded-atom method to liquid transition metals,” Physical Review B, vol. 32, no. 6, pp. 3409–3415, 1985.

R. E. Rozas and J. Horbach, “A comparison of EAM models of Ni close to the melting temperature,” Submitted

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.

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.

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.

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.

P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, “Bond-orientational order in liquids and glasses,” Physical Review B, vol. 28, no. 2, pp. 784–805, 1983.

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.

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.

H. Wadell, “Volume, shape, and roundness of quartz particles,” The Journal of Geology, vol. 43, no. 3, pp. 250–280, 1935.

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.

S. Gurevich, A. Karma, M. Plapp, and R. Trivedi, “Phase-field study of three-dimensional steady-state growth shapes in directional solidification,” Physical Review E, vol. 81, no. 1, Article ID 011603, 2010.