Predictive Microkinetic Model for Solid Oxide Fuel Cell Patterned

Subscriber access provided by BOSTON UNIV

Article

Predictive Micro-Kinetic Model for Solid Oxide Fuel Cell Patterned Anode: Based on Extensive Literature Survey and Exhaustive Simulations Shixue Liu, Arief Muhammad, Kazuya Mihara, Takayoshi Ishimoto, Tomofumi Tada, and Michihisa Koyama J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b05914 • Publication Date (Web): 02 Aug 2017 Downloaded from http://pubs.acs.org on August 8, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Predictive Micro-kinetic Model for Solid Oxide Fuel Cell Patterned Anode: Based on Extensive Literature Survey and Exhaustive Simulations Shixue Liu, *,† Arief Muhammad,∥ Kazuya Mihara,∥ Takayoshi Ishimoto,† Tomofumi Tada,§ and Michihisa Koyama*,†,∥ †

INAMORI Frontier Research Center, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan

∥

Department of Hydrogen Energy Systems, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan

§

Materials Research Center for Element Strategy, Tokyo Institute of Technology, 4259 Nagatsuta, Midori-ku, Yokohama 226-8503, Japan

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 37

ABSTRACT: A one-dimensional micro-kinetic model combining H and O migration mechanisms is used to simulate the surface diffusion, chemical and charge-transfer reactions near the triple phase boundary of Ni/YSZ patterned anode. A number of parameter sets are exhaustively examined in the micro-kinetic schemes to obtain a proper set for the explanation of experimental observations. We find a set of parameters free from unphysical fitting parameters, and can explain a large activity gap of patterned anodes reported in literature. From our simulation, we found that the lower activity patterned anode is kinetically governed by both H and O migration across the triple phase boundary, while the higher activity patterned anodes are governed by the O migration.


Page 3 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. INTRODUCTION A composite of Ni and 8 mol% yttria stabilized zirconia (Ni-YSZ) is the most widely used anode material for solid oxide fuel cells (SOFCs). The practical Ni-YSZ anode is synthesized as porous material in which the specific length of triple phase boundary (TPB) as the active site is three-dimensionally extended. To elucidate the intrinsic local activity per TPB length without the influence of the complex microstructure, thin and dense Ni anode with a pattern configuration prepared on a single crystalline YSZ substrate was investigated.1-4 The TPBs are well-defined as the straight lines of dense Ni patterns in contact with the single crystalline YSZ. The simplified structure of the patterned anodes allows us to directly relate the local activity per TPB length and working conditions such as gas pressure, temperature, and current density. It makes easier to discuss the coupled kinetic processes of physical transport, chemical and charge-transfer reactions at TPB. Mizusaki et al. studied the mechanism based on current density-overpotential curves of patterned anodes and concluded that the rate-limiting step appears on the Ni surface where the hydrogen oxidation by the surface adsorbed oxygen species.1 Bieberle et al. proposed that the rate limiting process is either the hydrogen adsorption/desorption on the Ni surface or the removal of O2- from the YSZ surface in a wet fuel gas atmosphere based on the patterned anode experiments.2 Utz et al. regarded that the charge-transfer reactions at TPB is rate-limiting while it is difficult to determine from the gas pressure dependency which mechanism (i.e., H or O migration mechanism) is major for the patterned anode.3 Yao et al. argued that the charge-transfer reactions coupled with H2O-related processes as the ratelimiting step.4 That is, many different mechanisms have been proposed from experimental



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 37

work of the patterned anodes. The activity of the patterned anodes is another controversial issue in addition to the reaction mechanisms. Figure 1 shows the current densityoverpotential curves reported by different authors, which are normalized by the TPB length and considering the gas pressure dependencies. It is clearly noticed that the patterned anode reported by Mizusaki et al.1 has current density ca. three orders of magnitude lower than those reported by Bieberle et al.2 and Yao et al.4

Figure 1. The normalized current density as a function of overpotential based on the experimental data for the patterned anode.

Bessler and Vogler et al.5,6 and Goodwin et al.7 carefully picked up the possible charge-transfer reaction paths and compared these paths using one-dimensional microkinetic modeling approach. The H, O and OH migration (or spillover) mechanisms were investigated in their micro-kinetic simulations, assuming a single mechanism in their micro-kinetic modeling. In other words, these competitive processes were not coupled in


Page 5 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


their models. Vogler et al.6 compared the simulation results of H, O, and OH migration mechanisms and experimental observations, and concluded that the H migration mechanism is plausible for the patterned anode by Bieberle et al.2 Goodwin et al.7 proposed that the activity of patterned anode by Mizusaki et al.1 can be explained by H migration mechanism. According to the transition state theory, we have two kinetic parameters, i.e. the preexponential factor and the activation energy for each elementary process. Typical preexponential factor of the reaction rate constant is known to be ca. kbT/h (2×1013 s-1 at 973 K), where kb is the Boltzmann constant, T is the temperature and h is the Planck constant. Carefully examining the parameters used in literature, one can find much deviated values as results of using some of kinetic parameters as fitting parameters so as to make the simulation result match with the experimental observations.6,7 The use of such fitting parameters may be useful in speculating the mechanism behind the experimental observations, however, it is obvious that the artificially fitted parameters do not necessarily retain the physical meaning and thus the simulation based on the fitting is no more predictive. To propose a guiding principle for SOFC anode from the simulation results, the unphysical parameters should not be used in the micro-kinetic modeling. To the best of our knowledge, no experimental or theoretical study explains the origin of the large activity gap between patterned anodes of Mizusaki et al.1 and of Bieberle et al.2 or Yao et al.4 In the theoretical point of views, the gap has to be explained in the context of differences in physical or chemical elementary processes without unphysical fitting parameters, when we admit that each experiment was carefully conducted with a



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 37

reasonable accuracy and the physics are appropriately considered in the simulation. To guarantee the predictivity of the simulation, it is generally desirable to adopt parameters derived from the first-principles method, while it is often difficult to obtain reliable parameters for high temperature. More concretely, one can find a wide range of firstprinciples calculation results in literature for the same elementary process, which will result in the large difference in the kinetic parameter values under the SOFC operation conditions. In the current work, parameters in literature are extensively re-surveyed for the surface diffusion, the surface reactions on Ni and YSZ surfaces, and charge-transfer reactions at TPB. Then, we have conducted an exhaustive number of simulations to explore the physically reasonable parameter set by using an in-house finite element code that incorporates both the H and O migration mechanisms into a single model8. In the analysis, we assumed that a plausible mechanism based on reasonable physical parameters must be able to explain each of patterned anode characteristics in literature and the difference in the mechanism should explain the gap between the patterned anode experiments in literature.

2. SIMULATION DETAILS Assuming that the Ni patterns and the current density through the TPB lines are ideally periodic on the single crystalline YSZ substrate and that the gas pressure are uniform, the two-dimensional cross sectional structure of the patterned anode can be simplified to a one-dimensional model as shown in Figure 2. Following the preceding studies5-8, one-dimensional model is adopted in this study. The simulated range includes half of the periodic YSZ surface and half of periodic Ni surface.


Page 7 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In our micro-kinetic model, the physical diffusion processes of the adsorbed species are simulated by Darcy’s law on the one-dimensional YSZ or Ni surfaces. The chemical surface reactions between the adsorbed species, also including adsorption/desorption reactions, are treated as sink and source terms in Darcy’s law for reactants and products. The electrochemical charge-transfer reactions are treated as flux type boundary condition at TPB for every migrated species.

Figure 2. Schematic of two-dimensional cross section and the simplified one-dimensional model of patterned anode.

2.1. Governing equations and boundary conditions. The surface diffusion processes of chemical species adsorbed on YSZ and Ni surfaces are governed by Darcy’s law as

∂Ci ∂  ∂C  +  − Di i  = si ∂t ∂x  ∂x 

,

(1)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 37

where Ci, Di and si are the surface concentration (mol/m2), diffusion coefficient (m2/s) and source/sink term (mol/m2s) of chemical species i, respectively. The diffusion coefficient of adsorbed species on surface can be expressed as

 E  D = D 0 exp − act   RT 

,

(2)

where D0 and Eact are respectively the pre-exponential factor (m2/s) and diffusion barrier (J/mol); R and T are gas constant and working temperature, respectively. Three types of adsorbed species, O2-, H2O and OH-, are considered on YSZ surface, and four types of adsorbed species, H, O, OH and H2O, are consider on Ni surface. The D0 and Eact of adsorbed species found in literature are summarized in Table S1, and their error bar plots with average, minimum and maximum values are shown in Figure S1. The average values of D0 and Eact shown in Table 1 are used to calculate the reference values of diffusion coefficient at the working temperature of 973 K.

Table 1. The average pre-exponential factor of diffusion coefficient and diffusion barrier from literature survey, and the calculated diffusion coefficient at 973 K. Species YSZ surface

Ni surface

D0[m2/s] -14

O2-

3.0×10

H2O

3.6×10

OHH O OH H2O

5.1×10 2.1×10-7 5.9×10-7 6.0×10-7 6.7×10-7

-7 -8

Eact[kJ/mol] D(973K) [m2/s] -18

77.4

2.1×10

33.2

6.0×10

55.0 15.6 48.5 27.3 18.7

5.6×10 3.1×10-8 1.5×10-9 2.1×10-8 6.6×10-8

-9 -11


Page 9 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The source/sink term si corresponding to surface reactions is expressed as si =

no.rxns

∑ j =1

SR SP  −ν ν ν j,i  k f , j ∏ C l j ,l − k b , j ∏ C l j ,l l =1 l =1 

  

,

(3)

where νj,i, kf,j and kb,j are respectively the stoichiometric factor (unitless) of species i, the forward and backward rate constants in reaction j, and SR and SP are the numbers of reactants and products. The units of the rate constants depend on the surface reactions. The surface reaction rate constant can be expressed by the Arrhenius equation as

 E act  kf , j = kf0, j exp − f , j   RT    k b, j = k

0 f,j

 E fact,j − ∆G j exp −  RT 

and

(4-1)

,

(4-2)

   

଴ ௔௖௧ , ‫ܧ‬୤,௝ , and ∆Gj are the pre-exponential factor (same unit as kf,j), reaction where ݇୤,௝

barrier (J/mol), and Gibbs free energy difference (J/mol) of reaction j, respectively. As shown in Figure 2, the one-dimensional model is defined in between the two end points, A and B. The symmetry boundary condition is adopted for the diffusion flux at the two end points. At the TPB point, i.e. the boundary between YSZ and Ni, flux boundary condition is used for the migrated H and O species involved in the charge-transfer reactions. The flux per unit TPB length (mol/(m s)) is calculated based on the chargetransfer reactions as Ji =

no. CT rxns

∑

j =1



SR



l =1

ν j , i  k fCT,j ∏ θ l

−ν

SP

j ,l

− k bCT,j ∏ θ l l =1

ν

j ,l

  

,

(5)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

where ߥ௝,௜ ,

Page 10 of 37

CT kfCT , j and k b , j are the stoichiometric factor (unitless) of species i, the forward

and backward reaction rate constant in the charge-transfer reaction j, and θl is the coverage of species l on YSZ or Ni surface (unitless). The units of the rate constants depend on the charge-transfer reactions. The rate constants of the charge-transfer reactions include the effect of transferred charges between YSZ and Ni and are related with the working potential ∆φwork (V) as

k

k

CT b, j

=k

CT,0 f, j

CT f,j

=k

CT,0 f,j

 Efact,j   zF   expα exp − ∆φwork     RT   RT

and

 Efact,j − ∆G j    exp − (1 − α ) zF ∆φwork  exp −   RT RT   

(6-1)

,

(6-2)

where α, z and F are the symmetry factor, the transferred charge number in the chargetransfer reaction and the Faraday constant, respectively. According to the electric circuit in the half cell, the working potential ∆φwork can be expressed as

∆φwork = η + ∆φref −

iρelytdelyt 2A

,

(7)

where η, ∆φref, i, ρelyt, delyt, and A are the overpotential (V), the reference potential (V), current density through the electrolyte (A/m2), Ohmic resistivity (Ω m), the thickness and the cross-sectional area of the electrolyte (m2), respectively. The cross-sectional area is equal to (lNi/2+lYSZ/2)×1 m2. The reference potential can be calculated at the equilibrium condition, and it is defined by the gas pressure and the temperature.


Page 11 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.2. Chemical and charge-transfer reactions. In this study, three reactions on the YSZ surface (R1~R3), five reactions on Ni surface (R4~R8) and three charge-transfer reactions (H migration (H1, H2) and O migration (Om)) are considered (see Table S2 for more details). We surveyed reaction parameters for the surface and charge-transfer reactions in the published first-principles theoretical studies; the details are summarized in Tables S3S5. Because there is no first-principles study reporting a complete set of parameters required for all the elementary reactions considered in the micro-kinetic modeling, one is obliged to refer different theoretical studies for the considered elementary reactions. As for YSZ surface reactions, we picked up the references as follows: i) the H2O desorption (reverse reaction of R1) and disproportionation of hydroxyl (reverse reaction of R2) reported by Gorski et al.,9 ii) the H2O adsorption (R1) and reverse reaction of disproportionation of hydroxyl on YSZ (111) surface (R2) reported by Chaopradith et al.,10 iii) the oxide ion transfer from YSZ bulk to YSZ surface (R3) reported by Ammal et al.11 For reactions on the Ni surface, we selected three theoretical reports: iv) hydrogen dissociation (R5) and two-step oxidation (R6, R8) on Ni(111) surface by Blaylock et al.,13 v) parameters for R7 on Ni(111) surface by Wang et al.,14 as well as vi) another set of choices for the Ni surface reactions (R6-R8) developed by Hecht et al. for the methane steam reforming reaction considering 42 elementary reactions,19 which is also cited in micro-kinetic modeling by Vogler et al.6 The charge-transfer reactions on the Ni-YSZ-gas TPB have been theoretically studied adopting different models in size and surface/interface orientations:11,12,15-18 vii) Ammal et al. reported H1 and H2 reaction barriers of 124, 128 kJ/mol and Om reaction



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 37

barrier of 79 kJ/mol,11 viii) Cucinotta et al. reported Om reaction barrier of 162 kJ/mol,17 ix) Fu et al. reported H1 and H2 reaction barriers of 132 kJ/mol, and Om reaction barrier of 131 kJ/mol,18 x) Liu et al. calculated Om reaction barrier from YSZ surface over the yttrium site at TPB to the Ni surface, which is as large as 262 kJ/mol.12 In this work, the H migration reaction parameters from Ammal et al.11 were coupled with Om reaction parameters of Ammal et al.11, Cucinotta et al.17 and Liu et al.12

2.3. Finite element method (FEM). The coupled physics including surface diffusion, chemical and charge-transfer reactions are solved by an in-house finite element code.8 Due to the fact that the typical time constant of the surface reaction is much shorter than that of the surface diffusion, the governing equation (Eq. 1) needs to be discretized implicitly in time for the stable computation. For spatial discretization, a conventional finite element technique is employed. The weak form of Eq. 1 is given by

∫

Ω

W(

C i( n +1) − C i( n ) − s i( n +1) )dx + ∆t

∫

Ω

Di

∂C i  dW dC i( n +1)  dx − WD i dx dx ∂ x  boundary 

=0

,

(8)

of Ω

where the superscript n represents the time step, W is the weight function, and Ω is the continuous YSZ or Ni surface region considered in the analysis. The surface concentration and the weight function are defined on the basis of a standard Galerkin approach using piecewise linear interpolation functions. The resultant nonlinear system of discretized equations is solved by employing a Newton-Raphson method.

2.4. Simulation conditions and parameter sets. Table 2 lists the calculation conditions


Page 13 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for pattern geometry and working conditions, which were extracted from the experimental setups by Mizusaki et al.1, Bieberle et al.2 and Yao et al.4

Table 2. The YSZ and Ni widths in pattern geometry and partial pressures used in simulations. Partial pressure [Pa] Width[µm] Specific TPB length 2 lYSZ lNi lTPB[m/m ] pH2 pH2O a 4 10000/3000/1000 1700/850/400 Mizusaki 25 10 3.03×10 b 4 Bieberle 20 20 15000 50 3.7×10 c 4 100 100 10000 3000 Yao 1.0×10 a b c Taken from ref 1. Taken from ref 2. Taken from ref 4.

To calculate the current density-overpotential relations, the reference potential ∆φref is firstly calculated under the chemical equilibrium condition, then the overpotential is varied from negative to positive conditions; the current density (A/m2) is calculated at each overpotential by summing up the contribution from all charge-transfer reaction fluxes as itot =

no. CT rxns

∑

j =1

SR SP  −ν j ,l ν j ,l CT zF ν j , i lTPB  k fCT θ − k ,j∏ l b, j ∏ θ l l =1 l =1 

  

,

(9)

where lTPB is the specific TPB length (m/m2) of the patterned anode. To obtain the normalized current density per TPB length with the unit of A/m, we need divide the current density by the TPB length. From the data selected in 2.2, we can prepare two sets of YSZ surface reactions: i) R1 and R2 by Gorski et al.9 and R3 by Ammal et al.11 and ii) R1 and R2 by Chaopradith et al.10 and R3 by Ammal et al.,11 two sets of Ni surface reactions: i) R5, R6, and R8 by



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 37

Blaylock et al.13 and R7 by Wang et al.14 and ii) R5-R8 by Hecht et al.,19 and four sets of charge-transfer reactions: i) H1, H2, and Om by Ammal et al.11, ii) H1 and H2 by Ammal et al.11 and Om by Cucinotta et al.,17 iii) H1 and H2 by Ammal et al.11 and Om by Liu et al.12 and iv) H1, H2, and Om by Fu et al.18 Thus, total sixteen reaction parameter sets as shown in Table S6 were prepared. As for the diffusion parameters of the seven adsorbed species on YSZ and Ni surfaces, two orders of magnitude smaller and larger diffusion coefficients (i.e., 0.01D/100D) were also evaluated in addition to the reference values listed in Table 1 considering the wide range of diffusion coefficient values found in literature. This means 2,187 (=37) combinations of diffusion coefficient values are possible, giving a total of 34,992 combinations (16×2187). In addition, 60 different simulations were conducted for each combination changing the overpotential and gas composition. In total, we conducted 2,099,520 simulations (i.e., 16×2187×60) to exhaustively explore the parameter sets reasonably explaining the experimental observations by Mizusaki et al. 1, Bieberle et al.2 and Yao et al.4 All simulations are conducted until the steady state is reached.

3. RESULTS AND DISCUSSION 3.1. Comparison of reaction parameter sets. We compared all the calculated results on the current density-overpotential curves with experimental data from Mizusaki et al.1 as shown in Figure S2 for the sixteen reaction parameter sets; every plot includes 2,187 calculated curves. The best correspondence was obtained in the parameter set 3 in Table S6: YSZ surface reaction parameters from Gorski et al.9 with Ammal et al.11, Ni surface reaction parameters from Blaylock et al.13 with Wang et al.14, H1/H2 parameters from


Page 15 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Ammal et al. 11, and Om parameters from Liu et al.12; the best parameters for surface diffusion will be discussed in the subsequent section.

3.2. Simulation of patterned anodes of Mizusaki et al.1 3.2.1. Current density-overpotential characteristic. Although the reaction energy of Om reaction from DFT calculation by Liu et al.12 is 79 kJ/mol, when the reaction energy of Om is adjusted from 79 to 50 kJ/mol, the calculated curves fit well with curves of Mizusaki et al.1 as shown in Figure 3. Because the reported reaction energy of Om reaction ranges from 43 to 123 kJ/mol depending on the local structures and element of cation atom at TPB,11,12,17,18 the value of 50 kJ/mol used here can be regarded as a reasonable parameter. Other reaction parameters are listed in Table 3. The calculated current density shows a good agreement with that by Mizusaki et al.1 for both anodic and cathodic currents, while it is admitted that the current density dependence on pH2O and pH2 are not well reproduced in cathodic side. The water and hydrogen reaction orders for anodic current corresponding to the fuel cell operation mode will be discussed more in detail in a later section based on the specific study on reaction orders by Mizusaki et al.20



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 37

(b)

(a)

Figure 3. The current density-overpotential curves under different partial pressures of (a) water, and (b) hydrogen. The open symbols are the experimental data from Mizusaki et al.1

Table 3. The reaction parameter set for Ni/YSZ patterned anode. The unit of energy is kJ/mol. Eact 0 164 0.96 0

∆H -42 123 -41 -18

∆G 77 103 -41 63

4

-92

34

98

33

15

63

24

63

2.1×10 m2/(mol s)

42

9

-47

k fCT ,j

Eact

∆H

∆G

124

-27

-53

128

96

50

262

50

15

Surface reactions R1 R2 R3 R4

kf,j 7 2.1×10 m3/(mol s) 1.6×1018 m2/(mol s) 1.6×1018 m2/(mol s)

R5

2.4×10 m5/(mol2 s)

R6 R7 R8

6

3.8×10 m3/(mol s) 11

17

1.4×10 m2/(mol s) 18

1.4×10 m2/(mol s) 15

Charge-transfer reactions 4

H1

3.0×10 m2/(mol s)

H2

1.4×10 m2/(mol s)

Om

4 4

3.0×10 m2/(mol s)


Page 17 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2.2. Dominant mechanism. To see the contribution of H1, H2 and Om reactions to the total current density, each contribution is shown in Figure 4 as a function of overpotential. The pressures of H2 and H2O are 3000 Pa and 850 Pa, respectively. The curve of Om current density is more symmetric than that of H1/H2. The current densities from H1 and H2 are almost the same. In the anodic side, the contribution from the H1/H2 are dominant, and the Om contribution is limited by the high reaction barrier of 262 kJ/mol. While in the cathodic side, the dominant mechanism gradually shifts from H1/H2 to Om as shown in Figure 4(b).

(a)

(b)

Figure 4. (a)The total current density, current density contributed by H1, H2 and Om migration reactions, and (b) the ratio of current density contributed by Om migration reaction to the total current density as a function of overpotential. The pressures of H2 and H2O are 3000 Pa and 850 Pa, respectively.

The coverage distributions of adsorbed species on the polarized patterned anode are shown in Figure 5, with the pressure of H2 and H2O at 3000 Pa and 850 Pa. There is no



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 37

obvious change in coverages of OYSZ, H2OYSZ, HNi, ONi, OHNi and H2ONi along the surfaces, whereas only the coverage of OHYSZ has obvious deviation on YSZ surface due to the polarization.

(b)

(a)

Figure 5. The coverage distribution of adsorbed species on YSZ (left side) and Ni (right side) surfaces at overpotential of (a) -0.2 V and (b) 0.2 V. The pressures of H2 and H2O are 3000 Pa and 850 Pa, respectively.

3.2.3. Rate-limiting step. The normalized fluxes of surface reactions and charge-transfer reactions at overpotential of 0.2 V are shown in Figure 6. There are four competitive paths in our model to produce current, as follows: 1) R3′&R5H2R2′R1′, 2) R3′&R5H2H1R1′,

3),

R3′&R5OmR6R8R4′

and

4)

R3′&R5OmR6R7′R4′, where the single quote mark means the backward direction (i.e., dotted arrows in Figure 6) of corresponding reactions in Table S2. Figure 6 indicated the main path for H2 oxidation is the second path shown above: (I) Oxide ions migrate from YSZ bulk to YSZ surface (R3′) and diffuse to TPB; H2 molecules dissociate on Ni surface


Page 19 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(R5) and diffuse to TPB; (II) H atoms migrate to YSZ side and form OH-YSZ with O2-YSZ (H2); (III) H atoms migrate to YSZ side and form H2OYSZ with OH-YSZ (H1); (IV) H2OYSZ molecules desorb to gas phase (R1′). The forward and backward fluxes of charge-transfer reaction H1 are the smallest in this path, and thus the rate-limiting reaction step is H1 reaction.

Figure 6. The forward reaction fluxes under overpotential of 0.2 V. The net fluxes, i.e. the forward reaction flux subtracted by the backward reaction flux, are shown in brackets. The fluxes are normalized by the flux of the total transferred electrons in all charge transfer reactions. The pressures of H2 and H2O are 3000 Pa and 850 Pa, respectively.

3.2.4. Pressure dependency. Mizusaki et al. reported the dependence of exchange



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 37

current density on the pressure of H2 and H2O from the measurement of the polarization resistance, which is inversely proportional to the exchange current density.20 The reported reaction orders of H2 and H2O at 973 K are about 0 and 0.6, respectively. In the microkinetic modeling, the fluxes of forward and backward charge-transfer reactions can be calculated, and they are in equilibrium at open circuit condition. Thus, the flux at open circuit condition can be used to define the exchange “local” current density of different charge transfer reactions. For H1 reaction, the exchange “local” current density can be expressed as 12

 Efact,H1 − αz∆GH1  kf ,1k b,2kf ,3kb,5 CVbulk CH2O(g )  0 0 CT  θH θ - , (10)  i0,local = zFlTPBkf,H1 exp − Ni OHYSZ  k k k k C C  RT   b,1 f ,2 b,3 f ,5 Obulk H2 (g ) 

where

CVbulk , CObulk , CH2O(g)

and

CH2 (g)

are bulk vacancy concentration, bulk oxygen

concentration, H2O concentration and H2 concentration, respectively. Similar equations can be derived for H2 and Om reactions. The calculated exchange “local” current density of H1, H2 and Om are shown in Figure 7 as functions of gas pressure. Each slope corresponds to the reaction order of the corresponding reaction. The exchange “local” current density of H2 is much larger than that of H1, and that is why H1 is the rate-limiting step in the H migration path. The reaction orders are determined by H1 and Om reactions. The H migration shows positive reaction orders on both pressure of H2 and H2O, while Om shows negative reaction order of H2 and positive order of H2O. The zero of the experimental H2 reaction order can be recognized as an intermediate between the reaction orders of H1 and Om which are 0.23 and -0.52, respectively; the H2O reaction order is also a middle between that of H1 and Om. Thus, these results indicate the patterned anode working as fuel cell is 20 ACS Paragon Plus Environment

Page 21 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


governed by mixed mechanism of H1 and Om.

(b)

(a)

Figure 7. The exchange “local” current density dependence on (a) H2 pressure with fixed H2O pressure of 850 Pa, (b) H2O pressure with fixed H2 pressure of 3000 Pa. The slope value marked for each migration reaction is the reaction order.

3.2.5. Tafel slopes. The Tafel plots can show simple linear relationship between logarithm of the current density and the overpotential at high overpotential regions. The comparison of simulated Tafel slopes with experimental curves is a good way of mechanistic discussion. The exchange “local” current density of H2 reaction is much larger than H1, and thus H2 reaction can be approximated as being in equilibrium even under polarization.

With

the

k f,CTH2 θ H Ni θ O 2 - YSZ = k b,CTH2 θ V Ni θ OH - YSZ

equilibrium

equation

for

H2

migration

reaction

at TPB site, and considering the independence of HNi, VNi and

O2-YSZ on the overpotential, the equation

0 log( θ OH - YSZ ) ≈ F η / 2 . 3 RT + log( θ OH ) YSZ

can be

0 derived at TPB site, where θ OH - YSZ corresponds to the coverage of OH-YSZ at the open



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

circuit condition. This linear dependence of log( θ OH

-

YSZ

)

Page 22 of 37

on η due to the fast H2 migration

reaction becomes important in the following discussion. The dependence or independence of coverage of all adsorbed species and surface vacancies on the overpotential can be found in Figure S3(a). The current density at high positive overpotential contributed by H1 migration reaction can be expressed as

log i = log i0 + log(

θ H θ OH Ni

0 θ H0 θ OH Ni

-

-

YSZ

)+

YSZ

αzFη ln(10 ) RT

≈ log i0′ +

(1 + α ) zFη . (11) ln(10) RT

In this Tafel equation the slope factor, (1+α)z, is equal to 1.5 (α=0.5, z=1). For the cathodic side, the O migration is dominant at high negative potential, the Tafel equation can be expressed as


θO θV θ O0 θ V0 Ni

YSZ

Ni

YSZ

)−

(1 − α ) zF η (1 − α ) zF η ≈ log i0′ − ln(10 ) RT ln(10 ) RT

. (12)

The slope factor -(1-α)z is equal to -1.0 (α=0.5, z=2). The analytical slopes of Tafel plots in anodic and cathodic sides quantitatively correspond well to the micro-kinetic model and experimental results by Mizusaki et al.1, as shown in Figure 8; the slopes with unit of F/ln(10)RT in anodic and cathodic sides are respectively 1.5 and -1. The value obtained by micro-kinetic simulation was 1.46 and -0.61, and Mizusaki et al.’s observation was 1.73 and -0.87. This means the slope in experimental curve can be explained by the dominant mechanism in the competitive H and O migration reactions. Bazant et al. have theoretically proved that, for a single step charge-transfer reaction, the symmetry factor α of 0.5 should be accepted to approximate the current


Page 23 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


density-overpotential relationship generalized by Marcus theory to a simplified relationship, which is same as equations expressed in Eqs. 5 and 6.21 If the symmetry factor α is shifted to other values, the index of the adsorbed species coverage, θl, in Eq. 5 should also be shifted from the stoichiometric factor, ߥ௝,௜ . Goodwin et al.7 obtained slopes corresponding well to Mizusaki’s experimental curve1 at anodic current side with H migration mechanism using the symmetry factor α of 0.5, while discrepancy persists in the cathodic current. Although these fitting symmetry factors result in slopes matching the experimental curves, the availability of the model is decreased.

Figure 8. Comparison of micro-kinetic plot and estimated Tafel slopes with experimental current density-overpotential curve of patterned anode of Mizusaki et al.1 The pressures of H2 and H2O are 3000 Pa and 850 Pa, respectively.

3.3. Simulation of patterned anodes of Bieberle et al.2 and Yao et al.4 3.3.1. The current density-overpotential characteristic. The activation barrier for oxide



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 37

ion migration from YSZ surface to Ni surface via yttrium TPB site in Om mechanism is about 262 kJ/mol according to the calculation by Liu et al.12 As other results in similar DFT calculations, smaller barriers in the range from 79 to 162 kJ/mol are also reported; the barrier height sensitively depends on a local atomistic structure of TPB sites.11,17,18 We thus chose a reaction barrier of 205 kJ/mol for Om migration reaction, and kept all other reaction parameters the same with those shown in Table 3. In addition, we investigated the influence of diffusion coefficients. The calculated curves correspond to the experimental curves well as shown in Figure 9 when we increase diffusion coefficient of ONi from 1.5×10-9 to 1.5×107

m2/s; the original diffusion coefficient of 1.5×10-9 m2/s for ONi was obtained by Eq. 2 and

the average value of D0 and Eact as shown in Table 1, and one hundred times larger D for the oxygen diffusion coefficient (See in Table S2) can be expected by looking at the variation of surveyed data, the maximum D0 and minimum Eact. If there is no notification, the diffusion coefficient of 1.5×10-7 m2/s is used for ONi in the following discussion of models for Bieberle et al.2 and Yao et al.4 When the overpotential is larger than 0.29 V, the curve trend of Yao et al.4 drastically changed, while this was not observed in the curve of Bieberle et al.2 The reason could be the high partial pressure of H2O in the experiment of Yao et al.4 leading to a high surface coverage of ONi on the Ni surface, as seen in Figure S3. This means part of the Ni surface was oxidized at the high overpotential, and the reaction mechanism on that part should be corrected depending on the NiO surface to reproduce that trend.


Page 25 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

Figure 9. The theoretical curves compared with experimental results by (a) Bieberle et al.2, and (b) Yao et al.4 The YSZ and Ni widths in pattern geometry and partial pressures used in the simulations are shown in Table 2.

3.3.2. Dominant mechanism. In our micro-kinetic model for Yao et al.4, the contributions of H1, H2 and Om reactions to the total current density are shown in Figure 10(a). The Om reaction is dominant in both cathodic and anodic sides as shown in Figure 10(b). Compared to the simulation results obtained for the comparison with Mizusaki et al.1, the dominant mechanism changes due to the smaller Om reaction barrier. The surface coverage distributions of adsorbed species under polarization are shown in Figure 11. The coverages of OH-YSZ, ONi and OHNi are sensitive to the polarization.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 37

(b)

(a)

Figure 10. (a)The total current density, current density contributed by H1, H2 and Om reactions, and (b) the ratio of current density contributed by Om migration reaction to the total current density as a function of overpotential. The pressures of H2 and H2O are 10000 Pa and 3000 Pa, respectively.

(b)

(a)

Figure 11. The coverage distributions of adsorbed species on YSZ (left side) and Ni (right side) surface at overpotential of (a) -0.2 V and (b) 0.2 V. The pressures of H2 and H2O are 10000 Pa and 3000 Pa, respectively.

3.3.3. Tafel slopes. The typical species coverage-overpotential relationship in simulated results for Bieberle et al.2 and Yao et al.4 are shown in Figures S3(b) and (c). The coverage


Page 27 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of OH-YSZ at TPB site still follows the linear dependence on overpotential as 0 log( θ OH - YSZ ) ≈ F η / ln( 10 ) RT + log( θ OH ). YSZ

The coverage of ONi has similar linear dependence

at high negative overpotential, which is also represented as

log(θ O Ni ) ≈ Fη / ln(10 ) RT + log(θ O0 Ni ) .

For anodic side dominated by Om reaction, the current density at high positive overpotential can be expressed as

θV θO


Ni

θ θ 0 VNi

2-

YSZ

0 O2 - YSZ

)+

αzFη ln(10) RT

≈ log i0′ +

αzFη ln(10) RT

.

(13)

Independence of coverage of VNi and O2-YSZ from the overpotential under anodic current is assumed to derive the final approximated form in Eq. 13. Here the slope factor, αz, is equal to 1 (α=0.5, z=2). The comparison of Tafel plots was shown in Figure 12. The slopes observed in our micro-kinetic modeling for Bieberle et al.2 and Yao et al.4 are 1.03 and 1.02, respectively. The experimental values of Bieberle et al.2 and Yao et al.4 are 0.69 and 1.04, respectively. Vogler et al.6 reported that the patterned anode of Bieberle et al.2 was dominated by H1&H2 migration reactions. Their corresponding slope factor of the H1&H2 mechanism is 1.5 without active water dissociation reaction on YSZ surface (R2), and decreases to 0.92 with active R2 reaction. Both O migration dominant mechanism in our model and H migration dominant mechanism in Vogler et al.6 cannot well corresponds to the slope of experiment by Bieberle et al.2 The Tafel slope of anodic side by Eq. 13 well corresponds to the slope of Yao et al.4 as shown in Figure 12(b).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

Page 28 of 37

(b)

Figure 12. Comparison of micro-kinetic plots and estimated Tafel slopes with experimental curves of (a) Bieberle et al.2 (Black open squares are data cited by Vogler et al.6) and (b) Yao et al.4 The YSZ and Ni widths in pattern geometry and partial pressures used in the simulations are shown in Table 2.

For the cathodic current side, the current density at high negative overpotenital can be expressed as


θO θV θ O0 θ V0 Ni

YSZ

Ni

YSZ

)−

(1 − α ) zF η (1 − 2α ) F η ≈ log i0′ − ln(10 ) RT ln(10 ) RT

(14)

The linear dependence on overpotential of log( θ O ) ≈ F η / ln(10 ) RT + log( θ O0 ) and Ni

Ni

independence of coverage of VYSZ on overpotential are considered in the approximation in Eq. 14. The analytical slope factor, -(1-2α), in the cathodic side is equal to 0 (α=0.5). The limiting current density at the cathodic side is due to the limited ONi diffusion coefficient on the Ni surface, similar as gas diffusion limit in the porous anode. Assuming a fast ONi surface diffusion, the coverage of ONi at TPB can be nearly independent of the overpotential, and then the slope factor is approaching to -(1-α), which is equal to -0.5


Page 29 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(α=0.5), smaller than the experimental slope factor of Yao et al.4 of ca. -0.29.

3.4. Sensitivity analysis. The sensitivity of diffusion coefficient values are analyzed as shown in Figure 13. For the parameter set for the patterned anode by Mizusaki et al.1, the diffusion coefficients of ONi and H2OYSZ are the most sensitive as can been seen from Figures 13(a) and (b). For the parameter sets for the patterned anodes by Yao et al.4, the most sensitive diffusion coefficient is that of ONi from Figures 13(c) and (d). Compared to the current density in the anodic current side shown in Figures 13(b) and (d), the current density in the cathodic current side shown in Figures 13(a) and (c) is more sensitive to the diffusion coefficient.

(a)

(b)

(c)

(d)

Figure 13. Sensitivity analysis using parameter set for patterned anodes of Mizusaki et al.1 at overpotential of (a) -0.2 V and (b) 0.2 V, and those for patterned anodes of Yao et al.4 at overpotential of (c) -0.2 V and (d) 0.2 V.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 37

The high sensitivity on the ONi diffusion coefficient means the effective width of Ni surface is sensitive to the ONi diffusion characteristics. The net flux of HNi oxidation reaction on Ni surface (R6) is plotted in Figure 14(a) to estimate the effective width of Ni surface starting from the TPB line. To estimate the effective width, deff, threshold of 99.99% of flux contribution from TPB is used. The dependence order of effective active width on the oxygen diffusion coefficient is about 0.4 as shown in Figure 14(b).

(a)

(b)

Figure 14. (a) Net flux of R6 of model for patterned anode of Yao et al.4 (b) Effective width of Ni surface from TPB as a function of diffusion coefficient of ONi. The overpotential is at -0.2 V.

3.5. Discussion of activity decrease. In the preceding sections, the large activity gap between patterned anodes of Mizusaki et al.1 and of Bieberle et al.2 or Yao et al.4 is addressed by a key process of Om reaction barrier from YSZ side to Ni side. The high


Page 31 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


barrier of 262 kJ/mol was used in the parameter set for the lower activity patterned anodes of Mizusaki et al.1 The value corresponds to oxygen migration over the yttrium site at TPB reported by Liu et al.12 Smaller barriers are reported for oxygen migration over the zirconium site at TPB, and the range from 79 to 162 kJ/mol depends on a local atomistic structure of TPB sites.11,17,18 Therefore, we can expect that oxygen must go over a higher barrier if we assume the most of the cations at TPB are Y atoms, due to the segregation indicated by Hughes et al. during the annealing in air.22 Of course, the segregation of some other elements such as Si, Ca or Na22-24 can influence the TPB activities and diffusion properties, resulting in a large overpotential difference as observed in our simulation results.

4. CONCLUSIONS By extensive parameter survey and exhaustive simulations, physically reasonable parameter sets were found to explain of the considerable activity gap observed in different experimental patterned anodes. The key factor governing the activity gap of the patterned anodes is the charge-transfer reaction barrier of oxide ion from YSZ to Ni. In the simulation using a parameter set for patterned anode of Mizusaki et al.1, the dominant mechanism at the cathodic current side and the anodic current side shifts from O migration to H migration. In the simulation using parameter sets for patterned anodes of Bieberle et al.2 and Yao et al.4, O migration is dominant in both anodic and cathodic current sides, and the rate limiting step for the cathodic current side is oxygen diffusion process on the Ni surface.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 37

AUTHOR INFORMATION Corresponding Authors * E-mail: [email protected]. (S. Liu) * E-mail: [email protected] (M. Koyama) Notes The authors declare no competing financial interests. ACKNOWLEDGMENTS This work was supported by JST-CREST (JPMJCR11C2). Activities of INAMORI Frontier Research Center are supported by KYOCERA Corporation. S. L. acknowledges very fruitful discussion with Prof. Dayadeep S. Monder from Indian Institute of Technology Bombay.

Supporting Information Surveyed data for surface diffusion, surface reactions and charge-transfer reactions. Current density-overpotential results of 16 parameter sets. Coverage at TPB site of different models under polarization.

REFERENCES (1) Mizusaki, J.; Tagawa, H.; Saito, T.; Yamamura, T. Kinetic studies of the reaction at the nickel pattern electrode on YSZ in H2-H2O atmospheres. Solid State Ionics 1994, 70/71, 52-58. 32 ACS Paragon Plus Environment

Page 33 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(2) Bieberle, A.; Meier, L. P.; Gauckler, L. J. The electrochemistry of Ni pattern anodes used as solid oxide fuel cell model electrodes. J. Electrochem. Soc. 2001, 148, A646A656. (3) Utz, A.; Störmer, H.; Leonide, A.; Weber, A.; Ivers-Tiffée, E. Degradation and relaxation effects of Ni patterned anodes in H2-H2O atmosphere. J. Electrochem. Soc. 2010, 157, B920-B930. (4) Yao, W.; Croiset, E. Stability and electrochemical performance of Ni/YSZ pattern anodes in H2/H2O atmosphere. Can. J. Chem. Eng. 2015, 93, 2157-2167. (5) Bessler, W. G.; Vogler, M.; Störmer, H.; Gerthsen, D.; Utz, A.; Weberd, A.; IversTiffée, E. Model anodes and anode models for understanding the mechanism of hydrogen oxidation in solid oxide fuel cells. Phys. Chem. Chem. Phys. 2010, 12, 13888-13903. (6) Vogler, M.; Bieberle-Hütter, A.; Gauckler, L.; Warnatz, J.; Bessler, W. G. Modelling study of surface reactions, diffusion, and spillover at a Ni/YSZ patterned anode. J. Electrochem. Soc. 2009, 156, B663-B672. (7) Goodwin, D. G.; Zhu, H.; Colclasure, A. M.; Kee, R. J. Modeling electrochemical oxidation of hydrogen on Ni–YSZ pattern anodes. J. Electrochem. Soc. 2009, 156, B1004–B1021.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 37

(8) Kohno, H.; Liu, S.; Ogura, T.; Ishimoto, T.; Monder, D. S.; Karan, K.; Koyama, M. Detailed transport-reaction models for SOFC Ni-YSZ patterned anodes: A critical inquiry. ECS Trans. 2013, 57(1), 2821-2830. (9) Gorski, A.; Yurkiv, C.; Starukhin, D.; Colpp, H. R. H2O chemisorption and H2 oxidation on yttria-stabilized zirconia: Density functional theory and temperatureprogrammed desorption studies. J. Power Sources 2011, 196, 7188-7194. (10) Chaopradith, D. T.; Scanlon, D. O.; Catlow, C. R. A. Adsorption of water on yttriastabilized zirconia. J. Phys. Chem. C 2015, 119, 22526-22533. (11) Ammal, S. C.; Heyden, A. Combined DFT and microkinetic modeling study of hydrogen oxidation at the Ni/YSZ anode of solid oxide fuel cells. J. Phys. Chem. Lett. 2012, 3, 2767-2772. (12) Liu, S.; Ishimoto, T.; Monder, D. S.; Koyama, M. First-principles study of oxygen transfer and hydrogen oxidation processes at the Ni-YSZ-gas triple phase boundaries in a solid oxide fuel cell anode. J. Phys. Chem. C 2015, 119, 27603-27608. (13) Blaylock, D. W.; Ogura, T.; Green, W. H.; Beran, G. J. O. Computational investigation of thermochemistry and kinetics of steam methane reforming on Ni(111) under realistic conditions. J. Phys. Chem. C 2009, 113, 4898-4908. (14) Wang, G. C.; Tao, S. X.; Bu, X. H. A systematic theoretical study of water dissociation on clean and oxygen-preadsorbed transition metals. J. Catal. 2006, 244, 10-16.


Page 35 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(15) Shishkin, M.; Ziegler, T. Hydrogen oxidation at the Ni/yttria-stabilized zirconia interface: A study based on density functional theory. J. Phys. Chem. C. 2010, 114 11209-11214. (16) Shishkin, M.; Ziegler, T. Direct modeling of the electrochemistry in the three-phase boundary of solid oxide fuel cell anodes by density functional theory: a critical overview. Phys. Chem. Chem. Phys. 2014, 16, 1798-1808. (17) Cucinotta, C. S.; Bernasconi, M.; Parrinello, M. Hydrogen oxidation reaction at the Ni-YSZ anode of solid oxide fuel cells from first principles. Phys. Rev. Lett. 2011, 107, 206103. (18) Fu, Z.; Wang, M.; Zuo, P. X.; Yang, Z.; Wu, R. Importance of oxygen spillover for fuel oxidation on Ni/YSZ anodes in solid oxide fuel cells. Phys. Chem. Chem. Phys. 2014, 16, 8536-8540. (19) Hecht, E. S.; Gupta, G. K.; Zhu, H.; Dean, A. M.; Kee, R. J.; Maier, L.; Deutschmann, O. Methane reforming kinetics within a Ni–YSZ SOFC anode support. Appl. Catal., A 2005, 295, 40-51. (20) Mizusaki, J.; Tagawa, H.; Saito, T.; Kamitani, K.; Ymamura, T.; Hirano, K.; Ehara, S.; Takagi, T.; Hikita, T.; Ippommatsu, M., et al. Preparation of Nickel Pattern Electrodes on YSZ and Their Electrochemical Properties in H2‐H2O Atmospheres. J. Electrochem. Soc. 1994, 141, 2129-2134.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 37

(21) Bazant; M. Z. Theory of chemical kinetics and charge transfer based on nonequilibrium thermodynamics. Acc. Chem. Res. 2013, 46, 1144–1160. (22) Hughes, A.E. Segregation in single-crystal fully stabilized yttria-zirconia. J. Am. Ceram. Soc. 1995, 78, 369-378. (23) Hansen, K. V.; Norrman, K.; Mogensen, M. H2-H2O-Ni-YSZ electrode performance effect of segregation to the interface. J. Electrochem. Soc. 2004, 151, A1436-A1444. (24) Liu, Y.-L.; Primdahl, S.; Mogensen, M. Effects of impurities on microstructure in Ni/YSZ-YSZ half-cells for SOFC. Solid State Ionics 2003, 161, 1-10.


Page 37 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Graphic


Predictive Microkinetic Model for Solid Oxide Fuel Cell Patterned

Recommend Documents