Optimization of electrode-electrolyte interface structure for solid oxide fuel cell cathode

Journal of Power Sources 449 (2020) 227565

Contents lists available at ScienceDirect

Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour

Optimization of electrode-electrolyte interface structure for solid oxide fuel cell cathode An He *, Junya Onishi, Naoki Shikazono Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo, 153-8505, Japan

H I G H L I G H T S

G R A P H I C A L A B S T R A C T

� Adjoint method is applied to optimize the electrode-electrolyte interface. � Interface TPB reaction is considered during optimization. � Optimized LSCF cathode structure dras­ tically enhances the reaction current. � Reducing cross-section area of pillar-like electrolyte region is preferred. � Effect of ionic conductivities and reac­ tion rates on optimization is investigated.

A R T I C L E I N F O

A B S T R A C T

Keywords: Solid oxide fuel cell Numerical optimization Electrode-electrolyte interface Reaction rate Ionic conductivity

A numerical method is developed to optimize the shape of electrode-electrolyte interface for maximizing the electrochemical performance of solid oxide fuel cells (SOFCs). The electrode-electrolyte interface structure is chosen as the design variable to be optimized, and the total amount of reaction current in the electrode is maximized as the objective function. Adjoint method is used to compute the sensitivity of the design parameter on the objective function, i.e. the sensitivity of the interface shape on the total reaction current is investigated. In this study, pure La0.6Sr0.4Co0.2Fe0.8O3 (LSCF) cathode and Gd0.1Ce0.9O1.95 (GDC) electrolyte is chosen for the design target. Both LSCF/pore double phase boundary (DPB) reaction and LSCF/GDC/pore triple phase boundary (TPB) reaction are considered. It is found from the computational results that the cathode with optimized electrode-electrolyte interface structure largely enhances the reaction current than the flat electrode-electrolyte interface. Finally, the sensitivities of the optimal sturcture on the conductivities and reaction rates are investigated.

1. Introduction Solid oxide fuel cell (SOFC) is attracting great attention for its high energy conversion efficiency and flexibility for the choice of the fuels [1, 2]. It is known that the SOFC electrochemical performance is closely related with the electrode material and its microstructure. At high

operation temperatures (800 � C-1000 � C), the choices of composing materials are significantly limited [1]. On the other hand, at interme­ diate temperatures (500 � C-800 � C), the SOFC suffers from the de­ teriorations of ionic conductivities and polarization characteristics. Therefore, it is of great importance to develop electrode materials which are effective under such intermediate operation temperatures. Mixed ionic and electronic conductors (MIECs) such as La0.6Sr0.4Co0.2Fe0.8O3

* Corresponding author. E-mail address: [email protected] (A. He). https://doi.org/10.1016/j.jpowsour.2019.227565 Received 13 October 2019; Received in revised form 19 November 2019; Accepted 3 December 2019 Available online 16 December 2019 0378-7753/© 2019 Elsevier B.V. All rights reserved.

A. He et al.

Journal of Power Sources 449 (2020) 227565

Nomenclature

i*reac

η*act

Electrode domain Electrolyte domain Electrode-electrolyte interface Level-set function (m) Ionic conductivity oxide ion (S/m) ;electrolyte Ionic conductivity of the electrolyte (S/m)

Ωelectrode Ωelectrolyte Γinterface Φ

σ O2 σ O2

σeff O2

;electrode

μ~O2

F H δI Hsmooth δI; ​ smooth

ρDPB ρTPB

i0;DPB i0;TPB T

α1 ηact pO2 Jc μ~*O2

L S I R

τI

Effective ionic conductivity of the electrode (S/m)

τ

Oxide ion electrochemical potential (J/mol) Faraday constant (C/mol) Heaviside function Dirac delta function Smoothed Heaviside function Smoothed Dirac delta function DPB density (m2 =m3 ) TPB density (m=m2 ) Exchange current density for DPB reaction (A= m2 ) Exchange current density for TPB reaction (A= m) Operation temperature (� C ) α2 Transfer coefficients Activation overpotential (V) Oxygen partial pressure (Pa) Total current (A)

ξ k λDPB λTPB L lelectrode lelectrolyte lreac; ​ TPB lreac; ​ DPB Δx

Adjoint variable for the reaction current Adjoint variable for the activation overpotential Lagrangian cost function Sensitivity (δL δΦ ) Current (A) Resistance (Ω) Discrete fictitious time increment for updating the level-set function Tortuosity factor Characteristic length of the transition layer (m) Dimensionless parameter for ionic conductivities σeff2 (σ O2 ;electrode ) O ;electrolyte Dimensionless parameter for DPB reaction (lelectrode =lreac; ​ TPB ) Dimensionless parameter for TPB reaction (lelectrode =lreac; ​ DPB ) Dimensionless parameter for cathode thickness (lelectrolyte =lelectrode ) Cathode domain thickness (m) Electrolyte domain thickness (m) Characteristic length for TPB reaction (m) Characteristic length for DPB reaction (m) Grid size (m)

Adjoint variable for the oxide ion electrochemical potential

(LSCF) [3,4] have now become a major cathode material because of its high conductivity and electro-catalytic property even at low operation temperatures. In addition to the material developments, microstructures of the electrode are modified in order to improve the cell performance and reliability. Not only the porous microstructures, but also the elec­ trode electrolyte interface shape has been investigated [5–12]. Konno et al. [5] demonstrated that the wavy shaped electrode-electrolyte interface structure can improve the electrochemical performance. Ber­ tei et al. [6] computed the electrochemical performance of the Ni-yttria stabilized zirconia (Ni-YSZ) anode with different designs of anode-electrolyte interfaces. Chueh et al. [7] also computed the per­ formance of Ni-YSZ anode with various designs of anode-electrolyte interfaces, and the mechanical strengths of the representative cases were investigated as well. In addition to the numerical studies, various electrode-electrolyte interface structures have been realized with different fabrication techniques, and the effects of such interface struc­ tures are investigated experimentally [10–12]. Okabe et al. [10] used UV-nanoimprint lithography to fabricate a ceramic electrolyte with micro cylinder pillars with the pillar height of 10.8 μm and pillar diameter of 7.76 μm. Shimura et al. [11] used laser ablation to fabricate the YSZ pillars on the surface of the YSZ electrolyte substrate, and the YSZ pillar is characterized with a pitch of 7.5 μm and a height of 5–10 μm. Zhi et al. [12] used electrospinning to grow YSZ nanofiber scaffold with diameter of 200 nm onto a YSZ substrate, and lanthanum strontium manganite (LSM) particles were infiltrated into the scaffold. Numerical optimization is helpful to find the optimal design pa­ rameters for an electrochemical system which can maximize the per­ formance [15]. Recently, several simulation models have been proposed to optimize the electrode-electrolyte interface structure [13,14,16,17]. For instance, Iwai et al. [16] optimized the electrode-electrolyte inter­ face structure using the adjoint method, and a wavy shaped interface structure is obtained as the optimal structure considering the trade-off between the enhancement of the oxide ion conduction and the deteri­ oration of gas diffusion. Onishi et al. [17] also used adjoint method to optimize the interface shape between Ni-YSZ anode and YSZ electrolyte considering the trade-off between the enhancement of effective ionic

Fig. 1. Reaction sites in (a) pure LSCF and (b) LSCF-GDC composite cathodes with GDC electrolyte.

conductivity and the decrease in the electrochemically reactive region. In their work [16,17], the electrochemical reaction is only considered inside the electrode region, but the reaction at the electrode-electrolyte 2

A. He et al.

Journal of Power Sources 449 (2020) 227565

interface is not taken into account. However, it is not clear whether it is reasonable or not to ignore the three phase boundary reaction at the electrode-electrolyte interface. As shown in Fig. 1 (a), for pure LSCF cathode with GDC electrolyte, LSCF/pore double phase boundary (DPB) reaction takes place inside the electrode, while LSCF/GDC/pore triple phase boundary (TPB) reaction might take place at the electrode-electrolyte interface if LSCF/GDC/pore TPB reaction is active. It is therefore not appropriate to ignore the electrochemical reaction at the electrode-electrolyte interface. For the LSCF-GDC composite cathode with GDC electrolyte as shown in Fig. 1 (b), both TPB and DPB reactions occur inside the electrode region, while only TPB reaction takes place at the electrode-electrolyte interface. Thus, interface and electrode bulk regions should be treated with numerically different electrochemical reaction combinations. It is known that the electrode materials, microstructures and oper­ ation conditions affect the effective conductivities and total reaction kinetics [18–21]. Therefore, the optimal microstructure should depend on materials and operating conditions. In order to give further guidance for the electrode design, it is important to know how the optimal interface shape will be affected by the ionic conductivities and reaction rates of the electrode. In the present study, a numerical model to optimize the shape of the electrode-electrolyte interface is developed. Adjoint method is applied for the optimization of pure LSCF cathode and GDC electrolyte interface structure, and the total amount of the electrochemical reaction current is chosen as the objective function to be maximized. The DPB and TPB reactions are assumed to take place inside the electrode region and at the electrode-electrolyte interface, respectively. The sensitivities of the optimal shape on ionic conductivities and reaction rates are evaluated. 2. Physical model for SOFC cathode 2.1. Computational domain In the present study, the computational domain is divided into two regions, i.e. the pure LSCF electrode region Ωelectrode and the dense GDC electrolyte region Ωelectrolyte . These two regions are divided by the electrode-electrolyte interface Γinterface . In the electrode region, the porous microstructure is not explicitly resolved, but simply treated as a uniform continuum phase. For the pure LSCF cathode, homogenized effective ionic conductivity and DPB area density are applied throughout the electrode region. LSCF/GDC/pore TPB reaction takes place at the electrode-electrolyte interface, and LSCF/pore DPB reaction takes place inside the electrode region. The method to obtain cathode microstructural parameters such as effective ionic conductivity and DPB area density from FIB-SEM reconstructed 3D microstructures obtained in Ref. [22]. As shown in Fig. 2 (a), a level-set function ΦðxÞ is used to quantify the interface between the electrode and the electrolyte regions, where x represents the spatial position. The level set function ΦðxÞ is defined as the distance between the position x and the interface, i.e. the level set function becomes ΦðxÞ ¼ 0 at the interface, ΦðxÞ < 0 inside the elec­ trode, and ΦðxÞ > 0 inside the electrolyte: 8 � < ΦðxÞ > 0 x 2 ​ Ωelectrolyte� ; (1) ΦðxÞ ¼ 0 x 2 ​ Γinterface ; : ΦðxÞ < 0 ðx 2 ​ Ωelectrode Þ; in this study, the electrode-electrolyte interface shape is chosen as the design variable to be optimized, and the interface shape is deformed by modifying the spatial distribution of the level-set function ΦðxÞ.

Fig. 2. (a) Schematics of the computational domain and the level-set function, (b) distributions of ionic conductivity with diffusive interface, (c) distributions of DPB electrochemical reaction current with diffusive interface, and (d) dis­ tributions of TPB electrochemical reaction current with diffusive interface.

2.2. Governing equations The current simulation is conducted at an operation temperature of T ¼ 700� C in pure oxygen. The electronic conductivity of the electrode 3

Journal of Power Sources 449 (2020) 227565

A. He et al.

Table 1 Microstructural parameters and ionic conductivities used in the present computation. VLSCF

τLSCF

σO2

;LSCF (S/m)

0:681

1:40

0.474

σ O2

� � � α1;TPB F ηact ireac;TPB ¼ i0; ​ TPB ρTPB exp RT

HðΦÞ in Eq. (6) is the Heaviside function expressed in Eq. (4); δI ðΦÞ in Eq. (7) is the Dirac delta function: � δI ðΦÞ ¼ 0; Φ 6¼ 0 : (9) δI ðΦÞ ¼ þ∞; Φ ¼ 0

~O2 is the electrochemical potential of the oxide ion, σO2 is the where μ ionic conductivity and ireac is the reaction current.

Following condition should be satisfied for Dirac delta function: Z

σ eff O2 ;electrode

;electrolyte

⋅ HðΦÞ þ σ eff O2

;electrode

⋅f1

To facilitate the adjoint optimization, linearized Butler-Volmer equations are applied in the present study. It is reported that the error which comes from the linearization of the Butler-Volmer equation re­ mains within 4% when ηact < 0:1 V ​ [25].

(3)

HðΦÞg;

is the ionic conductivity of the electrolyte,

ireac;DPB ðΦÞ ¼ i ​ 0;DPB ρ ​ DPB ðα1; ​ DPB

α2; ​ DPB Þ

ireac;TPB ðΦÞ ¼ i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

α2; ​ TPB Þ

is the effective ionic conductivity of the electrode, and HðΦÞ

is the Heaviside function which is expressed as: � HðΦÞ ¼ 1; Φ � 0 HðΦÞ ¼ 0; Φ < 0

electrolyte ( Φ > 0) and σO2 ¼

σ eff O2 ;electrode

;electrolyte

in the

(Φ < 0) in the electrode.

Since pure LSCF cathode and GDC electrolyte are adopted in this study, LSCF the conductivities can be written as σeff ¼ VτLSCF σO2 ;LSCF and O2 ;electrode ;electrolyte

¼ σ O2

;GDC .

F η ⋅ð1 RT act

HðΦÞÞ;

F η ⋅δI ðΦÞjrΦj: RT act

(11) (12)

The value of ireac;DPB in Eq. (11) is 0 in the electrolyte (Φ > 0) and F becomes ireac;DPB ¼ i ​ 0;DPB ρ ​ DPB ðα1; ​ DPB α2; ​ DPB Þ RT ηact in the electrode (Φ < 0). The value of ireac;TPB in Eq. (12) is 0 in the electrolyte and electrode (Φ 6¼ 0) and becomes ireac;TPB ðΦÞ ¼ i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB α2; ​ TPB Þ F η ⋅δ ðΦÞjrΦj only at the interface (Φ ¼ 0). RT act I Exchange current densities and transfer coefficients are obtained from Ref. [8]. Total DPB area and TPB length are obtained from the reconstructed 3D microstructure by the marching cubes method [26] and the centroid method [27], respectively. DPB density is the total DPB area divided by the volume of the reconstructed structure. TPB density is the total TPB area divided by the cathode-electrolyte interface area of the reconstructed structure. The values for i ​ 0;DPB , ρ ​ DPB , α1; ​ DPB , α2; ​ DPB , i ​ 0;TPB , ρ ​ TPB , α1; ​ TPB and α2; ​ TPB are shown in Table 2.

(4)

the value of σ O2 ðΦÞ in Eq. (3) changes as σ O2 ¼ σO2

σ O2

(10)

∞

The ionic conductivity in the whole computational domain with respective to the level-set function Φ is written as:

where σ O2

þ∞

δI ðΦÞdΦ ¼ 1

2.3. Ionic conductivity

;electrolyte

�� ⋅δI ðΦÞjrΦj;

where i0;DPB is the exchange current density for DPB reaction (A=m2 ), ​ i0;TPB is the exchange current density for TPB reaction (A=m), ρDPB is the DPB area density (m2 =m3 ), ρTPB is the TPB area density (m=m2 ), α1 and α2 are the transfer coefficients. ηact is the local activation over­ potential [24] which is expressed as: � � �� 1 1 p 2~ ηact ¼ μe μ~O2 þ RTlog O2o;gas : (8) 2F 2 p

(2)

r~ μO2 ¼ ireac ðΦÞ;

σ O2 ðΦÞ ¼ σO2

RT

ηact

(7)

5.43

σO2 ðΦÞ 2F

α2;TPB F

;GDC (S/m)

is several orders of magnitude larger than the ionic conductivity of the electrolyte material [22]. Therefore, the electronic electrochemical po­ tential can be assumed to be uniform over the electrode region. In addition, according to Refs. [22,23,29], the oxygen gas partial pressure is almost uniform throughout the electrode region when thin electrode is applied. Thus, gas diffusion is not solved in the present study. Accord­ ingly, only the transport equation of oxide ion is solved: r

� exp

VLSCF and τLSCF are the volume fraction and

tortuosity factor of the LSCF phase, which are computed according to the 3D microstructure obtained in Ref. [22]. σ O2 ;LSCF and σ O2 ;GDC are the conductivities of LSCF and GDC, which are obtained from Ref. [8]. The values for VLSCF , τLSCF , σ O2 ;LSCF and σO2 ;GDC are shown in Table 1.

2.4. Electrochemical reaction

2.5. Boundary conditions

For pure LSCF cathode, LSCF/pore DPB reaction is considered inside the electrode region, and LSCF/GDC/pore TPB reaction is considered at the electrode electrolyte interface. Butler-Volmer equations are used to express the DPB and TPB reaction currents, which are given as:

A Dirichlet boundary condition is imposed at the anode-electrolyte interface:

μ~O2 ¼ μ~O2

(5)

ireac ¼ ireac;DPB þ ireac;TPB ; � � � α1;DPB F ireac;DPB ¼ i0; ​ DPB ρDPB exp ηact RT

�� α2;DPB F ηact ⋅ð1 exp RT

(13)

;B ;

~O2 ;B is given as 0:1 � 2F (J/mol). At the interface between the where μ cathode and the current collector, a zero-flux boundary condition is applied:

�

HðΦÞÞ; (6)

∂μ~O2 ¼ 0; ∂n

(14)

Table 2 Parameters for the reaction current calculations. i0;TPB (A/m)

i0;DPB (A/m2)

7:95 � 10 6

40:99

ρTPB (m/m2) 6

1:22 � 10

ρDPB (m2/m3) 6

2:84 � 10

4

α1; ​ DPB 1.0

α2; ​ DPB 1.2

α1;TPB 2.0

α2; ​ TPB 2.0

A. He et al.

Journal of Power Sources 449 (2020) 227565

4. Numerical optimization The target of the present study is to find the optimal electrode electrolyte interface structure which can maximize the electrochemical performance. The total amount of reaction current is chosen as the objective function to be maximized: Z ~ O2 Þ ¼ ðireac;DPB ðΦÞ þ ireac;TPB ðΦÞÞdV: (21) Jc ðΦ; μ an iterative method proposed in Ref. [17] is applied to obtain an optimal solution for ΦðxÞ. The flowchart for this iterative method is shown in Fig. 3. To start with, an initial interface structure is given, i.e. the initial value of the level-set function Φðm¼0Þ is calculated according to the given initial interface structure, where the superscript m represents the ~ðmÞ number of updates. Then, the distribution of μ ðxÞ is obtained by O2 solving Eq. (2), this step is called forward analysis. After the forward

~*;ðmÞ analysis, the distribution of adjoint variable, μ ðxÞ, is obtained by O2 solving following equation. This step is called adjoint analysis. r

where n is the direction from electrolyte to current collector. Symmetric boundary conditions are applied at other boundaries.

In the numerical simulation, the Heaviside function and the Dirac delta function are smoothed in order to facilitate the optimization: 1 � Φ� Hsmooth ðΦÞ ¼ 1 þ tanh ; (15) 2 ξ tanh2

Φ� ; ξ

(16)

;electrode

ireac;DPB ðΦÞ ¼ i0; ​ DPB ρDPB ðα1; ​ DPB

α2; ​ DPB Þ

ireac;TPB ðΦÞ ¼ i0; ​ TPB ρTPB ðα1; ​ TPB

α2; ​ TPB Þ

⋅ð1

F η ⋅ð1 RT act

Hsmooth ðΦÞÞ; Hsmooth ðΦÞÞ;

F η ⋅δI;smooth ðΦÞjrΦj; RT act

F ⋅ ð1 RT

i*reac;DPB ¼ i ​ 0;DPB ρDPB ðα1; ​ DPB

α2;DPB Þ

i*reac;TPB ¼ i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

α2; ​ TPB Þ

Hsmooth ðΦÞÞ⋅η*act ;

F ⋅ δI;smooth ðΦÞjrΦj⋅η*act ; RT

(24) (25)

(26)

where ΔτI is a discrete increment of a fictitious time, and SðmÞ is the sensitivity, the procedure for obtaining the sensitivity is introduced in Appendix D, and it can be regarded as the combination of the effects from transport Strans , DPB reaction SDPB , and the TPB reaction STPB , which are expressed as follows:

therefore, the choice of δI;smooth matches with the requirements for Dirac delta function given by Eq. (10). The ionic conductivity and the reaction currents are expressed using smoothed Hsmooth and δI;smooth as: ⋅ Hsmooth ðΦÞ þ σeff O2

(23)

Φðmþ1Þ ¼ ΦðmÞ þ ΔτI ⋅SðmÞ ;

∞

;electrolyte

(22)

where η*act ¼ ð~ μ*O2 1Þ=2F. Then, the interface shape is updated according to the following equation:

where ξ is the characteristic length of the interface layer. In this study, δI;smooth is defined as the derivative of Hsmooth with respective to the levelset function Φ, which satisfies the following condition: Z þ∞ (17) δI;smooth ðΦÞdΦ ¼ Hsmooth ðþ∞Þ Hsmooth ð ∞Þ ¼ 1:

σ O2 ðΦÞ ¼ σO2

r~ μ*O2 ¼ i*reac ;

i*reac ¼ i*reac;DPB þ i*reac;TPB ;

3. Diffusive interface

δHsmooth 1� 1 ¼ 2ξ δΦ

2F

the derivation of Eq. (22) is shown in Appendix C. In Eq. (22), ionic conductivity σ O2 is given by Eq. (18). The source term i*reac , which represents the electrochemical reaction for adjoint variable, can be ~*O2 ðxÞ: written as a function of adjoint electrochemical potential μ

Fig. 3. Flowchart of the optimization process.

δI;smooth ¼

σO2 ðΦÞ

(27)

S ¼ Strans þ SDPB þ STPB ; �

(18) (19)

Strans ¼

σ O2

SDPB ¼

1

STPB ¼

(20)

þ 1

Fig. 2 (b)–(d) shows the distributions of σ O2 ðΦÞ, ireac;DPB ðΦÞ, and ireac;TPB ðΦÞ with diffusive interface. As can be seen from the figure, dis­ tributions of these electrochemical parameters change smoothly across the interface between the electrode and electrolyte regions. In addition, the validation of the numerical simulation for solving Eq. (2) with diffusive interface is shown in Appendix A.

;electrolyte

σ eff O2

� δH ;electrode

smooth

δΦ

�

μ~*O2 i ​ 0;DPB ρ ​ DPB ðα1; ​ DPB

i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

α2; ​ TPB Þ

�

μ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

α2; ​ TPB Þ

2F

α2; ​ DPB Þ

F h r 1 RT

α2; ​ TPB Þ F 1 RT

μ~O2

r

μ*O2 ; r~

~O2 δHsmooth F μ ; RT 2F δΦ

μ~*O2

(28) (29)

i rΦ �μ ~ O2 δI;smooth ðΦÞ 2F jrΦj

~O2 δδI;smooth ðΦÞ F μ jrΦj RT 2F δΦ

μ~*O2

�μ ~ O2 rΦ : δI;smooth ðΦÞr 2F jrΦj (30)

rΦ it is reported that local curvature can be expressed as rjrΦj [28].

Therefore, the above set of equations mathematically demonstrates that

5

A. He et al.

Journal of Power Sources 449 (2020) 227565

Jflat Fig. 4. (a) Transformation of the electrode-electrolyte interface shape during optimization, (b) Changes in the improvement factor (J Jflat ) against iteration steps, (c) Oxide ion electrochemical potential distributions and reaction current distributions of flat interface (left) and optimized interface (right).

6

A. He et al.

Journal of Power Sources 449 (2020) 227565

0:04, k ¼ 0:0426, λDPB ¼ 2:87, and λTPB ¼ 0:0501. Uniform computa­

tional grid Δx ¼ 100nm is applied in the present study, where 250. The thickness of the interface is chosen as A GDC pillar with a width of

linit;width lelectrolyte

ξ Δx

¼ 3.

¼ 4 and a length of

lelectrode Δx

linit;length lelectrolyte

¼

¼ 10

on a flat GDC surface is given as the initial interface structure. Fig. 4 (a) shows the transformation of the electrode-electrolyte interface shape during optimization. Fig. 4 (b) shows the change in the reaction current improvement during optimization. The reaction current improvement Jflat during optimization with respective to a flat interface (J Jflat ) is moni­

tored. From 0 to 10,000 steps, according to Fig. 4 (a), the pillar-like interface extended into the electrode region, while it remains rela­ tively smooth. At the 10,000th step, the improvement factor increased from 0.09 to 0.37. From 10,000 to 100,000 steps, the electrodeelectrolyte interface became rough. At this stage, the improvement factor increased from 0.37 to 0.70. At the 100,000th step, the optimized structure shows improvement factor as high as 0:70. The low ionic po­ tential region broadens for the final structure as shown in Fig. 4 (c). It is observed in Fig. 4 (a) that the cross-sectional area of the pillarlike electrolyte region of the optimized structure reduces along the cathode thickness direction. This phenomenon is consistent with the observation by Chueh et al. [7], in which reducing cross-sectional area of the pillar-like electrolyte region along the thickness direction more effectively improves the performance than the constant cross-sectional area case. Chueh et al. attributed this phenomenon to the effect of gas transport. However, in the present study, uniform oxygen pressure is assumed, while only the oxide ion transport is solved. This indicates that the reduction of the cross-sectional area of the pillar-like electrolyte region is attributed to the effects of ionic transport. Fig. 5 shows inter­ face structures with different cross-sectional areas along the thickness direction. In the figure, A, B and C represent different positions in the computational domain. For Fig. 5 (a), the cross-sectional area in the A-B region is larger than that in B–C. Thus, the resistance of electrolyte re­ gion in A-B region (R ΔR 2 ) is smaller than the resistance in B–C region (R þ ΔR 2 ). For Fig. 5 (b), the resistances of A-B and B–C are kept constant. For Fig. 5 (c), the resistance of A-B region (R þ ΔR 2 ) is larger than the resistance of B–C region (R ΔR 2 ). Since the ionic current decreases along the thickness direction due to electrochemical reaction, ionic current in A-B region (I) is larger than that in B–C region (I ΔI). Thus, the elec­ trochemical potential difference (Δ~ μO2 ) from A to C for the above three cases can be expressed as follows:

Fig. 5. Schematics of the resistances and ionic currents inside the pillar-like electrolyte regions with (a) reducing cross-sectional area, (b) constant crosssectional area and (c) increasing cross-sectional area along the thick­ ness direction.

the sensitivity STPB is related with the local curvature of the interface which appears in the last term of Eq. (30).

Δ~ μAC;reducing ¼ Δ~μAB;reducing þ Δ~μBC;reducing � � � � ΔR ΔR þ ðI ΔIÞ ⋅ R þ ¼ 2IR ¼I⋅ R 2 2

5. Computational results 5.1. Optimization of the interface between pure LSCF cathode and GDC electrolyte

(32)

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffieff σ O2 ;electrode RT λDPB ¼ lelectrode ; i ​ 0;DPB ρ ​ DPB ðα1; ​ DPB α2;DPB ÞF

(33)

, λTPB ¼ lelectrode

σeff O2

;electrode

i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

RT

α2; ​ TPB ÞF

;

ΔR ΔI⋅ ; 2

(35)

Δ~ μAC;constant ¼ Δ~μAB;constant þ Δ~μBC;constant ¼ I ⋅ R þ ðI

As indicated in Appendix B, for a flat interface shown in Fig. A.1(b), the electrochemical potential distribution is governed by four dimen­ sionless numbers, L, k, λDPB , and λTPB , which are defined as: � L ¼ lelectrolyte lelectrode (31)

σ eff O2 ;electrode k¼ ; σ O2 ;electrolyte

R ⋅ ΔI

ΔIÞ ⋅ R ¼ 2IR

(36)

R⋅ΔI;

Δ~ μAC;increasing ¼ Δ~μAB;increasing þ Δ~μBC;increasing � � � � ΔR ΔR þ ðI ΔIÞ ⋅ R ¼ 2IR ¼ I⋅ Rþ 2 2

ΔR R ⋅ ΔI þ ΔI⋅ : 2

(37)

Among these three cases, Δ~ μAC;reducing shown by Eq. (35) is the smallest. Thus, pillar-like electrolyte region with reducing crosssectional area is preferable. Furthermore, since the properties inside the electrode region are uniform in the present study, the gas diffusion path from the current collector side of the optimized structure will become shorter than the flat interface case due to the intrusion of the interface into the electrode. In addition, with the increase in effective ionic conductivity and with the resultant increase in reactive volume, it is expected that the local reac­ tion current at TPB can be reduced, which leads to the reduction of local overpotential. This may contribute to enhance the durability of the electrode, since the local oxygen potential gradient becomes more mild.

(34)

for the present pure LSCF cathode, the above parameters become L ¼ 7

A. He et al.

Journal of Power Sources 449 (2020) 227565

Fig. 6. Optimal interface shapes for (a) different k values (different electrolyte conductivities), (b) different λDPB values (different DPB electrochemical reaction rates), and (c) different λTPB values (different TPB electrochemical reaction rates).

8

A. He et al.

Journal of Power Sources 449 (2020) 227565

On the other hand, it is reported that the gas diffusion in the cathode becomes very important especially under the interconnector rib where the gas channel is blocked by the rib [30]. Gas diffusion should be considered for such cases, which is under development as an ongoing work. 5.2. Effects of material properties and microstructural parameters In practice, if different material, microstructure or operation condi­ tions are applied, the conductivities and reactive site densities will vary. According to Appendix B, the electrode performance is governed by four dimensionless parameters, L, k, λDPB , and λTPB . In this study, lelectrolyte and lelectrode are fixed as 1 μm and 25 μm. For this condition, only three dimensionless parameters are to be investigated, i.e., k, λDPB and λTPB . The sensitivities of the electrode-electrolyte interface structure on these dimensionless parameters are investigated by changing these parame­ ters independently. In this study, σ O2 ;electrolyte , i0;DPB ρDPB and i0;TPB ρTPB are changed in order to vary these dimensionless parameters independently.

5.2.1. Influence of electrolyte ionic conductivity The LSCF cathode with GDC electrolyte in Section 5.1 is regarded as the reference case. In this section, numerical optimization are conducted under k ¼ 0:0213, 0.0426, and 0:0852, which are characterized with 0.5, 1, and 2 times of the reference k value, respectively. The parameter k can be changed independently by varying the electrolyte ionic con­ ductivity (σO2 ;electrolyte ) which depends on the material and operation condition [18]. On the other hand, λTPB and λDPB are fixed as λTPB ¼ 0:0501 and λDPB ¼ 2:87. Fig. 6 (a) shows the optimized microstructures. Fig. 7 (a) shows the cross sectional area Aelectrolyte =Adomain , perimeter Pinterface =Pdomain , and Joptimized Jflat ) Jflat

improvement factor (

for different k values. Here, Aelectrolyte

and Pinterface are the cross-sectional area and the perimeter of the pillarlike electrolyte region, and Jflat is the total reaction current for the flat interface case with the corresponding electrochemical parameters. As shown in Fig. 7 (a), the deformation of the electrode-electrolyte inter­ face becomes smaller as k gets smaller (larger σO2 ;electrolyte ), which in­ dicates that if high ionic conductivity electrolyte is used, smaller deformation of electrode-electrolyte interface is required. This is because when electrolyte material is more conductive, small pillar-like electrolyte region is enough to attain sufficient effective conductivity. As shown in Fig. 7 (a), the improvement factor increases with the decrease in k (increase in σO2 ;electrolyte ), which indicates that if high ionic

conductivity electrolyte material is applied, the electrochemical per­ formance can be more effectively improved by the deformation of the interface. This is because the electrolyte material with larger σ O2 ;electrolyte requires small deformation which do not sacrifice the re­ gion which was originally electrochemically reactive.

5.2.2. Influence of DPB reaction rate Parameter λDPB can be changed independently by varying the DPB reaction rate, i0;DPB ρDPB . In this section, calculations are conducted for λDPB ¼ 2:03, 2.87 and 4:06, which are characterized with 0.5, 1, and 2 times of the reference i0;DPB ρDPB value, respectively. It is mentioned that i0;DPB ρDPB can be changed in reality by applying different electrode materials [20], different electrode particle sizes [21] or changing the operation conditions [19]. On the other hand, λTPB and k are fixed as λTPB ¼ 0:0501 and k ¼ 0:0426. Fig. 7 (b) shows the cross sectional area Aelectrolyte =Adomain , perimeter Pinterface =Pdomain , and improvement factor for different λDPB values. As shown in Fig. 7 (b), the improvement factor of the optimal structure increases with λDPB (increase in the DPB reaction rate), which indicates that if larger DPB reaction rate is realized, the electrochemical performance can be improved by the deformed inter­ face. According to Fig. 7 (b), for the microstructure with small λDPB

Fig. 7. Cross sectional areas Aelectrolyte =Adomain , perimeters Pinterface = Pdomain and performance improvements of optimal structures with (a) different k values (different electrolyte conductivities), (b) different λDPB values (different DPB electrochemical reaction rates), and (c) different λTPB values (different TPB electrochemical reaction rates).

9

A. He et al.

Journal of Power Sources 449 (2020) 227565

(small i0;DPB ρDPB ), the size of the deformed electrolyte-electrolyte interface is small. This is because small DPB reaction rate, i.e. rela­ tively large ionic conductivity, results in thick reactive layer. For thick reactive layer, only a small amount of interface deformation is required to provide enough ionic flux. Further deformation of the interface re­ places the electrode region which was originally reactive. On the other hand, according to Fig. 7 (b), for the microstructure with large λDPB (large i0;DPB ρDPB ), the size of the interface deformation is large. Large i0;DPB ρDPB results in thin reactive layer. For thin reactive layer, large interface deformation is required in order to extend the reactive layer.

electrolyte interface to maximize the total amount of reaction current in the whole computational domain. The proposed model is applied to optimize the interface shape between the pure LSCF cathode and the GDC electrolyte. The reaction current of the cathode with optimal electrode-electrolyte interface structure is largely enhanced compared with that of the flat electrode-electrolyte interface. The optimal elec­ trolyte region structure is characterized with a pillar shape with reducing cross-section area along the thickness direction. In addition, the sensitivities of the optimal interface shape on the electrolyte ionic conductivity and reaction rates are investigated. The electrochemical performance can be more effectively improved by the interface defor­ mation for cases with higher ionic conductive electrolyte, larger DPB reaction rate in the electrode bulk, and larger TPB reaction rate at the interface. Optimal structure has larger interface deformation for less ionic conductive electrolyte or larger DPB reaction rates. And optimal interface becomes rough for large TPB reaction electrodes. The current approch is expected to give helpful guidance for the SOFC electrode design and manufacturing.

5.2.3. Influence of TPB reaction rate Parameter λTPB can be changed independently by varying the TPB reaction rate, i0;TPB ρTPB . In this section, numerical simulations are con­ ducted for λTPB ¼ 0:0, 0:0501 and 0:1002, which are characterized with 0, 1, and 2 times of the reference i0;TPB ρTPB value, respectively. In reality, TPB reaction rate can be improved by reducing the particle size [21], by introducing more reactive materials [20] or by changing the operation conditions [19]. The λTPB and k are fixed as λDPB ¼ 2:87 and k ¼ 0:0426. Fig. 6 (c) shows the optimal structures for different λTPB values. Fig. 7 (c) shows Aelectrolyte =Adomain , Pinterface =Pdomain , and improvement factors. According to Fig. 7 (c), cross sectional area Aelectrolyte is the same for different TPB reaction rates, which implies that the amount of interface deformation is independent from the TPB reaction rate. While the interface surface area increases with the TPB reaction rate, which in­ dicates that rougher electrode-electrolyte interface is preferred for the enhanced TPB reaction rate electrodes. It is observed in Fig. 7 (c) that the improvement factor of the optimized structure increased with the in­ crease of TPB reaction rate, which implies that rough interface defor­ mation is more effective for enhanced TPB reaction rate electrodes.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This research was partly supported by MEXT “Priority Issue on PostK computer” project and by NEDO “Technology Development for Pro­ motion of Applications of Solid Oxide Fuel Cells” project. This research used computational resources of the K computer provided by the RIKEN Center for Computational Science through the HPCI System Research project (Project ID: hp150280, hp160226, hp170253, hp180187, hp190178).

6. Conclusions A numerical method is developed to optimize the shape of electrode-

Appendix A. Validation of the numerical model Figure A.1 (c) shows the experimental and computational results for pure LSCF cathode at 700� C with pure oxygen. The measurements were carried out in our previous study [22]. The computational result is obtained by solving transport equation, Eq. (2), with a flat interface shown in Fig. A.1(b). The structural parameters for the computational result are shown in Tables 1 and 2, which are obtained from the same sample used for the experiment [22]. At the anode-electrolyte interface, Dirichlet boundary is applied, and the ionic potential is varied in order to correspond to the experimental cathode overpotentials (ηcathode ), which is given as: � � 1 1 pO ;cc RTln 2o ; ηcathode ¼ μ~O2 ;cathode=electrolyte 2~μe ;cathode=cc (A.1) 2F 2 p where the subscript “cathode/electrolyte” represents the interface between the cathode microstructure and the electrolyte layer, while the subscript “cathode/cc” represents the interface between the cathode microstructure and the current collector, po ¼ 1atm is the reference pressure. At other boundaries, Newmann boundary with zero-flux is applied. The reaction current ireac can be obtained from the numerical simulation, and the current density is calculated from the following equation [31]: Z lelectrode I¼ ireac dx (A.2) 0

According to Fig. A.1 (c), there is only a slight difference between the simulated and measured current densities, which reveals the validity of the present simulation.

10

A. He et al.

Journal of Power Sources 449 (2020) 227565

Fig. A.1. (a) Reconstructed structure from FIB-SEM, (b) schematic of the computational volume with flat electrode-electrolyte interface, and (c) comparison between the predicted current density and the experimental result.

Appendix B. Dimensionless parameters for the electrochemical reaction in an electrode/electrolyte system with a flat interface A flat interface between an electrolyte and an electrode as shown in Fig. A.1(b) is considered. In this case, the ionic conduction equation can be reduced into one dimension, which is expressed as follows: ~O2 ;electrolyte d2 μ ¼0 dx2

~O2 ;electrode d2 μ 1 ¼ 2 μ~O2 dx2 lreac;DPB ~ O2 where μ

(B.1)

ðx < 0Þ;

;electrolyte

;electrode

~O2 and μ

(B.2)

ðx > 0Þ:

;electrode

are the electrochemical potential of oxide ion in electrolyte and electrode regions, x is the cathode thickness

direction, and x ¼ 0 indicates the electrode-electrolyte interface. The quantity lreac;DPB is the characteristic length scale for electrode region elec­ trochemical reaction which takes place inside the electrode bulk region: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kσ O2 ;electrolyte RT lreac;DPB ¼ ; (B.3) i ​ 0;DPB ρ ​ DPB ðα1; ​ DPB α2; ​ DPB ÞF where k is defined as: k¼

σ eff O2 ;electrode : σ O2 ;electrolyte

(B.4)

The electrochemical potential must follow the boundary conditions as explained in Eqs. (15) and (16). For one-dimension problem, the boundary condition is expressed as:

11

A. He et al.

μ~O2

;electrolyte

d~ μ O2

;electrode

dx

Journal of Power Sources 449 (2020) 227565

~ O2 ¼μ

;B

at x ¼

(B.5)

lelectrolyte ;

(B.6)

¼ 0 at x ¼ lelectrode ;

~O2 ;B is the electrochemical potential at the electrolyte boundary. In addition, the boundary condition for the interface (x ¼ 0) between the where μ electrolyte and the anode can be expressed as:

μ~O2

;electrolyte

~ O2 ¼μ

;electrode

~ O2 ¼μ

;I

(B.7)

at x ¼ 0;

~O2 ;I is the electrochemical potential at the electrode-electrolyte interface, which is to be determined later. where μ Because of the existence of the interface TPB reaction, the following relationship holds: d~ μ O2

;electrolyte

dx

¼k

d~ μ O2

k μ~ 2 ; ðx ¼ 0Þ; lreac;TPB O ;I

;electrode

dx

(B.8)

where lreac;TPB is the characteristic length scale for the interface TPB reaction which is expressed as: lreac;TPB ¼

σeff O2

;electrode

i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

RT

α2; ​ TPB ÞF

(B.9)

:

Equation (B.1) can be solved with the boundary conditions Eq. (B.5) and (B.8): � � μ~O2 ;electrolyte ðxÞ μ~O2 ;B x ¼1 1 : lelectrolyte μ~O2 ;I μ~O2 ;I

(B.10)

Equation (B.2) can be solved with the boundary conditions Eq. (B.6) and (B.8)

μ~O2

;electrode ðxÞ

μ~O2

¼

expð

;I

λDPB ÞexpðλDPB x=lelectrode Þ þ expðλDPB Þexpð expðλDPB Þ þ expð λDPB Þ

λDPB x=lelectrode Þ

;

where λDPB is the dimensionless parameter defined as: � λDPB ¼ lelectrode lreac;DPB : ~O2 ;electrode ðxÞ and μ ~ O2 Taking the derivatives of μ ~O2 ;I , can be expressed as: at the interface, μ

μ~O2 μ~O2

;B

;electrolyte ðxÞ

(B.11)

(B.12) with respect to x, and substituting the results into Eq.(B.8), the electrochemical potential

(B.13)

¼ L ⋅ kλDPB tanhðλDPB Þ þ L ⋅ kλTPB þ 1;

;I

where L and λTPB is expressed as � L ¼ lelectrolyte lelectrode ;

(B.14)

� λTPB ¼ lelectrode lreac;TPB :

(B.15)

The above analysis shows that the ionic potential distribution for flat interface is governed by four dimensionless parameters: L, k, λDPB and λTPB which are expressed in Eqs. (B.14), (B.4), (B.12) and (B.15), respectively. In this study, lelectrolyte and lelectrode are fixed as 1 μm and 25 μm, respectively. Thus, in the present condition, only three dimensionless parameters are to be investigated, i.e., k, λDPB and λTPB . Appendix C. Derivation of governing equations for adjoint analysis In the present study, the objective function, Eq. (21), is maximized with respect to the spatial distribution of the level set function ΦðxÞ. Any change ~O2 from the governing in the interface structure expressed by ΦðxÞ will cause changes in the distribution of σO2 , which results in the change in μ equation, Eq. (2). Therefore, the process for seeking the optimal structure can be regarded as a constrained optimization problem. Here, the Lagrange multiplier method is employed to solve the present constrained optimization problem. The Lagrangian is introduced as follows: Z � ~ O2 ; μ ~*O2 ¼ Jc ðΦ; μ ~ O2 Þ þ L Φ; μ μ~*O2 cðΦ; μ~O2 ÞdV; (C.1) � σ ðΦÞ ~*O2 is the Lagrange multiplier and c Φ; μ ~O2 ¼ r O22F r~ where Jc is the objective function in Eq. (21), μ μO2 ireac ðΦÞ is the operator which comes from the constraint of transport equation for oxide ion, Eq. (2). It is noted that the values of L and Jc are identical as far as the constraint, Eq. (2), is � � ~ O2 ; μ ~*O2 is related satisfied over the entire domain. Therefore, the Lagrangian L is to be maximized as well as the objective function Jc . Since L Φ; μ ~O2 and μ ~*O2 , the optimized result is obtained if the following equations are satisfied: with Φ; μ δL ~O2 Þ ¼ 0; ¼ cðΦ; μ δ~ μ*O2 12

(C.2)

A. He et al.

Journal of Power Sources 449 (2020) 227565

~ O2 Þ ~ O2 Þ δL δJc ðΦ; μ δcðΦ; μ ~*O2 ¼ þμ ¼ 0; δ~ δ~ δ~ μ O2 μ O2 μ O2

(C.3)

~ O2 Þ ~ O2 Þ δL δJc ðΦ; μ δcðΦ; μ ~*O2 ¼ þμ ¼ 0: δΦ δΦ δΦ

(C.4)

Equations (C.2)-(C.4) cannot be solved in a direct way. We apply an iterative method shown in Fig. 3 to solve these equations. The initialization and forward analysis are introduced in Section 4. ~*O2 by solving the Eq. (C.3). However, Eq. (C.3) should be transformed in order to facilitate the Adjoint analysis is to obtain the distribution of μ computation. Lagrangian function can be rearranged as the combination of transport effect L trans , DPB electrode region reaction effect L DPB , and the TPB electrode-electrolyte interface reaction effect L TPB : � ~ O2 ; μ ~*O2 ¼ L trans þ L DPB þ L TPB ; L Φ; μ (C.5) Z h

� μ ~ 2 �i μ~*O2 r σr O dV;

¼

L

trans

L

DPB ¼

L

TPB

(C.6)

2F

Z �

μ~*O2 i ​ 0;DPB ρ ​ DPB ðα1; ​ DPB

1

μ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

Z � ¼

�

1

�

According to Ref. [17],

δL trans δ~ μO2

� Hsmooth ðΦÞÞ dV;

(C.7)

� ~ O2 F μ ⋅ δI; ​ smooth ðΦÞjrΦj dV: RT 2F

(C.8)

α2; ​ DPB Þ α2; ​ TPB Þ

~ O2 F μ ⋅ ð1 RT 2F

and δL DPB can be transformed into: δ~ μO2

δL trans 1 ¼ r⋅σ r~ μ*O2 ; 2F δ~ μ O2 δL DPB 1 1 ¼ 2F δ~ μ O2

(C.9)

�

μ~*O2 i ​ 0;DPB ρ ​ DPB ðα1; ​ DPB

α2;DPB Þ

F ⋅ð1 RT

(C.10)

Hsmooth ðΦÞÞ:

It should be noted that Eq. (C.9) can be satisfied if the following boundary conditions are applied [17]. At the electrolyte side, Dirichlet condition ~*O2 is applied, for μ

μ~*O2 ¼ 0;

(C.11)

~*O2 is applied, and for the other boundaries except for the electrolyte side, Neumann condition for μ � d~ μ*O2 dn ¼ 0:

(C.12)

δL TPB δ~ μO2

is obtained as the same methodology of

δL TPB 1 1 ¼ 2F δ~ μ O2 Thus,

δL μO2 δ~

�

μ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

α2; ​ TPB Þ

δL DPB : δ~ μO2

F ⋅δI;smooth ðΦÞjrΦj: RT

(C.13)

is expressed as the combination of Eqs. (C.5), (C.6), and (C.9). Accordingly, Eq. (26) is transformed into the following equation which

can be solved with the same method as the Eq. (2): r

σO2 ðΦÞ 2F

r~ μ*O2 ¼ i*reac ;

(C.14)

where the conductivity σO2 ðΦÞ is given by the same function as that for the forward analysis step. The source term i*reac , which is the analogy of the reaction current for the adjoint variable, can be written as a function of the adjoint electrochemical potential. i*reac ¼ i ​ 0;DPB ρ ​ DPB ðα1; ​ DPB

α2; ​ DPB Þ

þi ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

α2; ​ TPB Þ

where η*act ¼ ð~ μ*O2

1Þ=2F.

F ⋅ ð1 RT

Hsmooth ðΦÞÞ⋅η*act

F ⋅ δI;smooth ðΦÞjrΦj⋅η*act ; RT

(C.15)

Appendix D. Derivation of sensitivity For the shape update in Section 4, the sensitivity S is required to be computed. In this study, S is defined as the derivative of the Lagrangian function L which is defined in Eq. (C.1) with respective to level-set function Φ: S¼

δL : δΦ

(D.1)

In this section, the transformation of δL is introduced in order to facilitate the computation. According to Eq. (C.5) to (C.8), the Lagrangian function δΦ 13

A. He et al.

Journal of Power Sources 449 (2020) 227565

can be rearranged as the combination of transport effect L trans , electrode region reaction effect L DPB , and the electrode-electrolyte interface reaction effect L TPB ; therefore, the sensitivity can also be rearranged as the combination of transport effect Strans , electrode region reaction effect SDPB , and the electrode-electrolyte interface reaction effect STPB , which can be obtained using the calculus of variation: � δL trans Strans ¼ (D.2) ¼ ðL trans ðΦ þ δΦÞ L trans ðΦÞÞ δΦ; δΦ �

SDPB ¼

δL DPB ¼ ðL δΦ

DPB ðΦ þ δΦÞ

STPB ¼

δL TPB ¼ ðL δΦ

TPB ðΦ þ δΦÞ

δΦ;

(D.3)

δΦ:

(D.4)

DPB ðΦÞÞ

L

� TPB ðΦÞÞ

L

can be simplified into the following equation if the boundary conditions of Eqs. (C.11) and (C.12) are held: � δσ μ~ 2 �i L trans ðΦ þ δΦÞ L trans ðΦÞ ¼ μ~*O2 r δΦr O dV δΦ 2F Z h � Z h �i i δσ μ~ 2 δσ μ~ 2 ¼ r⋅ μ*O2 r O μ*O2 dV δΦ~ dV δΦr O r~ δΦ 2F δΦ 2F Z h i δσ μ~O2 * μO2 dV: (D.5) ¼ δΦr r~ δΦ 2F The ðL

trans ðΦ þδΦÞ

L

trans ðΦÞ Þ

Z h

Therefore, Strans can be expressed as: Strans ¼

μ~ δσ δΦr O2F2 δΦ

δL trans ¼ δΦ

r~ μ*O2

δΦ

� δH

� ¼

σ O2

;electrolyte

σ eff O2

smooth

;electrode

δΦ

μ~O2

r

2F

μ*O2 : r~

(D.6)

SDPB in Eq. (D.3) can be obtained in a direct way: SDPB ¼

δL DPB ¼ δΦ

The L

�

μ~*O2 i ​ 0;DPB ρ ​ DPB ðα1; ​ DPB

1

TPB ðΦ þδΦÞ

1

and L

TPB ðΦÞ

TPB ðΦÞ ¼

L

TPB ðΦ þ δΦÞ ¼

1

~O2 δHsmooth F μ : RT 2F δΦ

(D.7)

in Eq. (D.4) can be expressed with the following equations:

�

μ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

L

α2; ​ DPB Þ

α2; ​ TPB Þ

�

μ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

~ O2 F μ ⋅δI;smooth ðΦÞjrΦj; RT 2F

α2; ​ TPB Þ

(D.8)

~ O2 F μ ⋅δI;smooth ðΦ þ δΦÞjrðΦ þ δΦÞj: RT 2F

(D.9)

where δI ðΦ þδΦÞ is expressed by δI ðΦ þ δΦÞ ¼ δI ðΦÞ þ

δδI;smooth δΦ; δΦ

(D.10)

and jrðΦ þδΦÞj in Eq. (D.9) can be solved with Taylor expansion: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jrðΦ þ δΦÞj ¼ jrΦ þ rδΦj ¼ jrΦ þ rδΦj2 ¼ jrΦj2 þ 2rΦ⋅rδΦ ! sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2rΦ⋅rδΦ rΦ⋅rδΦ rΦ⋅rδΦ ¼ jrΦj þ ¼ jrΦj 1 þ : ¼ jrΦj ⋅ 1 þ jrΦj jrΦj2 jrΦj2

(D.11)

Thus, Eq. (D.9) can be transformed to: L

TPB ðΦ þ δΦÞ

�

¼ 1

μ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

¼ 1

μ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

�

Then, L

TPB ðΦ þδΦÞ

L

� � ~ O2 � F μ δδI;smooth � rΦ⋅rδΦ ⋅ δI;smooth ðΦÞ þ δΦ jrΦj þ RT 2F δΦ jrΦj � � ~ O2 F μ rΦ⋅rδΦ δδI;smooth δδI;smooth rΦ⋅rδΦ α2; ​ TPB Þ þ : δI;smooth ðΦÞjrΦj þ δI;smooth ðΦÞ δΦjrΦj þ δΦ RT 2F jrΦj jrΦj δΦ δΦ

α2; ​ TPB Þ

TPB ðΦÞ

L TPB ðΦÞ TPB ðΦ þ δΦÞ � ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB ¼ 1 μ

can be expressed as:

L

� � ~ O2 F μ rΦ⋅rδΦ δδI;smooth þ δI;smooth ðΦÞ δΦjrΦj RT 2F δΦ jrΦj � � � ~ O2 F μ δδI;smooth rΦ⋅rδΦ * ~O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB α2; ​ TPB Þ þ 1 μ δI;smooth ðΦÞ δΦ RT 2F δΦ jrΦj � � � ~ O2 F μ rΦ⋅rδΦ δδI;smooth * ~O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB α2; ​ TPB Þ ¼ 1 μ þ δI;smooth ðΦÞ δΦjrΦj RT 2F jrΦj δΦ � � � ~ O2 F μ rΦ * ~O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB α2; ​ TPB Þ δΦ ¼r 1 μ δI;smooth ðΦÞ RT 2F jrΦj

α2; ​ TPB Þ

14

(D.12)

A. He et al.

� δΦr 1

Journal of Power Sources 449 (2020) 227565

�

�

~ O2 F μ rΦ δI;smooth ðΦÞ RT 2F jrΦj � ~O2 δδI;smooth F μ ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB α2; ​ TPB Þ þ 1 μ δΦjrΦj RT 2F δΦ � � � ~ O2 F μ rΦ ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB α2; ​ TPB Þ ¼ δΦr 1 μ δI;smooth ðΦÞ RT 2F jrΦj � ~O2 δδI;smooth F μ * ~O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB α2; ​ TPB Þ þ 1 μ δΦjrΦj : RT 2F δΦ

μ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

α2; ​ TPB Þ

(D.13)

Thus, Eq. (D.4) can be expressed as: i rΦ �μ ~ O2 δI;smooth ðΦÞ 2F jrΦj �μ ~ O2 F rΦ * ~ O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB α2; ​ TPB Þ 1 μ δI;smooth ðΦÞr RT jrΦj 2F � ~ F μ δδ 2 I;smooth O ~*O2 i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB α2; ​ TPB Þ þ 1 μ jrΦj: RT 2F δΦ

STPB ¼

δL TPB ¼ δΦ

i ​ 0;TPB ρ ​ TPB ðα1; ​ TPB

α2; ​ TPB Þ

F h r 1 RT

μ~*O2

(D.14)

Accordingly, S is obtained by combining Eq. (D.6), (D.7) and (D.14).

References [17]

[1] S.C. Singhal, K. Kendall, High Temperature Solid Oxide Fuel Cell, Elsevier, 2002. [2] Y. Suzue, N. Shikazono, N. Kasagi, Micro modeling of solid oxide fuel cell anode based on stochastic reconstruction, J. Power Sources 184 (2008) 52–59, https:// doi.org/10.1016/j.jpowsour.2008.06.029. [3] C. Ding, T. Hashida, High performance anode-supported solid oxide fuel cell based on thin-film electrolyte and nanostructured cathode, Energy Environ. Sci. 3 (2010) 1729–1731, https://doi.org/10.1039/c0ee00255k. [4] M.E. Lynch, L. Yang, W. Qin, J.J. Choi, M. Liu, K. Blinn, M. Liu, Enhancement of La0.6Sr0.4Co0.2Fe0.8O3-δ durability and surface electrocatalytic activity by La0.85Sr0.15MnO3�δ investigated using a new test electrode platform, Energy Environ. Sci. 4 (2011) 2249–2258, https://doi.org/10.1039/c1ee01188j. [5] A. Konno, H. Iwai, M. Saito, H. Yoshida, A corrugated mesoscale structure on electrode-electrolyte interface for enhancing cell performance in anode-supported SOFC, J. Power Sources 196 (2011) 7442–7449, https://doi.org/10.1016/j. jpowsour.2011.04.051. [6] A. Bertei, F. Tariq, V. Yufit, E. Ruiz-Trejo, N.P. Brandon, Guidelines for the rational design and engineering of 3D manufactured solid oxide fuel cell composite electrodes, J. Electrochem. Soc. 164 (2017) F89–F98, https://doi.org/10.1149/ 2.0501702jes. [7] C.-C. Chueh, A. Bertei, C. Nicolella, Design guidelines for the manufacturing of the electrode-electrolyte interface of solid oxide fuel cells, J. Power Sources 437 (2019), 226888, https://doi.org/10.1016/j.jpowsour.2019.226888. [8] A. He, T. Shimura, J. Gong, N. Shikazono, Numerical simulation of La0.6Sr 0.4Co0.2Fe0.8O3 - Gd0.1Ce0.9O1.95 composite cathodes with micro pillars, Int. J. Hydrogen Energy 44 (2019) 6871–6885, https://doi.org/10.1016/j. ijhydene.2019.01.171. [9] T. Shimura, K. Nagato, N. Shikazono, Evaluation of electrochemical performance of solid-oxide fuel cell anode with pillar-based electrolyte structures, Int. J. Hydrogen Energy 44 (2019) 12043–12056, https://doi.org/10.1016/j.ijhydene.2019.03.112. [10] T. Okabe, Y. Kim, Z. Jiao, N. Shikazono, J. Taniguchi, Fabrication process for micropatterned ceramics via UV-nanoimprint lithography using UV-curable binder, Jpn. J. Appl. Phys. 57 (2018), https://doi.org/10.7567/jjap.57.106501. [11] T. Shimura, K. Nagato, N. Shikazono, Evaluation of electrochemical performance of solid-oxide fuel cell anode with pillar-based electrolyte structures, Int. J. Hydrogen Energy 44 (2019) 12043–12056, https://doi.org/10.1016/j.ijhydene.2019.03.112. [12] M. Zhi, S. Lee, N. Miller, N.H. Menzler, N. Wu, An intermediate-temperature solid oxide fuel cell with electrospun nanofiber cathode, Energy Environ. Sci. 5 (5) (2012) 7066–7071, https://doi.org/10.1039/c2ee02619h. [13] F. Delloro, M. Viviani, Simulation study about the geometry of electrodeelectrolyte contact in a SOFC, J. Electroceram. 29 (2012) 216–224, https://doi. org/10.1007/s10832-012-9766-8. [14] M. Geagea, J. Ouyang, B. Chi, F. Delloro, A. Chesnaud, A. Ringued� e, M. Cassir, A. Thorel, Architectured interfaces and electrochemical modelling in an anode supported SOFC, ECS Trans 68 (2015) 2961–2969, https://doi.org/10.1149/ 06801.2961ecst. [15] M. Secanell, R. Songprakorp, A. Suleman, N. Djilali, Multi-objective optimization of a polymer electrolyte fuel cell membrane electrode assembly, Energy Environ. Sci. 1 (2008) 378–388, https://doi.org/10.1039/b804654a. [16] H. Iwai, A. Kuroyanagi, M. Saito, A. Konno, H. Yoshida, T. Yamada, S. Nishiwaki, Power generation enhancement of solid oxide fuel cell by cathode-electrolyte interface modification in mesoscale assisted by level set-based optimization

[18] [19] [20] [21] [22]

[23]

[24]

[25] [26] [27]

[28] [29] [30]

[31]

15

calculation, J. Power Sources 196 (2011) 3485–3495, https://doi.org/10.1016/j. jpowsour.2010.12.024. J. Onishi, Y. Kametani, Y. Hasegawa, N. Shikazono, Topology optimization of electrolyte-electrode interfaces of solid oxide fuel cells based on the adjoint method, J. Electrochem. Soc. 166 (2019) F876–F888, https://doi.org/10.1149/ 2.0031913jes. D.J. Brett, A. Atkinson, N.P. Brandon, S.J. Skinner, Intermediate temperature solid oxide fuel cells, Chem. Soc. Rev. 37 (8) (2008) 1568–1578, https://doi.org/ 10.1039/b612060c. A. Esquirol, J. Kilner, N. Brandon, Oxygen transport in La0.6Sr0.4Co0.2Fe0.8O3 δCe0.8Ge0.2O2 x composite cathode for IT-SOFCs, Solid State Ion. 175 (1–4) (2004) 63–67, https://doi.org/10.1016/j.ssi.2004.09.013. E.V. Tsipis, V.V. Kharton, Electrode materials and reaction mechanisms in solid oxide fuel cells: a brief review, J. Solid State Electrochem. 12 (2008) 1367–1391, https://doi.org/10.1007/s10008-008-0611-6. C. Ding, T. Hashida, High performance anode-supported solid oxide fuel cell based on thin-film electrolyte and nanostructured cathode, Energy Environ. Sci. 3 (2010) 1729–1731, https://doi.org/10.1039/c0ee00255k. Y.T. Kim, Z. Jiao, N. Shikazono, Evaluation of La0.6Sr0.4Co0.2Fe0.8O3 Gd0.1Ce0.9O1.95 composite cathode with three dimensional microstructure reconstruction, J. Power Sources 342 (2017) 787–795, https://doi.org/10.1016/j. jpowsour.2016.12.113. M. Kishimoto, H. Iwai, K. Miyawaki, M. Saito, H. Yoshida, Improvement of the subgrid-scale model designed for 3D numerical simulation of solid oxide fuel cell electrodes using an adaptive power index, J. Power Sources 223 (2013) 268–276, https://doi.org/10.1016/j.jpowsour.2012.09.077. K. Matsuzaki, N. Shikazono, N. Kasagi, Three-dimensional numerical analysis of mixed ionic and electronic conducting cathode reconstructed by focused ion beam scanning electron microscope, J. Power Sources 196 (2011) 3073–3082, https:// doi.org/10.1016/j.jpowsour.2010.11.142. P. Costamagna, P. Costa, V. Antonucci, Micro-modelling of solid oxide fuel cell electrodes, Electrochim. Acta 43 (1998) 375–394, https://doi.org/10.1016/S00134686(97)00063-7. W.E. Lorensen, H.E. Cline, Marching cubes: a high resolution 3D surface, Construction Algorithm. Comput. Graphics 21 (1987) 163–169. N. Shikazono, D. Kanno, K. Matsuzaki, H. Teshima, S. Sumino, N. Kasagi, Numerical assessment of SOFC anode polarization based on three-dimensional model microstructure reconstructed from FIB-SEM images, J. Electrochem. Soc. 157 (2010) B665, https://doi.org/10.1149/1.3330568. S.J. Osher, Fronts propagating with curvature dependent speed, Comput. Phys. 79 (1988) 1–5, https://doi.org/10.1007/s13398-014-0173-7.2. Q. Cai, C.S. Adjiman, N.P. Brandon, Investigation of the active thickness of solid oxide fuel cell electrodes using a 3D microstructure model, Electrochim. Acta 56 (2011) 10809–10819, https://doi.org/10.1016/j.electacta.2011.06.105. H. Geisler, A. Kromp, A. Weber, E. Ivers-Tiffee, Stationary FEM model for performance evaluation of planar solid oxide fuel cells connected by metal interconnectors, J. Electrochem. Soc. 161 (2014) F778–F788, https://doi.org/ 10.1149/2.079406jes. K. Miyawaki, M. Kishimoto, H. Iwai, M. Saito, H. Yoshida, Comprehensive understanding of the active thickness in solid oxide fuel cell anodes using experimental, numerical and semi-analytical approach, J. Power Sources 267 (2014) 503–514, https://doi.org/10.1016/j.jpowsour.2014.05.112.