Introduction
Pure zirconia (ZrO2) has a monoclinic structure up to a temperature of about 1446 K, where it changes to the tetragonal modification. For temperatures higher than 2643 K zirconia adopts a cubic fluorite structure [1]. However, the high temperature phases can be partially or completely stabilized at room temperature by doping with aliovalent oxides such as yttria (Y2O3), calcia (CaO), or magnesia (MgO). Yttria is the most commonly used dopant for stabilizing the cubic phase of zirconia; a fully (cubic) stabilized zirconia is obtained with a Y2O3-content of >7 mol% [2], while a Y2O3-content of about 2-6 mol% gives a partially stabilized zirconia. Cubic stabilized zirconia has improved mechanical and thermal properties such as high strength, toughness, and thermal-shock resistance. Furthermore the addition of substitutional cation (e.g., Ca2+, Mg2+, Y3+), which have lower valency than zirconium ion (Zr4+), induces the generation of oxygen vacancies for charge compensation. For example, the substitution of Zr4+ with Y3+ causes the negative net charge in the lattice; for every mole of yttria incorporated into the zirconia lattice, the charge neutrality condition is kept by forming an oxygen vacancy. In a Kroger-Vink notation, it is written as follows
(1)Where $Y^{'}_{Zr}$ means Y in the Zr site with the apparent negative charge, and $\ddot{V}_{O}$ is the vacancy in the oxygen site with double positive charge. $O^x_O$ is the lattice oxygen, i.e., oxygen in the oxygen site with net charge of zero. The existence of vacancies on the oxygen site gives rise to the high ionic conductivity of YSZ. Oxygen is transported by hopping through its vacancy sites (vacancy diffusion mechanism). The concentration of oxygen vacancy is determined by the concentration of the dopant.
Areas of Application of YSZ
One major application of YSZ is in solid fuel cell (SOFC) technology where it is used as electrolyte material due to its oxygen ion conductivity, high chemical and cryastallographic stability, and low electronic conductivity [3]. The SOFC comprises of two porous electrodes, an anode and a cathode, separated by a dense electrolyte and it is usually operated at high temperatures (~1000oC). The function of the electrolyte is to transport the oxygen ions from the cathode to the anode where oxidation of the fuel by the ions occurs (as depicted figure 1), and also to block the electrons produced at the anode from passing through the cell to the cathode [4]. The typical thickness of the electrolyte material is 100 -150 micron. For an effective operation of a fuel cell, with the electrolyte within this range of thickness, the ionic conductivity of the electrolyte material must be higher than 0.1 Scm-1. The ionic conductivity determines, to a large extent, the available power and operating temperature of a SOFC. The ionic conductivity of 8 mol% YSZ at 1000oC is 0.18 Scm-1 [5]
Figure 1. (a) Scanning electron microscopy (SEM) image of a typical SOFC showing the anode, electrolyte, and cathode layers; (b) Schematics of a SOFC showing the mode of operation.
YSZ is also widely used as thermal barrier coating (TBC) layer material for gas turbine engine (see figure 2) because of its low thermal conductivity and relatively high strength, toughness, thermal shock resistanc [6]. The low thermal conductivity of bulk YSZ results from the low intrinsic thermal conductivity of zirconia and phonon scattering defects (vacancies) introduced by the addition of yttria. Both the yttrium cations and the oxygen vacancies are effective phonon scattering sites the thermal conductivity is decreased as the yttria content is increased. In practice, a yttria concentration in the range of 6 to 8 wt.% is generally used.
Figure 2. Cross section SEM micrograph of a TBC superimposed onto a schematic showing the temperature reduction provided by the TBC. (b) SEM image a typical showing YSZ TBC
However, the high oxygen ion conductivity of YSZ is not desirable in its application as a TBC material [6]. Thermally grown oxide layer (TGO) usually forms at the boundary between the TBC layer and the bond-coat as a result of high diffusivity of ionic oxygen through the YSZ layer. This leads to spallation failure of the TBC when the TGO grows to a certain thickness.
The Anomaly of Oxygen Ion Conductivity in YSZ
Experimental studies have shown that the ionic conductivity of YSZ increases with increasing dopant concentration due to the increase of the ionic charge carrier (oxygen vacancy). However, the conductivity does not increase monotonically; it decreases after reaching the maximum at a certain dopant concentration as shown in figure 3a. The maximum conductivity is obtained at ~10 mol%. Furthermore, this trend in ionic conductivity/diffusivity is commonly observed in a wide range of oxide materials having the same fluorite structure [9]. A detailed understanding of oxygen mobility in stabilized zirconia is essential for predicting and designing material systems best suited for the applications mentioned above.
Figure 3. (a) Oxygen self-diffusivity (Do) as a function of Y2O3 mole fraction in YSZ; (a) YSZ cubic fluorite type structure.
Simulation of Oxygen Ion Diffusion in YSZ
Several theoretical studies have been carried out by different researchers in order to understand the underlying mechanism of oxygen-ion diffusion in YSZ. Murray et al. used lattice statistics to study the effect of dopant content on ionic conductivity in yttria-doped ceria [8]. They considered different cation distributions (i.e., different numbers of host and dopant ions) among the nearest neighbor cation sites of an exchanging oxygen vacancy-oxygen ion pair. While they observed the expected variation of ionic conductivity with dopant content, they attributed the considerable discrepancy in their calculated values of the ionic conductivity and experiment to vacancy-dopant association effects. In another study, Shimojo and Okazaki conducted a comprehensive molecular dynamics investigation to address this problem [10]. They analyzed their simulation data (for YSZ) in terms of the oxygen density and potential energy distributions in tetrahedra, the comers of which are filled by different combinations of the dopant and host cations. Their results show that oxygen migration preferentially takes place across tetrahedral edges that have two Zr ions at their corners (Zr-Zr edge) and almost no oxygen migration takes place across Y-Y edges. Meyer and Nicoloso have also used percolation model [10] to investigate this phenomenon. In their study, they assume simple model potentials for the dopant-vacancy interaction that mimic the three different possibilities that have been reported in the literature, namely, dopant-vacancy nearest neighbor attraction, dopant-vacancy nearest neighbor repulsion, and a modified barrier for oxygen vacancy hops across edges of the cation tetrahedra that contain dopant ions. In agreement with the Shimojo et al. conclusions, they found that only the dopant-affected migration barrier model reproduced experimental results for the ionic conductivity of oxides with the fluorite structure.
Nevertheless, a major limitation of the above mentioned methods in simulating the oxygen diffusivity in YSZ is that they do not reproduce the correct relative energies for the three different polymorphs of pure zirconia [11]. In contrast, density functional theory (DFT) methods are able to reproduce the experimentally observed relative energetics of the different phases of pure zirconia, making DFT the most credible theoretical methods for the simulation of ionic conductivity in YSZ. Since the DFT approach requires no empirical information, it does not suffer from some of the limitations inherent in using less reliable empirical potentials (fit to limited data) to conduct molecular dynamics or atomistic Monte Carlo simulations of these oxides.
DFT Approach to Oxygen Ion Diffusion Simulation in YSZ
Oxygen ion diffusion simulation using the DFT methods involves the calculation of the oxygen migration energies for different configurations, using first-principles methods and then using these calculated migration barriers in kinetic Monte Carlo (KMC) model to calculate the long-time self-diffusivity of oxygen as a function of temperature and dopant concentration. This approach has an advantage over molecular dynamics simulations, which are restricted to relatively short times and, therefore, insufficiently sample the statistical distribution of local environments [11]. However, because of the relatively large computational expense associated with the first-principles calculations, certain assumptions are made in the calculation of the oxygen migration barrier in YSZ, in order to reduce the computational load. Due to these assumptions, some possible interactions which may affect the migration energy are excluded in the DFT calculations. For example, none of the DFT approaches includes vacancy-vacancy interactions, and the interactions with extended neighbors beyond the two adjacent tetrahedra. While the significance of these excluded interactions is not known a priori, the good match between calculated results for oxygen diffusivity (using DFT) and experiment suggests that these assumptions are reasonable.
Calculating Oxygen Ion Migration Barriers using DFT Methods
Several approaches have been employed by various researchers to calculate the oxygen migration energy in YSZ using DFT. Three of such approaches will be discussed extensively in this write-up. For easy comparison, only the approaches that made use of Vienna ab initio simulation package (VASP) to perform their calculations, adopting the generalized gradient approximation (GGA), is reported here.
Oxygen Ion Migration Barriers along Different Crystallographic Directions
Eichler et al., in their approach, considered the diffusion barriers for oxygen-ion vacancies at several locations in the supercell, which provided a connected path for vacancy migration to positions equivalent to the initial position. The barriers between several vacancy positions were calculated until the vacancy arrives finally at a position equivalent to the starting point. The energy barriers for vacancy diffusion are expected to be different for vacancies occupying different positions. Vacancies have different local environments depending on the type of cation (Zr4+ or Y3+) occupying the corners of the two adjacent tetrahedra containing the oxygen ion – vacancy pair. The dissimilarity in the barriers, on one side, is due to the greater ionic radius of yttrium and on the other side, due to the smaller positive charge on the Y3+ (than on Zr4+) which is less attractive to the doubly negative charged oxygen ions.
Figure 4. Sketch of the optimized YSZ structures showing the distortion effect of a vacancy. Arrows indicate the direction of distortion; circles and squares denote the positions of Y dopant and O vacancies, respectively. The numbers describe the distortion (in Ǻ) with respect to the ideal bulk.
Two possible pathways for the vacancy diffusion in the supercell were considered: within the (x,y) plane [011], or perpendicular to it [001] (see figure 5). In the first case the pathway is longer—the vacancy encounters six positions until it is again in a position equivalent to the starting configuration, compared to only four inequivalent positions during diffusion in [001] direction(see table 1). However, since the pathway parallel to the [110] direction is symmetrical, only four positions are inequivalent and the computational effort remains the same.
Table I. Binding energies Eb (eV) for the vacancy positions that are encountered along [110] and [001] diffusion pathways in Y-doped zirconia (compare figure 5).
Fgure 5. Interaction energies for vacancy diffusion in Y-doped zirconia along the [011] solid line) and the [001] Direction (dashed line) as a function of a reaction path coordinated s. Calculated values are indicated by points, smooth lines are drawn on the basis of these points and corresponding tangents.
On the basis of the configurations along the pathways , a smaller barrier is expected for diffusion within the (x,y) plane, since in the second case the vacancy has to come very close to the Y ions - a configuration which is energetically unfavorable. Diffusion along [001] is dominated by a huge barrier between position B and C, where the oxygen ion (the vacancy) has to pass between the two Y ions. This high barrier was attributed to the bigger size of the Y3+ ionic radius so that the gap between the two ions is smaller. The migration energy for oxygen diffusion was then assumed to be the difference between the highest energy along the path and the lowest energy. From their computations, oxygen diffusion in YSZ (2.9 mol% Y2O3) proceed along [110] direction with an activation energy of 0.87 eV.
Oxygen Ion Migration Barriers Considering Cations at Two Tetrahedra Corners
The “crystallographic direction” approach used by Eichler et al. can serve as an approximation only because the migration energy is obtained from only a specified path. In an attempt to obtain more accurate oxygen migration energies, Krishnamurthy et al. [11] have employed a different approach wherein the energy barriers for oxygen migration across a particular type of cation edge (Zr-Zr, Zr-Y, or Y-Y, as shown in figure 6) were calculated for different configurations using DFT. The absolute energies corresponding to two different configurations were calculated. These include a configuration where all cations and oxygen ions and a vacancy occupy equilibrium lattice positions and a configuration corresponding to the transition state for the migration of a single oxygen vacancy. The difference between these energies is the activation energy for oxygen migration.
Figure 6. Schematics (a) YSZ structure showing an outline of two adjacent tetrahedral; (b) two adjacent tetrahedral showing the relative positions of cations, diffusing oxygen ion, oxygen vacancy, and the saddle-point plane.
Figure 7. Illustration of the migration energy barrier. The middle point is the energy corresponding to the saddle point where three atoms (diffusing oxygen ion and two cations at tetrahedral corner) align in the same plane.
To calculate the energy barrier for oxygen migration across a particular type of cation edge, attention was focused on one particular vacancy-oxygen ion pair and the edge of the cation tetrahedra through which the oxygen ion and vacancy interchange. Y ions, other than those that constitute this cation edge, and oxygen vacancies, other than the one considered, are placed as far away from the hopping oxygen vacancy as possible to minimize their effect on the migration energy. The saddle point in the energy landscape is approximated by the split-vacancy configuration corresponding to the oxygen ion occupying a location midway between the equilibrium positions of the oxygen vacancy pair, as shown in figure 7. This configuration was realized through the use of a constraint via the method of Lagrange multipliers. Car-Parrinello molecular dynamics was used for efficient structural relaxation of YSZ. The energy barrier for oxygen migration is calculated as the total energy difference between a fully relaxed equilibrium configuration and a saddle point configuration arranged from the relaxed geometry using the constraint described above. This calculation is repeated for oxygen vacancy-ion exchange across the three different types of cation edges that exist in YSZ, namely, Zr-Zr, Zr-Y, and Y-Y. The corresponding migration energies are listed in Table II.
Table II. Activation Energies, for oxygen migration across
A-cation/B-cation edges, AB
A simulation cell with 2 X 2 X 2 conventional fluorite unit cells was used for the modeling. This simulation cell contained 64 oxygen lattice sites and 32 cation lattice sites. All cation sites were initially filled with Zr ions and anion sites with O ions. Two, four, or six of the Zr ions were replaced with Y ions and. correspondingly, one, two, or three oxygen ions were removed from the anion sublattice, leaving vacancies. These correspond to approximate yttria mole tractions of 3%, 6%, or 9%, respectively. All the oxygen ion migration barriers were calculated using the cubic structure containing ~9 mol% yttria. This restriction implies that the activation barrier for oxygen migration was unaffected by the yttria concentration and yttrium ion distribution.
Oxygen Ion Migration Barriers Considering all Cations at the Tetrahedra Corners
Another approach for obtaining the oxygen ion migration barriers in YSZ using GGA-based DFT simulations was reported by Pornprasertsuk, et al. [12]. Their approach is very closely related to that used by Krishnamurthy et al., but a broader range of diffusion barriers in the vicinity of a diffusing oxygen ion and an oxygen vacancy was accounted for. They considered the influence of all six cations in two adjacent tetrahedra containing diffusing oxygen ion and oxygen vacancy (see figure 8) to the migration energy barrier, unlike in the approach used by Krishnamurthy et al., where only two cations at the saddle-point plane were considered. With DFT calculations, they were able to obtain the binding energy of $\ddot{V}_{O}Y^{'}_{Zr}$ and $Y^{'}_{Zr}Y^{'}_{Zr}$ for different cation configurations (as shown in table III), which quantified the influence of defect interactions on the oxygen ion diffusion mechanism.
Table III. List of calculated binding energies of $\ddot{V}_{O}Y^{'}_{Zr}$ and $Y^{'}_{Zr}Y^{'}_{Zr}$ from the first nearest-neighbor (1st NN) to the fourth nearest-neighbor (4th NN) interactions.
The supercell used in all the calculations consisted of 30 Zr, 6 Y, and 69 O atoms corresponding to 8.3 mol % YSZ. The saddle point was assumed to be located on the plane that has the following properties: (i) it is perpendicular to the shortest path and (ii) it contains the other two cations (Figure 8). To calculate the energy at the saddle point, all three ions (two cations and one oxygen ion) were allowed to relax only within this plane while the other atoms were allowed to relax in all directions. The corresponding migration energy is the energy difference between the saddle-point energy and the initial state energy (Figure 7). It was assumed that the barrier depends largely on the sites’ first nearest neighbors to the diffusion center.
Even though this last approach accounted for broader range of diffusion barriers in the vicinity of a diffusing oxygen ion and an oxygen vacancy, yet the interactions of oxygen ion vacancies and the other interactions beyond the components in the two tetrahedra were neglected in order to reduce the computational load. To further reduce the computational complexity, they assumed that there were no more than three Y ions in the six cation sites.
Performing Kinetic Monte Carlo Using the Calculated Migration Barriers
In order to determine the long-time self-diffusivity of oxygen as a function of temperature and dopant concentration, kinetic Monte Carlo (KMC) simulation was performed (in each case) based on the DFT calculated barriers. The KMC technique attempts to capture the effects of rare atomic processes that directly contribute to changes in macroscopic properties. It reflects the macroscopic behavior of the system by statistically averaging the atomic processes. KMC was used to simulate an activated random-walk process in a randomly distributed landscape of vacancy and Y atoms.
Y cations and oxygen vacancies are introduced randomly in the periodic supercells used in the KMC simulation, depending on the concentration of the Y ion. Oxygen ions are assumed to move on the anion sublattice by exchanging with oxygen vacancies on nearest neighbor sites on the anion sublattice. For an activated random-walk of vacancies in the oxygen sublattice, the root mean square displacement,$\langle r^{2}(t)\rangle$ can be determined. It varies linearly with time as demonstrated in the example in figure 9.
Figure 9. Mean square displacement of the oxygen vacancies is shown as a function of time in (ZrO2,)1-x(Y2O3)x/2( VO)x/2 for different yttrium ion fractions, x, at I400K [11].
(2)Where $D_{v}$ is vacancy diffusivity; and $t$ is time taken to reach equilibrium.
(3)and
(5)Where $v^{o}$ is frequency factor representing the lattice vibrational frequency which was assumed to be a constant at all dopant concentrations; $\Delta{E}_{M}$ is the first-principles calculated energy barrier for oxygen ion diffusion in YSZ; ${k}_{B}$ is the Boltzmann’s constant; and $T$ is temperature. The $\exp\left[ {-\Delta{E}_{M}\over{k}_{B}T}\right]$ term is the probability that an oxygen ion jumps from its lattice site onto an adjacent vacancy.
The KMC procedure consists of a series of moves. Each move consists of the following steps:
- Identify all possible events from the current configuration. In the YSZ structure (fluorite-type oxide), there are six possible events for each oxygen vacancy (i.e., six nearest-neighbor oxygen sites where the vacancy can hop into)), except in the case where any one of the first nearest neighbors is also a vacancy.
- Obtain the jump rates ( $v_{i}$) for each of the events; rate is proportional to $\exp\left[ {-\Delta{E}_{M}\over{k}_{B}T}\right]$ .
- Generate a pseudorandom number$\gamma$ between 0 and 1.
- Advance the time of each step ($\Delta{t}$) by ${\ln(\gamma)\over{\Sigma{v_{i}}}$
- Choose one of the events depending on the random number $\gamma$, consistent with the relative rate of all the events.
- Reconfigure the system according to the chosen event.
- Update and record the new position of the vacancy and time.
The vacancy diffusivity was then calculated using the expression
(6)The oxygen self-diffusivity, Do is obtained from Dv by considering the jump balance, hence
(7)Where cv is the vacancy concentration.
Oxygen diffusivity is an Arrhenius function of temperature for all yttrium ion concentrations. Therefore it can be written as follows
(8)Thus, the activation energy$Q$ for oxygen diffusion in YSZ can be determined from the slopes of the plot of logarithm of the oxygen diffusivity in YSZ versus the inverse temperature for several Y2O3, mole fractions, as shown in figure 10.
figure 10. Arrhenius plot of the log of oxygen diffusivity vs temperature inverse for different yttria mole fractions [11].
Figure 11. Activation energy for oxygen self-diffusion as a function of yttria mole fraction. Obtained using DFT calculated oxygen ion migration barriers from the (a) second approach: Krishnamurthy et al. [12]; (b) third approach: Pornprasertsuk, et al. [13].
Summary
Ionic conductivity in YSZ can be studied using DFT-calculated migration energy barriers in a kinetic Monte Carlo model to simulate the long-time self-diffusivity of oxygen as a function of temperature and dopant concentration. The inability to include all possible interactions in the DFT calculations leads to under estimation of the activation energy for oxygen diffusion in YSZ. For example, none of the DFT approaches included vacancy-vacancy interactions, and the interactions with extended neighbors beyond the two adjacent tetrahedra. Oxygen vacancies diffuse by jumping across tetrahedral edges; the types of cations bounding the edge determine the magnitude of the hopping barrier. Oxygen diffusivity increases with increasing yttria mole fraction, before attaining a maximum value at ~10 mol% yttria, after which it decreases. The DFT/KMC simulation, employing the approaches examined here, provides results with qualitative agreement relative to the experimental observations, thereby improving the understanding of the mechanism for oxygen diffusion in YSZ.
References
[1] A. Eichler, “Tetragonal Y-doped zirconia: Structure and ion conductivity,” Phy. Rev. B Physics (2001)
[2] C. Pascual, J.R. Jurado, and P. Duran, “Electrical behaviour of doped-yttria stabilized zirconia ceramic materials,” J. Mater. Sci. (1983)
[3] T. Kawada, H. Yokokawa, “Materials and Characterization of Solid Oxide Fuel Cell,” Key Eng. Mater., (1997).
[4] J. Goodenough, Y. Huang. J. Power Sources, (2007).
[5] C. Chena, C. Varanasia, J. P. Fellner, Electrical properties of heterogeneously doped yttria stabilized zirconia, J. Power Sources, (2005).
[6] N.P. Padture, M. Gell, E.H. Jordan, “Thermal barrier coatings for gas-turbine engine applications,” Science, (2002).Science, (2002).
[7] D, W. Stickler and W. G, Carlson, J. Am. Ceram. Soc. (1964)
[8] A. D. Murray, G.E. Murch, C.R.A Catlow, "A new hybrid scheme of computer simulation based on Hades and Monte Carlo: Application to ionic conductivity in Y3+ doped CeO2 " Solid State Ionics, (1986)
[9] F. Shimojo and H. Okazaki, J. Phys. Soc. Jpn., (1992)
[10] M. Meyer, N. Nicoloso, V. Jaenisch, “Percolation model for the anomalous conductivity of fluorite-related oxides" Phys. Rev. B,(1997)
[11] R. Krishnamurthy, Y.-G. Yoon, D. J. Srolovitz, and R. Car, "Oxygen Diffusion in Yttria-Stabilized Zirconia: A New Simulation Model" J. Am. Ceram. Soc., (2004).
[12] R. Pornprasertsuk, P. Ramanarayanan, C. B. Musgrave, F. B. Prinz, "Predicting ionic conductivity of solid oxide fuel cell electrolyte from first principles" J. Appl. Phys., (2005)
