Article A Microscale Modeling Tool for the Design and Optimization of Solid Oxide Fuel Cells

A two dimensional numerical model of a solid oxide fuel cell (SOFC) with electrode functional layers is presented. The model incorporates the partial differential equations for mass transport, electric conduction and electrochemical reactions in the electrode functional layers, the anode support layer, the cathode current collection layer and at the electrode/electrolyte interfaces. A dusty gas model is used in modeling the gas transport in porous electrodes. The model is capable of providing results in good agreement with the experimental I-V relationship. Numerical examples are presented to illustrate the applications of this numerical model as a tool for the design and optimization of SOFCs. For a stack assembly of a pitch width of 2 mm and an interconnect-electrode contact resistance of 0.025 Omega cm(2), a typical SOFC stack cell should consist of a rib width of 0.9 mm, a cathode current collection layer thickness of 200-300 mu m, a cathode functional layer thickness of 20-40 mu m, and an anode functional layer thickness of 10-20 mu m in order to achieve optimal performance.


Introduction
Solid oxide fuel cells (SOFCs) are clean and efficient chemical to electrical energy transfer devices.A traditional SOFC includes three layers: porous anode, dense electrolyte and porous cathode [1,2].The fuel and air are transported from the gas channels to the three phase boundary (TPB) near the electrode/electrolyte interfaces through the porous electrodes and the electrochemical reactions take place in the TPBs.To improve the fuel cell performance by increasing the TPB length to reduce the activation polarization while maintaining the concentration polarization at an acceptable level [3][4][5][6][7], the new design of an anode-supported SOFC cell often consists of five layers: an anode support layer (ASL), an anode functional layer (AFL), an electrolyte layer, a cathode functional layer (CFL) and a cathode current collector layer (CCCL), as shown schematically in Figure 1.The thin functional layers, CFL and AFL, are designed to increase the TPB length with low porosity and small particle sizes.The thick layers, CCCL and ASL, are of high porosity and large particle sizes to enable easy gas transport to the functional layers.In order to achieve high performance for the fuel cells, the cell geometries should be carefully optimized.The thickness of CFL and CCCL are usually the focus of structural optimization in experimental and theoretical works due to their strong influences on the cathode activation polarization and concentration polarization.The thickness of ASL is also an important parameter for the performance optimization [8].Experimentally, Zhao et al. have carried out a systematic study on the influence of the microstructures of cell components on the cell properties [3] and provided the following optimized cell parameters: an ASL thickness of 0.5 mm and porosity of 57%, an AFL thickness of 20 μm, a CFL thickness of 20 μm and a CCCL thickness of 50 μm.The experiments by Haanappel et al. [4] showed that a CFL thickness of at least 10 μm and a CCCL thickness of at least 45-50 μm were required for the optimal cell performance.As experiments are expensive and time consuming, however, it is difficult for the experiments to explore a large number of design parameters and different material choices.Numerical models that incorporate the proper physical mechanisms are more versatile and flexible and are highly desirable for the guiding analyses of a large engineering community.
Most existing theoretical and modeling works are concerned with the traditional SOFCs without the functional layers [8][9][10][11][12][13][14][15][16].For SOFCs with the functional layers, the electrochemical active areas are extended from the electrode/electrolyte interfaces to the functional layers and gaseous species are consumed and produced within the functional layers instead of only at the interfaces.The coupling of the governing equations is not only by the boundary conditions but also by the source or sink terms and the multi-physics theoretical models are complicated.Consequently, only recently have there been a few theoretical models considering SOFCs with the functional layers [17][18][19][20].
In the present work, a two dimensional stack cell model is developed for SOFCs with the electrode functional layers.Electrochemical reactions in the functional layers and on the electrode/electrolyte interfaces are considered in detail.A dusty gas model is used to describe the gas transfer process in porous layers and electric conducting equations are used to describe electron and oxygen ion conducting processes.In fitting the experimental data for model parameters, the effects of the interconnect ribs and the interconnect-electrode contact resistance on the fuel cell performance are also considered.The model is then used to optimize several geometry parameters for stack cell.

Method
The 2D cross-section geometry model for the SOFC stack is shown in Figure 2.Here d channel and d rib are one half of the channel width and one half of the interconnect rib width, respectively.The typical geometric, material and operational parameters are shown in Table 1.The air and the fuel are supplied to the reaction sites in the functional layers (CFL and AFL) from the gas channels through the outside porous layers (CCCL and ASL).Two important electrochemical reactions occur in the functional layers when cell operates with hydrogen fuel and air oxidant as illustrated in Figure 2. At the CFL reaction site, each half of an oxygen molecule gets two electrons conducted from the nearest cathodic interconnect rib and is reduced to an oxygen ion.The oxygen ion is then transported from the CFL reaction site to the AFL reaction site through the electrolyte (YSZ).Finally at the AFL reaction site, the oxygen ion reacts with a hydrogen molecule and produces a water molecule and two electrons and the electrons are conducted to the nearest anodic interconnect rib.
The overall cell performance depends on the operating cell potential and the output current.The operating cell potential (V cell ) can be formally expressed as [8] where 0 E is the Nernst potential, η conc;a and η conc;c are the concentration overpotentials in the anode and cathode, respectively, η act;a and η act;c the anode and cathode activation overpotentials, η ohm;e;a , η ohm;e;c the electronic ohmic overpotentials in the anode and cathode, η ohm;i;a , η ohm;el and η ohm;i;c the ionic ohmic overpotentials in the anode, electrolyte and cathode, η ASR;a and η ASR;c the anode-rib and cathode-rib interface overpotentials due to the contact resistance at the material boundaries.Combining the expressions for the Nernst potential and the concentration overpotentials, where subscripts f, air, AFL and CFL refer to the values in fuel channel, air channel, AFL area and CFL area, respectively.
The details of computing the other overpotential terms in Equation 1 will be described later.Notice, however, the concentration balance overpotentials ( 0 a  , 0 c  ), the activation overpotentials ( act;a act;c ,

  )
and the electrode ohmic overpotentials ( ohm;e;a ohm;i;a ohm;i;c ohm;e;c , , , ) in Equation 1a and Equation 2are dependent on the given active sites in the AFL and CFL layers that are fined meshed to reveal their spatial variations.The average output current density of the SOFC stack can be calculated by where d pitch is the pitch width and d pitch = d channel + d rib .The x and y axis are along the horizontal and vertical directions in the 2D model, respectively.j y is the y component of the current density flux vector.

Electrochemical Reactions in AFL and CFL
In the present model, empirical cathodic (anodic) Butler-Volmer equations are used to describe the electrochemical reaction rate functions in CFL (AFL) and on the CFL/electrolyte (AFL/electrolyte) interface.The transfer current densities per TPB length (A/m 1 ) in the AFL ( TPB,a i ) and CFL ( TPB,c i ) are expressed [17] as: where η is the activation polarization, p i the partial pressure of species i.
For the reaction occurs in the functional layers, the current source (A/m 3 ) can be calculated as V TPB TPB i  .Here V TPB  is the volume specific TPB length (mm -3 ) and can be expressed [17,21] V TPB 6 s i n 2 where d p is the particle diameter, ) the number density of all particles, j  the volume fraction of particle type j (el: electron conducting particle, io: ionic conducting particle or YSZ), j ) the probability of particle type j to form percolated or globally continuous clusters,  the angle of particle contact.
For the reaction occurs on the electrode/electrolyte interface, the current flow (A/m 2 ) can be calculated as A TPB TPB i  .Here A TPB  is the area specific TPB length (mm -2 ) and may be expressed [17] The activation polarizations, η act,a and η act,c , are related to the electric and balance potentials by [8] 0 act,a e,AFL i,AFL a where V e,AFL and V i,AFL are respectively the electronic and ionic potentials at the electrochemically active site in the AFL, V i,CFL and V e,CFL the ionic and electronic potentials at the electrochemically active site in the CFL.

Dusty gas model
The dusty gas model is used here as it is more accurate than the Darcy's law and Maxwell-Stefen model to describe the mass transport in porous electrodes.The original dusty gas model can be expressed as [10,14,22,23 where N i is the total molar flux of species i, x i ( / ) ) the binary diffusion coefficient, ε and τ the porosity and tortuosity, c i , ν i and M i the molar concentration, diffusion volume and molar mass of species i, respectively [24,25].The required parameters are shown in Table 1.
For binary gas, Equation 7 may be reduced to where the effective diffusion coefficients of the species, D i , are given by The effective molar flow velocity u is given by  

Governing equations
The mass transport processes in porous electrodes are governed by the mass diffusion and convective equations.For species i, the transport equation can be expressed by [26] i i N R    (11) where R i is the reaction rate of species i (mol m -3 s -1 ) in the functional layers and 0 i R  in ASL and CCCL because there is no reaction,

Electrical Conduction
The electronic charge transfer in the electrodes (ASL, AFL, CFL and CCCL) are governed by: e e e ( ) where e eff  is the effective electronic conductivity of the electrode, V e the electric potential in the electrode.The ionic charge transfer in the functional layers (AFL, CFL) and electrolyte are governed by where i eff  is the effective ionic conductivity of the functional layers or the electrolyte, V i the electric potential.
i i eff V    is the flux vector of the ionic current density.i Q is the oxygen ionic current source (A/m 3 ) and i 0 Q  in ASL and CCCL as there is no electrochemical reaction, V i T P B , a T P B , a in CFL.The ionic potential differences along the ionic current flux paths yield the ohmic overpotentials, η ohm;i;a , η ohm;i;c and η ohm;el .
The effective conductivity can be expressed [17]: ( ) where σ j is the conductivity of species j, m Bruggeman factor considering the effects of tortuous conduction paths and constricted contact areas between particles.The electric potential loss inside the interconnect plate is assumed to be negligible due to the high conductivity of the metallic material.The local current densities cross the interconnect/anode ( I a j  ) and the cathode/interconnect ( c I j  ) interfaces are determined by the associated electric potential changes, or the interface overpotentials: e,I/a e,a/I ASR;a I a contact contact where V e,I/a and V e,a/I are respectively the interconnect and ASL electric potentials at the anodeinterconnect boundary, V e,c/I and V e,I/c the CCCL and interconnect electric potentials at the cathodeinterconnect boundary, ASR contact the area specific contact resistance.

Boundary Conditions (BCs)
The boundary settings for the mass transport equations are shown in Table 2.The molar concentrations at the channel/electrode interface are related to the molar fractions by the ideal gas equation of state, and can be calculated by 0 0 atm / i i c x P RT  .
The boundary settings for the electronic and ionic charge transfer equations are shown in Table 3.For the electronic charge transfer equation the current flow on the electrode/Electrolyte interface is inward/outward current boundary, while for the ionic charge transfer equation the current flow is interior current source.The contact resistance is set on the interface between interconnect ribs and the electrodes.

Numerical Method
The finite element commercial software COMSOL MULTIPHSICS ® Version 3.4 [26] was used in the present study to solve the required partial differential equations (PDEs) such as the mass diffusion and convective equations, the electronic and ionic charge transfer equations.The Butler-Volmer equations are integrated into the PDEs as sources or boundary settings.The COMSOL stationary nonlinear solver uses an affine invariant form of the damped Newton method [26] to solve the discretized PDEs with a relative tolerance of 1 × 10 -6 .Free triangle meshes with a maximum element size of 5 × 10 -6 m for electrolyte and AFL areas were used.
It is worth mentioning a subtle feature in the numerical solution process.According to Equation 4b, some current may be generated by the electrochemical reaction even when 2 O P is very small and may result in some small negative value for

Model Fitting of the Experimental I-V Curve
As shown in [3], silver meshes were pressed against CCCL and ASL for single cell testing.The meshes may be regarded as interconnect ribs with small pitch width.Here a pitch width of 0.3 mm and a rib width of 0.05 mm corresponding to the experimental description were used for the model fitting.When the contact resistance is 0.008 Ωcm 2 and the boundary hydrogen/oxygen molar fraction is 0.97/0.21, the theoretical curve shown in Figure 3 agrees well with the experimental data.It should be noted here that when the pitch width is small, the current through the electrode/rib interface I a j  and c I j  are uniform [8], then the effective contact resistance is approximately . The effective contact resistance is an effective value assuming that the contact resistance was distributed over the entire pitch width instead of the rib width.Accordingly, eff contact ASR is found here to be about 0.096 Ωcm 2 , in good agreement with the experimental estimation that the total ohmic loss is 0.19 Ωcm 2 at 700 °C and half the total ohmic loss is due to the contact resistance [3].The agreement between the theoretical and experimental results validates the numerical model described here.Moreover, it is striking to see that the interconnect-electrode contact resistance is magnified by a factor of 2d pitch /d rib and even a very small contact resistance may substantially affect the fuel cell testing result.The ionic current density distribution is quite uniform along the x direction and may be viewed as one-dimensional along the y direction in the model for the testing single cell.Figure 4 shows the ionic current density in the CFL, electrolyte and AFL at x = 0 when operating cell potential is 0.7 V.The maximum current density across the electrolyte is 0.64 A/cm 2 .The current densities generated by the functional layer body reaction and the functional layer/electrolyte boundary reaction are respectively 0.591 A/cm 2 and 0.049 A/cm 2 in the cathode side, and 0.439 A/cm 2 and 0.201 A/cm 2 in the anode side.This means that the body reaction is the dominant contributor for the fuel cell current generation and the model picture with only boundary reactions may be too simplified.As shown in Figure 4, the ionic current density in AFL decreases rapidly with the distance from the AFL/electrolyte interface, and almost all current is generated in the 4 μm region around the interface.The electrochemical reaction occurs in the whole area of CFL, but a layer thickness of 11 μm near the interface contributes about 90% of the total current density.This indicates that the effective layer thicknesses for the AFL and CFL electrochemistry reactions are quite small.As the cathode activation overpotential for a given current density is higher than the anode counterpart, the cathode side requires an effective area larger than the anode side in order to generate the same amount of the total current.Therefore, the thickness of the effective cathode active layer is higher than that of the anode side.

Distributions of Physical Quantities in Stack Cell
The output current density is affected by the cathode layer thickness and a thicker layer benefits to the oxygen transport and the stack performance [8].Two CCCL thickness l CCCL , 50 μm and 200 μm, are used to show the physical quantities in stack cell.The pitch width is 2 mm and the rib width is 0.8 mm.The boundary hydrogen and oxygen molar fraction are 0.7 and 0.21 respectively.The area specific contact resistance on the anode/interconnect and the cathode/interconnect boundaries are uniform (0.025 Ωcm 2 ). Figure 5 shows the physical quantity distributions for l CCCL of 50 μm.The average output current density is 0.374 A/cm 2 .Figure 6 shows the distributions for l CCCL of 200 μm and the average output current density is 0.409 A/cm 2 .The hydrogen and oxygen distributions for a CCCL layer of 50 μm are shown in Figure 5a and are similar to the distributions in a stack cell without the functional layers [8].The oxygen distribution in the area under the rib is easily exhausted and the area is hardly active.With a thicker CCCL layer of 200 μm (Figure 6a), the oxygen transfer is more effective and can supply oxygen to the area under the rib, leading to less concentration polarization and more active area.The existence of thin functional layers almost has no effect on the gas transfer.
The total electronic current densities and the current flux vectors are shown Figures 5b and 6b.The total current density in the CFL is obviously less than that in the CCCL.The reason is that the electronic conductivity of CCCL (2.99 × 10 4 s m -1 ) is much largely than that of CFL (5307 s m -1 ), and the current flux in CCCL has a much larger x-component than that in CFL to carry the current to the interconnect rib.Understandably, the current density distribution is more uniform and with a less ohmic polarization for a l CCCL of 200 μm than that for a l CCCL of 50 μm.
Figures 5c and 6c show the y component of the electronic current density flux.The values for AFL and CFL other than at the AFL/electrolyte or CFL/electrolyte interfaces are approximately the local current densities produced by the bulk and boundary electrochemical reactions.At the ASL/AFL interface, the value decreases substantially from the area under gas channel to the area under the rib for l CCCL = 50 μm (Figure 5c) while it is quite uniform for l CCCL = 200 μm (Figure 6c), also indicating that the latter is more favorable.

Optimization of the CCCL thickness
As described in Section 3.2, different CCCL layer thicknesses for otherwise identical cells may have different performances.Figure 7 shows the variations of the average output current density with l CCCL .The current density increases by 9.4% when l CCCL increases from 50 to 200 μm, but the difference is less than 0.4% when the l CCCL changes between 200 and 400 μm.Considering the performance benefit and material cost, the optimized CCCL layer thickness should be 200-300 μm.Unless stated otherwise, the layer thickness of 200 μm is used for optimizing other parameters described below.

Optimization of the rib width
The optimal rib width is the result of balancing the polarizations of the current collection and the gas transport.The most important factor affecting the optimized rib width is the contact resistance [8]. Figure 8 shows the optimized rib width and the optimal output current density versus the area specific contact resistance.Larger contact resistance requires larger optimized rib width to balance the increased ohmic loss, η ASR;a and η ASR;c .The optimal current density decreases rapidly with the increase of the contact resistance for two main reasons: a larger contact resistance leads to a higher ohmic loss and the increased optimal rib width leads to a smaller active area.For the reference contact resistance of 0.025 Ωcm 2 , the optimized rib width is 0.9 mm. Figure 8.The optimized rib width and the optimal output current density vs the area specific contact resistance.

Optimization of the CFL thickness
Limited by the catalytic ability of the cathode materials, the cathode activation polarization is often the most important component of the total fuel cell polarization and a CFL layer with small particle size and small porosity is desirable for enlarging the TPB length and reducing the activation polarization.For a given microstructure, the CFL thickness is also an important factor affecting the cell performance.Figure 9 shows the variation of the average output current density with the CFL thickness (l CFL ).Significant improvement is observed when l CFL is increased from 5 μm to 20 μm and may be attributed to the increase of CFL active area as indicated in Figure 4.The output current density is quite similar when l CFL is increased from 20 μm to 40 μm as neither the active CFL thickness is extended nor the gas transport is affected appreciably.As shown in Figure 9, further increase of l CFL may cause reduced output current density due to the increased concentration polarization.Therefore, the optimal l CFL is between 20 μm and 40 μm.

Optimization of the AFL thickness
The effect of AFL thickness on the output current density is also studied and shown in Figure 10.As shown in Figure 4, the active AFL thickness is only 4 μm.It is not surprising that the output current density in Figure 10 decreases with l AFL .However, the decrease of the output current density is quite small and is less than 0.68% when l AFL varies from 5 μm to 20 μm as the concentration polarization is not significant for such thin AFL layer.An AFL thickness of less than 20 μm may be considered to be optimal.

Conclusions
A 2D finite element model for anode-supported planar SOFC stack cells with functional layers is demonstrated.The model considers the microstructures of electrode material, the gas transfer process in the porous electrode layers, the electronic conducting process in the electrode layers, the ionic conducting process in the functional layers and electrolyte, the bulk electrochemical reactions in the functional layers and the boundary electrochemical reactions at the functional layer/electrolyte interfaces.The model was validated by comparison with experimental results on button cell I-V data.Distributions of physical quantities in SOFC stack cells with functional layers are illustrated.The numerical model is also used to optimize geometry parameters of representative stack cells, demonstrating its capabilities as a design and modeling tool for the development of SOFC technology.

Figure 1 .
Figure 1.Schematic diagram of an anode-supported planar SOFC stack cell unit.

Figure 2 .
Figure 2. Schematic cross-section model of half repeating unit of SOFC stack.
molar fraction, k the permeability, µ the viscosity, p the total gas pressure,


for the water steam production in AFL and 2 V O T P B , c T P B / 4 R i F    for the oxygen consumption in CFL.
the flux vector of the electronic current density.e Q is the electronic current source (A/m 3 ) and e 0 Q  in ASL and CCCL because there is no electrochemical reaction, The electronic potential differences along the electrical current flux paths yield the ohmic overpotentials, η ohm;e;a and η ohm;e;c.

Table 3 .
Boundary settings for the electronic and ionic charge transfer equations.

2 OP
due to the oxygen consumption by the current production.This in turn may cause numerical failure for Equation 4b.To overcome the numerical problem, oxygen molar fraction was set as with the physical requirement.

Figure 3 .
Figure 3.Comparison of theoretical and experimental I-V curves.

Figure 4 .
Figure 4.The ionic current density distribution in the CFL, electrolyte and AFL at x = 0.

Figure 5 .Figure 6 .
Figure 5.The distributions of physical quantities when the CCCL layer thickness is 50 μm: (a) hydrogen and oxygen molar fractions; (b) total electronic current density and current flux vector; (c) the y component of the current flux vector.

Figure 7 .
Figure 7. Variation of the average output current density with the CCCL layer thickness.

Figure 9 .
Figure 9.The variation of the average output current density with the CFL thickness.

Figure 10 .
Figure 10.The variation of the average output current density with the AFL thickness.

Table 1 .
Geometric, material and operational parameters for SOFC with functional layers.

Table 2 .
Boundary settings for mass transports in electrodes.