Molecular dynamics simulation study on nucleation mechanisms of Cu$_3$Au superalloy

(1)

Molecular dynamics simulation study on nucleation mechanisms of Cu

₃

Au superalloy

FATIH AHMET CELIK¹^,∗and EBRU TANBO ˘GA KORKMAZ²

1Faculty of Arts and Sciences, Physics Department, Bitlis Eren University, 13000 Bitlis, Turkey

2Institute of Science, Bitlis Eren University, 13000 Bitlis, Turkey

∗Corresponding author. E-mail: facelik@beu.edu.tr

MS received 23 April 2021; revised 16 August 2021; accepted 7 October 2021

Abstract. The nucleation kinetics of ordered intermetallic alloy of Cu3Au based on L12superlattice are obtained by molecular dynamics simulation utilising Sutton–Chen potential under different pressures. The nucleation mechanism for the system is examined with classical nucleation theory at a slow cooling rate. The crystal-type bonded pairs described by indexes of the Honeycutt–Andersen method are considered as embryos in the system.

The critical nucleus radius, the interfacial free energy and the nucleation rate are determined at given temperatures under 0, 5, 10, 15 and 20 GPa pressures. The structural properties at different pressures of the system are analysed with radial distribution function. As a result, while the critical radius and the interfacial free energies of the system decrease, the nucleation rates increase with increasing pressure.

Keywords. Molecular dynamics; classical nucleation theory; Cu–Au alloys; Honeycutt–Andersen method.

PACS Nos 02.70.Ns; 36.40.–c; 68.55.A–

1. Introduction

Many binary intermetallic alloys display specific physical and mechanical characteristics such as low corrosion resistance, phase equilibria diagrams and good oxida- tion [1–4]. Among them, Cu–Au superalloy is a popular material due to order–disorder phase transition and the character of phase equilibria [5,6]. Recently, a lot of computational investigations have been done on these alloys to improve their different physical properties [7–

12]. Lekkaet al studied the structural and vibrational properties using an effective potential model based on molecular dynamics (MD) simulation [13]. Kart et al reported a comprehensive study on pure Cu, Au metal and Cu–Au intermetallic alloys using Sutton–

Chen (SC) and quantum SC potentials [14]. Although these investigations are important, there are still several shortcomings such as the nucleation mechanisms at the atomic level. Recently, the nucleation and crystallisation processes have been monitored by experimental studies or microscopic methods [4,15–17]. However, there are several experimental difficulties in understanding homogeneous nucleation processes. Hence, there is not much knowledge about the details of the nucleation processes on the microscopic level [18,19].

Alternative computational methods have been developed for understanding the phenomena of homogeneous nucleation mechanisms [20–22]. Molecular dynamics (MD) simulation methods provide an important oppor- tunity to determine the physical properties of materials [4]. Zhu and Chen [18] performed the nucleation of molten potassium bromide (KBr) clusters by using MD simulation. They calculated the nucleation rates and interfacial free energies of sodium bromide clusters comprising different ion pairs [19]. Wolde and Frenkel [23] numerically predicted the crystal nucleation rate and compared it with other simulation and experimental studies. Shibuta et al [24] studied the behaviour of the nucleation process of a monoatomic system at the atomic level. On the other hand, the interatomic interactions are important for reliable results of MD calculations. In this context, the EAM potentials are widely used in the MD simulations for binary and ternary metal- lic systems [25,26].

In this study, we employed the nucleation mechanism by using the MD simulation method utilising the SC-EAM potential function for the Cu3Au superlattice.

We also provided the details of critical nucleus radius, interfacial free energy and nucleation rates under different pressures of ordered Cu3Au alloy with classical

0123456789().: V,-vol

(2)

nucleation theory (CNT) by using MD. We analysed the structural development of the system using radial distribution function (RDF) analysis, the total energies and the melting temperatures at given temperatures under different pressures. The results were compared with other simulation studies in the literature. The effect of pressure on nucleation rate and critical nucleus radius was investigated by using the computational method.

2. Materials and methods

2.1 Simulation details

In our study, the potential function is employed as the interatomic interaction developed by the SC version of EAM [27]. The SC potential consists of two contribu- tions to the total energy (Etot)for the entire system with N atoms. The total energy can be expressed as [28]

i

((F(ρ¯˙I)+1 2

j=i

(φ(r_{i j})), (1)

where Fi is the embedding energy as a function of electron density.φi j is the pair potential and ri j is the distance between atoms. EAM parameters of the Cu3Au system is described in detail in [14]. The probability of finding an atom in a circle of radius ofr around a given centre atom enables to probe the structural evolution of a system [29,30]. This structural development can be defined using the radial distribution function (RDF).

The RDF,g(r), is defined as g(r)= V

N²

˙In_i(r) 4πr²r

, (2)

whereV is the volume of the MD box,Nis the number of atoms andn(r) represents the number of atoms around a central atom within the circle [29].

MD simulation of Cu3Au alloy system is performed using Parrinello–Rahman algorithm [31]. The present simulation is done with periodic boundary conditions by using a time step of 3.32 fs under the NPT ensemble.

The MD box is constructed by an L1₂ supercell with a cubic MD box with 4000 atoms. In our MD simulation, the cubic box is formed by L12superalloy for the Cu3Au system known experimentally [6] as seen in figure1. The atoms are located regularly on the superlattice points as seen in figure1. Firstly, the system waits for 50,000 MD steps (166 ps) at 300 K and then is heated from 300 to 2500 K with increments of 50 K with a total time of 1626 ps. Secondly, the pressure is applied from 0 to 20 GPa (0, 5, 10, 15 and 20 GPa) and the melting points of the system are determined for different pressures. The system under different pressures is equilibrated for 500,000

MD steps (1660 ps) at the determined melting temperatures. Finally, by applying a cooling rate, the liquid systems are cooled from melting points to 300 K with the decrements of 100 K. At each temperature, 166 ps are carried out for equilibrium with a cooling rate of 6.02×10¹¹ K/s. Thus, the total MD times during the cooling are obtained as 1662, 2160, 2493, 2825 and 3324 ps for 0, 5, 10, 15 and 20 GPa pressures, respectively.

2.2 Homogeneous nucleation kinetics

CNT is used as an important tool for understanding homogeneous nucleation kinetics [4]. According to CNT, the critical radius of a nucleus is estimated as r_c= 2σsl

G_v, (3)

whereσslis the interfacial free energy expressed as σsl = k_TH_fus

(V²NA)¹^/³, (4)

whereNAis the Avogadro number,V represents molar volume andHfusdefines the heat of fusion of the system. When kT = 0.32 the interfacial free energy is calculated [18,19]. In line with CNT, the nucleation rate (J) can be written as

ln Nn

N0

= −J V_c(t−t₀) , (5) whereNn is the number of unfrozen clusters,N0is the number of clusters,t₀is the time lag,tis the nucleation time andVcis the effective volume of the cluster [18,19].

3. Results and discussion

Figure2shows the total energy during the heating process. In general, the total energy of a system increases as the temperature increases and then displays a sharp increase at a given temperature value which is the melting point of the system. As seen from the figure, the melting temperature is about 1100 ±10 K at 0 GPa pressure. This value is compatible with the experimental value of 1250 K [6,14] and is also in good agreement with the reported simulation results [14]. Besides, the effect of pressure on melting temperature is seen from the figure. The melting temperature of the system increases with increasing pressure. The values of melting temperature are obtained as 1450±10, 1700±10, 1900±10 and 2050K±10 for 5, 10, 15 and 20 GPa respectively. In brief, the effect of pressure on melting point is observed during the melting.

Figure3shows the total energy of the system during the cooling process under different pressures. As shown

(3)

Figure 1. The initial structure of the MD simulation box based on the L12superlattice.

200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 T(K)

-340 -330 -320 -310 -300 -290 -280 -270

Etot (kJ/mol)

0 GPa 5 GPa 10 GPa 15 GPa 20 GPa

Figure 2. The total energy of the system vs. temperature for different values of pressure during heating.

in figure3, the energy of the system is dependent on the pressure. The energy during the cooling process displays sharp decreases, indicating that the system transforms to the crystal phase. It can be seen from the figure that the crystallisation begins at higher temperatures with the pressure increment.

The RDF is a useful tool to describe the structural changes of a system [32]. Figure4shows the RDF curves at 0, 5, 10, 15 and 20 GPa pressures at 300 K with a slow cooling rate at different MD times. One can see that the RDF peaks become intense and new crystal peaks emerge at farther atomic distances. At the same time, the height of RDF peaks increases as the pressure increases.

While the peaks of RDF curves are very few sharp peaks at 0 GPa, these curve peaks are sharper at 20 GPa. The

200 400 600 800 1000 1200 1400

T(K) -330

-320 -310 -300 -290

Etot(kJ/mol)

Figure 3. The total energy of the system vs. temperature for different values of pressure during cooling.

final structure is close to the ideal fcc crystal lattice at 20 GPa.

The RDF analysis enables the estimation of structural development during the phase transformation. How- ever, this method is not favoured for predicting the behaviour of atomic clusters in the crystal or amorphous matrix at an atomic level [33]. Honeycutt–Andersen (HA) proposed the common neighbour analysis technique (CNAT) to characterise the local order of an atomic cluster [34]. Many works have already employed this technique [30,34–36]. In the CNAT, the local order of a cluster forming in the crystal and amorphous matrix represents a set of four indices (ijkl). Ifaandbatoms are

(4)

2 3 4 5 6 7 r (Å)

0 1 2 3 4 5

g(r)

2 3 4 5 6 7

r (Å) 0

1 2 3 4 5

g(r)

2 3 4 5 6 7

r (Å) 0

1 2 3 4 5 6

g(r)

2 3 4 5 6 7

r (Å) 0

1 2 3 4 5 6 7 8

g(r)

10 GPa

15 GPa 20 GPa

0 GPaand 5 GPa

(a) (b)

(c) (d)

(1662and 2160 ps) (2493 ps)

(2825 ps) (3324 ps)

Figure 4. The variation of RDF curves of the system at 300 K temperature and different MD times with slow cooling rate:

(a) 0 and 5 GPa, (b) 10 GPa, (c) 15 Gpa and (d) 20 GPa (red line represents the curve at 0 GPa pressure).

close neighbours, called bonded pairs,iequals 1.jindi- cates the number of nearest-neighbour atoms (common neighbours) bonded with a and b atoms. k represents the number of bonds among the common neighbours under al is the distinguishing parameter. The crystal- type bonded pairs are 1421 (fcc), 1422 (hcp), 1441 and 1661 (bcc). For instance, the 1421 bonded pairs form an ideal fcc structure, namely, there are 100% 1421 bonded pairs in the fcc structure. 1551, 1541, 1431 and 1231 bonded pairs generally occur in the disorder structures [30,36]. Some bonded pairs are shown in figure5[30].

The total number of crystal-type bonded pairs at 300 K under a slow cooling rate is illustrated in figure6. The number of crystal-type bonded pairs increases from 28%

to 80% when the pressure increases from 0 to 5 GPa at 300 K. The number of crystal-type bonded pairs increases linearly with increasing pressure and reaches 98% under 20 GPa pressure at 300 K. We can infer that the 1421 bonded pairs are dominant, and the final solidi- fied structure is an almost fully perfect fcc crystal. These

Figure 5. The systematic drawing of the related different bonded pairs [30].

results show good agreement in RDF curves and explain the reason for stronger peaks on RDF curves at 20 GPa as seen in figure4.

(5)

0 5 10 15 20 25 Pressure (GPa)

20 40 60 80 100 120

Thetotalnumberofcrystal-typebondedpairs(%) ^{300 K}

Figure 6. The total number of crystal-type bonded pairs at 300 K under different values of pressure.

500 600 700 800 900 1000 1100 1200

T(K) 0

1 2 3 4 5 6 7

rc(Å)

Figure 7. The values of critical nucleus radius during the solidification at different values of pressures.

The understanding of the nucleation mechanism of materials is very important to explain the phase transition phenomenon. Hence various works on nucleation kinetics have been employed to examine the phase transition process in previous studies [18,19]. Figure 7 shows the radius of the critical nucleus in the temperature range under different pressures. The value of the critical radius decreases with decreasing temperature.

0 200 400 600 800

t (ps) -2.4

-2.2 -2 -1.8 -1.6 -1.4

ln(Ncry/N0)

Figure 8. The plots ln (Ncry/N0) vs. time at a range of temperatures under different pressures.

Table 1. The nucleation radius, the interfacial free energies and nucleation rates of the system with different pressures (Hfusis obtained from figure2).

Pressure (GPa)

Temperature (K) Interval of temperature (K)

300σsl(mJ s⁻²) 300 rc(Å)

J(×10⁻³⁵,m⁻³s⁻¹)

0 25.00 1.06 7.67 (1200–800 K)

5 11.21 1.09 11.75 (1500–1100 K)

10 10.05 1.11 12.02 (1700–1300 K)

15 9.32 1.13 12.44 (1900–1500 K)

20 9.01 1.15 16.27 (2200–1800 K)

The effect of pressure on the critical radius is seen clearly from the figure. High pressure causes a reduction of the critical radius at high temperatures. The interfacial free energies (Hfus) of the system at 300 K under different pressures are listed in table1(H_fus is obtained from the range of energy discontinuity as seen in figure2). As pressure increases, the interfacial free energy decreases.

These results also reveal that by coming together, many atoms easily overcome an energy barrier under the effect of pressure. The changes in interfacial energy contribute to the free energy across two existing phases in the solid/liquid interface. The early works declared that interfacial free energy tends to increase with increasing temperature [18,19]. Our results are consistent with previous studies and the homogeneous nucleation theory [37].

(6)

Figure 9. Schematic illustration of atomic distributions in the MD box for system (a) at 300 K and 0 GPa (blue dots represent crystal atoms and red dots represent unknown structures), (b) at 300 K and 20 GPa (blue dots represent Cu atoms and red dots represent Au atoms).

The estimation of nucleation rate is a central subject to understand the crystal growth phenomena of a system [4,18,19]. Under different pressures, the nucleation rates obtained by MD calculations according to eq. (5) are plotted in figure8 and summarised in table1. The nucleation rates can be derived from the slopes of the curves in the figure at a range of temperatures as seen from table1. Here,N_cryrepresents the total number of crystal-type bonded pairs andN0is the number of total clusters at a given temperature. Table 1 clearly shows that high pressure leads to fast nucleation. Therefore, it is concluded that the pressure attributes to the ear- lier crystallisation formation with increasing nucleation rate.

We present snapshots of a deeper investigation of the effect of pressure on the structure in a final configura- tion for the system by using the OVITO visualisation program [38]. Figure 9a illustrates the MD box at 300 and 0 GPa, where the crystal lattices occur, and unknown structures or non-crystalline atomic groups are available in the crystal phase. It is found that the number of fcc atoms does not reach the maximum value. Fig- ure 9b illustrates the MD box at 300 K and 20 GPa.

The unknown structures are almost absent in the crystal phase and the fcc-type crystal clusters are dominant in the system. Hence, the 1421 bonded pairs form the cubic fcc lattice and the crystal planes appear clearly in the MD box.

(7)

4. Conclusions

In this study, we have presented the nucleation processes based on CNT of Cu3Au superalloy by molecular dynamic simulation under different pressures. The potential energy is used to Sutton–Chen version of EAM. The nucleation kinetics including critical radius, interfacial free energies and nucleation rates are calculated by the computational method. The pressure dependence of nucleation kinetics for the model system is important to understand the nucleation. Moreover, the increment of pressure facilitates the tendency of crystallisation during the slow cooling process. On the other hand, the structural analysis with RDF and the pres- ence of crystal-type bonded pairs in the system show that the number of perfect fcc unit cells increases with increasing pressure. Especially, the increase in nucleation rate with applied pressure affects the completion of crystallisation. We emphasise that the determination of nucleation parameters using MD simulation provides reliable results compared to the experimental and computational studies.

References

[1] D Hu and R Kovacevic,Int. J. Mach. Tools Manuf.43, 51 (2003)

[2] J Campbell, Complete casting handbook: Metal cast- ing processes, metallurgy, techniques and design (Butterworth-Heinemann, 2015)

[3] D Uhlenhaut, J Kradolfer, W Püttgen, J Löffler and P Uggowitzer,Acta Mater.54, 2727 (2006)

[4] A Mahata, M A Zaeem and M I Baskes,Model. Simul.

Mater. Sci. Eng.26, 025007 (2018)

[5] Y Zhang, F Luo, Q Bai, C Zhang, B Yu and Z Zhang,J.

Phys. Chem. Solids150, 109879 (2021)

[6] P P Fedorov and S N Volkov,Russ. J. Inorg. Chem.61, 772 (2016)

[7] I Marković, V Grekulović, M R Vujasinović and S Mladenović,J. Alloys Compd.831, 154726 (2020) [8] A Y Tse, E K Karasz and K Sieradzki,Scr. Mater.176,

112 (2020)

[9] Z Wang, R Liu, Z Feng, L Lin, R Xie, J Li, H Liu, F Huang and Z Zheng,Cryst. Eng. Commun.22, 113 (2020)

[10] Q K. Li, Y Zhang and W Y Chu,Comput. Mater. Sci.

25, 510 (2002)

[11] N Artrith and A M Kolpak,Comput. Mater. Sci.110, 20 (2015)

[12] G Chen, C Wang and P Zhang,Comput. Mater. Sci.112, 80 (2016)

[13] C E Lekka, N I Papanicolaou and G A Evangelakis,Surf.

Sci.479, 287 (2001)

[14] H H Kart, M Tomak and T Ça˘gin,Model. Simul. Mater.

Sci. Eng.13, 657 (2005)

[15] J W P Schmelzer,J. Non-Cryst. Solids356, 2901 (2010) [16] H Vehkamäki, (Springer Science & Business Media,

2006)

[17] S Sen and T Mukerji, J. Non-Cryst. Solids 246, 229 (1999)

[18] X Zhu and K Chen, J. Phys. Chem. Solids 66, 1732 (2005)

[19] X Zhu,J. Mol. Struct. Theochem.680, 137 (2004) [20] S Sengul, M Celtek and U Domekeli,Comput. Mater.

Sci.172, 109327 (2020)

[21] F A Celik and A K Yildiz,J. Non-Cryst. Solids415, 36 (2015)

[22] Y Wang and C Dellago, J. Phys. Chem. B 107, 9214 (2003)

[23] P R ten Wolde and D Frenkel,Phys. Chem. Chem. Phys.

1, 2191 (1999)

[24] Y Shibuta, S Sakane, T Takaki and M Ohno,Acta Mater.

105, 328 (2016)

[25] L Hui, B Xiufang and Z Jingxiang,Mater. Sci. Eng. A 271, 116 (1999)

[26] S Ozgen and E Duruk,Mater. Lett.58, 1071 (2004) [27] M S Daw and M I Baskes,Phys. Rev. B29, 6443 (1984) [28] A P Sutton and J Chen,Philos. Mag. Lett.61, 139 (1990) [29] H Tanaka,J. Non-Cryst. Solids351, 678 (2005) [30] Y Q Yuan, X G Zeng, H Y Chen, A L Yao and Y F Hu,

J. Korean Phys. Soc.62, 1645 (2013)

[31] M Parrinello and A Rahman,Phys. Rev. Lett.45, 1196 (1980)

[32] R Liu, K Dong, J Li, A Yu and R Zou,J. Non-Cryst.

Solids351, 612 (2005)

[33] B Yu, Y Liang, Z Tian, R Liu, T Gao, Q Xie and Y Mo, J. Crystal Growth535, 125533 (2020)

[34] J D Honeycutt and H C Andersen,J. Phys. Chem.91, 4950 (1987)

[35] Z A Tian, R S Liu, H R Liu, C X Zheng, Z Y Hou and P Peng,J. Non-Cryst. Solids354, 3705 (2008)

[36] L Qi, S Feng, N Xu, M Ma, Q Jing, G Li, R Liu, L Qi, S Feng, N Xu, M Ma, Q Jing, G Li and R Liu,Mater.

Res.18, 78 (2015)

[37] W D Callister and D G Rethwisch,Materials science and engineering: An introduction(John Wiley & Sons, Inc., 2007) p. 311

[38] A Stukowski,Model. Simul. Mater. Sci. Eng.18, 015012 (2009)