Analysis of Pd-Ni Nanobelts Melting Process Using Molecular Dynamics Simulation

Abstract

The melting process of Pd-Ni alloy nanobelts with different Ni atom content has been simulated by molecular dynamic (MD) method. The radial distribution function, the Lindemann index, and pair analysis method were used to characterize Pd-Ni nanobelt models in simulation. The results indicate that the melting temperature of Pd-Ni nanobelt with composition far from pure metal was lower than that of other models, and the breaking point of the nanobelt can be illustrated by the Lindemann index. Pair analysis indicates that the number of FCC pairs will decrease and almost disappear at melting point with increasing temperature. The melting points of Pd-Ni alloy nanobelts were also calculated by thermodynamic method, and the results were close to that obtained by MD simulation.

1. Introduction

After carbon nanotube has been reported in 1991 [1], one-dimensional (1D) materials have received much attention due to their excellent physical property and have been expected to be used in a wide range of application fields [2–6]. Because of their dramatically different behavior from bulk materials, much research work has been done for the structure and properties of nanowires [7–10]. Thermal stability and melting behavior of nanomaterial are another subject that has been studied theoretically and experimentally for a rather long time. The melting process of nanomaterial determined by experiment is very difficult. Nevertheless, computer simulations provide an excellent method to study these melting processes at an atomic level. Up to date, the melting process and thermal properties of many kinds of 1D metal nanomaterial, including Au [8], ZnSe [11], GaN [12], Zr [13], Ag [14], Cu [15], Ni [16, 17], Pd [18], and Cu₃Au [19] nanowires, have been investigated by molecular dynamic (MD) or other computer calculation methods. However, few efforts have been focused on the thermal stability and melting behavior of the 1D binary alloy nanomaterial with different component.

Pd-Ni alloy can be used as catalyst in many reactions [20, 21] and hydrogen sensors [22]. However, its catalytic properties and performance of sensor relate to Pd or Ni concentration in the Pd-Ni alloy. Therefore, the structural stability and melting characteristics of Pd-Ni 1D alloy with different nickel concentration are studied by MD in the present work. And the relationship of melting points and fusion enthalpy of 1D Pd-Ni alloy models for MD process has also been investigated by using the model introduced in the literature [23].

2. Method of Simulation

The parallel code lrge-scale atomic/molecular massively parallel simulator (LAMMPS) [24] was used for all MD simulations. And visualization of MD simulation was performed using visual molecular dynamics (VMD) [25] and open visualization tool (Ovito) [26]. The embedded-atom method (EAM), a set of 𝑛-body potentials have been proved that it can accurately describe various dynamic properties of transition and noble metals. The interactions of atoms were modeled using the EAM potential, which is a convenient model to study metal clusters with FCC structure because the local density effects are included in its parameterization. In this work, the EAM potential data was obtained from the literature [27].

The geometry of nanowire is constructed as a regular FCC lattice with initial surface orientations of [100], [010], and [001] in the 𝑋, 𝑌, and 𝑍 directions of simulation box, respectively. An infinitely long nickel nanowire is modeled by applying periodic boundary conditions in the [001] direction. The other two Cartesian directions are also simulated with periodic boundary conditions, and the lengths of 𝑋, 𝑌 directions in the simulation box were large enough to prevent interaction between nanobelt models.

The initial configuration of nanobelt model with the size of 4.5𝑎∗4.5𝑎∗19.5𝑎 shown in Figure 1 has 2000 atoms, where the 𝑎 is the lattice parameter of crystal. In the present study, the Pd-Ni alloy nanobelt models were considered as the ideal solution, so the 𝑎 of Pd-Ni alloy nanobelt can be calculated from the lattice constant of nickel (𝑎_Ni= 0.352 nm) and palladium (𝑎_Pd= 0.389 nm) by using the expression of 𝑎=𝑥_Ni⋅𝑎_Ni+(1−𝑥_Ni)⋅𝑎_Pd, where 𝑥_Ni was the atomic concentration of nickel in the nanobelt model.

Figure 1 

The configuration of nanobelt model.

As the ideal binary mixture, positions of Ni and Pd atoms can exchange randomly. Through replacing the Pd atoms with Ni atoms randomly, the models with different Ni concentration can be made, and the content of Ni in the Pd-Ni alloy nanobelt was shown in Table 1. And the integer random distribution number (RND_i) was generated by RND_𝑖𝑅=INT^𝑖4,𝑅_𝑖=mod5∗𝑅_𝑖−1,4∗𝑀,𝑅₀≥1,𝑀=2^𝑘,𝑘=log₂𝑆+1. (1)

Table 1 

Content of nickel of every 1D Pd-Ni alloy model.

In order to remove the interactions which might lead to local structural distortion and result in the unstable simulation, energy minimization of the structure has been performed before starting molecular dynamics simulation. The MD simulations were carried out in the canonical ensemble with a constant number of atoms 𝑁 and volume 𝑉(𝑁𝑉𝑇). The equations of motion were integrated using the Verlet leapfrog algorithm with a time step of 5 × 10⁻⁴ ps. The nanobelt was heated to 1800 K in increments of 50 K, and the heating rate was 0.5 K/ps. Near the melting point (in the region of 900–1300 K), the temperature increments were reduced to 20 K to account for the large temperature fluctuations. At each temperature 10⁵ MD trajectory steps were propagated, which is sufficient to ensure the validity of the results presented here.

The radial distribution function (RDF) shown in (2), being regarded as one of the most important parameters, is used to describe the structure characterization of solid, amorphous, and liquid states [11]. Consider 𝑔(𝑟)=Ω𝑁²⋅∑^𝑁𝑖_𝑖=1𝑛_𝑖4𝜋𝑟²,Δ𝑟 (2) where 𝑔(𝑟) is the probability of finding an atom in a distance ranging from 𝑟 to Δ𝑟+𝑟 (Δ𝑟 is the step of calculation); Ω is the simulated volume of unit cell; 𝑁 is the number of atoms in the system; and 𝑁_𝑖 is the averaged number of atoms around the 𝑖th atom in the sphere shell ranging from 𝑟 to 𝑟+Δ𝑟.

The Lindemann index was recently applied to characterize nanoparticles, and it was proposed that melting of a nanoparticle at a critical value in the range of 0.03–0.05 [28]. The Lindemann indices of each atom in the models were calculated at each temperature as 𝛿_𝑖=1𝑁−1⋅_{𝑗()≠𝑖}𝑟²_𝑖𝑗𝑇−𝑟_𝑖𝑗²_𝑇𝑟_𝑖𝑗𝑇,1𝛿=𝑁⋅_𝑖𝛿_𝑖, (3)(4) where 𝛿_𝑖 and 𝛿 are the Lindemann indices of the 𝑖th atom and the whole model, respectively; ⟨⋯⟩_𝑇 in (3) denotes the thermal average at temperature 𝑇; 𝑟_𝑖𝑗 is the distance between the 𝑖th and 𝑗th atoms.

Pair analysis (PA) technique, which can effectively describe the characteristics of the geometric structural evolvement, has been used for analysing the geometric features of the atomic cluster [29]. In the present study, PA method was used to analyse the structural changes accompanying the melting process of Pd-Ni 1D models. Based on the regulation of bond pair, if two atoms are within a specified cutoff distance, they are called a bonded pair of models [29, 30]. In RDF curves, 𝑔(𝑟_𝑐) is the first minimum value; then 𝑟_𝑐 was defined as cutoff distance for PA. The four-index number 𝑖𝑗𝑘𝑙 is used to express bonded pairs of atomic clusters [29, 30]. If any atomic pair A-B formed a bond, 𝑖=1 and otherwise 𝑖=2; 𝑗 refers to the number of near neighbors which form bonds with both atom A and atom B; 𝑘 stands for the number of pairs among the neighboring atoms forming bonds; 𝑙 is a special distinguished index parameter. Based on the PA technique, the 1201 and 1311 bond pairs represent the rhombus symmetrical features of short-range order. The FCC structure has the type of 1421 bond pairs, whereas the HCP crystal has the equal number of 1421 and 1422 bond pairs. The difference between 1421 and 1422 bond pairs is the topological arrangement of the two bonds between the four neighbors. The bond pair 1551, 1541, and 1431, corresponding to a pentagonal bipyramid, is the characteristic of icosahedral order.

3. Results and Discussion

Figure 2 shows the total energy of nanobelts with different Ni content versus temperature (TE curve) during heating process. Total energies increase linearly with temperature in the early stage. When close to the melting transition temperature, simple jumps in total energy can be easily observed. The melting transition temperature of model can be estimated according to the TE curves. However, being different from that of cluster and bulk system, in the TE curve of nanobelt besides sharp increasing region, there is also a decreasing part after melting takes place. Therefore, the derivative curve of TE versus temperature shown in Figure 2, which can be considered as heat capacity [11], should have sharp peak and valley region. The melting points of the nanobelts can be determined by the position of the peak because of the heat absorption of melting. For reducing the total energy, especially surface energy, liquid 1D nanomaterial will change to spherical-like material due to the little interaction of atoms in the model at high temperature. That is to say, the nanobelt will break for decreasing the TE after 1D nanomaterial melting. From the derivative curve, the breaking temperature of 1D Pd-Ni can be estimated by the position of the valley. The melting temperature and breaking temperature of all the nanobelt models can be found in Figure 2. There are two sharp increasing and decreasing regions in the TE curve of PdNi1200 model. The reasons of emerging the first increasing and decreasing region were the melting of the model and decreasing of surface energy. In Figure 3(a), part of crystal structure can be found in the model, which results in PdNi1200 model holding 1D nanomaterial. With the temperature increasing, the melting of whole model (Figure 3(b)) and the collapse of nanobelt should induce the second increasing and decreasing region in the derivative curve.

Figure 2 

Total energy versus temperature curves for Pd-Ni nanobelt models during heating process.

(a)

(b)

Figure 3 

Snapshots of PdNi1200 nanobelt model at 1000 K and 1160 K.

Figure 4 shows the temperature dependence upon the Lindemann index of the nanobelts during heating. It is clear from Figure 4 that the Lindemann index increases slowly and linearly with temperature at low temperature stage because of the linear increase in atomic kinetic energy with temperature. The value of the Lindemann index during this stage is very small since most of the atoms do not have large amplitude motion but merely vibrate around their original lattice positions. Following the first stage, the Lindemann index increased and decreased rapidly. It is evident that the melting point and breaking temperature of the nanobelt can be estimated from the derivative curve of the Lindemann index versus temperature. At the temperature corresponding to the peak of derivative curve, lots of atoms in the nanobelt model exhibit large amplitude diffusion, which should be considered as melting of nanobelt. In the curve, there is also a sudden collapse of the 1D structure at the temperature corresponding to the valley. The melting temperature and breaking temperature estimated from derivative curves were the same as those shown in Figure 2. It can be found that with Ni concentration increase in the model, the fusing temperature decreases first and then increases. However, the breaking temperature of nanobelt with different Ni concentration did not obey this rule due to relating with the atom position and diffusion process of models. The Lindemann index of every atom in models can be used to determine the melting point and breaking temperature more precisely. From Figure 5, it can be found that the Lindemann index of atoms in PdNi1600 model was bigger than 0.05 at 1080 K which indicates that melting process occurs at this temperature [28]. At 1200 K, though there was no breaking, the Lindemann index of some atoms was very high, which suggests that the nanobelt will break at little higher temperature.

Figure 4 

The Lindemann index versus temperature curves of Pd-Ni nanobelt models.

Figure 5 

The atomic Lindemann index of PdNi1600 nanobelt model along 𝑍 direction.

Structural properties and changes of nanobelt during heating are of interest in understanding mechanical and other properties of materials. Figure 6 shows that there are only 1421, 1311, and 1211 pairs in the 1D Pd-Ni alloy models at very low temperature. It can be found that the number of 1421 pairs will decrease obviously with increasing temperature. Due to the collapse of 1421 pairs, atoms in the FCC structure should form other pairs, which can be identified by the elevation of the 1201, 1301, and 1311 pair numbers. At melting temperature, the number of 1421 and 1311 pairs in the nanobelt drops abruptly due to melting transition of nanobelt, and the 1421 pairs almost disappear after nanobelt melting. It can also be found that the number of 1201 and 1321 bondtypes increased obviously at melting temperature, which was caused by the collapse of 1311 and 1421 pairs.

Figure 6 

Relative amount of various atomic-bonded pairs in PdNi1600 nanobelt model at different temperature during heating process.

Through pure element thermodynamic data, the classical thermodynamics method can be used to predict many thermodynamic parameters of bulk alloy, such as enthalpy and melting point. Therefore, by building the relationship of thermodynamic properties of nanomaterial and bulk material, thermodynamic parameters of nanomaterial can also be calculated by thermodynamics method. In this study, the fusion enthalpy of 1D Pd-Ni alloy model was calculated by (5) which was reported by Guisbiers and Buchaillot [31]. Consider Δ𝐻_𝑚Δ𝐻_𝑚,∞=𝑇^∗_𝑚𝑇_𝑚,∞,Δ𝑆_𝑚Δ𝑆_𝑚,∞≈1, (5) where 𝑇^∗_𝑚 and Δ𝐻_𝑚 are the size-dependent melting point and fusion enthalpy, respectively; 𝑇_𝑚,∞ and Δ𝐻_𝑚,∞ are the melting point and fusion enthalpy of bulk material, respectively. Here, 𝑇^∗_𝑚 value can use the results of MD simulation. Thermodynamic parameters of bulk Pd-Ni alloy can be calculated from the pure Ni and Pd thermodynamic data [32] which has been shown in the Appendix. For example, Δ𝐻_𝑚,∞, Δ𝑆_𝑚,∞, and 𝑇_𝑚,∞ can be calculated by (6) [33] and (7), respectively. The expression of 𝐺^𝐿 and 𝐺^𝑆, which were Gibbs energy of liquid and solid phases of bulk material, used in (6) and (7) can be illustrated by (8). Then, using the thermodynamic data of bulk material calculated previously, Δ𝐻_𝑚 and Δ𝑆_𝑚 will be obtained by (5). Consider Δ𝐻_𝑚,∞=𝐺^𝐿−𝑇⋅𝑑𝐺^𝐿𝑑𝑇_{𝑇=𝑇𝐿,∞}−𝐺^𝑆−𝑇⋅𝑑𝐺^𝑆𝑑𝑇_{𝑇=𝑇𝑆,∞},Δ𝑆_𝑚,∞=−𝑑𝐺^𝐿𝑑𝑇_{𝑇=𝑇𝐿,∞}−−𝑑𝐺^𝑆𝑑𝑇_{𝑇=𝑇𝑆,∞},𝜕𝐺^𝑆𝜕𝑥_𝑥𝑆𝑖=𝜕𝐺^𝐿𝜕𝑥_𝑥𝐿𝑖,𝜕𝐺^𝑆𝜕𝑥_𝑥𝑆⋅𝑥^𝐿_𝑖−𝑥^𝑆_𝑖=𝐺^𝐿𝑥^𝐿_𝑖−𝐺^𝑆𝑥^𝑆_𝑖,𝐺=𝑥_𝑖×𝐺⁰_𝑖+𝑥_𝑗×𝐺⁰_𝑗+𝑅𝑇𝑥_𝑖ln𝑥_𝑖+𝑥_𝑗ln𝑥_𝑗+𝐺^ex,_ex𝐺7=_𝑖≠𝑗𝑥_𝑖⋅𝑥_𝑗⋅𝐿⁰+𝑥_𝑖−𝑥_𝑗⋅𝐿¹+𝑥_𝑖−𝑥_𝑗²⋅𝐿²+⋯. (6)(7)(8)

All the calculated results have been shown in Table 2. It can be found that the relation between melting point of nanobelt model and Ni concentration is the same as that of bulk phase Pd-Ni alloy. However, because the specific surface area of nanomaterial is very great, the Pd-Ni nanobelt melt temperature is smaller than that of bulk phase. The melting points of nanobelts can be reproduced by equation of 𝑇_𝑚^∗∗=Δ𝐻_𝑚/Δ𝑆_𝑚. The results shown in Table 2 indicate that the melting point obtained by MD simulation is close to 𝑇_𝑚^∗∗, and the error is below 5%. Particularly, the difference of melt point between MD and thermodynamic method for PdNi800 and PdNi1200 was only 3-4 K. The nanobelt models used in this work were considered as the ideal solution, so the error may be different when the other mix model was used, such as aggregation state model.

Table 2 

Thermodynamic parameters of 1D Pd-Ni alloy models and bulk Pd-Ni alloy.

4. Conclusion

Summarily, the heating process of Pd-Ni alloy nanobelts has been simulated by MD and the melting temperature was also calculated by thermodynamic method. The melting and breaking of Pd-Ni alloy nanobelts can be indicated by increasing and decreasing regions of total energy relating to temperature curve. With increasing temperature, number of FCC pairs will decrease and almost disappear at melting point. The melting point is lower than bulk phase Pd-Ni alloy, and the melt point obtained by MD is close to that of calculated by thermodynamic method.

Appendix

FCC-Ni𝐺^0solid_Ni=−5179.159+117.854𝑇−22.096𝑇ln(𝑇)−4.8407×10⁻³𝑇²()𝐺298.15<𝑇<1728_0solidNi=−27840.620+279.134977𝑇−43.1𝑇ln(𝑇)+1127.54×10²⁸𝑇⁻⁹(1728<𝑇<3000). (A.1)

Liquid-Ni𝐺^0liquid_Ni(𝑇)=11235.527+108.457𝑇−22.096𝑇ln−4.8407×10⁻³𝑇²−382.318×10⁻²³𝑇⁷()𝐺298.15<𝑇<1728_0liquidNi(𝑇).=−9549.817+268.597977𝑇−43.1𝑇ln)(1728<𝑇<3000 (A.2)

FCC-Pd𝐺^0solid_Pd=−10204.027+176.076315𝑇−32.211𝑇ln(𝑇)+7.120975×10⁻³𝑇²−1.919875×10⁻⁶𝑇³+168687𝑇⁻¹()𝐺298.15<𝑇<900_0solidPd=917.062+49.659892𝑇−13.5708𝑇ln(𝑇)−7.17522×10⁻³𝑇²+0.191115×10⁻⁶𝑇³−1112465𝑇⁻¹(900<𝑇<1828). (A.3)

Liquid-Pd𝐺^0liquid_Pd=1302.731+170.964153𝑇−32.211𝑇ln(𝑇)+7.120975×10⁻³𝑇²−1.919875𝐸−6𝑇³+168687𝑇⁻¹()𝐺298.15<𝑇<600_0liquidPd(𝑇)=23405.778−116.918419𝑇+10.8922031𝑇ln−27.266568×10⁻³𝑇²+2.430675×10⁻⁶𝑇³−1853674𝑇⁻¹()𝐺600<𝑇<1828_0liquidPd(𝑇)𝐿=−12373.637+251.416903𝑇−41.17𝑇ln)(1828<𝑇<4000_0liquidNi,Pd𝐿=9180.1080−17.2493𝑇_1liquidNi,Pd𝐿=424.2768+3.6219𝑇_2liquidNi,Pd𝐿=−5794.1283+6.9498𝑇_0solidNi,Pd𝐿=−15866.3670+5.6399𝑇_1solidNi,Pd𝐿=5276.7800+2.0358𝑇_2solidNi,Pd=12038.2371−8.3220𝑇. (A.4)

Acknowledgments

The authors acknowledge the support of Project supported by the National Natural Science Foundation of China (Grant nos. 51207031, 51177022), the Promotive Research Fund for Excellent Young and Middle-Aged Scientists of Shandong Province (Grant no. BS2011NJ002), China Postdoctoral Science Foundation (Grant no. 2013M541368), the Science Foundation of HIT at Weihai (Grant no. IMJQ10080026), and the Natural Scientific Research Innovation Foundation in HIT (Grant no. IMXK57080026).