Modelling the 3D microstructure and performance of solid oxide fuel cell electrodes: Computational parameters

Electrochimica Acta 56 (2011) 5804–5814

Contents lists available at ScienceDirect

Electrochimica Acta journal homepage: www.elsevier.com/locate/electacta

Modelling the 3D microstructure and performance of solid oxide fuel cell electrodes: Computational parameters Qiong Cai a,∗ , Claire S. Adjiman b , Nigel P. Brandon a a b

Department of Earth Science and Engineering, Imperial College London, London SW7 2AZ, United Kingdom Department of Chemical Engineering, Centre for Process Systems Engineering, Imperial College London, London SW7 2AZ, United Kingdom

a r t i c l e

i n f o

Article history: Received 22 December 2010 Received in revised form 12 April 2011 Accepted 16 April 2011 Available online 25 April 2011 Keywords: Solid oxide fuel cell electrode 3D microstructure model Microstructure discretization Sample size

a b s t r a c t In this paper, the computational parameters for a 3D model for solid oxide fuel cell (SOFC) electrodes developed to link the microstructure of the electrode to its performance are investigated. The 3D microstructure model, which is based on Monte Carlo packing of spherical particles of different types, can be used to handle different particle sizes and generate a heterogeneous network of the composite materials. Once formed, the synthetic electrodes are discretized into voxels (small cubes) of equal sizes from which a range of microstructural properties can be calculated, including phase volume fraction, percolation and three-phase boundary (TPB) length. Transport phenomena and electrochemical reactions taking place within the electrode are modelled so that the performance of the synthetic electrode can be predicted. The degree of microstructure discretization required to obtain reliable microstructural analysis is found to be related to the particle sizes used for generating the structure; the particle diameter should be at least 20–40 times greater than the edge length of a voxel. The structure should also contain at least 253 discrete volumes which are called volume-of-fluid (VOF) units for the purpose of transport and electrochemical modelling. To adequately represent the electrode microstructure, the characterized volume of the electrode should be equivalent to a cube having a minimum length of 7.5 times the particle diameter. Using the modelling approach, the impacts of microstructural parameters on the electrochemical performance of the electrodes are illustrated on synthetic electrodes. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Solid oxide fuel cells (SOFCs) are electrochemical devices, converting chemical energy into useful electrical energy in a clean and efficient manner. The electrodes of SOFCs, where the electrochemical reactions actually take place, are composites of ion and electron conducting materials and have a porous and tortuous microstructure, to allow the transport of ions, electrons and gases within the electrodes. In porous composite SOFC electrodes, electrochemical reactions are widely assumed to occur at percolated triple phase boundaries (TPBs) – the meeting points of the ionic, electronic and pore phases. Specific requirements of the electrodes are to facilitate the transport of reactant and product gases, provide a large TPB length in the electrode bulk, along with sufficient pathways for the transport of ions and electrons through the electrode. The abundance of TPBs is one primary factor determining an electrode’s performance, but the length of TPB alone is not sufficient to characterize an electrode microstructure and a comprehensive understanding of the percolated solid and pore net-

∗ Corresponding author. Tel.: +44 (0) 207 594 7366. E-mail address: [email protected] (Q. Cai). 0013-4686/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.electacta.2011.04.065

works is required to better model the transport and electrochemical phenomena in these porous composites. Access to this microstructural information is obtained through two primary routes: one is through the development of detailed three-dimensional (3D) electrodes by constructing representations of real electrodes using methods such as focused ion beam-scanning electron microscopy (FIB-SEM) [1–5] or X-ray computed tomography (XCT) [6,7]; the other is through a variety of modelling methodologies which provide numerically generated microstructural analogues. 3D synthetic microstructures are commonly generated using Monte Carlo techniques for random packing or distributing of spherical particles of different types (representing the ionically and electronically conducting phases, and pore phase in some cases) into either a fixed bed or cubic lattice [8–12]. The discrete element method (DEM) [11,13] has also been used to generate synthetic microstructures: starting with a set of non-contacting spherical ionically and electronically conducting particles with random initial positions and velocities in a periodic box, particle re-arrangement and numerical sintering through mechanical interactions between particles are used to mimic the densification process in real-life electrodes. The 3D synthetic microstructure provides the opportunity to manipulate and observe a wide range of design parameters – such as particle size, electrode thickness and

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

composition – that may be very difficult to alter and monitor in practice. Some researchers have used a correlated resistor network model to represent the electrodes by discretizing the particle packing into a resistance network, through which the resistance and the effective conductivity of the electrodes have been calculated [13,14]. 3D microstructures can also be generated using stochastic correlation reconstruction schemes [15], which are based upon statistical spatial correlations derived from 2D experimental images. This approach is useful in investigating microstructure and performance changes in the electrode during operation, although its application as a design tool is limited. These electrode microstructure models can be used to evaluate important microstructural parameters such as phase volume, phase surface and interface areas, phase connectivity and tortuosity, and TPB length. Transport and electrochemical modelling within a discretized electrode microstructure can also be modelled directly, providing the capability to model the fundamental effects of electrode geometry on mass transport and electrochemical performance of porous SOFC electrodes, without the need to input empirical parameters such as porosity, tortuosity and mean pore diameter, which can be difficult to measure. However, only a few studies can be found in the literature linking detailed microstructure to transport and electrochemical modelling. Modelling multi-species gas transport in a solid oxide fuel cell (SOFC) anode was first done using the lattice Boltzmann method (LBM) at a two-dimensional (2D) level, based on a 2D micrograph of the anode [16]; later the incorporation of a six-step kinetic electrochemical oxidation mechanism into this 2D LBM method model enabled the study of electrochemical performance [17]. The LBM method models the continuum transport of gases with different molecular weights, by considering a species-by-species collision process which is broken down into self-collisions, representative of self diffusion and bulk-viscous fluid motion, and inter-species collisions that are representative of molecular diffusion processes. The need to extend the 2D model to 3D is obvious as percolation is inherently three-dimensional in nature. Suzue et al. [15] used the LBM method to model the transport phenomena and electrochemical reactions in 3D SOFC anode microstructures derived using stochastic correlation reconstruction, through which 3D distributions of potential and current were obtained. Gas diffusion was assumed to be due to the combination of binary Fickian diffusion of H2 and H2 O and Knudsen diffusion; and the Butler–Volmer equation was used to describe the charge-transfer kinetics. Golbert et al. [9] employed the volume-of-fluid (VOF) method to analyze 3D SOFC electrode microstructures, modelling the transport of different species (electron, ion and gases) in corresponding phases, coupled with electrochemical reaction modelling. The modelling framework developed by Golbert et al. [9] provides a platform to link electrode design parameters to microstructural properties and in turn the performance of the electrode. The microstructures used in the 3D model can be generated by Monte Carlo packing of spherical particles of different types, followed by simulated sintering, or measured experimentally [1,2,18]. Monte Carlo generation is useful to study the effect of particle size and to generate a heterogeneous network of the composite materials from which a range of microstructural properties can be calculated, including phase volume fraction, percolation and TPB length. The structure is then divided into discrete volumes, and each volume i receives a vector of volume fractions, fi , whose elements fi,k denote the fraction of the volume occupied by each material k.  The sum of the fractions is such that k fi,k = 1. For convenience, these volumes are referred to as volume-of-fluid (VOF) units hereafter. Any VOF unit for which at least one phase fraction value is strictly between zero and one contains a phase interface. The VOF method, which captures exact interface information through discrete volume data, can represent complex multiphase structures

5805

and continuous phase boundaries which are well suited for diffusion/conduction. The focus of this paper is primarily to present our investigation of the space of computational parameters within which the model developed by Golbert et al. [9] can be used. The effect of the choice of resolution/discretization scheme for microstructure analysis and transport/electrochemical modelling, which is important to the accuracy and reliability of the simulation results, is examined. The required sample sizes for the electrodes are investigated, linking the simulated electrode volume to the experimental characterization of SOFC electrodes. An important aspect in reporting the volumetric properties of a heterogeneous structure is the recognition that a large enough volume needs to be considered if the true structure of the sample is to be correctly represented. There are two approaches to ensure the accuracy and validity of the measurement: one is to make measurements on several independent but neighbouring regions of the sample; the other is to make sure that the measured sample is of sufficient size that the quantities derived from it are representative of the electrode. Iwai et al. [3] reconstructed 3D microstructures from FIB-SEM characterization of SOFC anodes and reported variation in the Ni and YSZ volume fractions for samples from the same electrode, which implies that a large sample size is preferred for quantitative analysis. Shearing et al. [19] tested the homogeneity of electrode tomography data based on X-ray computed tomography for the study of Ni–YSZ composite electrode microstructures, showing the need to characterize sample volumes of a minimum size of 10 ␮m3 for the microstructures considered. However, there is uncertainty in the applicability of the 10 ␮m3 sample volume to the characterization of the SOFC anode microstructures of any origin, as the measurements by Shearing et al. [19] were based on only one electrode. In this work we use the model to investigate the effects of the sample size on the extracted characteristics of the sample, linking design parameters such as particle size to the sample size of the electrode, and in turn its performance.

2. The electrode model 2.1. Particle packing A Monte Carlo approach has been adopted to generate microstructural analogues of SOFC electrodes by packing three types of spherical particles into a fixed 3D box. The packed structure seeks to qualitatively mimic the characteristics of a physical three-phase electrode, which is commonly fabricated by sintering a mixture of solid particles, commonly nickel oxide and yttria stabilized zirconia (YSZ) for SOFC anodes. Typically a distribution of particle sizes is used to introduce porosity to the real-life microstructure. In NiO derived anodes, the porosity is also enhanced by the volume change associated with the reduction of NiO to Ni in the fuel environment of the working device. In the Monte Carlo structure generation, pore former particles are used to introduce porosity in the packed structure. At later stages, the pore former particles are removed, leaving pore space for gas diffusion. During the computer-based generation of electrodes, pore former size can be treated as an independent parameter so that particle sizes and volume fractions can be adjusted independently. This is not the case in practice, but the ability to manipulate these parameters freely in-silico provides insight into desirable microstructural characteristics, and can help identify promising avenues of research to improve electrode fabrication. The first step in generating the particle structure is to populate the system with constituent particles (electronic, ionic and pore phases) in a random fashion. An initial particle is randomly placed in a fixed packing box, and then each new particle is positioned

5806

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

such that it is in contact with some existing particles but does not overlap with others. The random selection of the particle type is weighted in order to obtain the desired volume fractions. 2.2. Sintering During fabrication, an electrode is sintered at high temperature to improve the contact between the different particles. The expansion of the particles leads to “necks” through which transport of electrons and ions can occur. Numerically, the initial particle network has singular contact points, and the system is ‘sintered’ by enlarging the particles’ radii. This leads to overlap with a corresponding neck radius, which directly affects the all-important TPB length. Moreover, it affects the resistance (ionic or electronic) between adjacent particles. There are different ways to simulate sintering. Ali et al. [8] simply modelled the sintering behaviour by enlarging the size of the particles to obtain a certain degree of overlap between particles. Schneider et al. [13] used a numerical approach to model the sintering process. During numerical sintering, grain boundary and surface diffusion are introduced to account for the formation of a sintered contact between particles. The particle centres approach each other in the progress of sintering. Kenney et al. [10] captured the effect of sintering by varying the minimum allowable distance between two contacting particles until the desired structure porosity was obtained. Metcalfe et al. [12] modelled sintering by adding a sintering neck to particle intersection regions to represent the result of diffusion of solid materials to the particle intersection regions during sintering, which is driven by the high surface energy at sharp particle–particle interfaces. For all these sintering models, the geometric features (contact size and particle position) of the sintered packings are independent of composition. A geometric approach similar to that of Ali et al. [8] is used to represent sintering in this model. The particles expand upon heating; the centres of the particles and the dimensions of the packing bed are kept fixed in the sintering process. Since this is a simple sintering model, we neglect forces between two expanding particles and assume geometric expansion of perfect spheres. In order to ensure conservation of mass, we have assumed that the particle radius expands at a uniform rate proportional to its radius. The ratio of the final radius of a particle after expansion to the initial radius of the same particle before expansion thus defines the extent of sintering. Upon the expansion of the particles, the adjacent particles may overlap depending on the extent of sintering. The overlap will be taken into account later when determining actual volume fractions, but does not influence the expansion at this stage. In our approach, we have accounted for the possibility of new contact points forming as particles expand. This is done by keeping track of the distances between all the particles and their radii during the sintering and updating the contact points accordingly. After sintering, a complex network of the three phases is formed, as shown in Fig. 1.

Fig. 1. Illustration of the electrode model after sintering, with blue representing Ni phase, red representing YSZ phase and empty representing pore phase.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.).

Particles of different types have different percolating directions. For example, the ion conducting particles should percolate from the electrode/electrolyte interface, whilst the electronic conducting particles and the pore phase should percolate from the current collector/electrode interface. In order to identify the percolation of each single phase, three matrices, each containing only the information on one phase (i.e. type 1, 2, or 3), are extracted from the 3D matrix of the electrode structure. A three-dimensional mask is applied to find the percolating clusters: matrix elements are assigned position indices; each matrix column is processed to check the neighbours of each index. This function is applied to search for the percolating network of each phase independently, effectively yielding three separate cluster-labelled matrices. Sets of four voxels sharing a corner where all three phases are present are then identified as a TPB point. All TPB points are recorded in sparse matrix format as non-zero values and their positions in the matrix are

2.3. Microstructure analysis Once the spherical particles have been sintered, the discretized structure is initialized, by dividing the packing box along its three dimensions into cells of equal size, as shown in Fig. 2. Each cell contains a small volume of the packed structure and is called a voxel, with each voxel given a phase fraction value corresponding to the dominant phase, based on the positioning and sizes of the spheres: for example, 1 indicates a YSZ phase, 2 a Ni phase and 3 a pore phase (i.e. gap between solid spheres). The result of this analysis is a three-dimensional labelled matrix, in which each element represents a labelled voxel. Using this 3D matrix, the percolation and TPB length can be analyzed and calculated.

Fig. 2. Illustration of discretization of the sintered electrode into voxels (as indicated by the small grey cubes). The bigger red cube shown in the figure indicates a VOF unit which contains eight voxels. The electrode microstructure can also be taken as a 3D matrix of VOF units, providing a framework for transport and electrochemical modelling, as described in Section 2.4.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

a 4.0

0

3

5

3

10

Number of VOF units

3

15

3

20

3

25

3

30

3

35

3

40

3

45

3

50

b

3

55

1.0

3

Number of VOF units 3 3 3 3 3 3 3 20 25 303 35 40 45 50 55

3

10 15

0.8

3.0 -2

0.7

TPB density / µm

-2

5

0.9

3.5

TPB density / µm

3

0

5807

2.5 2.0 Vox25 Vox50 Vox100 Vox200 Vox400 Vox500

1.5 1.0 0.5

0.6 0.5 Vox25 Vox50 Vox100 Vox200 Vox400 Vox500

0.4 0.3 0.2 0.1

0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5

0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5

D / L_VOF

D / L_VOF

c

Number of VOF units 4.0

0

3

5

3

10

3

15

3

20

3

25

3

3

35

30

3

40

3

45

3

3

50 55

3.5

TPB density /µm

-2

3.0 2.5 2.0 Vox25 Vox50 Vox100 Vox200 Vox400 Vox500

1.5 1.0 0.5 0.0 0

1

2

3

4

5

6

7

8

9

10

11

D / L_VOF Fig. 3. Variations in the percolated TPB density with the dimensions of the voxels and VOF units for microstructure (a) L10P1, (b) L20P2 and (c) L5P1. D/L VOF is the ratio of the particle diameter to the size of the VOF unit. See Table 4 for further descriptions of each run.

noted. The TPB points on the percolating clusters of all the three phases are identified as percolated TPBs. Note that although the TPB points can be distributed within the overall electrode structure, only those on the percolating path of the three phases are active TPB points (or percolated TPBs) where the electrochemical reactions can actually occur. The TPB length is calculated by joining the neighbouring TPB points, which is dependent on the resolution of the structure discretization, i.e. the original voxel size. The resolution of the discretization required to achieve reliable results is examined and reported in Section 3.

or more variables (e.g. materials interfaces or interfaces between fluid and deformable structures) are tracked and located. Details of the model including the model equations, boundary conditions and input parameters can be found in our recently published paper [2]. The phase fractions in a VOF unit i are represented by fi,pore , fi,io , and fi,el , for the pore, ionic and electronic phases, respectively, with their sum equal to 1. The characteristics of a VOF unit i, for example its conductivity i , are assumed to be dependent on the characteristic values in the pure phases (e.g., no , n = pore, io, el) and the phase fractions in that VOF unit:

2.4. Transport and electrochemical modelling

i =



no fi,n

(1)

n

Transport phenomena in porous electrodes are intimately linked to electrode microstructure; specifically to the porosity, the pore size and the tortuosity factor. In order to capture this behaviour, the transport phenomena are modelled for each phase separately and linked at the percolated TPB, where charge-transfer occurs. Fickian diffusion is assumed for the transport of species in the microstructure. The volume-of-fluid (VOF) method [20,21] is used to solve gas transport, oxygen ion transport, and electron transport in the three phases. Using the VOF methodology, the heterogeneous network is represented by the phase fractions in a VOF unit (which contains a number of voxels) as shown in Fig. 3. The free boundaries on which discontinuities exist in one

Note that materials with mixed conductivity can be handled in this framework (for example, anodes using CGO rather than YSZ), since each material can have a different conductivity for each charged species. In this paper, we assume that each phase conducts only one species. Assuming steady state and ideal gases in the pores, and considering only the diffusion of H2 and H2 O, we model the ionic and electronic potentials in an SOFC anode using the following system of equations:

∇ (io fio ∇ Vio ) = −TPB j

(2)

∇ (el fel ∇ Vel ) = 2TPB j

(3)

5808

 ∇

 ∇

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

DH2 O fpore RT DH2 fpore RT

 ∇ PH2 O

 ∇ PH2

= TPB

j 2F

(4)

j = −TPB 2F

(5)

where TPB is the local TPB density (in ␮m), i.e. the total TPB length in one local VOF unit, and j is the charge-transfer current per unit TPB length (in A ␮m−1 ),  io and  el are the ionic conductivity and the electronic conductivity (in S ␮m−1 ) respectively, DH2 O and DH2 are the diffusion coefficients of H2 O and H2 (in m2 s−1 ) respectively, fio , fel , and fpore are the volume fractions of ion conducting, electron conducting and pore phases, Vio and Vel are the ionic potential and the electronic potentials (in V) respectively, PH2 O and PH2 are the partial pressures of H2 O and H2 , respectively, F is the Faraday constant (in s A mol−1 ), R is the gas constant (in J mol−1 K−1 ), and T is the temperature (in K). All quantities except for F and R depend on the position (x, y, z). In our model, charge-transfer is assumed to be the rate determining electrode process, and thus the Butler–Volmer equation is used to calculate the current density as shown by Eq. (6):



j = jo



exp



nˇF  RT



− exp



−

n(1 − ˇ)F  RT

(6)

where jo is the exchange current per TPB length in A ␮m−1 , n is the number of electrons transferred, ˇ is the transfer coefficient (ˇ = 0.5) and  is the over-potential. For the SOFC anode, the potential terms are related as follows: Veq −  = Vio − Vel

(7)

where Vio the ionic potential, Vel the electronic potential, and Veq is the potential at equilibrium, defined by Eq. (8). Veq

RT = Voe − ln nF



PH2 O



PH2

 = Veq − (Vio − Vel ) = Voe + Vel − Vio +

(9)

RT ln nF



PH2 PH2 O



(10)

Substituting Eq. (10) into Eq. (6) results in a definition of the charge-transfer current as a function of the state variables:

− exp

 P H2

(Voe + Vel − Vio )

2RT

 −nF 2RT

Vio = [0, 1.23] ∂Vel /∂x = 0 ∂PH2 O /∂x = 0 ∂PH2 /∂x = 0

ranges from 0 to 1.23 V, depending on the operating temperature. The electronic conductivity of Ni ( el ) is assumed to be high enough such that its electronic resistance can be neglected, so the electrode potential, Vel , is assumed to be zero throughout the electrode. The boundary conditions as given in Table 1 are as follows. It is assumed that the current collector/electrode interface (ion insulating) is at x = 0 and the electrode/electrolyte interface (electron insulating and impermeable to gas) is at x = L. The symbols PH2 O , ∞ and PH2 , ∞ indicate the partial pressure of H2 O and H2 at the inlet. All the boundary planes perpendicular to the y and z axes are assumed to be insulating. The system of equations represented by equations (2)–(5), (10), (11) and the boundary conditions in Table 1 is solved by developing mass/charge balances for each voxel surrounding the discretization points, using the finite volume method [23]. The details of the solver methodology can be found elsewhere [9]. The input parameters used are given in Table 2. The electronic conductivity is taken from experimental measurement of the conductivity of Ni. The ionic conductivity, calculated using Eq. (12), is similar to the experimentally measured conductivity of YSZ which is 4.0 × 10−6 S ␮m−1 [24]. The exchange current per unit TPB length is derived from an experimental SOFC anode [2]. For the study reported in this paper, the partial pressures of H2 and H2 O at the inlet (i.e. x = 0) are set equal.

3.1. Resolution required for microstructure discretization

Combining Eqs. (7) and (8) gives an expression for the overpotential;

 nF

@ x=L

∂Vio /∂x = 0 Vel = 0 PH2 O = P, ∞ PH2 = PH2 , ∞

3. Results and discussion

Voe = 1.23 − 2.304 × 10−4 (T − 298.15)



@ x=0

(8)

where Voe is the standard potential which is determined as a function of the Gibbs free energy change in the reaction, taking the partial pressures of H2 , H2 O and O2 to be at standard pressure. Voe is calculated using

j = jo exp

Table 1 Electrochemical simulation boundary conditions.



(Voe + Vel − Vio )

PH2 O PH2 O PH2

(11)

Three electrode microstructures are used to examine the resolutions required for discretizing different microstructures: L10P1, L20P2 and L5P1. Details of the three microstructures are given in Table 3, showing the overall dimension, the diameter of the particles for generating the structure, and the phase volume fractions of each microstructure. Different resolutions for microstructure discretization are evaluated by determining the percolated TPB density of the microstructures when different sizes of voxels and VOF units are applied (Note that both voxels and VOF units are cubes.). Details of the microstructure discretization scheme including the number of voxels, voxel size and the ratio of particle diameter to voxel size are given in Table 4. For each microstructure, six cases of voxelization (i.e. Vox25, Vox50, Vox100, Vox200, Vox400 and Vox500) are investigated; for each case of voxelization, different sizes of the VOF units (which give different number of VOF units) are examined. The variations in the percolated TPB density with the dimensions of voxels and VOF units are shown in Fig. 3. As shown in Fig. 3(a) and (b), microstructure L10P1 and L20P2, which are generated with different overall dimensions and differ-

The ionic conductivity,  io , is estimated by using Eq. (12) [22]: io =

0 io

T

 E io

exp −

RT

(12)

0 = 36 S ␮m−1 is the pre-exponential factor (or referwhere io ence ionic conductivity) and Eio = 8.0 × 104 J mol−1 is the activation energy for the ionic transport. These phenomena are modelled in each VOF unit i, and volumeaveraged values are obtained for each VOF unit. Eqs. (2)–(5) are discretized through finite differences, taking a step size of 10−5 . Vio

Table 2 Electrochemical simulation input parameters. Temperature (K) Ionic conductivity (S ␮m−1 ) Electronic conductivity (S ␮m−1 ) [25] DH2 (m2 s−1 ) [26] DH2 O (m2 s−1 ) [26] Voe (V) Exchange current per unit TPB length (A ␮m−1 ) [2] PH2 ,∞ = P (atm)

1073 4.28 × 10−6 2.4 8.28 × 10−4 8.28 × 10−4 1.05 9.4 × 10−11 0.5

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

5809

Table 3 Details of the three microstructures, including the overall dimension, the diameter of the particles for generating the structure, the volume fractions of Ni, YSZ and pore phases. Microstructure

Overall dimension (␮m3 )

Particle diameter (␮m)

VNi

VYSZ

Vpore

L10P1 L20P2 L5P1

10 × 10 × 10 20 × 20 × 20 5×5×5

1.0 2.0 0.5

0.43 0.43 0.46

0.35 0.36 0.32

0.22 0.21 0.22

ent particle diameters but have the same overall length to particle diameter ratio, show similar development of the percolated TPB density: at Vox25, Vox50 and Vox100, the percolated TPB density shows large variations with the ratio of particle diameter to the edge length of the VOF cube unit, D/L VOF, or the number of VOF units; whilst at Vox200, Vox400 and Vox500, the percolated TPB density stays constant over a wide range of D/L VOF values. For microstructure L5P1, variations in the percolated TPB over a wide range of D/L VOF are seen with all the different voxel dimensions, as shown in Fig. 3(c). For all three microstructures, the percolated TPB density increases with the increase in the number of voxels (or D/L vox); however, the changes in the percolated TPB density values become very small or even negligible at high numbers of voxels (or a high D/L vox value). For microstructure L10P1 and L20P2, a change of less than 5% in the percolated TPB density is observed when the number of voxels is greater than 2003 (corresponding to a D/L vox value of 20 as indicated in Table 4); a higher accuracy, with the change in the percolated TPB density less than 1%, can be found when the number of voxels is greater than 4003 (corresponding to a D/L vox value of 40 as indicated in Table 4). For microstructure L5P1, a change of less than 5% in the percolated TPB density is found with the number of voxels greater than 1003 (corresponding to a D/L vox value of 20) and a change in the percolated TPB density of less than 1% is seen when the number of voxels is greater than 2003 (corresponding to a D/L vox value of 40). We thus conclude that, to obtain a reliable TPB density, the structure discretization scheme should be based on a D/L vox value of at least 20, with 40 desirable for high accuracy of the simulation results. In terms of VOF units, the higher the number of VOF units, the more detail regarding the transport and electrochemical phenomena can be determined, as these are modelled at the level of VOF units. Nevertheless the choice of the number of VOF units needs to be balanced against the high computational demand arising from a high number of VOF units. For example, based on a MatlabTM implementation, a single simulation run with Vox400 and 503 VOF units for microstructure L10P1 occupies 1 h of CPU time on a HPZ800 work station of 24 GB memory; when the number of VOF units is reduced to 253 , the same simulation completes within 6 CPU minutes on the same machine, with other programmes running properly at the same time. Taking into consideration the observation that the percolated TPB density of different microstructures shows little variation when the number of VOF units reaches 253 , the optimum number of VOF units is chosen

as 253 . Depending on the specific microstructure, the size of the VOF units can be different for a fixed number of VOF units. For example, when the number of VOF units is 253 , Fig. 3(a) shows D/L VOF = 2.5 for microstructure L10P1, corresponding to a VOF unit size of L VOF = 0.4 ␮m as D = 1 ␮m in this case; Fig. 3(b) shows D/L VOF = 2.5, giving L VOF = 0.8 ␮m as D = 2 ␮m in this case; Fig. 3(c) shows D/L VOF = 5 ␮m for microstructure L5P1, corresponding to L VOF = 0.2 ␮m as D = 1 ␮m for this structure. 3.2. Sample size required for characterizing electrodes To study the effects of the sample size, 28 electrodes were generated by Monte Carlo packing of uniformly sized particles into cubic boxes. Four particle diameters, D = 0.5, 1.0, 2.0, and 4.0 ␮m, are tested; for each particle size, packing boxes of seven different dimensions are tested. The edge length of the packing box, L, changes in accordance with the diameter of the particles packed in that box, to achieve a given range of L/D ratios for a given particle size. Although each electrode is generated using uniformly sized particles for particles of different types (i.e. Ni, YSZ and pore former), it is recognized that in reality, SOFC electrodes are fabricated from powders of a wide range of particle sizes. Kenney et al. [10] and Metcalfe et al. [12] introduced a particle size distribution during their particle packing processes. By way of a simplified example, the seed powders for an electrode can be approximated by mono-disperse spheres; such an approximation is still useful in giving insights into the link between the particle size and the microstructure and performance of the electrode. Given the random nature of the Monte Carlo method used to generate electrodes, for a given set of manufacturing parameters, for example, one particle diameter and one box length, each instance of the particle generation process produces a different geometrical configuration. Hence a sufficient number of such realizations of electrode generation are needed to obtain results that are statistically invariant with respect to the number of realiza¯ over an ensemble of randomly tions. The average of a property, X, generated structures (e.g. X¯ =

1 N

N 

Xi ) can then be taken as rep-

i=1

resentative of the whole ensemble. The number of realizations required to achieve statistical significance may vary depending on the electrode size. Simulating a relatively small sample size requires a larger number of realizations. For the 28 electrodes inves-

Table 4 Details of the microstructure discretization scheme for three different electrode structures, including the number of voxels, voxel size and the ratio of particle diameter to voxel size. L10P1

Vox25 Vox50 Vox100 Vox200 Vox400 Vox500

L20P2

L5P1

n vox

L vox (␮m)

D/L vox

n vox

L vox (␮m)

D/L vox

n vox

L vox (␮m)

D/L vox

253 503 1003 2003 4003 5003

0.4 0.2 0.1 0.05 0.025 0.02

2.5 5 10 20 40 50

253 503 1003 2003 4003 5003

0.8 0.4 0.2 0.1 0.05 0.04

2.5 5 10 20 40 50

253 503 1003 2003 4003 5003

0.2 0.1 0.05 0.025 0.0125 0.01

5 10 20 40 80 100

n vox: the number of voxels used for discretizing the microstructures; L vox: the length of the voxel cube, i.e. voxel size; D/L vox: the ratio of particle diameter to voxel size.

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

tigated here, a number of realizations ranging from 20 to 400 have been used in order to obtain a representative ensemble; a standard deviation of less than 10% of the ensemble mean is used as an indicator of convergence. In the case of uniform particle sizes, the total number of particles packed in one electrode is expected to be independent of the volume fraction of each particle type. It is worth noting that the packing algorithm described in the electrode model is not limited by the number of particles simulated. However, larger samples require greater computational resources for simulation. The 28 electrodes generated all have the same composition (phase volume fractions) at 0.215 for porosity, 0.350 for YSZ and 0.435 for Ni, which corresponds to the composition of the experimental electrode, based on which the exchange current per unit TPB length was derived [2]. The volume fractions were chosen equal for all electrodes in order to eliminate the effects of volume fraction on other structural properties or performance indicators. For uniformly sized particle packing, the phase volume fractions of the generated electrodes are co-determined by the initial particle fraction and the sintering extent. A particle fraction of YSZ: Ni: pore former = 4.0:4.8:1.0 has been used, together with a sintering extent of 1.2 to achieve the target electrode composition. Fig. 4 shows that these 28 electrodes, although generated from particles of different sizes, achieve similar volume fractions in the final microstructures, very close to the desired values. Small variations in the phase vol-

0.50

0.40 YSZ

0.35 0.30 0.25

Pore

0.20 0.15

D=4 um D=2 um D=1 um D=0.5 um

0.10 0.05 0.00 0

2

4

6

8

10

12

14

16

L/D Fig. 4. Volume fraction of Ni, YSZ and pore phase, as a function of L/D.

ume fractions are observed for electrodes with a small L/D ratio (L/D < 6), which is mainly attributed to the more heterogeneous nature of the structures at small L/D ratios. Fig. 5 shows the percolated TPB density (i.e. percolated TPB length per unit volume) as a function of L/D ratio. The four curves

0.9

0.25

(b) D = 2 µm

(a) D = 4 µm 0.8

TPB density / µm

-2

-2

0.20

TPB density / µm

Ni

0.45

Volume fraction

5810

0.15

0.10

0.05

0.7 0.6 0.5 0.4 0.3

0.00 0

2

4

6

8

10

12

14

0

16

2

4

6

8

10

12

14

16

10

12

14

16

L/D

L/D 14

4.0 (c) D = 1.0 µm

(d) D = 0.5 µm

13

3.5 -2

TPB density / µm

TPB density / µm

-2

12 3.0 2.5 2.0 1.5

11 10 9 8 7 6

1.0 0

2

4

6

8

L/D

10

12

14

16

0

2

4

6

8

L/D

Fig. 5. Percolated TPB density as a function of L/D ratio when (a) D = 4 ␮m, (b) D = 2 ␮m, (c) D = 1 ␮m and (d) D = 0.5 ␮m.

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

0.8

14

L/D = 7.5 L/D = 10 L/D = 15 Fit curve

10

0.7 -2

8 6 4 2

D = 4 µm D = 2 µm D = 1 µm D = 0.5 µm

0.6

Current density / A cm

-2

12

TPB density / µm

5811

0.5 0.4 0.3 0.2 0.1

0

0.0

0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

0

4.5

2

4

6

Particle diameter / µm

in Fig. 5 can be expressed using a universal function as indicated by Eq. (13): TPB

L D

(13)

where f(L/D) = 1 − exp(pt − L/D) is a universal function which is independent of D. Fitting equation (13) to the data in Fig. 5 gives ˛ = 3.32 and pt = 1.85, which are both universal parameters for all the four curves. It can be seen clearly that the percolated TPB density increases with increasing L/D ratio, and then reaches a near constant value when L/D ≥ 7.5. This is observed for all the four cases of D = 0.5, 1.0, 2.0 and 4.0 ␮m, indicating that percolated networks are only fully formed when L/D ≥ 7.5. Together with the observation from Fig. 4 that phase volume fractions stay constant at L/D ≥ 7.5, a conclusion can be drawn that a representative sample of one electrode should be equivalent to a cube whose edge length is at least 7.5 times of the diameter of the particles forming that electrode. The sample volume required for characterizing the SOFC anode has been reported by others. Grew et al. [27] calculated the solid surface area of the 3D microstructure of an SOFC anode and concluded that for the case of uniform particle diameter of 1.0 ␮m, the representative sample should have a minimum volume of 5 ␮m3 . Metcalfe et al. [12] presented the total and connected TPB length per unit volume of the 3D electrode microstructure generated by introducing a particle size distribution for particle packing on the face-centred cubic (fcc) lattice; and concluded that a cubic sample having an edge length of approximately 14 times of the mean particle diameter is of sufficient size to represent an SOFC anode. In our future work, particle size distributions will be introduced in the electrode generation process of our model, to examine whether the particle size distribution impacts on the representative sample size of an SOFC anode. It is also clear from Fig. 5 that the particle diameter has a large influence on the percolated TPB density: the percolated TPB density increases dramatically with decreasing particle diameter. Fig. 6 shows interestingly the same decay behaviour of the percolated TPB density with particle diameter for all three L/D values, i.e. 7.5, 10 and 15. This is because that at L/D ≥ 7.5, the function f(L/D) in equation (13) becomes unity, making equation (13) a simple form as given in equation (14). TPB =

˛ D2

10

12

14

16

Fig. 7. Area current density as a function of L/D, at over-potential of 50 mV. The operating temperature is 1073 K; the inlet gas consists of 50 vol% H2 and 50 vol% H2 O, as shown in Table 2.

Fig. 6. Percolated TPB density as a function of particle diameter.

˛ = 2f D

8

L/D

(14)

where ˛ = 3.32. Eq. (14) indicates that at large L/D values the percolated TPB density decreases in the order of 1/D2 as particle diameter D increases.

The high percolated TPB density found with electrodes made of small particles is attributed to the increased number of contact points introduced by the reduced particle size. Although the advantage of using small particles to achieve high percolated TPB density is obvious, this may also induce higher concentration overpotentials in the electrode at high current densities. Fig. 7 shows current density per unit area (A cm−2 ) as a function of L/D at over-potential of 50 mV. For the same particle diameter D, the area current density keeps increasing as L increases. This behaviour is observed for all four particle diameters, and can be described using a universal function as given in Eq. (15). j=

˛ f D

L D

(15)



where f = ln(L/D) − pt is independent of D. Fitting equation (15) to the data in Fig. 7 finds ˛ = 1.62 and pt = 0.64, which are both universal parameters independent of D. The increase in current density with the increase in L/D is in accordance with the sum of the percolated TPB length. For the same particle diameter D, as L increases, the sum of the percolated TPB length increases as a function of L3 , although the percolated TPB density per unit volume tends to stay constant after L/D reaches 7.5. Fig. 7 also shows that electrodes fabricated using smaller particles generate more current than those using bigger particles, corresponding to the high percolated TPB density found with small particles, shown in Fig. 6. However, examining current density per unit volume reveals a different picture. Fig. 8 shows the volumetric current density (i.e. current per unit volume, in A cm−3 ) as a function of L/D. For the same particle diameter D, the volumetric current density increases as L increases at small L/D, then peaks at L/D = 5.0, before decreasing with further increases of L/D. This behaviour is observed for all the four particle diameters. A polynomial function stands for the four curves in Fig. 8. i=

1 f D2

L D

(16)

where f = a(L/D)4 + b(L/D)3 + c(L/D)2 + d(L/D) + e is the same function for all the curves, with a = −0.0176, b = −0.8723, c = −15.712, d = 111.57, and e = 32.68. Note that the value of f is very sensitive to the values of a, b, c, d and e. This reflects the fact that the total current generated increases as the volume of the electrode is enlarged, but that an optimum electrode volume exists to give the highest volumetric current density, as revealed by Fig. 8. At electrode volumes beyond this optimum,

5812

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

22 -2

Volumetric current density / A cm

Volumetric current density / A cm

-3

85 (a) D = 4 µm

20 18 16 14 12 10 0

2

4

6

8

10

12

14

(b) D = 2 µm

80 75 70 65 60 55 50 45 40

16

0

2

4

6

L/D 1300

350

Volumetric current density / A cm

-3

-3

Volumetric current density / A cm

300

250

200

150

100 2

4

6

8

10

12

14

16

10

12

14

16

(d) D = 0.5 µm

(c) D = 1 µm

0

8

L/D

10

12

14

1200 1100 1000 900 800 700 600 500 0

16

2

4

6

8

L/D

L/D

Fig. 8. Volumetric current density as a function of L/D when (a) D = 4 ␮m, (b) D = 2 ␮m, (c) D = 1 ␮m and (d) D = 0.5 ␮m. The operating over-potential is 50 mV; the operating temperature is 1073 K; the inlet gas consists of 50 vol% H2 and 50 vol% H2 O, as shown in Table 2.

further increases in electrode volume do not generate significant further current; hence the volumetric current density starts to fall. The L/D ratio for optimum volumetric current density appears to lie around 5–6 for all the particle diameters studied. In the analysis so far, the electrode material properties and performance have been shown as either the total or average value for

the whole electrode. The modelling platform developed is also able to give detailed information of the electrode’s material properties and performance, distributed over the electrode in three dimensions. Here we examine an electrode having a cube edge of 3.75 ␮m, and comprising particles of 0.5 ␮m diameter, to demonstrate this.

0.53 500 -3

H2O

Current development / A cm

Partial pressure

0.52 0.51 0.50 0.49 0.48 0.47 0.0

H2

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

Axial position / µm Fig. 9. Hydrogen and steam partial pressure along the electrode thickness (x = 0 represents the interface of current collector with the electrode, whilst x = 3.75 represents the interface of the electrode with the electrolyte). Black line corresponds to an over-potential of 50 mV; grey thin line corresponds to an over-potential of 200 mV; grey thick line corresponds to an over-potential of 400 mV.

450 400 350 300 250 200 150 100 50 0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

Axial position / µm Fig. 10. Current generation along the electrode thickness at over-potentials of 50 mV (light grey thick line), 200 mV (grey thin line) and 400 mV (black thin line (x = 0 ␮m represents the interface of current collector with the electrode, whilst x = 3.75 ␮m represents the interface of the electrode with the electrolyte).

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

5813

Fig. 11. 3D current generation map of the electrode when over-potential is (a) 50 mV, (b) 200 mV and (c) 400 mV. x is the electrode thickness dimension, y and z dimensions form the cross section of the electrode. At x = 0 ␮m is the interface of current collector with the electrode, and at x = 3.75 ␮m, the interface of the electrode with the electrolyte.

Fig. 9 shows the development of the hydrogen and steam partial pressure along the electrode thickness direction. As can be seen, the hydrogen partial pressure and steam partial pressure are both the same at the inlet, but develop differently along the electrode thickness, as hydrogen is consumed and steam generated. Fig. 9 also shows the response of the electrode to different overpotentials: operating the electrode at high over-potential of 400 mV gives the most significant hydrogen consumption and steam generation, compared to the operation at low over-potentials of 50 and 200 mV reflecting the higher current density generated with increased over-potential. Fig. 10 shows current generation (expressed as A cm−3 ) along the electrode thickness, operating at different over-potentials. To illustrate current generation along the electrode thickness, the electrode is segmented along its thickness into identical sections, each having the same cross-sectional area of 3.75 × 3.75 ␮m2 and a thickness of 0.15 ␮m; the current generated per unit volume is calculated by summing up the current generated in all the VOF units contained in one segment and then dividing it by the volume of the segment. As shown in Fig. 10, operating the electrode at higher over-potential produces higher current density. Fig. 10 also shows that, in the case of operating at a high over-potential of 400 mV, more current is generated in the vicinity of the electrode/electrolyte interface; whilst the current generation curves at lower over-potentials of 50 and 200 mV are relatively flat. Fig. 11 shows the 3D map of the current generation within the electrode at different over-potentials. At an over-potential of 50 mV, current generation throughout the whole electrode

is evenly distributed; at an over-potential of 200 mV, there is relatively more current generated in the vicinity of the electrode/electrolyte interface compared to the region close to the current collector/electrode interface; at an over-potential of 400 mV, it is obvious that significant current generation happens in the vicinity of the electrode/electrolyte interface. This reflects the fact that, as shown in Table 2, the electronic conductivity of the Ni phase is orders of magnitude larger than the ionic conductivity of the YSZ phase; hence the TPB points closer to the electrode–electrolyte interface are preferred reaction sites to those at the gas–electrode interface, such that current can be carried in the Ni phase for as long a distance through the electrode as possible. As such the effective reaction area is preferentially distributed in the vicinity of the electrode–electrolyte layer, especially at higher over-potentials where a greater driving force exists to enable rapid charge-transfer between the YSZ and Ni phases. Fig. 11 also shows the un-even distribution of current density at the electrode/electrolyte interface, which highlights the markedly heterogeneous nature of the microstructure in this area. Such heterogeneity may be important in establishing the durability of the electrode operation, as local variations in current density could be expected to impact durability. 4. Conclusions We have investigated the computational parameters of a 3D microstructure model for linking electrode design parameters to microstructural properties, and in turn the performance of SOFC

5814

Q. Cai et al. / Electrochimica Acta 56 (2011) 5804–5814

electrodes. The 3D modelling procedure starts with a microstructure generation process by Monte Carlo packing of spherical particles of different types, followed by simulated sintering to form a complex heterogeneous microstructure. The synthetic electrodes then go through a microstructure discretization process in which the structure is divided into elements of small volume. The 3D model contains two levels of microstructure discretization: at the smallest level of voxels, microstructure analysis is done to obtain a range of microstructural properties including phase volume fraction, percolation and TPB length; at the level of VOF units (which normally consist of one or more voxels) transport and electrochemical modelling are performed using the volume-of-fluid method to predict the performance of the synthetic electrode. The appropriate resolution for microstructure discretization at the level of voxels was found to be related to the particle sizes used for generating the structure; the particle diameter should be at least 20–40 times the cube length of the voxels. Furthermore, the structure should contain at least 253 VOF units for accurate transport and electrochemical modelling. Twenty-eight electrodes were generated by packing particles of different sizes into cubes of different dimensions. The TPB density was found to be largely dependent on particle size, displaying an exponential decay function with increasing particle diameter. To adequately represent the electrode structure, the characterized volume of the electrode should be equivalent to a cube having a minimum length of 7.5 times the particle diameter. The evolution of gas concentration and current generation within the synthetic electrodes were also predicted using the volume-of-fluid method. The results showed the heterogeneous nature of the synthetic microstructures, especially at the electrode/electrolyte interface where more current was generated, calling for special attention to be paid to design of these interfaces. Acknowledgement This work is supported by the seventh framework of European Commission under the project RELHY (“Innovative solid oxide electrolyser stacks for efficient and reliable hydrogen production”) (Grant No. 213009).

References [1] J.R. Wilson, W. Kobsiriphat, R. Mendoza, H.-Y. Chen, J.M. Hiller, D.J. Miller, K. Thornton, P.W. Voorhees, S.B. Adler, S.A. Barnett, Nature Materials 5 (2006) 541. [2] P.R. Shearing, Q. Cai, J.I. Golbert, V. Yufit, C.S. Adjiman, N.P. Brandon, Journal of Power Sources 195 (2010) 4804. [3] H. Iwai, N. Shikazono, T. Matsui, H. Teshima, M. Kishimoto, R. Kishida, D. Hayashi, K. Matsuzaki, D. Kanno, M. Saito, H. Muroyama, K. Eguchi, N. Kasagi, H. Yoshida, Jounal of Power Sources 195 (2010) 955. [4] D. Gostovic, J.R. Smith, D.P. Kundinger, K.S. Jones, E.D. Wachsman, Electrochemical and Solid-State Letters 10 (2007) B214. [5] N. Shikazono, D. Kanno, K. Matsuzaki, H. Teshima, S. Sumino, N. Kasagi, Journal of the Electrochemical Society 157 (2010) B665. [6] J.R. Izzo, A.S. Joshi, K.N. Grew, W.K.S. Chiu, A. Tkachuk, S.H. Wang, W. Yun, Journal of the Electrochemical Society 155 (2008) B504. [7] P.R. Shearing, J. Gelb, N.P. Brandon, Journal of the European Ceramic Society 30 (2010) 1809. [8] A. Ali, X. Wen, K. Nandakumar, J. Luo, K.T. Chuang, Journal of Power Sources 185 (2008) 961. [9] J. Golbert, C.S. Adjiman, N.P. Brandon, Industrial & Engineering Chemical Research 47 (2008) 7693. [10] B. Kenney, M. Valdmanis, C. Baker, J.G. Pharoaha, K. Karana, Journal of Power Sources 189 (2009) 1051. [11] J. Sanyal, G.M. Goldin, H. Zhu, R.J. Kee, Journal of Power Sources 195 (2010) 6671. [12] C. Metcalfe, O. Kesler, T. Rivard, F. Gitzhofer, N. Abatzoglou, Journal of the Electrochemical Society 157 (2010) B1326. [13] L.C.R. Schneider, C.L. Martin, Y. Bultel, D. Bouvard, E. Siebert, Electrochimica Acta 52 (2006) 314. [14] J. Abel, A.A. Kornyshev, W. Lehnert, Journal of Electrochemical Society 144 (1997) 4253. [15] Y. Suzue, N. Shikazono, N. Kasagi, Journal of Power Sources 184 (2008) 52. [16] A.S. Joshi, K.N. Grew, A.A. Peracchio, W.K.S. Chiu, Journal of Power Sources 164 (2007) 631. [17] K.N. Grew, A.S. Joshi, A.A. Peracchio, W.K.S. Chiu, Journal of Power Sources 195 (2010) 2331. [18] S.H. Chan, X.J. Chen, K.A. Khor, Journal of the Electrochemical Society 151 (2004) A164. [19] P.R. Shearing, J. Gelb, J. Yi, W.-K. Lee, M. Drakopolous, N.P. Brandon, Electrochemistry Communications 12 (2010) 1021. [20] C.W. Hirt, B.D. Nichols, Journal of Computational Physics 39 (1981) 201. [21] M. Sussman, E.G. Puckett, Journal of Computational Physics 162 (2000) 301. [22] H. Zhu, R.J. Kee, Journal of Power Sources 117 (2003) 61. [23] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, second ed., Pearson Educational Ltd., London, 2007. [24] K. Chen, X. Chen, Z. Lü, N. Ai, X. Huang, W. Su, Electrochimica Acta 53 (2008) 7825. [25] M. Yousuf, P.C. Sahu, K.G. Rajan, Physical Review B 34 (1986) 8086. [26] E. Hernández-Pacheco, D. Singh, P.N. Hutton, N. Patel, M.D. Mann, Journal of Power Sources 138 (2004) 174. [27] K.N. Grew, A.A. Peracchio, J.R. Izzo, W.K.S. Chiu, in: S.C. Singhal, H. Yokokawa (Eds.), ECS Transations, Solid Oxide Fuel Cells 11 (SOFC-XI) – Part 3, The Electrochemical Society, USA, 2009, p. 1861.