Molecular dynamics calculation of thermodynamic properties of iron solidification

The aim of this study is to identify the best available inter-atomic potentials for molecular dynamics (MD) calculation of solidification of iron and then to use the best potential to calculate thermodynamic properties such as equilibrium melting temperature, enthalpy, heat capacity and solid-liquid interfacial free energy. Our study reveals that embedded atom method (EAM) potential developed by Ackland et al. [2004 J. Phys.: Condens Matter. 16 S2629] appears to be the most accurate model for MD simulation of iron solidification. Simulations with the above EAM model predict the equilibrium melting temperature of iron is 1790K, the solid-liquid interfacial energy 214 mJ/m2. The difference with the experimental data is 1.2%, and 4.9% respectively.


Introduction
Thermodynamic properties of solid-liquid interface affect solidification structure and properties of the solidified components [1]. How to obtain accurate thermodynamic properties of metal solidification is crucial for both industry application and fundamental research. During the last decades, various experimental technologies [2,3] such as electromagnetic, electrostatic and acoustic levitation are employed to measure the thermodynamic data with high accuracy. However, measuring the thermodynamic properties is expensive and some data such as solid-liquid interfacial free energy are difficult to measure [4]. As a consequence, experimentally measured thermodynamic data at solidliquid interface is rather limited. To overcome the experimental difficulties, computer simulations based on atomistic models such as molecular dynamics (MD) simulations have been developed and becoming a powerful tool for gaining fundamental knowledge of materials processes and properties over the past two decades [1,5].
Atomistic simulation starts with establishing a potential model which describes atomic interactions in the material of interest [5]. The accuracy of the simulation results is determined by the potential models. Due to the industrial importance of iron and steel, in recent years, several semi-empirical potential models for iron [6][7][8][9] have been developed and used to calculate solid-liquid interface properties including solid-liquid interfacial free energy [10,11]. However, the calculated properties vary with the potentials selected for the calculation. For example, the calculated equilibrium melting temperatures using Finnis-Sinclair potential 6 and EAM potential developed by Mendelev et al. [9] are 2400K, and 1772K, respectively. The difference is significant compared to the experimental value 1811K [11]. The main purpose of the present work is to first identify an appropriate inter-atomic potential for MD simulations of iron solidification and then use the potential to calculate the thermodynamic properties of iron solidification. 2 To whom any correspondence should be addressed.

EAM potential model
The potential model is the basis of molecular dynamic simulation. Currently one of the most widely used potential formats for metallic systems is given by the EAM model proposed by Daw and Baskes [12] on the basis of the density functional theory (DFT). In this study, we will use an EAM potential recently developed by Ackland et al [13] to calculate the thermodynamic properties of iron solidification. In this model, the total energy of iron system is expressed as where is the embedding energy which includes many-body effects.
i is the electronic density at atom i due to the surrounding atoms: where H(x) is the Heaviside step function.
The second term of the EAM, ), ( ij ij r describes a two-body interaction potential was chosen as [13]  . In the above equations, the B i coefficients are determined by continuity of the potential and its derivative. A k and a k are fitted by using ab initio calculation and experimental data [13]. The quality of this EAM potential in describing the atomic interactions of Fe was checked by simulating radial distribution functions (RDF) and density of the system. The NPT (constant number, pressure, and temperature ensemble) algorithm in combination with periodic boundary conditions was applied in the process of simulations. The pressure is set to zero. The time step is 5fs.
The RDF is defined in equation (6) rN r Here ) , ( r r r n i is the average number of atoms surrounding the ith atom in a spherical shell between r and r r . N is the total number of atoms involved in the system and V is the volume of the system. The calculated RDFs of liquid iron based on different potentials at 1820 K are shown in figure 1 together with experimental data [9]. As indicated in the figure, the RDF predicted by current EAM model fits the experimental data best.

Calculation of equilibrium melting temperature
The equilibrium melting temperature is a crucial reference point for studying material processes and properties. The accuracy of calculated T m strongly affects the accuracy of results of other thermodynamic and kinetic parameters. During calculation of the T m , an orthorhombic box subject to the periodical boundary conditions was filled with 43200 atoms. Half of the atoms were fixed and the others were maintained at high temperature of 3000K employing a Nose-Hoover thermostat [17]. The simulation system will contain two solid-liquid interfaces. Keeping the solid atoms fixed, the liquid atoms were subsequently quenched to an estimated melting temperature and equilibrated over 100 ps. After that, the whole system was allowed to relax up to 500 ps with the temperature maintained at the estimated value of T m . The equilibrium melting temperature was then determined from the evolution of potential energy of the solid-liquid coexistence system, and the details will be described in next paragraph.  Figure 3 shows that the potential energies increase with time at temperature above 1790K, indicating that the system is melting; while the potential energies decrease with time at temperature below 1790K, which means the system is freezing. At 1790K, the potential energy of the system remains unchanged which means the solid and liquid atoms coexist in the system. The equilibrium melting temperature was estimated to be about 1 1790 K, which agrees well with the experimental value of 1811K [11]. The difference is only 1.2%. Comparing with previous potential models [6][7][8][9], the EAM potential model used in present study appears to be the most accurate one. Table 1 summarized the calculated T m by using different potential models together with the experimental data. It shall be noted that for temperatures at 1770K and 1760K, there are plateaus in the potential energy curves near the end of simulation, which means all liquid has solidified fully.

enthalpy, specific heat and latent heat
The enthalpies of solid and liquid phases were obtained from homogeneous solid and liquid phases where ) (T H is the enthalpy as a function of temperature T, E the internal energy; P is the pressure, and V is the volume of the system. (11) The latent heat was obtained from the enthalpy difference between solid and liquid phases at the melting temperature, that is, The calculated latent heat is 16.3 kJ/mol. For a comparison, the calculated data using different potentials and the experimentally measured value are summarized in Table 1.
The isobaric specific heat can be obtained by the differential of temperature dependent enthalpy, The calculated isobaric specific heat for the solid phase is   [16]. In this temperature range, our results slightly increase from 32.9J/mol/K to 33.6J/mol/K. The experimental value from Beutl, Pottlacher, and Jager's work is about 46.2J/mol/K which is measured under a pressure from atmospheric pressure up to 3800 bar [15]. Trevorton and Margrave reported a value of 42.7J/mol/K [18]. In contrast, our result is a little smaller than the reported value. The simulation provides a much wider temperature range for the enthalpy and specific heat. This can be used in the theoretical analysis of solidification and metastable liquid structure.

Solid-liquid interfacial free energy
The solid-liquid interfacial energy was calculated from Gibbs-Thomson coefficient in relation between critical nucleus sizes and the inverse of the critical undercooling temperature by treating spherical crystal [10,19]. The relation is defined by where is the Gibbs-Thomson coefficient; K is the curvature, and 1 r and 2 r are the principal radii of curvature. In the case of the curvature of a sphere, the Gibbs-Thomson equation (16) can be rewritten as where * r is the critical radius of a spherical nucleus. Therefore, the solid-liquid interfacial free energy can be calculated by m V SL T L (18) In order to obtain the critical nucleus size r * , a series of spherical crystal clusters were generated from a large BCC block of Fe, using various spherical cut-off radii r ranging from 1nm to 5nm, and embedded into a liquid matrix. Then the system was quenched to a desired temperature and relaxed for 500ps to investigate the structure evolution of the embedded solid sphere. According to classical nucleation theory (CNT), if the radius of the solid sphere r is smaller than the critical radius r * at a certain undercooling, the solid sphere dissolves in the liquid; if r > r * , the solid sphere grows; if r = r * , the solid sphere and the liquid can coexist [19]. Figure 5 shows the size dependence of structure evolution of a spherical BCC crystal Fe embedded in the undercooled liquid at T=1600K. The total number of atoms in the system is 250000. It is found that when the radius of the solid sphere is 1.89nm, the spherical crystal grows in the undercooled liquid until the whole system becomes solid. However, when the radius is 1.75nm, the sphere melts and dissolves into the liquid surroundings. So the critical radius of the spherical nucleus can approximately be taken as 1.82nm at 1600K. The error in estimating the critical radius is within The critical radii r * as a function of temperature are presented in figure 6. It clearly indicates that the critical radius has an inverse relationship with the undercooling of the system. A linear fitting of the calculated data between T and 1/r* produces * / 3542 r T as indicated in figure 6(b). The Gibbs-Thomson coefficient for spherical Fe nucleus defined in equation (16) was found to be about mK 10 8 . 1 7 , which agrees with the experimental value mK 10 9 . 1 7 of iron measured by Turnbull and previous results calculated for the Finnis-Sinclair (FS) potential [10]. Different sizes of simulation systems containing 54000, 128000 and 250000 atoms were used to test the size effect on our final results. The three systems yield similar results as shown in figure 6(a). The system size has no obvious effect on the critical radius, which has also been confirmed in previous research work [10,19]. The solid-liquid interfacial energy at melting temperature was calculated from equation (18). Our simulations give a value of 214mJ/m 2 , which is about 4.9% higher than the experimental value 204mJ/m 2 . The FS potential used in reference 10 generates a value of 170mJ/m 2 based on the same simulation method. It is about 16.7% lower than the experimental one [10], which can contribute to a relatively high overestimation (about 32.5%) of the melting temperature (2400K vs 1811K). A summary of spherically averaged values of solid-liquid interfacial energy at T m calculated by different potentials and methods is presented in table 2. The ABCH potential predicts the closest value to the experimental data by using capillary fluctuation method. However, considering the value was obtained at the melting temperature 2328.7K, which is 28.6% higher than experimental one, we think the value calculated by using current potential should be more reasonable.

Conclusions
Molecular dynamics simulations (MD) based on the embedded atom method (EAM) potential were used to calculate the density, enthalpy, specific heat and solid-liquid interfacial free energy. The calculated results show reasonable agreement with the experimental data. It indicates that the present EAM potential model developed by Ackland et al. appears to be the most accurate model for MD simulation of iron solidification. Simulations with the above EAM model predict the equilibrium melting temperature of iron is 1790K, the solid-liquid interfacial energy 214 mJ/m2. The difference with the experimental data is 1.2%, and 4.9% respectively.