Molecular Dynamics Study of Ni/YSZ Systems Based on Improved Interatomic Model

Porous nickel/yttria-stabilized zirconia (YSZ) composite is a common choice as an anode material for solid oxide fuel cells (SOFCs). Because of its complicated structure, there are several processes occurring in the anode that affect total performance of SOFCs. For instance, oxygen ion diffusion provides oxygen transport to triple phase boundaries (TPBs) where fuel is oxidized. Another important process is nickel sintering, which causes decreasing of the active TPB length and thus leads to degradation of SOFCs. Understanding the mechanism of each individual phenomenon in Ni/YSZ anodes can contribute to improving total performance and stability of SOFCs.

Atomistic simulations can provide a deep insight into the processes taking place in Ni/YSZ anodes. While the density functional theory method is one of the most rigorous atomistic approaches used in materials science, its applicability is limited to small systems of hundreds of atoms. Larger systems can be studied by molecular dynamics (MD) and Monte-Carlo simulations, which are based on empirical potentials. However, most of the potentials for YSZ are based on simple pairwise interactions and designed in the way to reproduce specific bulk properties. Unfortunately, their validity for studying surface and interface properties of Ni/YSZ systems is usually only presumed [1], which can lead to questionable results for such processes as nickel sintering and oxygen surface diffusion.

Thus, in the present paper we develop a DFT-based model that purposefully aims for studying a wide range of processes in Ni/YSZ systems, including oxygen diffusion in bulk and surface regions of YSZ, yttrium surface segregation, and nickel sintering. Our model combines three types of interatomic interactions. The first one is based on the dipole model [2] and was parameterized for bulk and surface YSZ regions in our previous work [3]. It reproduces with a good accuracy not only bulk but also surface properties of YSZ. The second describes nickel by the embedded atom method [4]. The third is for interactions between Ni and YSZ and is described by a pairwise Buckingham potential, which was parameterized individually for different Ni/YSZ interfaces.

Having confirmed that the model well reproduces the bulk properties of YSZ, we employed the developed model to demonstrate that it reproduces with good accuracy the surface structures and bulk conductivity of YSZ. This allowed us to use the model to further study oxygen ion diffusion near (111) and (110) YSZ surfaces where YSZ structure can be very much distorted [5], and thus oxygen diffusion can be affected. Obtained results suggest suppressed oxygen diffusion along (111) because of the lower oxygen vacancy concentration near the surface. In contrast, oxygen vacancies tend to migrate to the (110) surface, which, together with activation of jumps along <110> directions, facilitates diffusion along the surface. This indicates that it can be easier for oxygen ions to migrate to TPBs near the (111) surfaces from bulk regions than along the surface, while near the (110) surface ions prefers to migrate along the surface.

We also demonstrate how the model performs for the nickel sintering problem, based on the structures consisting of two Ni particles being 30 Å in diameter and placed on a YSZ slab. We studied the effect of binding energy between the particles and the slab by considering different interfaces at the contact areas and using corresponding interactions. Simulations at 1250 K indicate that a relatively low binding energy of 0.5 J/m², which corresponds to the Ni(111)/YSZ(111) interface, allows the nanoparticles to slide over the YSZ support and coalesce. When the binding energy is above 1.0 J/m², the sliding is restricted and particle migration may occur only by means of nickel atom diffusion on nickel surface.

[1] H. Lee et al., Acta Materialia 58 (2010) 2197.

[2] P. Tangney et al., J. Chem. Phys. 117 (2002) 8898.

[3] A. Iskandarov et al., J. Phys.: Condens. Matter 25 (2015) 015005.

[4] Y. Mishin et al., Phys. Rev. B 59 (1999) 3393.

[5] G. Ballabio et al., Phys. Rev. B. 70 (2004) 075417.