- Email: [email protected]
Molecular dynamics investigations of grain boundary phenomena in cubic zirconia
Computational Materials Science 14 (1999) 177±184
Molecular dynamics investigations of grain boundary phenomena in cubic zirconia Craig A.J. Fisher *, Hideaki Matsubara Synergy Ceramics Laboratory, Japan Fine Ceramics Center, 2-4-1 Mutsuno, Atsuta-ku, Nagoya 456-8587, Japan
Abstract Yttria stabilized zirconia (YSZ) is a fast oxide ion conducting ceramic with the cubic ¯uorite structure that is used in a number of applications, including solid oxide fuel cells (SOFCs). A molecular dynamics (MD) study has been performed on symmetrical tilt grain boundaries with the R5 (3 1 0)/[0 0 1] h 36.9° misorientation to investigate the structure and dynamics of interfaces in this technologically important material. Simulations were performed on systems of 1920 atoms at constant temperatures up to 2673 K. Atomic interactions were described by a simple pair potential model of the Buckingham form. Structural relaxation produced an open structure corresponding to the introduction of a row of Schottky defects adjacent and parallel to the interface. Oxygen diusion along the boundary was observed at high temperature, even without vacancies in the bulk being introduced explicitly by aliovalent doping. However, the diusion rate was lower than that in single crystals of 8 mol% Y2 O3 stabilized zirconia. Further simulations demonstrated that interfaces between perfect zirconia crystals are sources of resistance in these ionically conducting systems. Ó 1999 Elsevier Science B.V. All rights reserved. Keywords: Atomistic simulation; Molecular dynamics; Zirconia; YSZ; Grain boundary; Diusion; Oxide ion
1. Introduction As well as stabilizing the cubic form of zirconia (ZrO2 ) to lower temperatures, addition of lower valence dopants such as yttrium (Y3 ) introduces compensating oxygen vacancies which enable the ceramic to display fast oxide ion conductivity at elevated temperatures. This makes stabilized zirconia a favoured candidate for use as the electrolyte in oxygen gas separators and solid oxide fuel cells (SOFCs) [1]. Improvements in conductivity have been achieved over the past decade mainly by reducing the amount of impurities in the starting
* Correspondence author. Tel.: 81 52 871 8815; fax: 81 52 871 8826; e-mail: ®[email protected].
powders and thus the amount of insulating phases formed at the grain boundaries during sintering [2]. Further improvements in conductivity are desirable, however, particularly in order to lower the temperature of operation (and consequently costcompetitiveness) of these environmentally friendly technologies. The purpose of the present study is to apply computer simulation techniques to evaluate the eect of grain boundaries and hence grain size on the oxide ion conductivity of zirconia ceramics. Molecular dynamics (MD) has only recently been applied to the study of oxide ion diusion in zirconia systems [3±6]. These simulations have all been concerned with single crystals and have successfully reproduced the maximum in oxide ion conductivity observed at around 8±10 mol% Y2 O3
0927-0256/99/$ ± see front matter Ó 1999 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 9 8 ) 0 0 1 0 4 - 9
178
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
using only simple eective pair potentials. While two-body potentials are sucient for reproducing the high symmetry cubic and tetragonal phases, more sophisticated models are required to successfully reproduce phenomena such as relative phase stabilities and transformations between zirconia's various polymorphs [7,8]. However, the orders-of-magnitude increase in calculation time required by these models makes them impractical for long time-scale MD simulations at the present time. Static and dynamic simulation techniques have been used widely to model grain boundaries in metallic systems based on the coincidence site lattice (CSL) concept [9]. However, grain boundary studies in ceramic materials have largely been con®ned to static (0 K) simulations of simple rocksalt structured compounds [10]. An early attempt to model grain boundaries in cubic and tetragonal zirconia has been reported by Bingham et al. [2]. The structures of two twist boundaries and one tilt boundary were calculated at 0 K using two sets of empirically derived pair potentials. Calculations of the defect energies and segregation pro®les of oxide vacancies and yttrium dopant ions at the different interfaces led to an estimate of the ratio between grain boundary resistivity and bulk resistivity, qgb /qbulk , of 50±120. In this paper we report the ®rst use of the MD technique to investigate oxygen diusion along grain boundaries of zirconia at ®nite temperatures. The ®rst part of the study examines the grain boundary structure and oxygen motion along the R5 (3 1 0)/[0 0 1] grain boundary in pure zirconia. This is followed by an examination of the eect of the same boundary on the overall ionic conductivity of a bicrystal of 8 mol% YSZ. 2. Simulation method The MD code MOLDY used in this study utilizes a modi®cation of the Beeman algorithm to eciently integrate Newton's equations of motion [11]. The atomic interactions are described by a potential function which divides the forces into long-range interactions (described by Coulomb's Law and summated by the Ewald method) and
short-range interactions treated by a pairwise function of the Buckingham form ! qi qj ÿr Cij Aij exp 1 ÿ 6; /ij r r r qij where qi and qj are the charges of ions i and j, respectively, r is the distance between them and Aij , qij and Cij are the parameters particular to each ion±ion interaction. The exponential term corresponds to electron cloud overlap and the rÿ6 term any attractive dispersion or van der Waal's forces. As with most MD simulations, all ions were treated as rigid spheres, since the incorporation of electronic polarizability, for example via the Shell Model, would render the calculations prohibitively time-expensive. The potential parameters used in the present study were taken from Lewis and Catlow [12]. Of several potential sets found in the literature, these were found to most accurately reproduce the experimentally determined unit cell geometry. This gives us some degree of con®dence that the potential model accurately represents atomic interactions in the real material. Fig. 1 shows the (3 1 0) plane passing through three unit cells of ZrO2 with the ideal cubic ¯uorite structure. The initial con®guration for the grain boundary simulations was formed according to CSL theory by taking two ZrO2 crystals and tilting them through an angle of 36.9° until their (3 1 0) planes coincided. The resulting R5 (3 1 0)/[0 0 1] symmetrical tilt interface is illustrated in Fig. 2(a). This con®guration was contained in a simulation box of size 8 CSL unit cell lengths in the x direction, two in the y direction and four ¯uorite unit cell lengths in the z direction, giving a total of 1920 ions. Application of 3D periodic boundary conditions created a second boundary at the box edges parallel to the central boundary, but oriented in the opposite direction. For simulations of 8 mol% YSZ (8YSZ), 94 Y3 ions were substituted for Zr4 ions and 47 O2ÿ ions removed to maintain zero net charge. Simulations were initially performed on systems of pure ZrO2 without Y3 substitution in order to clearly distinguish between oxide ion diusion at the grain boundary and that due to extrinsic
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
179
new constraints at each temperature before commencing the simulation runs of 60 000 time steps. Particle positions were periodically recorded during the course of each simulation to allow calculation of transport properties, in particular the mean square displacement (msd). The msd of a species numbering N particles is de®ned as msd
Fig. 1. Ideal ¯uorite structure of cubic zirconia with the (3 1 0) plane shown passing through three unit cells.
vacancies introduced by doping. Although the cubic phase is only stable above 2640 K in pure zirconia, in our model we maintained this structure at lower temperatures for ready comparison with yttria-stabilized systems. Systems of grain boundaries between crystals of 8YSZ were simulated and the overall oxide ion diusion compared with that in single crystals of the same composition. The time step used in all simulations was 2 fs. The R5 interface was equilibrated commencing from the ideal CSL structure at a constant pressure of 101.3 kPa at temperatures between 1273 and 2673 K. The particle velocities were rescaled every ten time steps during equilibration to maintain constant system temperature. Once equilibrium was reached, the average box dimensions were calculated for a few thousand time steps. The simulation box was subsequently held ®xed at these values during the actual simulation runs so that the simulations were performed in the microcanonical (NVE) ensemble. A second equilibration was performed for 5000 time steps under these
N D E 1X 2 rn t0 t ÿ rn t0 ; N nÿ1
2
where the angular brackets are taken to indicate an ensemble average over all possible initial times, t0 , and ri (t) is the position of ion i at time t. The msd gives a measure of the diusion undergone by a particular species during the course of the simulation run, since according to the Einstein equation, the gradient of a plot of msd versus time is proportional to the diusion coecient. A second set of simulations was performed on systems of 8YSZ. A run of 60 000 time steps was performed on a single crystal containing 1873 atoms at 1273 K to calculate the msd and diusion coecient of this material. Simulations of the R5 interface in 8YSZ were then performed at 1273 K in a manner identical to that used for the pure ZrO2 simulations and the results compared to those for the single crystal.
3. Results and discussion 3.1. Grain boundary structure Equilibration of the initial R5 CSL interface produced the same structure at temperatures of 1273, 1673, 1873, 2073 and 2673 K. The boundary energies at each temperature were also very similar, varying between 2.7 and 2.9 J/m2 . Fig. 3 shows the average positions of ions around the grain boundary at 1273 K. The initial R5 structure appears to have been maintained, but in fact the grain boundary has expanded by one layer perpendicular to the interface in order to reduce the strong repulsive interactions between ions in close proximity across the boundary in Fig. 2(a). This expansion is equivalent to the introduction of a layer of Schottky defects parallel to the grain
180
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
Fig. 2. (a) R5 (3 1 0)/[0 0 1] grain boundary in ZrO2 constructed according to coincident site lattice (CSL) theory. (b) Schematic illustration of introduction of Schottky vacancy layer to initial con®guration during relaxation of the grain boundary.
boundary. At the same time, half the Zr4 and O2ÿ ions opposite the Schottky layer (Fig. 2(b)) move across the new boundary to occupy some of the vacancies. As well as reducing the density of sites at the grain boundary, the vacancies allow ions to move easily between corresponding sites either side of the interface. The ¯uorite structure of cubic zirconia can be considered to consist of an fcc network of Zr4 cations with the O2ÿ anions ®lling the tetrahedral interstices. It is therefore not surprising that the structure is similar to that predicted for fcc metals such as copper [13] (except for the introduction of the vacancy layer to reduce the repulsive forces in the ionic system). In contrast, rock-salt structured oxides such as NiO [14], where both the cation and anion subarrays are in fcc arrangements, have an even higher density of like-charged ions across the interface in the initial R5 con®guration. In these
systems the boundary relaxes by sliding one crystal parallel to the boundary to the other so that anions and cations face each other across the interface, so that the mirror symmetry is destroyed. Our predicted structure for the R5 interface is very similar to recently published high-resolution transmission electron microscopy (HRTEM) images of symmetric tilt boundaries in cubic zirconia [15]. 3.2. Grain boundary diusion vs. temperature Simulations of pure ZrO2 at temperatures between 1273 and 2673 K found noticeable oxide ion diusion occurring at the grain boundaries. Since no extrinsically introduced vacancies from doping were present, this must be due to the Schottky vacancies intrinsic in the relaxed R5 structure. In order to distinguish between ions at the grain boundary and those in the bulk, the labelling
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
181
Fig. 3. Average positions of ions at R5 (3 1 0)/[0 0 1] grain boundary in ZrO2 at 1273 K.
scheme shown in Fig. 4 was used for msd calculations. Since the grain boundary structure contains sites with dierent coordination environments (and thus energies), these calculations must necessarily be averages of the behaviour of ions of a particular species at the grain boundary. Although somewhat arbitrary, our labelling scheme appears intuitive since it re¯ects the boundary symmetry, and is useful for highlighting
the dierence between behaviour of ions within the grain boundary region from those in the bulk. At 1273 K, only some of the oxide ions at the boundary travelled far from their starting positions. Although movement between adjacent sites across the grain boundary was relatively easy, the ions tended to quickly return to their starting positions without contributing to net diusion. As shown in Fig. 5(a), however, after about 30 ps
182
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
Fig. 4. Labelling scheme used to distinguish ions at the grain boundary (II) and second layer from the grain boundary (I) from bulk ions for msd calculations.
Fig. 5. Mean square displacements (msds) of ions at various temperatures, T, according to the labelling scheme shown in Fig. 4. (a) T 1273 K, (b) T 1673 K, (c) T 1873 K, (d) T 2073 K.
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
there is a de®nite increase in the msd of O±II ions, indicating that some diusion has occurred. There is also a small amount of diusion in the O±I ions, showing that the vacancies initially at the grain boundary have penetrated into the second layer from the boundary. When the temperature is increased to 1673 K (Fig. 5(b)), the diusion of O± II, and to some extent O±I, ions becomes more distinct, and even after short time periods net diusion has taken place. Analysis of the components of the msds in each direction, however, shows that most of the diusion is due to motion across the boundary (the x direction) rather than along it. Fig. 5(c) shows the msds at 1873 K. The increase and subsequent decrease in the O±II msd shows that after diusing a certain distance, ions tended to return to their starting positions or remain where they were. Analysis of the trajectories of the ions con®rmed this conclusion. The ions appeared to be trapped within the immediate region about the grain boundary and were unable to penetrate further into the ¯uorite crystal or traverse the larger distances between sites along the interface. This suggests that there is a signi®cant energy barrier to movement from the grain boundary into the bulk and for long-scale diusion along the R5 grain boundary. A further increase in temperature to 2073 K enables the oxide ions to diuse rapidly in the grain boundary and also the bulk of the crystal, as shown in Fig. 5(d). Although most of the diusion is associated with ions starting at the grain boundary, some diusion of O±I and O (bulk) ions also becomes apparent. The msds of ions at 2673 K were very similar in form to those at 2073 K except that much more rapid diusion of O±II ions took place. In both cases, diusion of O±II ions in the x direction was slightly greater than in the y or z directions, indicating that diusion along the grain boundary is not as favourable as that across the boundary from one crystal to the other. Despite the very high temperatures simulated, no cation diusion was observed. This is not surprising considering the very short time scales involved, and is consistent with the low cation diusion coecients determined experimentally. This result can also explain why very long anneal
183
times are required before the equilibrium structure is reached in real materials. Nevertheless, the msd plots reveal that cations at the boundary oscillate around their mean positions with greater amplitude than those in the bulk. This is obviously a result of the greater amount of space between ions at the grain boundary. 3.3. Diusion in single crystals and bicrystals of YSZ A plot of the msds of a single crystal and bicrystal of 8YSZ is given in Fig. 6. This plot illustrates clearly the decrease in conductivity caused by introducing grain boundaries. The single crystal msd also shows the change in slope ®rst reported by Brinkman et above about 0.4 A al. [5], postulated as being due to a transition from neighbour±neighbour diusion to long-scale ion diusion. The bulk tracer diusion coecient, D, was calculated according to the Einstein relation from the slope of the single crystal msd curve to be 1.1 ´ 10ÿ11 m2 /s. From the Nernst±Einstein equation, this is equivalent to an oxide ion conductivity of 3.0 S/m. If we assume a Haven ratio similar to that in a single fcc crystal, the corresponding conductivity for the bicrystal containing two R5 grain boundaries is half this value. In other words, the presence of the grain boundaries has reduced the conductivity of the system by a factor of two. In real materials the volume of the bulk crystal is
Fig. 6. Mean square displacements (msds) of ions in a single crystal (SC) and R5 bicrystal (GB) of 8 mol% YSZ at 1273 K.
184
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
much larger than that of the grain boundary regions, so this eect is expected to be greatly reduced at typical grain sizes (>1 lm). However, in nano grain-sized zirconia the volume of grain boundaries will be comparable to that of the bulk crystals. It has recently been suggested that nanosized grains are a promising means of increasing the conductivity of 8YSZ [16], based on the expectation that oxygen diusion along the boundaries will be more rapid than in the bulk. However, our simulations suggest that this is not likely to be the case, concurring with the limited experimental evidence available so far [16]. Despite the lower density and hence more ``free space'' at the grain boundary, the decrease in symmetry and increase in disorder at the boundary are not as conducive to ion transport as the equi-energy sites in the bulk. 4. Conclusions MD simulations of pure ZrO2 bicrystals containing R5 (3 1 0)/[0 0 1] grain boundaries have shown that the interface contains a number of intrinsic vacancies which contribute to grain boundary diusion at elevated temperatures. On addition of yttrium ions and removal of the corresponding number of bulk oxide ions, however, the grain boundaries were found to decrease the overall oxide ion conductivity in 8 mol% YSZ. Although this is not expected to have a strong eect at typical grain sizes, it is expected to lead to lower conductivities in materials containing nano-sized grains, since diusion of oxide ion vacancies along the boundary is less rapid than that in the bulk.
Acknowledgements This work has been carried out as part of the Synergy Ceramics Project under the Industrial Science and Technology Frontier (ISTF) Program promoted by AIST, MITI, Japan. Under this program, part of the work has been funded through NEDO. The authors are members of the Joint Research Consortium of Synergy Ceramics.
References [1] N.Q. Minh, J. Am. Ceram. Soc. 76 (1993) 563. [2] D. Bingham, P.W. Tasker, A.N. Cormack, Phil. Mag. A 60 (1989) 1. [3] F. Shimojo, T. Okabe, F. Tachibana, M. Kobayashi, H. Okazaki, J. Phys. Soc. Japan 61 (1992) 2848. [4] X. Li, B. Hafskjold, J. Phys.: Condens. Matter 7 (1995) 1255. [5] H.W. Brinkman, W.J. Briels, H. Verweij, Chem. Phys. Lett. 247 (1995) 386. [6] H. Okazaki, H. Suzuki, K. Ihata, Phys. Lett. A 188 (1994) 291. [7] M. Wilson, U. Sch onberger, M.W. Finnis, Phys. Rev. B 54 (1996) 9147. [8] A.F. Kohan, G. Ceder, Comp. Mater. Sci. 8 (1997) 142. [9] H.F. Fischmeister, Journal de Physique C4 (1985) 3. [10] D.M. Duy, J. Phys. C: Solid State Phys. 19 (1986) 4393. [11] K.D. Refson, MOLDY 2.11, MOLecular DYnamics, Version 2.11, University of Oxford, 1996. [12] G.V. Lewis, C.R.A. Catlow, J. Phys. C: Solid State Phys. 18 (1985) 1149. [13] G.H. Bishop, Jr., R.J. Harrison, T. Kwok, S. Yip, J. Appl. Phys. 53 (1982) 5596. [14] D.M. Duy, P.W. Tasker, Phil. Mag. A 47 (1983) 817. [15] K.L. Merkle, G.-R. Bai, Z. Li, C.-Y. Song, L.J. Thompson, Phys. Stat. Sol. A 166 (1998) 73. [16] R. Vassen, D. St over, L.G.J. de Haart, M. Cappadonia, Brit. Ceram. Proc. 56 (1996) 35.
Molecular dynamics investigations of grain boundary phenomena in cubic zirconia Craig A.J. Fisher *, Hideaki Matsubara Synergy Ceramics Laboratory, Japan Fine Ceramics Center, 2-4-1 Mutsuno, Atsuta-ku, Nagoya 456-8587, Japan
Abstract Yttria stabilized zirconia (YSZ) is a fast oxide ion conducting ceramic with the cubic ¯uorite structure that is used in a number of applications, including solid oxide fuel cells (SOFCs). A molecular dynamics (MD) study has been performed on symmetrical tilt grain boundaries with the R5 (3 1 0)/[0 0 1] h 36.9° misorientation to investigate the structure and dynamics of interfaces in this technologically important material. Simulations were performed on systems of 1920 atoms at constant temperatures up to 2673 K. Atomic interactions were described by a simple pair potential model of the Buckingham form. Structural relaxation produced an open structure corresponding to the introduction of a row of Schottky defects adjacent and parallel to the interface. Oxygen diusion along the boundary was observed at high temperature, even without vacancies in the bulk being introduced explicitly by aliovalent doping. However, the diusion rate was lower than that in single crystals of 8 mol% Y2 O3 stabilized zirconia. Further simulations demonstrated that interfaces between perfect zirconia crystals are sources of resistance in these ionically conducting systems. Ó 1999 Elsevier Science B.V. All rights reserved. Keywords: Atomistic simulation; Molecular dynamics; Zirconia; YSZ; Grain boundary; Diusion; Oxide ion
1. Introduction As well as stabilizing the cubic form of zirconia (ZrO2 ) to lower temperatures, addition of lower valence dopants such as yttrium (Y3 ) introduces compensating oxygen vacancies which enable the ceramic to display fast oxide ion conductivity at elevated temperatures. This makes stabilized zirconia a favoured candidate for use as the electrolyte in oxygen gas separators and solid oxide fuel cells (SOFCs) [1]. Improvements in conductivity have been achieved over the past decade mainly by reducing the amount of impurities in the starting
* Correspondence author. Tel.: 81 52 871 8815; fax: 81 52 871 8826; e-mail: ®[email protected].
powders and thus the amount of insulating phases formed at the grain boundaries during sintering [2]. Further improvements in conductivity are desirable, however, particularly in order to lower the temperature of operation (and consequently costcompetitiveness) of these environmentally friendly technologies. The purpose of the present study is to apply computer simulation techniques to evaluate the eect of grain boundaries and hence grain size on the oxide ion conductivity of zirconia ceramics. Molecular dynamics (MD) has only recently been applied to the study of oxide ion diusion in zirconia systems [3±6]. These simulations have all been concerned with single crystals and have successfully reproduced the maximum in oxide ion conductivity observed at around 8±10 mol% Y2 O3
0927-0256/99/$ ± see front matter Ó 1999 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 9 8 ) 0 0 1 0 4 - 9
178
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
using only simple eective pair potentials. While two-body potentials are sucient for reproducing the high symmetry cubic and tetragonal phases, more sophisticated models are required to successfully reproduce phenomena such as relative phase stabilities and transformations between zirconia's various polymorphs [7,8]. However, the orders-of-magnitude increase in calculation time required by these models makes them impractical for long time-scale MD simulations at the present time. Static and dynamic simulation techniques have been used widely to model grain boundaries in metallic systems based on the coincidence site lattice (CSL) concept [9]. However, grain boundary studies in ceramic materials have largely been con®ned to static (0 K) simulations of simple rocksalt structured compounds [10]. An early attempt to model grain boundaries in cubic and tetragonal zirconia has been reported by Bingham et al. [2]. The structures of two twist boundaries and one tilt boundary were calculated at 0 K using two sets of empirically derived pair potentials. Calculations of the defect energies and segregation pro®les of oxide vacancies and yttrium dopant ions at the different interfaces led to an estimate of the ratio between grain boundary resistivity and bulk resistivity, qgb /qbulk , of 50±120. In this paper we report the ®rst use of the MD technique to investigate oxygen diusion along grain boundaries of zirconia at ®nite temperatures. The ®rst part of the study examines the grain boundary structure and oxygen motion along the R5 (3 1 0)/[0 0 1] grain boundary in pure zirconia. This is followed by an examination of the eect of the same boundary on the overall ionic conductivity of a bicrystal of 8 mol% YSZ. 2. Simulation method The MD code MOLDY used in this study utilizes a modi®cation of the Beeman algorithm to eciently integrate Newton's equations of motion [11]. The atomic interactions are described by a potential function which divides the forces into long-range interactions (described by Coulomb's Law and summated by the Ewald method) and
short-range interactions treated by a pairwise function of the Buckingham form ! qi qj ÿr Cij Aij exp 1 ÿ 6; /ij r r r qij where qi and qj are the charges of ions i and j, respectively, r is the distance between them and Aij , qij and Cij are the parameters particular to each ion±ion interaction. The exponential term corresponds to electron cloud overlap and the rÿ6 term any attractive dispersion or van der Waal's forces. As with most MD simulations, all ions were treated as rigid spheres, since the incorporation of electronic polarizability, for example via the Shell Model, would render the calculations prohibitively time-expensive. The potential parameters used in the present study were taken from Lewis and Catlow [12]. Of several potential sets found in the literature, these were found to most accurately reproduce the experimentally determined unit cell geometry. This gives us some degree of con®dence that the potential model accurately represents atomic interactions in the real material. Fig. 1 shows the (3 1 0) plane passing through three unit cells of ZrO2 with the ideal cubic ¯uorite structure. The initial con®guration for the grain boundary simulations was formed according to CSL theory by taking two ZrO2 crystals and tilting them through an angle of 36.9° until their (3 1 0) planes coincided. The resulting R5 (3 1 0)/[0 0 1] symmetrical tilt interface is illustrated in Fig. 2(a). This con®guration was contained in a simulation box of size 8 CSL unit cell lengths in the x direction, two in the y direction and four ¯uorite unit cell lengths in the z direction, giving a total of 1920 ions. Application of 3D periodic boundary conditions created a second boundary at the box edges parallel to the central boundary, but oriented in the opposite direction. For simulations of 8 mol% YSZ (8YSZ), 94 Y3 ions were substituted for Zr4 ions and 47 O2ÿ ions removed to maintain zero net charge. Simulations were initially performed on systems of pure ZrO2 without Y3 substitution in order to clearly distinguish between oxide ion diusion at the grain boundary and that due to extrinsic
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
179
new constraints at each temperature before commencing the simulation runs of 60 000 time steps. Particle positions were periodically recorded during the course of each simulation to allow calculation of transport properties, in particular the mean square displacement (msd). The msd of a species numbering N particles is de®ned as msd
Fig. 1. Ideal ¯uorite structure of cubic zirconia with the (3 1 0) plane shown passing through three unit cells.
vacancies introduced by doping. Although the cubic phase is only stable above 2640 K in pure zirconia, in our model we maintained this structure at lower temperatures for ready comparison with yttria-stabilized systems. Systems of grain boundaries between crystals of 8YSZ were simulated and the overall oxide ion diusion compared with that in single crystals of the same composition. The time step used in all simulations was 2 fs. The R5 interface was equilibrated commencing from the ideal CSL structure at a constant pressure of 101.3 kPa at temperatures between 1273 and 2673 K. The particle velocities were rescaled every ten time steps during equilibration to maintain constant system temperature. Once equilibrium was reached, the average box dimensions were calculated for a few thousand time steps. The simulation box was subsequently held ®xed at these values during the actual simulation runs so that the simulations were performed in the microcanonical (NVE) ensemble. A second equilibration was performed for 5000 time steps under these
N D E 1X 2 rn t0 t ÿ rn t0 ; N nÿ1
2
where the angular brackets are taken to indicate an ensemble average over all possible initial times, t0 , and ri (t) is the position of ion i at time t. The msd gives a measure of the diusion undergone by a particular species during the course of the simulation run, since according to the Einstein equation, the gradient of a plot of msd versus time is proportional to the diusion coecient. A second set of simulations was performed on systems of 8YSZ. A run of 60 000 time steps was performed on a single crystal containing 1873 atoms at 1273 K to calculate the msd and diusion coecient of this material. Simulations of the R5 interface in 8YSZ were then performed at 1273 K in a manner identical to that used for the pure ZrO2 simulations and the results compared to those for the single crystal.
3. Results and discussion 3.1. Grain boundary structure Equilibration of the initial R5 CSL interface produced the same structure at temperatures of 1273, 1673, 1873, 2073 and 2673 K. The boundary energies at each temperature were also very similar, varying between 2.7 and 2.9 J/m2 . Fig. 3 shows the average positions of ions around the grain boundary at 1273 K. The initial R5 structure appears to have been maintained, but in fact the grain boundary has expanded by one layer perpendicular to the interface in order to reduce the strong repulsive interactions between ions in close proximity across the boundary in Fig. 2(a). This expansion is equivalent to the introduction of a layer of Schottky defects parallel to the grain
180
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
Fig. 2. (a) R5 (3 1 0)/[0 0 1] grain boundary in ZrO2 constructed according to coincident site lattice (CSL) theory. (b) Schematic illustration of introduction of Schottky vacancy layer to initial con®guration during relaxation of the grain boundary.
boundary. At the same time, half the Zr4 and O2ÿ ions opposite the Schottky layer (Fig. 2(b)) move across the new boundary to occupy some of the vacancies. As well as reducing the density of sites at the grain boundary, the vacancies allow ions to move easily between corresponding sites either side of the interface. The ¯uorite structure of cubic zirconia can be considered to consist of an fcc network of Zr4 cations with the O2ÿ anions ®lling the tetrahedral interstices. It is therefore not surprising that the structure is similar to that predicted for fcc metals such as copper [13] (except for the introduction of the vacancy layer to reduce the repulsive forces in the ionic system). In contrast, rock-salt structured oxides such as NiO [14], where both the cation and anion subarrays are in fcc arrangements, have an even higher density of like-charged ions across the interface in the initial R5 con®guration. In these
systems the boundary relaxes by sliding one crystal parallel to the boundary to the other so that anions and cations face each other across the interface, so that the mirror symmetry is destroyed. Our predicted structure for the R5 interface is very similar to recently published high-resolution transmission electron microscopy (HRTEM) images of symmetric tilt boundaries in cubic zirconia [15]. 3.2. Grain boundary diusion vs. temperature Simulations of pure ZrO2 at temperatures between 1273 and 2673 K found noticeable oxide ion diusion occurring at the grain boundaries. Since no extrinsically introduced vacancies from doping were present, this must be due to the Schottky vacancies intrinsic in the relaxed R5 structure. In order to distinguish between ions at the grain boundary and those in the bulk, the labelling
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
181
Fig. 3. Average positions of ions at R5 (3 1 0)/[0 0 1] grain boundary in ZrO2 at 1273 K.
scheme shown in Fig. 4 was used for msd calculations. Since the grain boundary structure contains sites with dierent coordination environments (and thus energies), these calculations must necessarily be averages of the behaviour of ions of a particular species at the grain boundary. Although somewhat arbitrary, our labelling scheme appears intuitive since it re¯ects the boundary symmetry, and is useful for highlighting
the dierence between behaviour of ions within the grain boundary region from those in the bulk. At 1273 K, only some of the oxide ions at the boundary travelled far from their starting positions. Although movement between adjacent sites across the grain boundary was relatively easy, the ions tended to quickly return to their starting positions without contributing to net diusion. As shown in Fig. 5(a), however, after about 30 ps
182
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
Fig. 4. Labelling scheme used to distinguish ions at the grain boundary (II) and second layer from the grain boundary (I) from bulk ions for msd calculations.
Fig. 5. Mean square displacements (msds) of ions at various temperatures, T, according to the labelling scheme shown in Fig. 4. (a) T 1273 K, (b) T 1673 K, (c) T 1873 K, (d) T 2073 K.
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
there is a de®nite increase in the msd of O±II ions, indicating that some diusion has occurred. There is also a small amount of diusion in the O±I ions, showing that the vacancies initially at the grain boundary have penetrated into the second layer from the boundary. When the temperature is increased to 1673 K (Fig. 5(b)), the diusion of O± II, and to some extent O±I, ions becomes more distinct, and even after short time periods net diusion has taken place. Analysis of the components of the msds in each direction, however, shows that most of the diusion is due to motion across the boundary (the x direction) rather than along it. Fig. 5(c) shows the msds at 1873 K. The increase and subsequent decrease in the O±II msd shows that after diusing a certain distance, ions tended to return to their starting positions or remain where they were. Analysis of the trajectories of the ions con®rmed this conclusion. The ions appeared to be trapped within the immediate region about the grain boundary and were unable to penetrate further into the ¯uorite crystal or traverse the larger distances between sites along the interface. This suggests that there is a signi®cant energy barrier to movement from the grain boundary into the bulk and for long-scale diusion along the R5 grain boundary. A further increase in temperature to 2073 K enables the oxide ions to diuse rapidly in the grain boundary and also the bulk of the crystal, as shown in Fig. 5(d). Although most of the diusion is associated with ions starting at the grain boundary, some diusion of O±I and O (bulk) ions also becomes apparent. The msds of ions at 2673 K were very similar in form to those at 2073 K except that much more rapid diusion of O±II ions took place. In both cases, diusion of O±II ions in the x direction was slightly greater than in the y or z directions, indicating that diusion along the grain boundary is not as favourable as that across the boundary from one crystal to the other. Despite the very high temperatures simulated, no cation diusion was observed. This is not surprising considering the very short time scales involved, and is consistent with the low cation diusion coecients determined experimentally. This result can also explain why very long anneal
183
times are required before the equilibrium structure is reached in real materials. Nevertheless, the msd plots reveal that cations at the boundary oscillate around their mean positions with greater amplitude than those in the bulk. This is obviously a result of the greater amount of space between ions at the grain boundary. 3.3. Diusion in single crystals and bicrystals of YSZ A plot of the msds of a single crystal and bicrystal of 8YSZ is given in Fig. 6. This plot illustrates clearly the decrease in conductivity caused by introducing grain boundaries. The single crystal msd also shows the change in slope ®rst reported by Brinkman et above about 0.4 A al. [5], postulated as being due to a transition from neighbour±neighbour diusion to long-scale ion diusion. The bulk tracer diusion coecient, D, was calculated according to the Einstein relation from the slope of the single crystal msd curve to be 1.1 ´ 10ÿ11 m2 /s. From the Nernst±Einstein equation, this is equivalent to an oxide ion conductivity of 3.0 S/m. If we assume a Haven ratio similar to that in a single fcc crystal, the corresponding conductivity for the bicrystal containing two R5 grain boundaries is half this value. In other words, the presence of the grain boundaries has reduced the conductivity of the system by a factor of two. In real materials the volume of the bulk crystal is
Fig. 6. Mean square displacements (msds) of ions in a single crystal (SC) and R5 bicrystal (GB) of 8 mol% YSZ at 1273 K.
184
C.A.J. Fisher, H. Matsubara / Computational Materials Science 14 (1999) 177±184
much larger than that of the grain boundary regions, so this eect is expected to be greatly reduced at typical grain sizes (>1 lm). However, in nano grain-sized zirconia the volume of grain boundaries will be comparable to that of the bulk crystals. It has recently been suggested that nanosized grains are a promising means of increasing the conductivity of 8YSZ [16], based on the expectation that oxygen diusion along the boundaries will be more rapid than in the bulk. However, our simulations suggest that this is not likely to be the case, concurring with the limited experimental evidence available so far [16]. Despite the lower density and hence more ``free space'' at the grain boundary, the decrease in symmetry and increase in disorder at the boundary are not as conducive to ion transport as the equi-energy sites in the bulk. 4. Conclusions MD simulations of pure ZrO2 bicrystals containing R5 (3 1 0)/[0 0 1] grain boundaries have shown that the interface contains a number of intrinsic vacancies which contribute to grain boundary diusion at elevated temperatures. On addition of yttrium ions and removal of the corresponding number of bulk oxide ions, however, the grain boundaries were found to decrease the overall oxide ion conductivity in 8 mol% YSZ. Although this is not expected to have a strong eect at typical grain sizes, it is expected to lead to lower conductivities in materials containing nano-sized grains, since diusion of oxide ion vacancies along the boundary is less rapid than that in the bulk.
Acknowledgements This work has been carried out as part of the Synergy Ceramics Project under the Industrial Science and Technology Frontier (ISTF) Program promoted by AIST, MITI, Japan. Under this program, part of the work has been funded through NEDO. The authors are members of the Joint Research Consortium of Synergy Ceramics.
References [1] N.Q. Minh, J. Am. Ceram. Soc. 76 (1993) 563. [2] D. Bingham, P.W. Tasker, A.N. Cormack, Phil. Mag. A 60 (1989) 1. [3] F. Shimojo, T. Okabe, F. Tachibana, M. Kobayashi, H. Okazaki, J. Phys. Soc. Japan 61 (1992) 2848. [4] X. Li, B. Hafskjold, J. Phys.: Condens. Matter 7 (1995) 1255. [5] H.W. Brinkman, W.J. Briels, H. Verweij, Chem. Phys. Lett. 247 (1995) 386. [6] H. Okazaki, H. Suzuki, K. Ihata, Phys. Lett. A 188 (1994) 291. [7] M. Wilson, U. Sch onberger, M.W. Finnis, Phys. Rev. B 54 (1996) 9147. [8] A.F. Kohan, G. Ceder, Comp. Mater. Sci. 8 (1997) 142. [9] H.F. Fischmeister, Journal de Physique C4 (1985) 3. [10] D.M. Duy, J. Phys. C: Solid State Phys. 19 (1986) 4393. [11] K.D. Refson, MOLDY 2.11, MOLecular DYnamics, Version 2.11, University of Oxford, 1996. [12] G.V. Lewis, C.R.A. Catlow, J. Phys. C: Solid State Phys. 18 (1985) 1149. [13] G.H. Bishop, Jr., R.J. Harrison, T. Kwok, S. Yip, J. Appl. Phys. 53 (1982) 5596. [14] D.M. Duy, P.W. Tasker, Phil. Mag. A 47 (1983) 817. [15] K.L. Merkle, G.-R. Bai, Z. Li, C.-Y. Song, L.J. Thompson, Phys. Stat. Sol. A 166 (1998) 73. [16] R. Vassen, D. St over, L.G.J. de Haart, M. Cappadonia, Brit. Ceram. Proc. 56 (1996) 35.