Structural origin underlying the effect of cooling rate on solidification point*
Li Chen-Hui, Han Xiu-Jun, Luan Ying-Wei, Li Jian-Guo
Laboratory of Advanced Materials Solidification, School of Materials Science and Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

Corresponding author. E-mail:

*Project supported by the National Basic Research Program of China (Grant No. 2011CB012900), the National Natural Science Foundation of China (Grant No. 51171115), the Natural Science Foundation of Shanghai City, China (Grant No. 10ZR1415700), the Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20100073120008), the Program for New Century Excellent Talents in Universities of China. This work is partially supported by Alexander von Humboldt Foundation.


Solidification behaviors of liquid aluminum at different cooling rates were examined via classical molecular dynamics simulation with an embedded atom method potential. The results demonstrate that solidification point decreases with increasing cooling rate. To explain this phenomenon, solid-like cluster in liquid was analyzed by the structural analysis method of bond order parameters. The results reveal that the size of the largest solid-like cluster in deeply undercooled liquid decreases with the increase of cooling rate, which can provide a structural interpretation to the above phenomenon.

Keyword: 61.25.Mv; 64.70.D–; 61.20.Ja; liquid metal; undercooled solidification; local structure; molecular dynamics
1. Introduction

Solidification of liquid metal is a fundamental issue in material sciences since almost all kinds of metals we used have to undergo such a procedure at least once or more.[1] In experiments, solidification process is always influenced by impurity, even advanced levitation techniques are employed, [2, 3] which hinders our further understanding of the fundamental principle of solidification. To overcome the experimental problem, molecular dynamics (MD) simulation has been widely used since the problem of impurity does not exist. MD simulation is a powerful tool to explore the phase transformation and structure of condensed matter at the level of atomic resolution.[4] In MD simulation, aluminum (Al) is often selected as a model metal for investigation[58] due to the wide application of Al as a good engineering material. For example, Shen et al.[9] recently presented different solidification processes of liquid Al with the respective cooling rate, and they found that the high cooling rate leads to glass formation, while the low cooling rate leads to crystallization.

As we know, the cooling rate is an important factor during a solidification process, and has a profound effect on the final structures.[10, 11] However, the dependence of solidification behavior on cooling rate is still poorly understood. Shibuta et al.[12] investigated the effect of cooling rate on the solidification behavior of molybdenum nanoparticle, and they observed that the solidification point decreases with increasing cooling rate. A similar observation was found in bulk nanocrystalline Al.[13] These studies provide us with a deeper knowledge of solidification behavior, but one question coming is why the solidification point decreases with increasing cooling rate. It is well known that, solidification process stems from nucleation except glass transition. Kelton et al.[14] demonstrated that a correlation exists between nucleation barrier and icosahedral short-range order, which is a kind of liquid local structure. Inspired by their work, we think that the decrease of solidification point with increasing cooling rate might be associated with liquid local structure. Recently, liquid local structure has attracted great attention both in experiments[1518] and in MD simulations, [19, 20] showing more and more importance in material science. Therefore, in this work, we try to provide an explanation to the above question from the view of liquid local structure.

Firstly, the solidification of liquid Al is conducted using MD simulation with different cooling rate. The cooling rates we selected differ by an order of magnitude, allowing us to clearly observe the discrepancy caused by cooling rate. Secondly, solidification point at each cooling rate is summarized. At last, we present the result of liquid local structure. Based on the result and the classical nucleation theory, [1] the correlation between cooling rate and solidification point is well interpreted. Much attention has been directed to the study of the equilibrium liquid structure, while the liquid structure during a non-equilibrium process like solidification has been rarely reported. We believe this work could attract more attention to further explore the liquid structure at non-equilibrium states.

2. Simulation details

All simulations were performed under external constant pressure of zero using the open source code LAMMPS.[21] The Finnis– Sinclair potential was adopted to describe the atomic interactions. The original file of the potential for Al is provided by Mendelev.[22] The potential has already been applied to predict the solidification and devitrification of Al.[23] Within the framework of embedded atom method, the total energy Etot consists of two terms

where i, j denote different atoms, N the number of atoms, ri j the distance between atoms i and j, φ i j(ri j) the pair potential between atom i and j at a distance of ri j, Fi(ρ i) the embedded energy, and ρ i the background charge density at the site of atom i

with ψ i(ri j) being the contribution to the charge density at atom i due to atom j.

The system size (8.1 × 8.1 × 8.1  nm with 32000 atoms) we used in simulation is large enough to observe the liquid-solid transition and exhibit negligible size effect of results. During simulations, periodical boundary conditions (PBCs) are used in three directions. Isothermal-isobaric ensemble (NPT) was used and controlled by the Nosé – Hoover method.[24, 25] A time step of 1  fs was applied with Verlet algorithm[26] to integrate Newton’ s equation of motions. The solidification of liquid Al was performed as follows. Firstly, the whole system was heated up to 2050  K, much higher than the experimental melting point (933.47  K), and held there for 600000 steps to make sure that the system is in complete liquid state. Then, the system was cooled linearly in time from 2050 to 50  K at the cooling rate ranging from 2.0 × 1010 to 2.0 × 1013  K/s, i.e., T(t) = T0Qt, where Q is the cooling rate. Each cooling process was repeated 3 times. At each desired temperature, the configuration was recorded for further analysis. Table  1 lists the cooling rates and the corresponding numbers of steps.

Table 1. Cooling rates and correponding numbers of steps.
3. Results and discussion

Figure  1 shows the evolution of the total energy per atom (E) with decreasing temperature at the four cooling rates. At Q1, Q2, and Q3, there is an abrupt change of the linear temperature dependence of E, which means liquid Al crystallizes during the cooling. It is observed that, within the transition region, E drops sharply at Q1 and Q2, while smoothly at Q3. This means that the nulceation resistance becomes significantly higher at Q3 than that at Q1 and Q2. For Q1, Q2, and Q3, we define the solidification point as the temperature where the linear decrease in ET curve disappears. It is seen that the solidification points at the three cooling rates are all far lower than the simulated melting temperature (Tm), Tm = 912.5  K. The simulation details of Tm are illustrated in Appendix A. At Q4, E decreases continuously within the whole temperature range, and no abrupt change can be observed, which is a typical feature of glass transition. The glass transition temperature (Tg) is 350  K, which is determined by the Wendt– Abraham parameter[27] (seen Appendix B). For Q4, Tg is considered as the solidification point. From Fig.  1, we can conclude that the solidification point decreases with increasing cooling rate, which is consistent with the result of Ref.  [12].

Fig.  1. ET curves under various cooling rates. The inset is the summarized solidification point.

There is now a general consensus that liquid metals are composed of various kinds of local structure.[1520] One important kind of local structure is the solid-like cluster which is generated by liquid structure fluctuation. Solid-like cluster is directly related to the solidification process because one proper solid-like cluster in liquid, especially in deeply undercooled liquid, is enough to initiate nucleation and induce crystallization. In order to further examine the potential mechanism underlying the correlation between cooling rate and solidification point, we have calculated the maximum size of solid-like cluster at different temperature in liquid using the recorded configurations. To do this, we determine the solid-like cluster by using the method of bond-order parameters.[28, 29] The advantage of this method is that it can efficiently distinguish solid-like and liquid-like atoms and has little relation to specific type of crystal. This method has been successfully applied in previous studies by other researchers.[30, 31] For each atom i, the normalized order parameter is defined as

where A is a normalization factor, B(i) the number of neighbors of atom i, Ylm a spherical harmonic function, and i j the unit vector pointing from atom i to j. If the complex inner product is larger than 0.5, the neighboring atoms i and j are deemed to be ordered neighbors. On the basis of predecessor’ s work and the usual convention, atoms with more than 8 ordered neighbors (including 8) are identified as solid-like.[32] If an atom has less than 8 ordered neighbors, it is considered to be liquid-like. The validity of the criterion is demonstrated in Ref.  [33]. Using this criterion, solid-like cluster formed by adjacent solid-like atoms can be well located. The size of the largest solid-like cluster, NMAX, as a function of temperature is presented in Fig.  2. The temperature range we focus on is from 650 to 1050  K, because no solid-like cluster is found at a temperature higher than 1050  K, while the impending nucleation is detrimental to the reasonability of the information we care about at temperatures lower than 650  K. The first minimum of pair distribution function g(r) is used as a cutoff distance for the determination of solid-like cluster. As shown in Fig.  2, solid-like clusters have already existed in superheated liquid. Solid-like cluster first appears under Q1 as temperature decreases. Obviously, the probability of finding solid-like clusters in undercooled liquid is much larger than that in superheated liquid. NMAX shows a general trend of increasing with the decrease of temperature, and increases rapidly when the temperature is close to crystallization. The change of NMAX with temperature can be approximated by an exponential expression, and the fitted results are presented in the small inserted figure. It is observed that, at Q1, Q2, and Q3, NMAX grows with decreasing temperature by orders of magnitude over a narrow temperature range near crystallization, and the growth rate becomes higher if the liquid is cooled with a slower cooling rate. At 650  K, NMAX in the undercooled liquid has a clear difference among the four cooling rates, and the slower cooling rate produces the larger NMAX. This observation may uncover why a slower cooling rate causes a higher solidification point. According to the classical nucleation theory, the Gibbs free energy difference Δ G to form a spherical nucleus in undercooled liquid consists of a volume term and an interfacial term,

where r is the nucleus radius, σ sl the solid/liquid interfacial energy, and Δ GV the volumetric free energy change. The maximum Δ G marked as Δ G* is received at a critical radius r*

If the size of solid-like cluster is smaller than r* , the cluster will dissolve. Otherwise, it will grow further and induce crystallization. From Fig.  2, it is most probable for Q1 to firstly form a solid-like cluster with a radius larger than r* . That is to say, it is easier for the liquid to firstly crystallize at Q1 than at the other three cooling rates, and consequently the solidification point at Q1 will be higher than that at Q2 and Q3. The similar discussion can be applied to Q2 and Q3, which can well explain why the solidification point decreases with increasing cooling rate. If the cooling rate (like Q4) is high enough to exceed the critical glass transition cooling rate, the solidification process is a glass transition rather than crystallization. Then the solidification point should be Tg which is considerably lower than the crystallization temperature. Because the solidification mechanism has changed, here we limit our discussion to the crystallization.

Fig.  2. The maximum size of solid-like cluster as a function of temperature. The inset is the fitting line. The size of solid-like cluster is characterized by the number of atoms it contains.

4. Conclusions

In summary, a series of MD simulations were performed to examine the cooling rate effect on the solidification behavior of liquid Al. The observed solidification points are all far lower than the simulated melting point, which indicates that the solidifications at the four cooling rates need a quite large undercooling. At the cooling rate less than or equal to 2.0 × 1012  K/s, crystallization takes place. At the cooling rate of 2.0 × 1013  K/s, glass transition is obtained and the corresponding glass transition temperature is 350  K. Following the increase of cooling rate, there is a clear drop of the solidification point. In order to elucidate that, solid-like clusters in undercooled and superheated liquid were analyzed. The results suggest that, within the temperature range near crystallization, the maximum size of solid-like cluster diminishes remarkably if the cooling rate increases a lot. This observation combined with the classical nucleation theory can reveal the structural origin underlying the effect of cooling rate on solidification point.

1 Kurz W and Fisher D J 1986 Fundamentals of Solidification1984 edition Switzerland Trans. Tech. Publications LTD [Cited within:2]
2 Notthoff C, Feuerbacher B, Franz H, Herlach D M and Holland -Moritz D 2001 Phys. Rev. Lett. 86 1038 DOI:10.1103/PhysRevLett.86.1038 [Cited within:1]
3 Willnecker R, Herlach D M and Feuerbacher B 1986 Appl. Phys. Lett. 49 1339 DOI:10.1063/1.97371 [Cited within:1]
4 Rapaport D C 2004 The Art of Molecular Dynamics Simulation2004 edition Cambridge Cambridge University Press [Cited within:1]
5 Kramer M J, Mendelev M I and Asta M 2014 Philos. Mag. 94 1876 DOI:10.1080/14786435.2014.886786 [Cited within:1]
6 Li F, Liu X J, Hou H Y, Chen G and Chen G L 2011 J. Appl. Phys. 110 013519 DOI:10.1063/1.3605510 [Cited within:1]
7 Mendelev M I, Schmalian J, Wang C Z, Morris J R and Ho K M 2006 Phys. Rev. B 74 104206 DOI:10.1103/PhysRevB.74.104206 [Cited within:1]
8 Li M Z, Wang C Z, Mendelev M I and Ho K M 2008 Phys. Rev. B 77 184202 DOI:10.1103/PhysRevB.77.184202 [Cited within:1]
9 Shen B, Liu C Y, Jia Y, Yue G Q, Ke F S, Zhao H B, Chen L Y, Wang S Y, Wang C Z and Ho K M 2014 J. Non-cryst. Solids 383 13 DOI:10.1016/j.jnoncrysol.2013.05.004 [Cited within:1]
10 Xia J C, Zhu Z G and Liu C S 1999 Chin. Phys. Lett. 16 850 [Cited within:1]
11 Liu C S, Zhu Z G, Xia J C and Sun D Y 2000 Chin. Phys. Lett. 17 34 [Cited within:1]
12 Shibuta Y and Suzuki T 2011 Chem. Phys. Lett. 502 82 DOI:10.1016/j.cplett.2010.12.020 [Cited within:2]
13 Hou Z Y, Tian Z, Mo Y F and Liu R S 2014 Comp. Mater. Sci. 92 199 DOI:10.1016/j.commatsci.2014.05.044 [Cited within:1]
14 Kelton K F, Lee G W, Gangopadhyay A K, Hyers R W, Rathz T J, Rogers J R, Robinson M B and Robinson D S 2003 Phys. Rev. Lett. 90 195504 DOI:10.1103/PhysRevLett.90.195504 [Cited within:1]
15 Mauro N A, Bendert J C, Vogt A J, Gewin J M and Kelton K F 2011 J. Chem. Phys. 135 044502 DOI:10.1063/1.3609925 [Cited within:2]
16 Hirata A, Guan P F, Fujita T, Hirotsu Y, Inoue A, Yavari A R, Sakurai T and Chen M W 2011 Nat. Mater. 10 28 DOI:10.1038/nmat2897 [Cited within:1]
17 Di Cicco A, Iesari F, De Panfili S, Celino M, Giusepponi S and Filipponi A 2014 Phys. Rev. B 89 060102 DOI:10.1103/PhysRevB.89.060102 [Cited within:1]
18 Mauro N A, Wessels V, Bendert J C, Klein S, Gangopadhyay A K, Kramer M J, Hao S G, Rustan G E, Kreyssig A, Goldman A I and Kelton K F 2011 Phys. Rev. B 83 184109 DOI:10.1103/PhysRevB.83.184109 [Cited within:1]
19 Ding J, Cheng Y Q and Ma E 2014 Acta Mater. 69 343 DOI:10.1016/j.actamat.2014.02.005 [Cited within:1]
20 Xiong L H, Lou H B, Wang X D, Debela T T, Cao Q P, Zhang D X, Wang S Y, Wang C Z and Jiang J Z 2014 Acta Mater. 68 1 DOI:10.1016/j.actamat.2014.01.004 [Cited within:2]
21 [Cited within:1]
22 Mendelev M I, Kramer M J, Becker C A and Asta M 2008 Philos. Mag. 88 1723 DOI:10.1080/14786430802206482 [Cited within:1]
23 Mendelev M I 2012 Modelling Simul. Mater. Sci. Eng. 20 045014 DOI:10.1088/0965-0393/20/4/045014 [Cited within:1]
24 Nosé S 1984 J. Chem. Phys. 81 511 DOI:10.1063/1.447334 [Cited within:1]
25 Hoover W G 1985 Phys. Rev. A 31 1695 DOI:10.1103/PhysRevA.31.1695 [Cited within:1]
26 Verlet L 1967 Phys. Rev. 159 98 DOI:10.1103/PhysRev.159.98 [Cited within:1]
27 Wendt H R and Abraham F F 1978 Phys. Rev. Lett. 41 1244 DOI:10.1103/PhysRevLett.41.1244 [Cited within:1]
28 Steinhardt P J, Nelson D R and Ronchetti M 1983 Phys. Rev. B 28 784 DOI:10.1103/PhysRevB.28.784 [Cited within:1]
29 ten Wolde P N, Ruiz-Montero M J and Frenkel D 1996 J. Chem. Phys. 104 9932 DOI:10.1063/1.471721 [Cited within:1]
30 Dullens R P A, Aarts D G A L and Kegel W K 2006 Phys. Rev. Lett. 97 228301 DOI:10.1103/PhysRevLett.97.228301 [Cited within:1]
31 Gasser U, Weeks E R, Schofield A, Pusey P N and Weitz D A 2001 Science 292 258 DOI:10.1126/science.1058457 [Cited within:1]
32 Guzmán J H and Weeks E R 2009 Proc. Nati. Acad. Sci. USA 106 15198 DOI:10.1073/pnas.0904682106 [Cited within:1]
33 Lechner W and Dellago C 2008 J. Chem. Phys. 129 114707 DOI:10.1063/1.2977970 [Cited within:1]
34 Han X J, Chen M and Guo Z Y 2004 J. Phys. : Condens. Matter 16 705 DOI:10.1088/0953-8984/16/6/001 [Cited within:1]
35 Han X J and Schober H R 2011 Phys. Rev. B 83 224201 DOI:10.1103/PhysRevB.83.224201 [Cited within:1]