**1. Introduction**

#### **1.1 Solid state ionics and SOFC**

Solid ionic compounds normally form insulating compounds because each ion usually has a closed-shell electronic configuration, which prevents delocalization of electrons and thus electronic conduction in the system. In an appropriate temperature regime, charge can be transported in these solid materials, however, by the motion of highly mobile ionic species (e.g. Li+, Na+, O2-, etc.) (Fig. 1 Top) (Chu et al., 2006). This phenomenon is termed "ionic conductivity". Such ionic compounds exhibit liquid-like conductivities whilst still in the solid state, i.e. at temperatures are well below their melting points (Hull, 2004). Exploiting this unusual property, considerable progress has been made in recent years in fuel cells and batteries, which are promising key technologies to meet our rising energy and environmental needs.

Using "superionic conductivity" (e.g. typically ~ 10-1 Scm-1) (Fig. 1 Top), high temperature fuel cells, e.g. solid oxide fuel cells (SOFC) operate between 600 and 1200 C. Compared to batteries which are electrochemical energy storage devices, SOFC's are energy conversion device that produces electricity by electrochemically combining fuel (e.g. H2, CH4, CO, hydrocarbon etc.) and oxidant (e.g. air, O2) chemical reactions at anode and cathode respectively (Fig. 1 Middle), with high thermodynamic efficiency and low environmental impact (Andersson et al., 2010; Etsell et al., 1970; Lashtabeg et al., 2006). A SOFC is a complex ceramic-solid-based electrochemical device consisting of three main components: an anode, electrolyte and cathode. The main component requiring fast ion transport and diffusion is the solid electrolyte. With external electron transfer, oxygen is reduced to oxygen ions (O2−) at the cathode side through the oxygen reduction reaction (ORR): ½ O2 +2e<sup>−</sup> O2−. Then the oxygen ions are incorporated into the oxygen ion conductor, the SOFC electrolyte. The oxygen ions (O2−) are transported through the electrolyte, but electrons are not. An ideal electrolyte in a SOFC is an electronic insulator, but an ionic conductor, and therefore permits only the oxygen ions to pass through to the anode. At the anode, the oxygen ions will combine with the fuel (e.g. H2) to form water: H2 + O2- H2O + 2e−.

Through this oxidation process, electrons are released and lead via an external circuit to the cathode where the ORR occurs again, and thus electrical power is generated (Andersson et

The Roles of Classical Molecular Dynamics Simulation in Solid Oxide Fuel Cells 343

al., 2010). In short, the overall electrochemical cell reaction is built upon the chemical potential difference between the cathode and anode sides, which is given by the Nernst

To maximize this superionic conductivity in a SOFC, alloying zirconia (ZrO2) with various metal oxides (e.g. Y2O3, Sc2O3, La2O3, MgO, CaO, etc) is a plausible approach (Lashtabeg et al., 2006). In certain regimes, alloying can stabilize the highly symmetric cubic-fluorite phase of zirconia, which consequently facilitates the ionic conductivity through the introduction of oxygen vacancies as the host Zr4+ cations are replaced by dopant aliovalent cations. Above all, yttria-stabilized zirconia, YSZ (i.e. the Zr1−*<sup>x</sup>*Y*x*O2−*x/*2 system, with *x/*2 being the Y2O3 dopant concentration), is the most common choice for the electrolyte (Fig. 1 Top), due to its good oxide ion conductivity over a wide range of oxygen partial pressures, its stability under oxidizing and reducing conditions, and its good high-temperature mechanical properties (Ralph et al., 2001). This versatility ultimately arises from atomic defects (i.e. oxide ion vacancies) in the cubic zirconia crystalline lattices. The vacancies, through the coupled interactions among the vacancies and ions in these different alloys can dramatically affect the structural, thermal, mechanical and electrical properties of the system. Furthermore, the optimum ionic conductivity of each alloy varies with synthesis route and sintering conditions due to the resultant diverse local morphologies and microstructures (Badwal et al., 1992; Butz et al., 2006; Chen et al., 2002; Fukui et al., 2004; Ioffe et al., 1978; Korte et al., 2008; Zhang et al., 2007; Zhu et al., 2005). Thus a variety of different morphologies and microstructures of YSZ can be found in experiments (Badwal et al., 1992; Butz et al., 2006; Chen et al., 2002; Cheng et al., 2011; Etsell et al., 1970; Fukui et al., 2004; Ioffe et al., 1978; Korte et al., 2008; Lashtabeg et al., 2006; Zhang et a., 2007; Zhu et al., 2005). Non-cubic crystalline phases, grain boundaries, and disordered lattices of amorphous features can commonly be found. To understand how the local microstructures and system morphologies affect the ionic conductivity and degradation of the electrolyte (YSZ) upon SOFC operation, theoretical simulation can provide detailed *in situ* atomistic information

To increase the basic understanding of the scientific phenomena that underlie current experimental findings, which could most dramatically affect engineering design, atomistic modeling based on quantum mechanics (e.g. first-principles or *ab initio* methods), molecular dynamics and Monte Carlo simulation (Fig. 2) are particularly useful. They provide relevant predictions of crystal structure, energetics, and vibrational frequencies through detailed atomistic descriptions that complement the standard analytical continuum model of ion diffusion at the macroscopic scale (Andersson et al., 2010; Cheng et al., 2011). Furthermore, various fundamental physical processes and chemical reactions can be described via quantum mechanics. The most accurate quantum mechanical atomistic descriptions, firstprinciples or *ab initio* methods, however are limited to small system size (less than about 1000 atoms) and short times (less than a few nanoseconds). For long-time-scales and larger length scales in atomistic modeling, recent sophisticated studies use kinetic Monte Carlo (KMC) (Gatewood et al., 2011; Lau et al., 2008, 2009; Pornprasertsuk et al., 2007; Turner et al., 2010; Wang et al., 2010, 2011) simulations. KMC can treat the longest time scales for chemical reactions, because reactants are not followed to products in time. Instead, chemical reactions are stochastically chosen to occur at the assigned rates. This is important for SOFC modeling, because charge-separating electrochemical reactions tend to have slower rates. KMC has modeled the ionic transport in a SOFC under various operating conditions

equation in the simplest approach.

that is difficult to obtain experimentally.

Fig. 1. (Top) Temperature dependence of the best solid ionic conductors of Ag+, Na+, Li+, and O2− ions (adapted from Chu et al., 2006). The ionic conductivities fall within ranges which are highlighted by the shaded areas. (Middle) A simple atomistic model of SOFC that consists of three basic components: anode, cathode and electrolyte. (Bottom) The basic properties of solid electrolyte of SOFC can be simulated using standard atomistic modeling.

Fig. 1. (Top) Temperature dependence of the best solid ionic conductors of Ag+, Na+, Li+, and O2− ions (adapted from Chu et al., 2006). The ionic conductivities fall within ranges which are highlighted by the shaded areas. (Middle) A simple atomistic model of SOFC that consists of three basic components: anode, cathode and electrolyte. (Bottom) The basic properties of solid electrolyte of SOFC can be simulated using standard atomistic modeling.

al., 2010). In short, the overall electrochemical cell reaction is built upon the chemical potential difference between the cathode and anode sides, which is given by the Nernst equation in the simplest approach.

To maximize this superionic conductivity in a SOFC, alloying zirconia (ZrO2) with various metal oxides (e.g. Y2O3, Sc2O3, La2O3, MgO, CaO, etc) is a plausible approach (Lashtabeg et al., 2006). In certain regimes, alloying can stabilize the highly symmetric cubic-fluorite phase of zirconia, which consequently facilitates the ionic conductivity through the introduction of oxygen vacancies as the host Zr4+ cations are replaced by dopant aliovalent cations. Above all, yttria-stabilized zirconia, YSZ (i.e. the Zr1−*<sup>x</sup>*Y*x*O2−*x/*2 system, with *x/*2 being the Y2O3 dopant concentration), is the most common choice for the electrolyte (Fig. 1 Top), due to its good oxide ion conductivity over a wide range of oxygen partial pressures, its stability under oxidizing and reducing conditions, and its good high-temperature mechanical properties (Ralph et al., 2001). This versatility ultimately arises from atomic defects (i.e. oxide ion vacancies) in the cubic zirconia crystalline lattices. The vacancies, through the coupled interactions among the vacancies and ions in these different alloys can dramatically affect the structural, thermal, mechanical and electrical properties of the system. Furthermore, the optimum ionic conductivity of each alloy varies with synthesis route and sintering conditions due to the resultant diverse local morphologies and microstructures (Badwal et al., 1992; Butz et al., 2006; Chen et al., 2002; Fukui et al., 2004; Ioffe et al., 1978; Korte et al., 2008; Zhang et al., 2007; Zhu et al., 2005). Thus a variety of different morphologies and microstructures of YSZ can be found in experiments (Badwal et al., 1992; Butz et al., 2006; Chen et al., 2002; Cheng et al., 2011; Etsell et al., 1970; Fukui et al., 2004; Ioffe et al., 1978; Korte et al., 2008; Lashtabeg et al., 2006; Zhang et a., 2007; Zhu et al., 2005). Non-cubic crystalline phases, grain boundaries, and disordered lattices of amorphous features can commonly be found. To understand how the local microstructures and system morphologies affect the ionic conductivity and degradation of the electrolyte (YSZ) upon SOFC operation, theoretical simulation can provide detailed *in situ* atomistic information that is difficult to obtain experimentally.

To increase the basic understanding of the scientific phenomena that underlie current experimental findings, which could most dramatically affect engineering design, atomistic modeling based on quantum mechanics (e.g. first-principles or *ab initio* methods), molecular dynamics and Monte Carlo simulation (Fig. 2) are particularly useful. They provide relevant predictions of crystal structure, energetics, and vibrational frequencies through detailed atomistic descriptions that complement the standard analytical continuum model of ion diffusion at the macroscopic scale (Andersson et al., 2010; Cheng et al., 2011). Furthermore, various fundamental physical processes and chemical reactions can be described via quantum mechanics. The most accurate quantum mechanical atomistic descriptions, firstprinciples or *ab initio* methods, however are limited to small system size (less than about 1000 atoms) and short times (less than a few nanoseconds). For long-time-scales and larger length scales in atomistic modeling, recent sophisticated studies use kinetic Monte Carlo (KMC) (Gatewood et al., 2011; Lau et al., 2008, 2009; Pornprasertsuk et al., 2007; Turner et al., 2010; Wang et al., 2010, 2011) simulations. KMC can treat the longest time scales for chemical reactions, because reactants are not followed to products in time. Instead, chemical reactions are stochastically chosen to occur at the assigned rates. This is important for SOFC modeling, because charge-separating electrochemical reactions tend to have slower rates. KMC has modeled the ionic transport in a SOFC under various operating conditions

The Roles of Classical Molecular Dynamics Simulation in Solid Oxide Fuel Cells 345

intensively studied using classical MD simulations (Devanathan et al., 2006; Fisher et al., 1998, 1999; Khan et al., 1998; Kilo et al., 2003; Lau et al., 2011; Li et al., 1995; Okazaki et al., 1994; Sawaguchi et al., 2000; Schelling et al., 2001; Shimojo et al., 1992; van Duin et al., 2008;

In crystalline solids, ionic conductivity is fundamentally different from electronic conductivity. Electronic conduction in a metal, for example, occurs on a three-dimensional array of ion cores whose excess valence electrons have dissociated to form a continuous "sea of free electrons" partially filling the electronic bands around the Fermi-level. Because the electron has a small mass, its de Broglie wavelength is large and therefore quantum mechanical effects force the electrons into those bands. As ions are much heavier than electrons, their motion is far less governed by quantum mechanics. Below the typical atomic vibrational frequencies (< 100 GHz), ionic motion is best described by thermally activated

The dynamics of mobile ions in a disordered inorganic ionic conductor (e.g. amorphous YSZ) is clearly a complex multi-particle problem. Unlike a perfect crystal, the potentialenergy landscape experienced by an ion in a disordered solid is irregular and contains a distribution of depths and barrier heights. The ions' interaction with the dynamic atomic lattice network is fundamental: first, because the lattice supplies a persistent disordered potential energy landscape for the mobile anions and second, because the local fluctuations of the lattice atoms promote anionic jumps. Additional multi-particle behavior stems from the interaction among the mobile ions and vacancies. All these distinct coupled interactions contribute to the complete theoretical description of ionic conductivity in amorphous solids (Dyre et al., 2009; Lammert et al., 2010). A complete analytical microscopic theory is not available by now and, due to the complexity of the problem, would be extremely difficult to formulate without some basic understanding at the atomistic level. The direct approach to ionic conduction in the YSZ electrolyte of the SOFC is via molecular dynamics (MD) simulations. From the time-evolved atomic trajectories detailed microscopic information about the underlying mechanisms is available. Understanding them would make theoretical

It is possible to model a few thousands to millions of particles in classical MD using phenomenological interatomic and intermolecular potentials. They are obtained by using the phenomenological approach of selecting a parameterized mathematical form for the interaction between atoms, and fitting its unknown parameters to various experimental or higher-level theoretical (e.g. *ab initio* quantum mechanics simulation) properties. In general, the flexibility, accuracy, transferability, and computational efficiency of the interatomic

For the YSZ electrolyte in a SOFC, the simplest, relevant interatomic potentials are those commonly used to describe rigid ionic compounds (Lewis et al., 1985). Under the rigid ion model, the potential energy is a simple function of the distance between the ions. It consists

potentials each have to be carefully considered (Frenkel & Smit, 1996).

hopping between (usually) charge-compensating sites (Dyre et al., 2009).

Yamamura et al., 1999) in the past few years.

**1.2 Basic justification of classical MD** 

prediction possible.

**2.1 Interatomic potentials** 

**2. The model of classical MD for YSZ** 

Fig. 2. Computational methods at different length and time scales to model a SOFC material adopted from Cheng et.al. 2011. (Adapted figure reproduced from Cheng et.al. 2011. Copyright 2011 RSC Publishing.)

(Gatewood et al., 2011; Lau et al., 2008, 2009; Pornprasertsuk et al., 2007; Turner et al., 2010; Wang et al., 2010, 2011). Specifically, KMC can probe SOFC performance by simultaneously capturing various reaction pathways of electrochemical and physicochemical reactions in the electrolyte and at the three-phase boundary (TPB), i.e. the interface where the gas reactants, electrolyte, and electrode meet. The electrical current through the YSZ is simulated in direct current (dc) and alternating current (ac), as electrochemical impedance spectra under various operating conditions using a minimal set of uniform chemical reaction rates on an assumed cubic YSZ lattice via KMC (Gatewood et al., 2011; Lau et al., 2008, 2009; Pornprasertsuk et al., 2007; Turner et al., 2010; Wang et al., 2010, 2011). Despite the robustness of KMC simulation, this approach is based on a rigid lattice gas model and can not predict the experimentally observed ionic conductivity maximum as a function of Y2O3 dopant concentration (Hull, 2009). Beyond the kinetics driven atomic motion as implemented in minimal KMC models, the complex dynamics of real lattices and the realtime multi-particle ion-vacancy interactions at a finite temperature can be computed 'on-thefly' in simulation. However simple or complex, the ionic conductivity of any atomistic model of solid YSZ is, the conductivity can accurately be derived from standard equilibrium classical MD simulation (Frenkel & Smit, 1996). Thus to explore and understand the nature of unique ionic conductivity in the solid electrolytes of SOFC's, YSZ solids have been intensively studied using classical MD simulations (Devanathan et al., 2006; Fisher et al., 1998, 1999; Khan et al., 1998; Kilo et al., 2003; Lau et al., 2011; Li et al., 1995; Okazaki et al., 1994; Sawaguchi et al., 2000; Schelling et al., 2001; Shimojo et al., 1992; van Duin et al., 2008; Yamamura et al., 1999) in the past few years.
