CC BY
65
УДК 538.913
Molecular Dynamics Studies of B2 NiAl Alloy Melting Point
N.T.H. Trung1-2, H.S.MPhuong1-2, M.D. Starostenkov1
'Polzunov Altai State Technical University (Barnaul, Russia) 2Nuclear Research Institute (Da Lat, Vietnam)
Исследование точки плавления интерметаллида NiAl сверхструктуры В2
Н.Ч.Х. Чунг1,2, Х.С.М. Фуонг1,2, М.Д. Старостенков1
'Алтайский государственный технический университет им. И.И. Ползунова (Барнаул, Россия)
2Институт ядерных исследований (Далат, Вьетнам)
In this paper, we perform a computer simulation to investigate the melting point at zero pressure of B2 NiAl intermetallic alloy by using LAMMPS with EAM potential developed by Mishin et al. Simulation box contains 20^20x20 unit cells with 16000 atoms, periodic boundary conditions are applied in three directions. To verify the quality of Mishin potential we first conduct several simulations to calculate defects formation energy, cohesive energy, equilibrium lattice constant and elastic constants of this alloy at absolute zero. The main simulation is performed by using one-phase method in NPT ensemble. Simulation results are analyzed and visualized by Ovito using radical distribution function and common neighbor analysis method. We observe B2 NiAl bulk alloy that begins to melt at 1840 K and crystallizes at 1153 K with critical cooling rate higher than those of almost other alloys. The good agreement between simulation results and experiment suggests that we should continue using Mishin potential for further work in B2 NiAl case study with more sophisticate simulation. Key word: molecular dynamics, lammps, eam potential, one-phase method, melting point, phase transition, critical cooling rate, crystallization, time-temperature-transformation diagram.
Применено компьютерное моделирование для исследования точки плавления при нулевом давлении интерметаллида №А1 со сверхструктурой В2 с использованием Lammps программного комплекса с потенциалом взаимодействия атомов, полученного по методу погруженного атома в работе Т Mishin и др. Рассчитанный блок кристалла содержал 20x20x20 элементарных ячеек, включающих 16000 атомов, к границам расчетного блока прикладывались периодические граничные условия в трех направлениях. Для верификации качества потенциала Т. Mishin и других были выполнены компьютерные эксперименты по вычислению энергий образования точечных дефектов, энергии когезии, равновесного параметра решетки и упругих констант сплава №А1 при абсолютном нуле. Основная процедура моделирования была проведена по методу однофазного PVT ансамбля. Результаты моделирования были проанализированы и визуализированы с использованием программного комплекса Ovito и с функции радиального распределения и анализа общего числа соседей атомов. Исследовался интер-металлид №А1 сверхструктуры В2 в объеме начиная с температуры плавления 1840 К и возникновения кристаллизации при температуре 1153 К с высокой скоростью охлаждения. Хорошее согласие результатов моделирования и экспериментальных данных подтверждает возможность применения потенциала Т Mishin и других для дальнейших работ по исследованию сплава №А1 сверхструктуры В2.
Ключевые слова: молекулярная динамика, 1аттр8, метод потенциала погруженного атома, точка плавления, фазовый переход, критическая скорость охлаждения, диаграмма температурно-временной трансформации.
DOI 10.14258/izvasu(2017)4-01
1. Introduction
NiAl possesses seven different ciystal phase: B2, B1, L10, L1p B32, "40" (NbP prototype) and B20 [1]. Among these structures, B2 (CsCl-type, Pm3m space group) is the most stable one both at ambient and very high pressure at room temperature [2]. The B2 NiAl compound is of particular interest from a mechanical point of view that motivate its extensive use in industry [3]. This class of intermetallic alloys is very promising materials for high temperature and pressure applications. Currently, it is being examined for use in diesel engine turbocharger rotors, high-temperature die and molds, hydro turbines, and cutting tools due to excellent oxidation resistance, high melting temperature, low density (5.9 g/cm3) and high thermal conductivity (4 to 8 times those of Ni3Al superalloys) [4-6].
Calculation of the thermodynamic melting point is not trivial for bulk metal, as has been shown by various study groups [7], in which several methods have been published. Generally, evaluating the melting point is not a straightforward approach because it requires the calculation of Gibbs free energies, which is a multistage process that needs several different molecular dynamics simulations. In case of B2 NiAl alloy many studies have conducted to find its melting point under various pressure, however the results are significant different between these research groups [1, 2, 4]. The discrepancies may come from different method and potential applied. Even in the simplest case when pressure equal zero, there are two different experiment values of B2 NiAl melting point 1911 K and 1995 K were
reported from literature [2, 4, 8]. The authors of this work would like to provide another perspective to this classical problem in thermodynamics.
Many potentials have been developed for Ni-Al system, including the Finnis-Sinclair potential of Yan and the EAM potentials by Foiles and Daw, Voter and Chen. However, in this work, to investigate the melting point for bulk B2 NiAl we chose the embedded-atom-method (EAM) potential developed by Mishin et al [1]. This potential was developed in 2009 using empirical fitting method. The version, which is compatible with LAMMPS, is available for download from http://www.ctcms.nist.gov/ potentials/.
To assess the qualities of this potential, at first we performed several simulations to determinate the defect formation energy, lattice constant, elastic constants and cohesive energy at absolutely zero. The results were given on Table 1, 2.
According to the table 1 the smallest formation energy is NiAl and then VA and VNi. Therefore in NiAl alloy, when Ni is rich, the most possible defect is antisite Ni on Al sublattice. But when Al is rich because the formation energy of AlNi is higher than that of VNi then VNi is most possible defect in the NiAl alloy. These predicts are good agreement with experiment result for defect behavior of Ni-rich and Al-rich NiAl alloy [3].The calculation results from table 1, 2 prove that, this potential accurately reproduces many properties on B2 phase, therefore, it is appropriate to use Mishin potential to perform more challenging molecular dynamic simulations.
Table 1: The defect formation energy from our simulation and other results
Type of defect This work Caro [9] Wang [10]
VAl 1.14 1.07 1.063
VM 1.76 1.63 1.335
Antisite AL. Nl 2.25 2.11 2.347
Antisite NiAl -0.32 -0.25 -0.150
Table 2: Lattice constant, cohesive energy and elastic constants of B2 NiAl
Parameter This work Experiment [1, 6]
Lattice constant a (A0) 2.83 2.88
Cohesive Energy E0(eV) -4.51 -4.50
Cu (GPa) 186 199
C12 (GPa) 132 137
C44 (GPa) 113 116
Bulk Modulus (GPa) 150 157
2. Calculation of melting point
In this paper we were performed a molecular dynamics simulation based on one-phase method [7] to calculate the melting point of B2 NiAl alloy at zero pressure.
This method is empirical, which based on the classical nucleation theory [2]. In this method, there is a sudden change in volume (potential energy) when solid is heating at superheating temperature Tsh far beyond
the melting temperature Tm, and when liquid is cooling at supercooling temperature Tsc far below Tm. The volume-temperature V-T (or potential energy-temperature) curve show the hysteresis effect of heating solid and cooling liquid. The melting point Tm can be estimated from Tsh and T by the following formula [2, 7, 11]
T =T
(1)
The simulation was performed within NPT thermodynamic ensemble using LAMMPS-an open classical molecular dynamics code [12]. In our lab LAMMPS was compiled with MPI and GPU package by NVIDIA CUDA toolkit to accelerate computing efficiency. Simulation box contains 20x20x20 unit cells with 16000 atoms, periodic boundary conditions
are applied in three directions. During simulation, thermodynamic properties such as energy, enthalpy, temperature, pressure, volume... were closely monitored.
The solid phase simulation began from a perfect B2 lattices The first 1 ns (106 timesteps) was using to equilibrate to the temperature of 300 K with Nose-Hoover themostat and barostat method [13] was applied to control temperature and pressure (0 bar). After that the equilibrium system was used by increasing temperature from 300 K to 2500 K far beyond estimate melting point in 3 ns with heating rate ~ 7.33 x 1011 K/s. The next 1 ns was used to archive equilibrium at 2500 K. This approach was used for the equilibrium liquid phase in which the temperature was decreased in 20.5 ns with cooling rate ~ 1.07X1011 K/s from 2500 K to 300 K, far below estimate Tm. The Potential Energy-Temperature curve (figure 1) was plotted from LAMMPS's data file after simulation finished.
Figure 1: Heating-cooling curve of B2 NiAl at zero pressure, the red line and the blue line are heating curve and cooling curve, respectively
From the figure 1 it is apparent that in heating curve there is a sudden change in potential energy at superheating temperature Tsh=2345 K and at supercooling temperature
Tsc =1106 K in cooling curve. Therefore, according to formula [2, 7, 11] the melting point at zero pressure of B2 NiAl is approximately 1840 ± 25 K.
Table 3: Melting point of B2 NiAl at zero pressure
Results This work Zhang [2] Experiment 1 [2] Experiment 2 [8]
Tm (K) 1840 1850 1911 1995
To verify the melting point we were using radical distribution function g(r) (RDF), it is one of the most powerful technique to analyze the structure of material, particularly for liquid. The RDF for B2 NiAl system can be calculated as follow:
gij(r) =
V
NN
"I
i i
(
nij( r -5г/2,г + ôr/2 )
Л
4nr 2ôr
(2)
where V is volume of the system, n^ (r -Sr/2, r + 5r/2) is average number of atom j surrounding atom i in a
spherical shell within r + Sr/2, Ni, N. are total number of atom i and j.
Figure 2 depicts the RDF curves of B2 NiAl from our simulation at 300 K, 1840 K, and 2500 K, respectively. At 300K curve depicting crystalline state with well-defined RDF peak up to 5.5 A0, indicate solid structure. With the temperature increasing up to 1840 K, these peaks begin to vanishing peak height with very broad RDF peaks, indicate liquid state. At 2500 K the second RDF peak is complete vanished means it was completely melted.
Figure 2: RDF curve of NiAl from melting simulation at 300 K, 1840 K, 2500 K
It should mention that we chose the simulation time for cooling process (20.5 ns in 20.5 million time steps) much longer than heating process (3 ns) because we have to achieve the cooling rate lower than critical cooling rate k of B2 NiAl so crystallization process can be happen. Without crystallization, the sudden change of potential energy in cooling curve will no longer exist. The temperature, at which the crystallization process begins to occur can be estimated from heating-cooling diagram as figure 4, or from RDF curve when the height
of second peak starting to increase, indicate ordered solid structure. However, in this article we used the common neighbor analysis method integrated in Ovito [14, 15] because this method gives more detail information about structure type of bulk NiAl during simulation. It was found that the crystallization of disordered liquid NiAl in cooling process is starting at approximately 1153 K. At the supercooling temperature Tsc =1106 K as visualized on figure 3 by Ovito, 52% atoms in simulation box have recovered to B2 ordered from disordered liquid NiAl.
Figure 3: Structure type of simulation box at 1106 K in cooling process
When metallic liquids solidify at the cooling rates higher than critical cooling rate k, nucleation can be suppressed, producing metallic glasses with an amorphous structure. Under such circumstances, the liquid structure with a short-range order remains in a solid amorphous state without long-range translation order. Each material
has its own critical cooling rate k, which can be estimated from the position of the "nose" on the time-temperature-transformation (TTT) diagram [16]. Most alloys require a critical cooling rate of about 106 K/s for the formation of an amorphous structure, while pure-element metals require critical cooling rates as high as 1012 K/s [17, 18].
In our simulation of cooling process, the temperature of simulation cell was decreased from 2500 K to 300 K in 20.5 ns. The default simulation time step is 1 fsmeans the cooling rate in our simulation is approximately
1.07*10" K/s. We concluded that the critical cooling rate of B2 NiAl have to be 105 times bigger than critical cooling rate for most alloys.
Figure 4: Time-Temperature-Transformation diagram for crystallization of the supercooling (undercooled) liquid NiAl alloy [19]
Figure 4 shows the Time Temperature Transformation diagram of NiAl intermetallic in the work of Levchenko et al [19]. This is the first time the TTT diagram for crystallization of stoichiometric intermetallic NiAl has been determined using MD simulation. From this diagram we can be estimate the critical cooling rate for the formation of amorphous NiAl. It was found to be very high, about 1011 - 1012 K.s-1, this value is agree well with the cooling rate from our MD simulation. Unlike almost others alloy, the critical cooling rate of liquid NiAl alloy is comparable with similar estimations for pure metal [18, 19].
Conclusion
The good agreement between the results from our work and experiment, simulation results from other
research group have proven that the Mishin potential well suited for simulation of NiAl behavior in B2 phase. It should be noted that the time-temperature-transformation diagram in the work by Levchenko et al [19] was obtained from molecular dynamic simulation with another version ofMishin EAM potential, which was particularly developed for B2 NiAl in 2002 [20]. However, this potential predicted the melting point of B2 NiAlby using two-phase method at 1520 K, 20% lower than experiment value of 1911 K.Therefore, it is expected to obtain more accurate time-temperature-transformation diagram for B2 NiAl alloy with the same potential using in our work. This problem will be discussed more detail in our future work.
References
1. Purja Pun G.P., Mishin Y. Development of an interatomic potential for the Ni-Al system // Philosophical Magazin. — 2009. — No.89.
2. Zhang W., Peng Yu., Liu Zh. Molecular dynamics simulations of the melting curve of NiAl alloy under pressure // AIP Advances. — 2014. — No.4.
3. Lozovoi A.Y., Mishin Y. Point defects in NiAl: The effect of lattice vibrations // Physical Review B. — 2003. — 68(18).
4. Darolia R., Walston W.S., Nathal M.V. NiAl Alloys for Turbine Airfoils. Superalloys // The Minerals, Metals & Materials Society. — 1996.
5. Wang Y., Liu Z.-K., Chen L.-Q. Thermodynamic properties of Al, Ni, NiAl, and Ni3Al from first-principles calculations // Acta Materialia. — 2004. V.52.
6. Meyer B. and Fahnle M. Atomic defects in the ordered compound B2-NiAl: A combination of ab initio electron theory
and statistical mechanics // Physical Review B. — 1999. — V.59. — No.9.
7. Epiphaniou N. Modeling of dynamic friction across solid material interfaces using molecular dynamics techniques, PhD Thesis, Cranfield University. 2009.
8. Yu Z., Hass D.D., Wadley H.N.G. NiAl bond coats made by a directed vapor deposition approach // Materials Science and Engineering A. — 2005. — V.394.
9. Caro J.A., Pedraza D.F. The stability of irradiation-induced defects in NiAl // Nuclear Instruments and Methods in Physics Research. — 1991. — No.59-60, Part 2.
10. Wang B., Wang Y Calculation of point defects NiAl alloy // Journal of Materials Science and Technology. — 1997. — No.13.
11. Luo S.N., Strachan A., Swift D.C. Nonequilibrium melting and crystallization of a model Lennard-Jones system // Journal of Chemical Physics. — 2004. —120(24).
12. Plimpton S. Fast Parallel Algorithms for Short-Range Molecular Dynamics // Journal of Computational Physics. — 1995. — V.117.
13. Hoover W.G. Canonic dynamics: Equilibrium phasespace distribution // Physical Review A. — 1985. — V.31.
14. Stukowski A. Structure identification methods for atomistic simulations of crystalline materials // Modelling and Simulation in Materials Science and Engineering. — 2012. — 20(4).
15. Stukowski A. Visualization and analysis of atomistic simulation data with OVITO — the Open Visualization Tool // Modelling and Simulation in Materials Science and Engineering. — 2010. — 18(1).
16. Schroers J.Wu., Busch Y., R.& Johnson W. L. Transition from nucleation controlled to growth controlled crystallization in Pd43Ni10Cu27P20 melts // Acta Materialia. — 2001. — V.49.
17. Poletaev G.M, Krasnov V.Yu., Starostenkov M.D, Medvedev N.N., The research of the structure of amorphous metals by molecular dynamics method // Journal of Physics: Conference Series. — 2008. — 98.
18. Han J.J., Wang C.P., Liu X.J., Wang Y., Liu Z.-K., Zhang T.-Y., Jiang J.Z. Abnormal correlation between phase transformation and cooling rate for pure metals // Scientific Reports. — 2016. — 6.
19. Levchenko E.V., Evteev A.V., Belova I.V., Murch G.E. Molecular dynamics determination of the time-temperature-transformation diagram for crystallization of an undercooled liquid Ni50Al50 alloy // Acta Materialia. — 2011. — No.59, Issue 16.
20. Mishin Y., Mehl M.J., Papaconstantopoulos D.A. Embedded-atom potential for B2-NiAl // Physical review B. — 2002. — V.65.
