Molecular Dynamics Simulation of Ni Nanoparticles Sintering Process

Article pubs.acs.org/JPCC

Molecular Dynamics Simulation of Ni Nanoparticles Sintering Process in Ni/YSZ Multi-Nanoparticle System Jingxiang Xu, Ryota Sakanoi, Yuji Higuchi, Nobuki Ozawa, Kazuhisa Sato, Toshiyuki Hashida, and Momoji Kubo* Fracture and Reliability Research Institute (FRRI), Graduate School of Engineering, Tohoku University, 6-6-11 Aoba, Aramaki, Aoba-ku, Sendai 980-8579, Japan ABSTRACT: We have developed a molecular dynamics (MD) simulation method to investigate the sintering of nickel nanoparticles in the nickel and yttria-stabilized zirconia (Ni/YSZ) anode of a solid oxide fuel cell (SOFC). The conventional sintering model consists of only two or three nickel nanoparticles. Therefore, it does not reflect the properties of the porous structure of the Ni/ YSZ anode or reproduce realistic sintering. Our Ni/YSZ multi-nanoparticle MD simulation method uses a multi-nanoparticle model based on the porosity and Ni/YSZ nanoparticle ratio of a realistic anode. The Ni and YSZ nanoparticles are packed randomly in the simulation cell, and compressed to achieve the correct porosity. Furthermore, because the reliable potential parameters for MD simulation between nickel and YSZ have not been reported, we determine reliable interatomic potential parameters between nickel and YSZ by using the nonlinear leastsquares method to fit the Morse potential function to interaction energies obtained by density functional theory. The sintering simulation using our Ni/YSZ multi-nanoparticle model and our potential parameters reveal that the YSZ nanoparticle framework suppresses the sintering of nickel nanoparticles by disrupting the growth of the neck between two nickel nanoparticles. The previously reported model of two nickel nanoparticles did not produce these results. Our multi-nanoparticle MD simulation method is effective for investigating the realistic sintering process in the porous structure of the Ni/YSZ anode and for designing durable anode structures for SOFCs.

1. INTRODUCTION

Molecular dynamics (MD) simulations have played an important role in elucidating the sintering of two nickel,16 two copper,17 two gold,18 and two titanium oxide nanoparticles19 on the atomic scale, by tracking the neck and shrinkage between the sintered nanoparticles. Sintering of three aluminum nanoparticles21 and three copper nanoparticles22 has also been investigated. In addition, Hawa et al. have performed pioneering works of sintering of 40-silicon nanoparticles in a straight-chain without a substrate by using MD simulation. They found that the sintering of 40 nanoparticles forms a cylindrical shape and then they aggregate into the sphere.48 These studies suggest that MD simulations are effective for investigating nanoparticle sintering. However, the majority of the studies focus on simple metal/metal or metal oxide/metal oxide systems, and metal/metal oxide systems have rarely been studied. Moreover, most MD sintering simulations have concentrated on systems that contain only two or three nanoparticles16−26 (Figure 1a). Recently, the pioneering works of sintering of 40silicon nanoparticles in a straight chain without a substrate were reported.48 However, the anode structure is porous and consists of Ni/YSZ nanoparticles (Figure 1b). The porosity of the anode has a significant effect on the electrocatalysis.28 The YSZ

Fuel cells, which convert chemical energy to electrical energy, are expected to become a widely used technology for generating electricity. Solid oxide fuel cells (SOFCs) have attracted considerable interest because of their efficient energy production and low pollutant emissions. They are also relatively inexpensive because they do not require a precious metal catalyst such as platinum.1−8 Currently, nickel and yttriastabilized zirconia (Ni/YSZ) cermet (a mixture of a ceramic and a metal) is the most widely used anode material for SOFCs,8−10 because of its good electrocatalytic oxidation properties.3 However, nickel has a relatively low melting temperature and a high surface mobility at the high operating temperatures in SOFCs, which results in nickel sintering in the anode.11 The sintering reduces the catalytic activity of the Ni/ YSZ anode, because it reduces the surface area of the nickel nanoparticles and thus the electrical conductivity.12−15 Hence, the sintering of the anode is a major obstacle to realizing commercial SOFCs. A detailed understanding of the sintering process in the Ni/YSZ anode is necessary for improving the anode performance. However, little is known about the atomicscale sintering dynamics and changes in the microstructure during the operation of SOFCs. Furthermore, the sintering of the nickel nanoparticles in the Ni/YSZ anode is not well understood, although atomic forces play a role in the sintering behavior. © 2013 American Chemical Society

Received: November 5, 2012 Revised: April 16, 2013 Published: April 16, 2013 9663

dx.doi.org/10.1021/jp310920d | J. Phys. Chem. C 2013, 117, 9663−9672

The Journal of Physical Chemistry C

Article

2. COMPUTATIONAL DETAILS A. Molecular Dynamics Simulation Method. The sintering of nickel nanoparticles in the Ni/YSZ anode was simulated by using the MD program package, New-Ryudo.31 The Verlet algorithm32 was employed to integrate the equation of motion. The MD simulations were executed with a 2.0 fs time step, and the NVT canonical ensemble was performed at a temperature of 1073 K, which is a common SOFC operating temperature. For the force field selection, the Born-Mayer-Huggins (BMH) potential was used to describe the ionic interactions of YSZ: U (rij) =

ZiZje 2 rij

⎛ ai + aj − rij ⎞ + f0 (bi + bj) × exp⎜⎜ ⎟⎟ ⎝ bi + bj ⎠

(1)

where Zi, Zj is the effective charge of ion i and j, e is the elementary electric charge, rij is the interatomic distance, and f 0 is a constant to adjust the units. The repulsion term parameters, ai and bi, are related to the size and stiffness, respectively, and were determined in our previous work (Table 1).33 Our

Figure 1. Schematic diagrams. (a) The two-nanoparticle model of sintering, and (b) the three-dimensional porous structure of Ni/YSZ.

Table 1. BMH Potential Parameters for YSZ33 29,30

nanoparticle framework allows oxide ion conduction, and the ratio of the nickel to YSZ nanoparticles affects the SOFC performance.10 The two or three nanoparticle models16−21 cannot take into account the effect of the YSZ nanoparticle framework on the sintering of the nickel nanoparticles. The pioneering works of sintering of 40-silicon nanoparticles in a straight-chain48 did not reflect the effect either. In addition, the previous sintering simulations of metal nanoparticles on metal oxide surface23 are unable to reflect the porosity and the ratio of nickel and YSZ. Thus, the sintering process of nickel particles in the realistic Ni/YSZ anode has not been clear. On the other hand, in the metal/metal oxide systems, it is difficult to determine the interatomic potential parameters for MD simulations of metal/metal oxide systems because it is impossible to obtain the physical properties of the metal/metal oxide interface experimentally.20 Therefore, in our previous MD simulation studies on the sintering and deposition of palladium and gold nanoparticles on a MgO surface,23−26 the potential parameters were determined by reproducing their structures. This revealed that the parameters for the metal and metal oxide were not sufficiently accurate. We established an effective method for accurately determining the interatomic potential parameters between the metal and the metal oxide by using density functional theory (DFT).27 These DFT-based potential parameters accurately reproduce the structure of a palladium particle on a MgO(001) surface, although a sintering simulation of metal/metal oxide systems has not yet been reported. In the present study, we investigated the sintering of nickel nanoparticles in a porous Ni/YSZ structure. We developed a multi-nanoparticle MD simulation method that reproduces the features of the actual anode structure, such as the porosity and the nanoparticle framework. The interatomic potential parameters in the simulation were calculated by the DFT method from the interaction energy between nickel and YSZ. We have successfully used our method to simulate the effect of the porous anode structure on the sintering observed in SOFCs.

atom

Zi

ai /Å

bi /Å

O Zr Y

−1.2 +2.4 +1.8

1.503 1.227 1.327

0.075 0.070 0.070

previous studies of the sintering and deposition of Pd and Au on a MgO surface23−26 showed that the Morse potential function is suitable for describing the metal/metal interaction. The Morse potential was used for reproducing the interaction between the nickel atoms: E(rij) = Dij{exp[−2βij(rij − rij*)] − 2 exp[−βij(rij − rij*)]} (2)

where rij is the interatomic distance, and the parameters Dij, βij, and r*ij are related to the bond energy, stiffness, and bond length, respectively. The previously reported potential parameters for nickel34 are presented in Table 2. Table 2. Morse Potential Parameters for Nickel34 i−j

Dij (kcal/mol)

βij (1/Å)

rij* (Å)

Ni−Ni

9.70

1.42

2.78

In the metal/metal oxide system, the studies with respect to the sintering and deposition of Pd and Au on the MgO surface23−26 have been performed in our group and show that the Morse potential function is suitable to describe the interaction between metal/metal oxide. Considering interactions between the nickel atoms and YSZ nanoparticles, we employ the Morse potential interatomic forces for describing the interaction of Ni−Zr, Ni−Y, and Ni−O bonds. In this study, the simulations for two nickel nanoparticle and multiparticle models are performed in two 6-core 3.46 GHz Intel Xeon Processors X5690 with 12 M cache. B. Parameterization Method for the MD Simulation. Accurate interatomic potential parameters between nickel and YSZ were determined by fitting the Morse potential function to the interaction energies obtained from the DFT calculations. 9664



Article

for the metal/metal oxide system to the interaction energies between nickel and YSZ calculated by DFT. Previous studies have shown that the {111} plane of the YSZ surface is more stable than other low-index planes39 and that a nine-layer structure is sufficient to represent the YSZ(111) surface.40 Thus, the nine-layer structure of the YSZ(111) surface was adopted as the substrate for the DFT calculations. In addition, the stability of the YSZ bulk depends on the positions of the yttrium atoms and the oxygen vacancies.41 Therefore, the YSZ(111) surface model was built on the basis of DFT calculations on the total energies of the YSZ(111) surface, depending on the positions of the yttrium atoms and oxygen vacancies, with DMol3.35,36 The total energies showed that the most stable YSZ(111) surface structure occurred where the oxygen vacancy occupies the nearest neighbor sites to the pair of yttrium atoms (Figure 2). The interaction energies between nickel and YSZ were calculated by using one nickel atom on the YSZ(111) surface as the model for the DFT calculations. The interaction energy is defined as

The interaction energies between one nickel atom and the YSZ surface at surface several sites and at different vertical distances from the surface were calculated with DFT. We used the nonlinear least-squares method to fit the Morse potential parameters describing the interaction energies between nickel and YSZ to the interaction energies obtained by DFT calculations. Spin-polarized DFT was employed to calculate the interaction energies between the nickel and the YSZ(111) surface with DMol3.35,36 In the DFT calculations, the effective core potentials were used to model the core electrons, and the double numerical plus polarization basis sets were used to describe the atomic orbitals. A generalized gradient approximation using the Perdew-Burke-Ernzerhof37,38 correlation functional was applied to the exchange-correlation term in all the calculations. One nickel atom on the YSZ(111) surface was used in the DFT calculation model (Figure 2). For the

E = E Ni/YSZ(111) − (E Ni atom + E YSZ(111) slab)

(3)

where E is the interaction energy of the nickel atom and the YSZ(111) surface, ENi/YSZ(111) is the total energy of the adsorbed nickel atom and the YSZ(111) surface, and ENi atom and EYSZ(111) slab are the total energy of the isolated nickel atom and the YSZ(111) surface, respectively. Moreover, the interaction energy is almost zero when the distance between the nickel atoms and the YSZ(111) surface is large. To fit the Morse potential function, the interaction energies were set to zero when the distance was greater than or equals to 8 Å. The interaction energies obtained from the DFT calculations were divided into three sets of interaction energies (Ni−Zr, Ni−O, and Ni−Y) for fitting the MD parameters. The interaction of the Ni/YSZ system was described by dividing the basic Morse potential function (eq 2) into three interaction energy terms, for Ni−O, Ni−Zr, and Ni−Y. This Morse potential function can be rewritten as

Figure 2. Front view of one nickel atom on the nine-layer YSZ(111) surface. The interaction energies were calculated for each distance, R, at each site. R is the vertical distance between the nickel atom and sites on the YSZ(111) surface, and was between 1.5 and 5.0 Å.

E=

∑ DNi− Zr {exp[−2βNi− Zr (rNi− Zr − rNi* − Zr)] Zr

* − Zr )]} − 2exp[−βNi − Zr (rNi − Zr − rNi

calculation of the interaction energies, the nickel atom was located at the zirconium, oxygen, and yttrium sites on the YSZ(111) surface, at a distance, R, perpendicular to the surface, where R was varied between 1.5 and 5.0 Å at intervals of 0.1 Å. The interaction energies at 36 points were calculated for each site. Thus, the interaction energies at a total of 108 points at the three different sites were calculated and were reproduced by the Morse potential function; fitting many interaction energies at different points based on DFT calculations realized the high accuracy of the interaction potentials.

+

∑ DNi− O{exp[−2βNi− O(rNi− O − rNi* − O)] O

* − O)]} − 2exp[−βNi − O(rNi − O − rNi +

∑ DNi− Y {exp[−2βNi− Y (rNi− Y − rNi* − Y )] Y

* − Y )]} − 2exp[−βNi − Y (rNi − Y − rNi

(4)

where the parameters DNi−A, βNi−A, and DNi−A * are related to the bond energy, stiffness, and the bond length between nickel and A, respectively. Here, A represents any atom in the YSZ(111) surface (Zr, O, or Y). The interaction energies at a total of 108 points for three sites were calculated and the interaction energies at 13, 17, and 21 points for sites at the O, Zr, and Y atoms on the YSZ(111) surface were obtained because the calculations at other points did not converge. The three sets of interaction energies, E, were obtained as a function of the vertical distance between the nickel atom and the YSZ(111) surface. The least-squares criterion for fitting was defined as

3. RESULTS AND DISCUSSION A. Determination of Ni/YSZ Interaction Potential Parameters. The potential parameters between the metal and metal oxide for MD simulation are difficult to determine, because it is impossible to obtain the physical properties of the metal/metal oxide interface experimentally. Thus, empirical potential parameters have been used in previous MD simulations of the metal/metal oxide system.23−26 We used accurate interatomic potential parameters for nickel and YSZ, which were determined by fitting the Morse potential function 9665



Article

min{∑ ∑ (Efit − E DFT)2 } site

(5)

R

where R is the vertical distance between the nickel atom and the sites on the YSZ(111) surface, and Efit and EDFT are the interaction energy values obtained by fitting the Morse potential and from the DFT calculations, respectively. The Levenberg−Marquardt (LM) method42−44 was used to solve the test function (eq 5) and determine the three sets of parameters, D, β, and r*, which are related to the bond energy, stiffness, and bond length of the Ni−O, Ni−Zr, and Ni−Y bonds in the rewritten Morse functions (eq 4). The fitting relied on simultaneously satisfying the minima of eq 5 with the interaction energies at 51 points between the nickel atom and the three sites on the YSZ(111) surface, while applying the periodic boundary condition within the cutoff value corresponding to half of the unit cell length. The interatomic potential parameters calculated for nickel and YSZ are presented in Table 3. The three sets of interaction Table 3. Morse Potential Parameters for Ni/YSZ i−j

Dij (kcal/mol)

βij (1/Å)

rij* (Å)

Ni−O Ni−Zr Ni−Y

31.92 0.56 0.84

2.35 1.31 1.08

1.77 3.80 4.17

energies, E, are shown as a function of the vertical distance between the nickel atom and the YSZ(111) surface in Figure 3. The interaction energies obtained from the DFT calculations and the fitted Morse potential curve of the nickel atom at the oxygen, zirconium, and yttrium sites are shown in Figure 3a−c. The fitted Morse potential curves accurately reproduced the DFT interaction energies at 51 points for sites at the O, Zr, and Y atoms on the YSZ(111) surface. The Morse potential function could easily be fitted to the interaction energies at one site. However, we used the LM method42−44 to fit and concurrently reproduce the interaction energies at 51 points for sites at the O, Zr, and Y atoms on the YSZ(111) surface in order to obtain highly accurate potential parameters for nickel and YSZ. The values of Dij, which is related to the bond energy, were 0.56, 31.92, and 0.84 kcal/mol for Ni−Zr, Ni−O, and Ni− Y, respectively. This indicates that the strongest interaction was between nickel and oxygen. The values of r*ij , which is related to bond length, were 3.80 and 4.17 Å for Ni−Zr and Ni−Y, respectively. Hence, the parameters related to the ionic bond length of Zr were smaller than those related to ionic bond length of Y, which is consistent with the experimentally determined ionic radii of Zr4+ and Y3+, 0.84 Å and 1.02 Å, respectively.45 Our fitting method successfully reproduced the results of the DFT calculations for the three sites on the YSZ(111) surface. Sequentially, to validate the potential parameters obtained by fitting the potential energy of a single nickel atom on an YSZ surface, we calculated the interaction energies between Ni clusters and YSZ(111) surface by both DFT and MD methods using our potential parameters. Ni4 and Ni7 are set on the surface as typical models of the Ni clusters because Ni4 is regular tetrahedron and shows spherical shape, and Ni7 is obtained by adding one layer to Ni4. The calculation models of Ni atom, Ni4 cluster, and Ni7 cluster on YSZ(111) surface are shown in Figure 4 and the obtained interaction energies by the DFT and MD methods are summarized in Table 4. The results

Figure 3. Interaction energies obtained from the DFT calculations and the fitted Morse potential curve of the nickel atom at the (a) oxygen, (b) zirconium, and (c) yttrium sites. The squares indicate the interaction energies obtained by the DFT calculation and the solid line corresponds to the potential curves determined by the least-squares fitting method.

Figure 4. Top view of (a) single nickel atom, (2) Ni4 cluster, and (c) Ni7 cluster on the YSZ(111) surfaces. The green, water color, red, and dark blue correspond to Y, Zr, O, and Ni.

show that MD simulations using our parameters are in excellent agreement with DFT calculations. The errors of the interaction energy between Ni4 cluster and YSZ surface and between Ni7 cluster and YSZ surface are less than 3%. This implies that our parametrization method by fitting the potential energy to DFT calculation results on the YSZ surface is sufficient to determine 9666



Article

the sintering mechanism. The diffusion of the initial surface and interior atoms of the nickel nanoparticles was tracked over time (Figure 5). The surface thickness was 2 Å, and the color of the initial surface atoms of the nickel nanoparticles in Figure 5 is different to distinguish them from the initial interior atoms. Initially, the two nickel nanoparticles were attracted to each other, with no external force. The attraction and deformation of the two nickel nanoparticles shows there was no change in the position of the initial surface atoms relative to the initial interior atoms of the nickel nanoparticles before 16 ps (Figure 5a,b). This indicates that the sintering took place after 16 ps. In addition, the four nickel atoms were given different color to track the trajectory of the atoms as shown in the Figure 5a,b. We observed that the left nanoparticle rotated clockwise and the right nanoparticle slightly rotated counterclockwise (Figure 5a,b). This implies that the two nanoparticles rotated to adjust the lattice plane in order to find the best contact point. This phenomenon has also been observed in previous MD simulations of the sintering of gold nanoparticles18 and of titanium oxide nanoparticles.20 The initial nickel nanoparticle surface atoms then moved to the contact area between the two nanoparticles to form a neck, and the gap between the two nanoparticles gradually filled, broadening the neck (Figure 5b− d). Figure 5c,d shows that some initial interior atoms migrated outward and became new surface atoms (black circles, Figure 5c,d). These new surface atoms were able to move to the neck region. Therefore, the diffusion of the initial surface atoms contributed to the growth of the neck region. The diffusion of nickel atoms in the nanoparticles was quantified and evaluated by measuring the mean square displacement (MSD) during the sintering. The MSD of nickel

Table 4. Interaction Energies between Nickel and YSZ E

Ni Atom

Ni4 Cluster

Ni7 Cluster

DFT (kcal/mol) MD (kcal/mol) Error (%)

−35.732a −35.972a 0.670b

−99.455a −102.334a 2.894b

−129.807a −133.064a 2.509b

E = ENi/YSZ(111) − (ENi atom + EYSZ(111) slab). bError = (EDFT − EMD)/ EDFT × 100%.

a

the highly accurate parameters between nickel and YSZ for describing the interaction in the Ni/YSZ system. B. Simulation of the Sintering of Two Nickel Nanoparticles on the YSZ(111) Surface. Initially, a simplified model consisting of two nickel nanoparticles on the YSZ(111) surface was used to investigate the fundamental characteristics of the sintering of nickel nanoparticles in the Ni/YSZ system. The simplified model contained 24 030 atoms, and the two nickel nanoparticles with a diameter of 40 Å were separated by 5 Å. The size of the YSZ(111) surface was 116.1 Å × 100.5 Å. The equilibrated Ni/YSZ(111) structure was obtained by stabilizing the simulation for 50 ps at 300 K in the NVT canonical ensemble before carrying out the sintering simulation. We observed that the MSD value does not change after 18 ps. Then, sintering simulations were performed for a total time of 1 ns. After 500 ps to 1 ns, the sintering process does not change qualitatively, which is shown by snapshots and mean square displacement. Furthermore, the nickel nanoparticles sintered within a few hundred picoseconds. Then, we showed the sintering simulation results within 500 ps in this paper. Figure 5 shows cross-sectional snapshots of the sintering of the nickel nanoparticles on the YSZ(111) surface at 0, 16, 60, and 500 ps. The diffusion of the nickel atoms played an important role in

Figure 5. Cross-sectional snapshots of the sintering of the nickel nanoparticles in the two nickel nanoparticle model on the YSZ(111) surface at 1073 K for (a) t = 0 ps, (b) t = 16 ps, (c) t = 60 ps, and (d) t = 500 ps. The four nickel atoms were given different colors to track easily the trajectory of the atoms as shown in Figure 4a,b. 9667



Article

nanoparticles in the Ni/YSZ nanoparticle system is calculated with the following expression: ⟨r 2(t )⟩ =

1 N

N

∑ |ri(t ) − ri(0)|2 i=0

(6)

where ⟨r (t)⟩ is the MSD, N is the total number of diffusing nickel atoms, and t is the time of the simulation. The diffusion mechanism during the sintering was investigated by evaluating the MSD of the initial surface and interior atoms of the nickel nanoparticles on the YSZ(111) surface (Figure 6). Before 2

Figure 7. Snapshots of the sintering of the nickel nanoparticles with a diameter of 30 (a)−(c) and 50 Å (d)−(f) in two nickel nanoparticle models on the YSZ(111) surface at 1073 K. (a) t = 0 ps, (b) t = 6 ps, (c) t = 500 ps, (d) t = 0 ps, (e) t = 21 ps, and (f) t = 500 ps.

previous studies investigated the effect of the nanoparticle size on the sintering in the two or three nanoparticle model on the plane surface by MD simulation.16,19,21,47 They showed that increasing in the nanoparticle size slows down the sintering process while the sintering mechanism does not change. Our result is in good agreement with previous results. According to these results, we conclude that the change of the nanoparticle size does not affect the sintering mechanism. C. Ni/YSZ Multi-nanoparticle Sintering Simulation. The sintering in the Ni/YSZ anode is closely related to the porosity of anode,28 the YSZ nanoparticle framework,29,30 and the ratio of nickel to YSZ nanoparticles.10 However, previous sintering simulations, which have focused on two or three nanoparticle models,16−24 cannot model the properties of porous SOFC anode materials. The pioneering works of sintering of 40-silicon nanoparticles in a straight chain without a substrate48 did not reveal the effect of porous structures either. Therefore, a multi-nanoparticle model which reflects the porosity and framework of YSZ nanoparticles in the actual anode structure is necessary for investigating sintering in the porous Ni/YSZ structure, but the simulation method has not been established yet. In this study, we developed a modeling methodology of Ni/YSZ multi-nanoparticle model. Because the typical porosity of anode materials currently used in SOFCs is 30−45%,28,8 we used a porosity of 40% in our simulation. Five nickel and five YSZ nanoparticles with a diameter of 40 Å were located randomly in a 100 Å × 100 Å × 100 Å simulation cell (Figure 8a). Thus, we compressed the initial Ni/YSZ nanoparticle model by using the MD simulation in the NPT canonical ensemble under a pressure of 5 GPa, until its porosity reached 40%, and the size of the simulation cell reached 81.3 Å × 84.1 Å × 81.9 Å (Figure 8b). The Ni/YSZ multi-nanoparticle model was equilibrated by a stabilization simulation for 50 ps at 300 K and we observed that the MSD value does not change after 30 ps. Then, a more realistic sintering simulation of the Ni/YSZ anode model was performed for 1 ns with a 2.0 fs time step in the NVT canonical ensemble at a temperature of 1073 K. Although we performed sintering simulations for a total time of 1 ns, we showed the sintering simulation results within 500 ps, since the nickel nanoparticles fully sintered within a few hundred picoseconds. Figure 9 shows snapshots of the sintering in the Ni/YSZ multi-nanoparticle model. The initial state of the model, which was stabilized at 300 K, is shown in Figure 9a. The nickel and YSZ nanoparticles were uniformly distributed around the nickel nanoparticles and formed the YSZ framework. The sintering

Figure 6. Mean square displacement of nickel atoms for the model of two nickel nanoparticles on the YSZ(111) surface at 1073 K.

sintering, the MSD values of the initial surface and interior atoms increased slowly because of the thermal diffusion. However, when the two nickel nanoparticles came into contact with each other, the MSD of the surface and interior nickel atoms showed a rapid, linear increase from 16 to 26 ps. The increase in the MSD of the surface atoms was greater than that of the interior atoms (Figure 6). This indicates that the surface atoms diffused more rapidly than the interior atoms, which is consistent with the MD simulation reported by Güvenç et al.46 As the surface atoms diffused, the interior atoms of the nickel nanoparticles rearranged themselves and migrated to the nanoparticle surface (Figure 5). This may have increased the MSD of the interior atoms before 26 ps. However, some of the interior atoms diffused to the neck region. When the interior atoms came into contact with each other in the neck region, they tend to be relatively diffusionless. In contrast to the previous stage, the interior atoms behaved like a solid, whereas the surface atoms remained liquid-like. This may explain why the MSD values of the surface atoms were larger than those of the interior atoms and why the difference between the MSD values of the surface and interior atoms increased after 26 ps (Figure 6). Therefore, surface diffusion was the dominant mechanism in the sintering of nickel nanoparticles on the YSZ(111) surface, and it contributed to the formation of the grain boundary and the growth of the neck. Moreover, we discuss the influence of the nanoparticle size on the sintering process of the two nickel nanoparticles on the YSZ surface. We calculate the different nanoparticle sizes with diameters of 30 and 50 Å in two nickel nanoparticle models to compare with the result of the diameter of 40 Å. Figure 7 shows snapshots of these results. First, two nickel nanoparticles are attracted each other. Then, they contact and form a neck. This indicates that the sintering mechanism is consistent with the result of a diameter of 40 Å. An increase in the size of nickel nanoparticles slows down the sintering processes. Several 9668



Article

Ni/YSZ anodes coalesce during the operation of SOFCs at 1073 K. Hence, our multi-nanoparticle model successfully reproduces experimentally observed sintering.12 The details of the sintering process in Ni/YSZ multinanoparticle model were determined by examining the crosssectional snapshots of the sintering process of the multinanoparticle model. Figure 10a−c shows the sintering process in the large pore of the YSZ framework at 0, 11, and 500 ps. Figure 10d−f shows the sintering process in the small pore of the YSZ framework at 0, 7, and 500 ps. For clarity, the region shown in the snapshots has been extended based on the periodic boundary condition. The two nickel nanoparticles in Figure 10a−c are defined as No. 1 and No. 2, and the nanoparticles in Figure 10d−f are defined as No. 3 and No. 4. Figure 10a,d shows the initial stage of the model. The two nickel nanoparticles in (a) glided over the spherical surface of the YSZ nanoparticle and approached each other through the large pore in the YSZ framework structure, because of the attractive force between the nanoparticles, with no external force. The two nickel nanoparticles made contact at 11 ps (Figure 10b), then a neck formed between them and broadened until the particles coalesced and the neck shape almost disappeared at 500 ps (Figure 10c). The two nickel nanoparticles finally formed a single, compact, spherical nanoparticle. Nanoparticle No. 3 approached nanoparticle No. 4 passed through the small pore in the YSZ framework and made contact at 7 ps (Figure 10e), forming a neck which then broadened (Figure 11f). However, the growth of the neck between nanoparticles No. 1 and No. 2 was different from that between nanoparticles No. 3 and No. 4. The neck between nanoparticles No. 3 and No. 4 could not disappear because the YSZ nanoparticles around them limited the broadening of the neck. Therefore, the small pores in the YSZ framework suppressed the neck growth and prevented the nickel nanoparticles from coalescing. The degree of sintering was quantified by calculating the distance between the mass of two nanoparticles as a function of time (Figure 11). The sintering processes between nickel nanoparticles in the large pore of YSZ framework (Figure 10a− c) and the small pore of YSZ framework (Figure 10d−f) are referred to as case 1 and case 2, respectively. The distance in case 1 reached a minimum more rapidly than that in case 2; the minimum distances for case 1 and case 2 were 30.7 and 38.3 Å, respectively. Therefore, the degree of sintering in case 2 was greater than that of case 1, because of the contribution of the framework formed by the YSZ nanoparticles. The two YSZ nanoparticles around nanoparticles No. 3 and No. 4 (Figure 10d−f) suppressed the growth of the neck, suggesting that they also prevented the nickel nanoparticles from coalescing. Furthermore, we performed the sintering simulations by using ten different initial configurations of Ni and YSZ nanoparticles in order to confirm the generality of the simulation. Figure 12a,b shows the cross-sectional snapshots of the Ni/YSZ multi-nanoparticle sintering process different from Figures 8−10. Here, case 1 and case 2 indicate the sintering process of nickel nanoparticles No. 1 and No. 3 in the large pore of YSZ framework and the sintering process of nickel nanoparticles No. 1 and No. 2 in the small pore of YSZ framework, respectively (Figure 12a,b). In this model, we observed that the neck of nickel nanoparticles No. 1 and No. 3 grew due to the sintering. However, the neck growth of nickel nanoparticles No.1 and No. 2 was suppressed by the YSZ framework. Figure 12c shows the change in the distance

Figure 8. Methodology for the Ni/YSZ multi-nanoparticle model. (a) Ni and YSZ nanoparticles were packed into a 100 Å × 100 Å × 100 Å simulation cell at random. (b) The Ni and YSZ multi-nanoparticle model was compressed until its porosity reached 40%, and the size of the simulation cell was 81.3 Å × 84.1 Å × 81.9 Å.

Figure 9. Time evolution of the Ni/YSZ multi-nanoparticle model. (a) The initial state of the Ni/YSZ multi-nanoparticles. (b) The coalescence of the nickel nanoparticles at 500 ps. Black circles indicate the coalescence of the nickel nanoparticles.

simulation at 1073 K was performed using this initial structure. The structure of the Ni/YSZ multi-nanoparticle model at 500 ps is shown in Figure 9b. The aggregation of the nickel nanoparticles was observed after the sintering simulation at 1073 K. This indicates that the nickel nanoparticles in porous 9669



Article

Figure 10. Cross-sectional snapshots of the Ni/YSZ multi-nanoparticle sintering simulation at 1073 K. Sintering processes in the large pore of YSZ framework at (a) t = 0 ps, (b) t = 11 ps, and (c) t = 500 ps, and in the small pore of YSZ framework at (d) t = 0 ps, (e) t = 7 ps, and (f) t = 500 ps. The simulation region is indicated by thick lines. For clarity, the region shown in the snapshots has been extended beyond the simulation region based on the periodic boundary condition.

realistic simulation. In the model, five nickel and five YSZ nanoparticles were located randomly in 100 Å × 100 Å × 100 Å simulation cells. The diameters of the nanoparticles were distributed from 10 to 30 Å, randomly. We also used the porosity of 40% in varying sized multi-nanoparticle sintering simulation. Figure 13a,b shows the cross-sectional snapshots of the sintering process in the large pore of YSZ framework of the varying sized Ni/YSZ multi-nanoparticle model at 0 and 500 ps, respectively. The cross sections of the sintering process in the small pore at 0 and 500 ps are shown in Figure 13c,d, respectively. Here, the sintering between nickel nanoparticles No. 1 and No. 2 in the large pore is referred to as case 1, and the sintering between No. 3 and No. 4 in the small pore is referred to as case 2. We observed that the sintering of the nickel nanoparticles No. 1 and No. 2 in the large pore occurred. The neck between them grew (Figure 13b). However, the neck between No. 3 and No. 4 in the small pore did not grow (Figure 13d). Figure 14 shows the change in the distance between the center of mass of the nickel nanoparticles for case 1 and case 2 in the varying sized multi-nanoparticle model. The decrease in the distance of case 1 was larger than that of case 2. It implied that the degree of the sintering between nickel nanoparticles No. 3 and No. 4 was smaller than that between nanoparticles No. 1 and No. 2. Thus, we conclude that the YSZ framework suppresses the sintering by disrupting the growth of the neck, which is consistent with the results of the equivalently sized Ni/YSZ multi-nanoparticle model. Therefore, our sintering simulation using the multi-nanoparticle model is effective for investigating the role of the YSZ framework in the sintering of nickel nanoparticles in the SOFC anode. Our Ni/YSZ multi-nanoparticle model produced more realistic results compared with those of the simplified model of two nanoparticles on the YSZ(111) surface used in section B. The method was successfully used to determine the effects of

Figure 11. Changes in the distance between the centers of mass of the nickel nanoparticles in the Figure 8 as a function of time. Case 1 indicates the distance between nickel nanoparticles in large pore (Figure 10a−c), and case 2 indicates the distance between nickel nanoparticles in small pore (Figure 10d−f).

between the center of mass of the nickel nanoparticles for case 1 and case 2 in the multi-nanoparticle model (Figure 12a,b). The decreases in the distance of case 1 and case 2 are 16.1 and 11.8 Å, respectively. It means that the sintering of No. 1 and No. 2 was inhibited by the small pore of the YSZ framework. The same phenomena are also observed in ten calculations. It shows that the repeatability is good in these calculations. Thus, we conclude that YSZ inhibits the sintering process of nickel nanoparticles. Hence, we propose that the YSZ nanoparticles form a supporting matrix in the anode, which prevents the nickel nanoparticles from coalescing. This was not observed in previous sintering simulations with two or three nanoparticles.16−24 Certainly, this phenomenon cannot be observed in 40-silicon nanoparticles in a straight chain without a substrate model.48 In addition, the sintering simulation of the varying sized multi-nanoparticle model was also performed for a more 9670



Article

Figure 14. Change in the distance between the center of the mass of the nickel nanoparticles for case 1 and case 2 in the varying sized multi-nanoparticle model. Here, the sintering between nickel nanoparticles No. 1 and No. 2 in the large pore is referred to as case 1, and the sintering between No. 3 and No. 4 in the small pore is referred to as case 2.

cannot be achieved in the previous sintering simulations on the conventional nanoparticles models.16−24,48

4. CONCLUSION We have developed an MD simulation method for investigating the sintering of nickel nanoparticles in the Ni/YSZ anode. The interatomic potential parameters for the sintering simulation represented the interaction between the nickel and YSZ framework well, and accurately reproduced the interaction energies at 108 points at the three sites on the YSZ(111) surface, when compared with the DFT calculations. The interatomic potential parameters were used in a sintering simulation of two nickel nanoparticles on a YSZ(111) surface, which showed that the two nickel nanoparticles approached and interpenetrated each other. However, previous sintering simulations cannot reproduce the experimentally observed results. Our multi-nanoparticle model for accurately simulating the anode structure reflected the effects of the porosity, the YSZ nanoparticle framework, and the ratio of nickel to YSZ nanoparticles in the Ni/YSZ SOFC anode, on the sintering of nickel nanoparticles. The Ni/YSZ multi-nanoparticle model was constructed by packing nickel and YSZ nanoparticles into a simulation cell, and the cell was compressed to achieve the desired porosity. The simulation indicated that the YSZ framework around the nickel nanoparticles prevented the nickel nanoparticles from coalescing. These results were not observed in the conventional sintering models. Therefore, our multi-nanoparticle model can effectively simulate the effects of the porosity, the YSZ nanoparticle framework, and the ratio of nickel to YSZ nanoparticles, on the sintering mechanism in the SOFC anode.

Figure 12. Cross-sectional snapshots of the second simulation of the Ni/YSZ multi-nanoparticle sintering simulation at 1073 K for (a) t = 0 ps and (b) t = 500 ps. The initial structure is different from that in Figure 8.

■

AUTHOR INFORMATION

Corresponding Author

*Tel: +81-22-795-6930. Fax: +81-22-795-6931. E-mail: [email protected].

Figure 13. Cross-sectional snapshots of the sintering process of the varying sized Ni/YSZ multi-nanoparticle model. Sintering processes in the large pore of YSZ framework at (a) 0 ps and (b) 500 ps, and in the small pore of YSZ framework at (c) 0 ps and (d) 500 ps, respectively.

Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Florian-Patrice, N. Electricity from Wood through the Combination of Gasification and Solid Oxide Fuel Cells; Ph.D. Thesis; ETH Zurich, 2008. (2) Larminie, J.; Dicks, A. Fuel Cell Systems Explained; Wiley: New York, 2000; pp 207−226.

the porosity of the SOFC anode, the YSZ nanoparticle framework, and the ratio of Ni to YSZ nanoparticles on the sintering of nickel nanoparticles, and could be used to predict the experimental properties of anodes. These advantages 9671



Article

(3) Jiang, S. P.; Chan, S. H. A Review of Anode Materials Development in Solid Oxide Fuel Cells. J. Mater. Sci. 2004, 39, 4405− 4439. (4) Singhal, S. C. Advances in Solid Oxide Fuel Cell Technology. Solid State Ionics 2000, 135, 305−313. (5) Ormerod, R. M. Solid Oxide Fuel Cells. Chem. Soc. Rev. 2003, 32, 17−28. (6) Singh, P.; Minh, N. Q. Solid Oxide Fuel Cells: Technology Status. Int. J. Appl. Ceram. Technol. 2004, 1, 5−15. (7) Singhal, S. C. Solid Oxide Fuel Cells. Electrochem. Soc. Interface 2007, 4, 41−44. (8) Zhu, W. Z.; Deevi, S. C. A Review on the Status of Anode Materials for Solid Oxide Fuel Cells. Mater. Sci. Eng., A 2003, 362, 228−239. (9) Atkinson, A.; Barnett, S.; Gorta, R. J.; Irvine, J. T. S.; McEvoy, A.; Mogensen, M.; Singhal, S. C.; Vohs, J. Advanced Anodes for HighTemperature Fuel Cells. Nat. Mater. 2004, 3, 17−27. (10) Jiang, S. P. Sintering Behavior of Ni/Y2O3-ZrO2 Cermet Electrodes of Solid Oxide Fuel Cells. J. Mater. Sci. 2003, 38, 3775− 3782. (11) Primdahl, S.; Sørensen, B. F.; Mogensen, M. Effect of Nickel Oxide/Yttria-Stabilized Zirconia Anode Precursor Sintering Temperature on the Properties of Solid Oxide Fuel Cells. J. Am. Ceram. Soc. 2000, 83, 489−494. (12) Simwonis, D.; Tietz, F.; Stöver, D. Nickel Coarsening in Annealed Ni/8YSZ Anode Substrates for Solid Oxide Fuel Cells. Solid State Ionics 2000, 132, 241−254. (13) Skarmoutsos, D.; Tietz, F.; Nikolopoulos, P. Structure Property Relationships of Ni/YSZ and Ni/(YSZ+TiO2) Cermets. Fuel Cells 2001, 1, 243−248. (14) Iwata, T. Characterization of Ni-YSZ Anode Degradation for Substrate-Type Solid Oxide Fuel Cells. J. Electrochem. Soc. 1996, 143, 1521−1525. (15) Minh, N. Q. Ceramic Fuel Cells. J. Am. Ceram. Soc. 1993, 76, 563−588. (16) Song, P.; Wen, D. Molecular Dynamics Simulation of the Sintering of Metallic Nanoparticles. J. Nanoparticle Res. 2010, 12, 823− 829. (17) Zhu, H.; Averback, R. S. Sintering Processes of Two Nanoparticles: A Study by Molecular Dynamics Simulations. Phil. Mag. Lett. 1996, 73, 27−33. (18) Arcidiacono, S.; Bieri, N. R.; Poulikakos, D.; Grigoropoulos, C. P. On the Coalescence of Gold Nanoparticles. Int. J. Multiphase Flow 2004, 30, 979−994. (19) Koparde, V. N.; Cummings, P. T. Molecular Dynamics Simulation of Titanium Dioxide Nanoparticle Sintering. J. Phys. Chem. B 2005, 109, 24280−24287. (20) Buesser, B.; Gröhn, A. J.; Pratsinis, S. E. Sintering Rate and Mechanism of TiO2 Nanoparticles by Molecular Dynamics. J. Phys. Chem. C 2011, 115, 11030−11035. (21) Raut, J. S.; Bhagat, R. B.; Fichthorn, K. A. Sintering of Aluminum Nanoparticles: A Molecular Dynamics Study. Nanostruct. Mater. 1998, 10, 837−851. (22) Zeng, P.; Zajac, S.; Clapp, P. C.; Rifkin, J. A. Nanoparticle Sintering Simulations. Mater. Sci. Eng., A 1998, 252, 301−308. (23) Miyamoto, A.; Yamauchi, R.; Kubo, M. Atomic Processes in the Deposition and Sintering of Ultrafine Metal Particles on MgO(001) as Investigated by Molecular Dynamics and Computer Graphics. Appl. Surf. Sci. 1994, 75, 51−57. (24) Kubo, M.; Stirling, A.; Miura, R.; Yamauchi, R.; Miyamoto, A. Molecular Dynamics Simulation for Ultrafine Gold Particles Deposited on Metal Oxides. Catal. Today 1997, 36, 143−151. (25) Kubo, M.; Miura, R.; Yamauchi, R.; Katagiri, M.; Vetrivel, R.; Broclawik, E.; Miyamoto, A. Atomic Control of Ultrafine Gold Particles on MgO(100) as Investigated by Molecular Dynamics and Computer Graphics. Jpn. J. Appl. Phys. 1995, 34, 6873−6877. (26) Kubo, M.; Oumi, Y.; Miura, R.; Stirling, A.; Miyamoto, A. Molecular Dynamics Study of Epitaxial Growth and Cluster Formation on MgO(001). AIChE J. 1997, 43, 2765−2772.

(27) Endou, A.; Teraishi, K.; Yajima, K.; Yoshizawa, K.; Ohashi, N.; Takami, S.; Kubo, M.; Miyamoto, A.; Broclawik, E. Potential Energy Surface and Dynamics of Pd/MgO(001) System as Investigated by Periodic Density Functional Calculations and Classical Molecular Dynamics Simulations. Jpn. J. Appl. Phys. 2000, 39, 4255−4260. (28) Brandon, N. P.; Brett, D. J. Engineering Porous Materials for Fuel Cell Applications. Philos. Trans. R. Soc. A 2006, 364, 147−159. (29) Jang, J. H.; Ryu, J. H.; Oh, S. M. Microstructure of Ni/YSZ Cermets According to Particle Size of Precursor Powders and Their Anodic Performances in SOFC. Ionics 2000, 6, 86−91. (30) Wang, F. H.; Guo, R. S.; Wei, Q. T.; Zhou, Y.; Li, H. L.; Li, S. L. Preparation and Properties of Ni/YSZ Anode by Coating Precipitation Method. Mater. Lett. 2004, 58, 3079−3083. (31) Yoshimoto, M.; Maeda, T.; Ohnishi, T.; Koinuma, H.; Ishiyama, O.; Shinohara, M.; Kubo, M.; Miura, R.; Miyamoto, A. Atomic-Scale Formation of Ultrasmooth Surfaces on Sapphire Substrates for HighQuality Thin-Film Fabrication. Appl. Phys. Lett. 1995, 67, 2615−2617. (32) Verlet, L. Computer ″Experiments″ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98−103. (33) Suzuki, K.; Kubo, M.; Oumi, Y.; Miura, R.; Takaba, H.; Fahmi, A.; Chatterjee, A.; Teraishi, K.; Miyamoto, A. Molecular Dynamics Simulation of Enhanced Oxygen Ion Diffusion in Strained YttriaStabilized Zirconia. Appl. Phys. Lett. 1998, 73, 1502−1504. (34) Girifalco, L. A.; Weizer, V. G. Application of the Morse Potential Function to Cubic Metals. Phys. Rev. 1959, 114, 687−690. (35) Delley, B. An All-Electron Numerical Method for Solving the Local Density Functional for Polyatomic Molecules. J. Chem. Phys. 1990, 92, 508−517. (36) Delley, B. From Molecules to Solids with the DMol3 Approach. J. Chem. Phys. 2000, 113, 7756−7764. (37) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (38) Perdew, J. P.; Burke, K.; Ernzerhof, M. Erratum: Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1997, 78, 1396−1396. (39) Ballabio, G.; Bernasconi, M.; Pietrucci, F.; Serra, S. Ab Initio Study of Yttria-Stabilized Cubic Zirconia Surfaces. Phys. Rev. B 2004, 70, 075417. (40) Grau-Crespo, R.; Hernandez, N. C.; Sanz, J. F.; De Leeuw, N. H. Theoretical Investigation of the Deposition of Cu, Ag, and Au Atoms on the ZrO2(111) Surface. J. Phys. Chem. C 2007, 111, 10448− 10454. (41) Xia, X.; Oldman, R.; Catlow, C. R. A. Computational Modeling Study of Bulk and Surface of Yttria-Stabilized Cubic Zirconia. Chem. Mater. 2009, 21, 3576−3585. (42) Levenberg, K. A Method for the Solution of Certain Non-Linear Problems in Least Squares. Q. Appl. Math. 1944, 2, 164−168. (43) Marquardt, D. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. SIAM J. Appl. Math. 1963, 11, 431−441. (44) Madsen, K.; Nielsen, H. B.; Tingleff, O. Methods for Non-Linear Least Squares Problems (2nd). [Online] Informatics and Mathematical Modelling, Technical University of Denmark, DTU, 2004; 4, pp 24− 29. http://www2.imm.dtu.dk/pubdb/views/publication_details. php?id=3215. (accessed September, 2011). (45) Lide, D. R. CRC Handbook of Chemistry and Physics, 90th ed.; CRC press: Boca Raton, 2009; pp 12-11−12. (46) Güvenç, Z. B.; Jellinek, J. Surface Melting in Ni55. Z. Phys. D 1993, 26, 304−306. (47) Koparde, V. N.; Cummings, P. T. Sintering of Titanium Dioxide Nanoparticles: A Comparison between Molecular Dynamics and Phenomenological Modeling. J. Nanoparticle. Res. 2008, 10, 1169− 1182. (48) Hawa, T.; Zachariah, M. R. Molecular Dynamics Simulation and Continuum Modeling of Straight-Chain Aggregate Sintering: Development of A Phenomenological Scaling Law. Phys. Rev. B 2007, 76, 054109.

9672


Molecular Dynamics Simulation of Ni Nanoparticles Sintering Process

Recommend Documents