Modeling of solid oxide fuel cells with optimized interconnect designs

International Journal of Heat and Mass Transfer 125 (2018) 506–514

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Modeling of solid oxide fuel cells with optimized interconnect designs Shumao Zeng a, Xiaoqiang Zhang a, Jun Song Chen a,b, Tingshuai Li a,b,⇑, Martin Andersson a,c a

School of Materials and Energy, University of Electronic Science and Technology of China, 2006 Xiyuan Ave, West Hi-Tech Zone, 611731 Chengdu, Sichuan, PR China Center for Applied Chemistry, University of Electronic Science and Technology of China, 2006 Xiyuan Ave, West Hi-Tech Zone, 611731 Chengdu, Sichuan, PR China c Department of Energy Sciences, Faculty of Engineering, Lund University, P.O. Box 118, SE-221 00 Lund, Sweden b

a r t i c l e

i n f o

Article history: Received 13 February 2018 Received in revised form 9 April 2018 Accepted 19 April 2018

Keywords: Solid oxide fuel cells Thermal stress Interconnect Interface Cathode

a b s t r a c t A 3D model is developed to investigate solid oxide fuel cells (SOFCs) contacting with optimized interconnect designs and the results indicate that the current density and thermal stress are closely related to both the shape of tip in interconnects and the depth of it in the cathode. The interconnect with triangular rips can yield the best electrochemical performance compared to those with tips of rectangle and trapezium, and the current densities increase with the depth of tips in cathodes, except the trapezoidal ribs, which shows a concaving change with the depth. The 1st principle stress reaches around 21.9 MPa and 16.6 MPa at the interfaces of electrodes and electrolytes, but it rises to 60 MPa and 18 MPa for the rectangular tips at the air and fuel inlets, respectively, which sharply decreases to nearly 25 MPa and 10 MPa with the depth in cathodes approaching 5 lm. The maximum shear stresses are found to reach 34.4 MPa and 32.1 MPa at the two interfaces, and the triangular tips will induce the most intensive stresses at electrolyte-cathode interface. The resulting conclusions are beneficial to optimize interconnect design to improve the efficiency of current collection and also reduce the risk of generation of remarkable thermal stresses. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction The solid oxide fuel cells (SOFC) is a promising power generator that can directly convert chemical energy in fuels to electricity [1,2], and the efficiency and stability depend mostly on the conductivities of materials [3,4], catalyst activity [5–7], the compositions of fuels [8,9] and the design of interconnects [10,11]. The interconnect is not only used to connect single solid oxide fuel cells in series in a stack, but it also acts as gas channels to transport air and fuels, which can thus strongly affect both the efficiency of current collection and the whole stack’s stability due to mechanical mismatches between components. The optimization of interconnect can be thus a significant part in research and development of SOFCs [12]. The patterns of interconnect contacting cathodes can strongly affect the power density and stability of a SOFC stack [13,14]. It is reported that the maximum power output increased from 0.07 W cm
2 to 0.1 W cm
2 with the tip of interconnect increasing from 40 to 60 cm2 [15]. Moreover, the tip area and its depth in cathode can strongly affect the power density and degradation [16], which ⇑ Corresponding author at: School of Materials and Energy, University of Electronic Science and Technology of China, 2006 Xiyuan Ave, West Hi-Tech Zone, 611731 Chengdu, Sichuan, PR China. E-mail address: [email protected] (T. Li). https://doi.org/10.1016/j.ijheatmasstransfer.2018.04.096 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved.

also confirms that the optimum tip contact area was 45% (accounting for the whole area of an interconnect) and its depth in the cathode was 105 lm for achieving the best performance and duration property. Besides, Liu et al. studied effect of the interconnect rib width on the cell performance [17], which presented that the rib surface contact contributed to the ohmic polarization and the optimal rib width was linear to the pitch width. Based on previous experimental results, it can be concluded that the shape and size of tip or rib in interconnects tend to affect the performance and durability of a stack. The contact area and the depth of tip or rib in cathodes are both expected to be large for more efficient current collection. However, there are mechanical mismatches between the metal interconnects and ceramic cathodes, which will induce thermal stresses in a single cell. Thermal stresses in single SOFCs were widely studied as functions of the operation conditions [18–20], cell configuration [21], porosity in electrodes [22] and sealing design [23], but there is lack of detailed investigations into stresses stemming from the patterns of the interconnect-electrode contact. We previously evaluated distributions of thermal stresses by changing the interconnect-electrode contact area and the results demonstrated that a larger contact area is beneficial to effectively reduce the 1st principle stress close to the corner of interconnectcathode contact region at the fuel inlet and a wider rib at the anode

S. Zeng et al. / International Journal of Heat and Mass Transfer 125 (2018) 506–514

507

Nomenclature i

r U V T V EOCV E0 R F p Eeq L iv i0 Ai Ea k

qg

! u ! F xi bn DTi Deff;ij Sj D Mij lij

X M r r m cp;g Qh DS DH C

current density, A m
2 electric/ionic conductivity, S m
1 electric/ionic potential, V volume fraction, – temperature, K working voltage, V open circuit voltage, V ideal voltage before partial pressure consideration, V gas constant, 8.314 J mol
1 K
1 Faraday constant, 96,485 C mol
1 gas partial pressure, Pa equilibrium potential for anode/cathode, V characteristic thickness, m volumetric current density, A m
3 exchange current density, A m
2 pre-exponential factor, S m
2 activation energy, J mol
1 permeability, m2 gas density, kg m
3 velocity vector, m s
1 volume force vector, N m
3 mole fraction for species i, – constant for viscosity calculation, – thermal diffusion coefficient, kg m
1 s
1 effective binary diffusion coefficient, m2 s
1 mass source term for species j, kg m
3 s
1 diffusion coefficient, m2 s
1 reduced molar mass, kg mol
1 characteristic length, Å diffusion collision integral, – mole mass, kg mol
1 radius of pores, m reaction rate, mol m
3 s
1 mole concentration, mol m
3 specific heat at constant pressure, J kg
1 K
1 heat source term, W m
3 entropy change, J mol
1 K
1 enthalpy change of reaction, J mol
1 elasticity matrix, Pa

Greek symbols g polarizations, V b transfer coefficient, – l dynamic viscocity, Pa S

e  w

x

s s

k

e a

porosity, – strain, – viscous stress tensor, Pa mass fraction, – tortuosity for gas transfer, – tortuosity for ion/electron transfer, – thermal conductivity, W m
1 K
1 stress, Pa thermal expansion coefficient, K
1

Subscribes l ion s electron eff effective a anode c cathode act activation con concentration b boundary between gas channel and electrode g gas re reversible p polarizations R reactions M methane W water th thermal ref reference Abbreviations SOFC solid oxide fuel cell TPB triple phase boundary ASL anode support layer CSL cathode support layer AAL anode active layer ASA active specific area MSR methane steam reforming WGSR water gas shift reaction PEN positive electrolyte negative FI fuel inlet AI air inlet FEM finite element method PDE partial differential equation OCV open circuit voltage EAi interface of electrolyte-anode interface of electrolyte-cathode ECi

side can lower thermal stresses [24]. Although the proper contact area and depth are required to enhance the cell performance, the effects of the tip shape and depth on thermal stresses are not revealed, which will be significantly related to the durability. Therefore, we develop a 3D model of a unit SOFC using the COMSOL Multiphysics (Version 5.3) to study the thermal stresses resulting from different patterns of the interconnect-cathode contact, aiming to optimize the interconnect design in lights of shape, size and depth of the tips that are embedded in cathodes.

listed in Table 1. Five groups of partial differential equations (PDE) are selected as governing equations that describe the ion/electron, mass, momentum and heat transport, and thermal expansion, which include:

2. Mathematic model

H2 þ O2
! H2 O þ 2e


ð1Þ

Three typical contact modes as shown in Fig. 1 are discussed in this study, in which Dr, Dm and De mean the contact depth for the rectangular, trapezoidal and triangular tips. We set depths of the tips in cathodes to zero as a standard case, whose outline can be seen in Fig. 2 and the corresponding geometry parameters are

1 O2 þ 2e
! O2
2

ð2Þ

(1) Ion and electron transport Hydrogen as fuel in the anode and oxygen as oxidant in cathode are considered in the electrochemical reaction, which is the source of electricity (in Eq. (1)) and the rink of electron (in Eq. (2)) [25]:

Ohm’s law is used as the governing equation of ion/electron transport:

508

S. Zeng et al. / International Journal of Heat and Mass Transfer 125 (2018) 506–514

Fig. 1. The shapes of ribs in interconnects and the possible patterns of contacting with cathode.

Fig. 2. The current density at the outermost boundary for the anode (a) and cathode (b).

Table 1 The main mechanical properties and direction definition in the model. Component

X (cm)

Y (mm)

Z (mm)

Young’s modulus (GPa)

Poisson’s ratio (1)

TEC (10
6 K
1)

Interconnect (cathode side) Cathode support layer Cathode active layer Electrolyte Anode active layer Anode support layer Interconnect (anode side) Air channel Fuel channel


10 to 0


650 to 0 0–50 50–70 70–80 80–95 95–495 495–1145
500 to 0 495–995

0–1.5

205 42.7 60.8 205 81.1 84.2 205 –

0.28 0.28 0.3 0.3 0.3 0.3 0.28

12.3 12.4 11.4 10.3 11.4 12.5 12.3

0.5–1.5

il ¼
reff;l rUl

ð3Þ

is ¼
reff;s rUs

ð4Þ

where i is the current density, r the corresponding conductivity, U the potential. The subscribe l means ion, s the electron and eff the abbreviation of effective. As for the conductivity, it is determined by the material and the microstructure of porous electrodes, which will affect the ionic and electric transport path [26], so the effective conductivity of electrodes is corrected as [27]:

reff;a;s ¼ rNi 

VNi;a

reff;a;l ¼ rYSZ 

ss;a VYSZ;a

reff;c;s ¼ rLSM 

sl;a VLSM;c

ss;c

reff;c;l ¼ rYSZ 

VYSZ;c sl ; c

ð8Þ

here V means the volume fraction and s the corresponding tortuosity for electron/ion transfer, which can be found in Ref. [26]. r is the conductivity for specific material [28,29]:

rNi ¼


 9:5  107 1150 exp
T T


rYSZ ¼ 33:4  103 exp
ð5Þ

rLSM ¼ ð6Þ

 10; 300 T


 4:2  107 1200 exp
T T

ð9Þ ð10Þ

ð11Þ

The working voltage is determined by the open circuit voltage (OCV) and various polarizations, which is set as 0.7 V [30,31]:

ð7Þ

V ¼ EOCV
gact
gcon
gohm

ð12Þ

S. Zeng et al. / International Journal of Heat and Mass Transfer 125 (2018) 506–514

EOCV ¼ E0


pH2 O RT ln pffiffiffiffiffiffiffiffi 2F pH2 pO2

! ð13Þ ð14Þ

where V means the working voltage, EOCV the open circuit voltage,

gi specific polarization, E0 the reference voltage for H2-O2 reaction,

R the gas constant (8.314 J mol
1 K
1), F the Faraday constant (96,485 C mol
1), pi the partial pressure of specific species, T the temperature (K). The activation polarization is defined as [32]:

gact;a ¼ Us
Ul
Eeq;a

ð15Þ

gact;c ¼ Us
Ul
Eeq;c

ð16Þ

here Eeq;x is the equilibrium potential, which is set as 0 V and EOCV for anode and cathode, respectively. The concentration polarization is expressed as:

gcon;c ¼

pH2 ;b  pH2 O;TPB RT ln 2F pH2; TPB  pH2 O;b pO2 ;b RT ln 2F pO2 ;TPB

!

r

ð18Þ

ð19Þ

where i is the current density, L the thickness of corresponding component, and r the specific conductivity. Butler-Volmer equation is used to restrict the current density [33,34]:



 b  2F  gact
ð1
bÞ  2F  gact iv ¼ ASA  i0  exp
exp RT RT

ð20Þ here iv means current in the specific volume, ASA the abbreviation of active specific area, which is set as 4:875  104 m
1 for the anode side and 4:875  105 m
1 for the cathode side in this work, i0 the exchange current density, which is expressed in Eq. (21) [35], and b is the transfer coefficient, assumed to be 0.5.

i0 ¼ Ai 


 RT Ea exp 2F RT

ð21Þ

In Eq. (21), Ai is the pre-exponential factor, which is 2:35  1011 S m
2 for the cathode and 6:54  1011 S m
2 for the anode, Ea the activation energy for the reaction in Eqs. (1) and (2), which is 137 kJ mol
1 for oxygen and 140 kJ mol
1 for hydrogen. (2) Momentum, Mass and Heat transport Navier-Stokes equation is used to restrict the momentum transport in this model:



  ! 1 2l ! ! ! wþ ¼ F þ qg  r  u  u
r
p þ ðr  u Þ k e 3

l

ð22Þ

here l is dynamic viscosity of gas mixture, which will be presented in Eqs. (23) and (24), k the permeability of porous electrodes (1:76  10
11 m2 ), qg the average density of gas, which can be calcu! ! lated by ideal-gas equation, u the velocity vector, F the volume force vector, e the porosity (0.3) and w viscous stress tensor.

l¼

X

li  xi

i

ð24Þ


   rp rT !  ðr u Þ
r DTi ðxj
xj Þ  p T ! X !
r qg xj Deff;ij  rxj þ qg u  rxj ¼ Sj

ð25Þ

n

where xj is the mole fraction of species j, xj the mass fraction,

!

iL

n T 1000

where li is the viscosity of specific species, which depends on the local temperature (T), and parameters determined by the category of species (bn ) [36], xi the mole fraction of species i. It should be noted that the fully developed laminar flow is assumed in the momentum transport part, and the velocity at the fuel inlet is 0.323 m s
1, the AI velocity 2.87 m s
1. The pressures at both the fuel and air outlets are kept at 1 atm. Moreover, the forward flow and no backward flow are assumed. Counter-flow is kept for solution in all cases. Eq. (25) is given to describe mass transfer in gas channels and porous electrodes [26]:

ð17Þ

In Eqs. (17) and (18), pi is still the partial pressure and the subscribe b means the boundary between gas channel and ASL/CSL, TPB in porous electrodes. The ohmic polarization is evaluated from the ohm’s law:

gohm ¼


bn 

n¼1

E0 ¼ 1:253
2:4516  10
4 T

gcon;a ¼

7 X

li ¼

509

ð23Þ

DTi the thermal diffusion coefficient, which is not considered in our model, n the number of species categories, Deff;ij the effective binary diffusion coefficient, which will be discussed in following part, Sj the source term, which is coupled with the electrochemical reaction in the active layers and reforming reactions in the anode active layer (AAL) and anode support layer (ASL). The binary diffusion coefficient in gas channels is calculated as [26,37]:

Dij ¼

Mij ¼

2:66  10
8 T1:5 pffiffiffiffiffiffiffi 2 p Mij  lij  X

1 Mi

2 þ M1j

ð26Þ

ð27Þ

here M means molecular mass (kg mol
1) lij , X the characteristic length and diffusion collision integral [38], respectively. While Knudsen diffusion cannot be neglected in the porous electrodes as the average pore radius (0.34 lm) is comparable to the mean free path of gas molecules (0.2–0.5 lm) [39]:

Dij;k

sffiffiffiffiffiffiffiffiffiffi 2 8RT ¼ r 3 pMij

ð28Þ

It can be seen that in the Knudsen diffusion model, the binary diffusion coefficient is related to the radius of pores (r), and the reduced molecular mass (Mij ) and temperature (T). After taking into consideration the microstructure of electrodes, the effective binary diffusion coefficient in electrodes can be expressed as [39]:

Deff;ij ¼

Dij  Dij;k e  DijþDij;k s

ð29Þ

where e means the porosity of electrodes and s represents tortuosity for gas transport. The methane steam reforming (MSR) and water gas shift reaction (WGSR) in the ASL and AAL are considered:

MSR : CH4 þ H2 O () CO þ 3H2

ð30Þ

WGSR : CO þ H2 O () H2 þ CO2

ð31Þ

The reaction rate is computed as [40]:





27; 063 mCH4 mH2 O
3:7 T

232:78 mCO m3H2  10
14 T4 exp T

rMSR ¼ 63:6T2 exp

ð32Þ

510

r WGSR

S. Zeng et al. / International Journal of Heat and Mass Transfer 125 (2018) 506–514



12; 509 ¼ 1199T2 exp mCO mH2 O
6:77 T

16; 909  104 T2 exp mCO2 mH2 T

3. Results and discussion

ð33Þ

where mi is the mole concentration of species i. Note that at the fuel inlet we set the mole fraction of each species as: xH2 : xCH4 : xCO : xCO2 : xH2 O ¼ 0:2626 : 0:171 : 0:0294 : 0:0436 : 0:4934, which has been defined by the International Energy Agency [41]. It is defined as ambient air at the air inlet, which is assumed to only contain oxygen and nitrogen. It is set as the convective flux at both outlets [42]. The governing equation for heat transport is [43]:

qg  cp;g  ! u  rT ¼ rðkeff rTÞ þ Q h

ð34Þ

where qg is the average density and cp;g the average heat capacity, which have been defined elsewhere [36], keff the effective thermal conductivity. Note that the thermal conductivity for gas is calculated as an average value of each species, while for porous electrodes it is related to thermal conductivity of conducting solid (km ) and that of gases (kg ), and the porosity (e): keff ¼ km ð1
eÞ þ kg e. Q h represents the heat source term, which contains both the reversible heat source (Q re ) and irreversible one and the heat generated/consumed by reforming reactions:

Q re ¼


i  ASA  DS  T 2F

Q h;p ¼
ASA  i  ðgact þ gcon Þ þ

ð35Þ X i2

Q h;R ¼ r MSR DHM þ rWGSR DHW

r

ð36Þ ð37Þ

where DS is the entropy change for reactions in Eqs. (1) and (2), which can be got in a previous report [44], DH the enthalpy change (J mol
1) for reforming reactions (MSR and WGSR) [45]:

DHM ¼
ð206205:5 þ 19:5175TÞ

ð38Þ

DHW ¼ 45; 063
10:28T

ð39Þ

Note that in the heat transport part the temperature at both inlets is assumed to be 1000 K.

The current density at outermost boundaries, distribution of thermal stresses and the maximum thermal stresses are discussed to study the cell performance and its mechanical stability with changes of the interconnect-cathode contact modes. Note that the standard case is defined as no tip or rib insertion into the CSL, and the definition of three directions is summarized in Table 1. Besides, since the distribution of thermal stress for different cases is similar, we just provide thermal stresses distribution in the standard case. More attentions have been paid to the variation of thermal stresses, which could be beneficial to obtain an optimized contact mode for practical applications. Fig. 2 presents distribution of the current density at the outermost boundaries of interconnects in the standard case, in which it can be seen that the maximum current density reaches 4541.7 A m
2 at the anode side and 4482.1 A m
2 at the cathode side, since the reduction of oxygen is slightly more difficult than the oxidation of hydrogen (by comparing the exchange current density). Meanwhile, the current density at the lower part of boundary is higher than that at the upper part, which is mainly because the lower part directly contacts with the tips while the upper part contacts with channels. The maximum current densities for different interconnectcathode contact modes and depths are displayed in Fig. 3, where the cell connected by interconnects with triangular tips has the most excellent current density, reaching a maximum of around 5060 A m
2 when the depth of tip in cathode is 1 lm, 11.4% higher than that for the standard case. Note that in this work the voltage remains constant as 0.7 V for all cases, so a higher current density means a better performance. Moreover, the maximum current density for the triangle-insertion mode almost linearly increases with the increase of the contact depth, which rises to 5080 A m
2 when the contact depth is enhanced to 10 lm. For the rectangleinsertion mode, the maximum current density shows a convex changing trend with the depth, but for the trapezium-insertion mode, the maximum current density shows a concave changing tendency with the depth. The triangle-insertion mode shows the best performance because it can shorten the direct transfer path and increase the effective area for electron collection, which results in a minimum active specific resistance. While the reason for similar performance of rectangle and trapezium insertion is that when contact depth obviously varies, the difference between these two shapes can be

(3) Thermal stress couplings The thermal stress is related to material mechanical properties, the reference temperature for thermal strain (Tref , set as 1000 K in this model, because the sintering temperature for ceramics will surpass 1000 °C and the cell needs to be heated to 1000 K before it works), and the temperature distribution, which is expressed as:

th ¼ C  ðe
aðT
Tref ÞÞ

ð40Þ

where C is a matrix regarding some material properties [22,46], e the total strain, a the thermal expansion coefficient, which is summarized in Table 1. Note that we use the symmetry condition to simplify the calculation in this work, and the outermost boundaries of interconnect are fixed along the main flow direction and normal to the PEN (positive-electrolyte-negative) components direction. The fixing assumption is based on the fact that the cells are sealed between metal interconnects in practice, and the stiffness of interconnects is much higher than that of the porous electrodes (seen in Table 1), while other boundaries are set as free for the thermal expansion.

Fig. 3. Variation of the Maximum current density under different embedded shapes and depths.

S. Zeng et al. / International Journal of Heat and Mass Transfer 125 (2018) 506–514

Fig. 4. The 1st principle stresses at different locations including the EAi (a), ECi (b), AI (c) and FI (d).

Fig. 5. The maximum 1st principle stresses at the EAi (a), ECi (b), AI (c) and FI (d) at different modes and depths.

511

512

S. Zeng et al. / International Journal of Heat and Mass Transfer 125 (2018) 506–514

negligible, because the length of bottom edges of the trapezium is not changed, i.e. the similar shape produces comparable current density when the contact or insertion depth is large enough. However, the direction for electron transporting along the rectangular and trapezoidal ribs could be different, which can account for the convex and concave changing trends. Since thermal stresses can lead to crack in the PEN region and also shorten the lifetime of a fuel cell [47], more attentions should be paid to the electrolyte-electrode interface. Fig. 4(a and b) displays the 1st principle stress at the interface of electrolyte-anode (EAi, parallel to the XZ plane) and the interface of electrolytecathode (ECi), from which it can be seen that EAi is subjected to a larger tensile stress than that at the ECi. Note that a serious tensile stress locates at the top of interconnect rib (the maximum 1st principle stress becomes 21.9 MPa), because of the orthogonal angle of the shape of interconnects. Meanwhile the 1st principle stress at the lower part of EAi is approximately 4 MPa higher than that of ECi, and this means a relatively larger likeliness of thermal crack than the cathode side at the lower part. Distributions of the 1st principle stress (MPa) at AI and FI for the standard case are demonstrated in Fig. 4(c and d). Both inlets suffer less from the high 1st principle stress (the maximum value of 8.39 MPa belongs to AI, while 6.41 MPa locates at FI). In addition, at the contact part between CSL and interconnect of AI a relatively high stress is found (around 4 MPa), but at FI this stress is concentrated on the corner of air channel. The 1st principle stress at AI is mainly observed at the corner of fuel channel. Besides, the distribution of the 1st principle stress is homogeneous at the contact area between the ASL and interconnects of FI. At the inlet of gas channel, the 1st principle stress at the corresponding electrode-interconnect contact area shows a uniform distribution, but at the outlet it concentrates at the corner of gas channel, which is mainly because at the inlet the temperature is kept at 1000 K, and the thermal conductivity of interconnect is high, leading to a uniform distribution of the temperature. Moreover, the electrochemical reaction is not serious at inlets, which will not result in thermal gradients at these locations. However, the temperatures at outlets are raised by electrochemical reactions and the current mostly passes the corner (because of the shortest path and smallest resistance), which induces a higher temperature at the corner of outlet, corresponding to concentration of the 1st principle stress. Besides, the upper part of electrolyte and the adjacent AAL are subjected to the strongest tensile stress, possibly stemming from the mismatch of mechanical properties between the electrolyte and AAL. The outer boundaries of the lower part of interconnects withstand a negative 1st principle stress, which can be a compressive stress [46], because we keep these boundaries fixed in the X and Y direction, so they may be compressed by the adjacent unit cells (in a SOFC stack) during the thermal expansion. Fig. 5 presents the 1st principle stress at different contactmodes and insertion depths. The variation of the 1st principle stress at EAi is shown in Fig. 5(a), where the maximum 1st principle stresses increases with the insertion depth. When the insertion depth is 1 lm, it is the trapezium-insertion mode that possesses the largest 1st principle stress of nearly 31 MPa. As the depth increases, the 1st principle stress for the rectangle-insertion reaches the maximum value (while the depth is between 3 and 7 lm), but if the depth is larger than 7 lm, and the maximum value for the trapezium-insertion mode is around 31.3 MPa with the depth of 10 lm. As for the variation of the 1st principle stress at EAi, although this triangle-insertion mode yields the best electrochemical performance, it has no the maximum value, which may be because that we merely insert the interconnects into CSL, so the thermal stresses at EAi will depend most on the local temperature. Since the rectangle-insertion mode enhances the electron

transportation, which corresponds to a higher efficiency of electrochemical reaction at the anode side, leading to a high temperature but a small thermal gradient, so the 1st principle stress in the rectangle-insertion mode is not that high. Fig. 5(b) displays variation of the 1st principle stress at the ECi, showing that the absolute value of stress is lower than that at the EAi because of a smaller Young’s modulus and comparable CTE at the cathode side. Besides, the triangle-insertion mode brings about a more serious thermal stress at the ECi at any insertion depth. It can be speculated that the triangle shape can enlarge the area for electron delivery and improve efficiency of the involved electrochemical reaction (non-homogenously inside the cathode), resulting in a greater thermal gradient at the contact area. Note that the 1st principle stress at trapezium-insertion mode decreases when the insertion depth rises to 3 lm, this variation corresponds to the tendency of current density. Moreover, a decreasing current density contributes to a lower temperature and no other special factors lead to a great thermal gradient in this mode. Besides, the rectangle-insertion mode slightly increases the thermal stress at ECi, and the 1st principle stress shows no an obvious increase when the insertion depth is larger than 5 lm, accounting for the sluggish increase of current density at the rectangle-insertion mode. Fig. 5(c and d) shows the evolution of the 1st principle stress at the AI and FI, respectively. Note that we use double Y axis for the two figures because of the large differences in the three curves, and the values should be identified from the axis based on the same color. The rectangle-insertion mode induces the highest 1st principle stress of 60 MPa at AI, and the corresponding 1st principle stress decreases to 20 MPa as the insertion depth is increased to 10 lm. Moreover, the 1st principle stress for the trapeziuminsertion mode decreases with the augment of insertion depth, while the stress of triangle-insertion mode slowly increases, which reaches around 13 MPa at the deepest insertion. However, the value of the stress at FI is lower than that at AI, which is consistent with some previous report [46]. The 1st principle stress becomes saturated (around 10 MPa) for the rectangle-insertion mode when the depth increases to 10 lm. While the triangle-insertion mode possesses a higher stress than that for the trapezium-insertion, and the difference for the two modes becomes larger as the insertion depth is raised, which is not same to the variation at AI. On the other hand, dislocation at the electrolyte-electrode interface due to mismatch of thermal expansion will take place, which will worsen the durability of a cell. The dislocation can be analyzed based on the maximum shear stress. The maximum shear stress is defined as the half of difference between the 1st principle stress and the 3rd principle stress, which is proportional to the thermal expansion. The thermal expansion parallel to the electrolyte-electrode interface is expressed by the black arrows as shown in Fig. 6. Note that the length and thickness are proportional to the level of thermal expansion, and the arrow direction is the same to thermal expansion. It can be found that the maximum shear stress of 34.4 MPa locates at the lower part of EAi for the standard case, where the most distinctive thermal expansion occurs due to the maximum current density causing the highest temperature. However, the maximum shear stress is lower at the ECi (D ¼ 2:3 MPa for the maximum value). Besides, there exists a shear stress of higher than 20 MPa at the upper part of EAi, which is 2–3 MPa higher than that at the ECi, indicating a larger probability of distortion in the electrolyte. Fig. 7 shows the maximum shear stress for different contact modes between the interconnect and cathode. In Fig. 7(a), it can be seen that the variation of the maximum shear stress at the EAi is similar to that in Fig. 5(a) because the maximum shear stress is determined by the 1st principle stress to some extent and another thing can be presumed that the insertion at the cathode

S. Zeng et al. / International Journal of Heat and Mass Transfer 125 (2018) 506–514

513

Fig. 6. The maximum shear stresses (MPa) at the EAi (a) and ECi (b) for the standard case.

Fig. 7. The maximum shear stresses for different modes and depths.

side induces no more changes to the 3rd principle stress at the anode side. The rectangle-insertion mode has the lowest value of 48.8 MPa at EAi when the insertion depth is below 1 lm. When the insertion depth reaches 10 lm, the trapezium-insertion mode has a dominant stress and the triangle-insertion gives the smallest stress. However, the maximum shear stress at the ECi is lower for all cases as shown in Fig. 7(b), and the trapezium-insertion mode produces the smallest stress (around 45.8 MPa) with the insertion depth of 10 lm, while the triangle and trapezium-insertion mode show the same magnitudes of shear stresses at ECi when the insertion depths are smaller than 1 lm. Besides, the effect of insertion on the maximum shear stress at ECi for the triangle-insertion mode cannot be neglected, which leads to the highest shear stress when depth varies between 0.5 and 4 lm, and the varying values as the depth exceeds 5 lm. As for the assembly of SOFC stacks, normally coarse LSM powers will be pasted onto the cathode surface to enhance the conductivity, and also acts as a buffer layer to relieve the mechanical contacting risk as we previously report [15], where the multiple layers including current collecting layer, active layer and sealing layer are assembled in series to obtain a short stack. The proposed contacting modes are thus expected to improve the stability of a short SOFC stack by reducing various thermal stresses.

1. The triangle-insertion mode is the most effective design to improve the electrochemical performance of SOFCs, which reduces the active specific resistance at the interface of interconnect-cathode, resulting in relatively higher temperature gradients and thermal stresses at ECi, but as for EAi, it causes no large thermal stresses because of the smaller thermal gradient. 2. The rectangle-insertion mode is beneficial to improvement of the performance, but the effectiveness is not that remarkable as the ribs are inserted into the CSL but it will introduce the largest thermal stresses at the gas inlets compared with the other insertion modes. 3. The trapezium-insertion mode has a complex effect on the cell performance, which firstly decreases and then increases to the same magnitude with that is resulted from the rectangleinsertion mode. The difference between the trapezium- and rectangle-insertion modes can be neglected when the insertion depth is distinctive, so similar results are obtained for the rectangle- and trapezium-insertion when the insertion depth increases large enough. Conflict of interest The authors declare that there is no conflict of interest.

4. Conclusion Acknowledgement A 3D model is developed to reveal the distribution of current density and thermal stresses by varying the insertion modes and depths, and conclusions can be drawn as follows:

This work is supported by National Natural Science Foundation of China (No. 51702039).

514

S. Zeng et al. / International Journal of Heat and Mass Transfer 125 (2018) 506–514

References [1] A. Atkinson, S. Barnett, R.J. Gorte, J.T.S. Irvine, A.J. McEvoy, M. Mogensen, S.C. Singhal, J. Vohs, Advanced anodes for high-temperature fuel cells, Nat. Mater. 3 (2004) 17–27. [2] S. Park, J.M. Vohs, R.J. Gorte, Direct oxidation of hydrocarbons in a solid-oxide fuel cell, Nature 404 (2000) 265–267. [3] B.C.H. Steele, Appraisal of Ce1
yGdyO2
y/2 electrolytes for IT-SOFC operation at 500 °C, Solid State Ion. 129 (1–4) (2000) 95–110. [4] L. Yang, S. Wang, K. Blinn, M. Liu, Z. Liu, Z. Cheng, M. Liu, Enhanced sulfur and coking tolerance of a mixed ion conductor for SOFCs: BaZr0.1Ce0.7Y0.2–xYbxO3–d, Science 326 (5949) (2009) 126–129. [5] J. Rossmeisl, W.G. Bessler, Trends in catalytic activity for SOFC anode materials, Solid State Ion. 178 (31–32) (2008) 1694–1700. [6] T. Takeguchi, R. Kikuchi, T. Yano, K. Eguchi KazutoshiMurata, Effect of precious metal addition to Ni-YSZ cermet on reforming of CH4 and electrochemical activity as SOFC anode, Catal. Today 84 (3–4) (2003) 217–222. [7] Y.-L. Lee, J. Kleis, J. Rossmeisl, Y. Shao-Horn, D. Morgan, Prediction of solid oxide fuel cell cathode activity with first-principles descriptors, Energy Environ. Sci. 4 (2011) 3966–3970. [8] M. Andersson, J. Yuan, B. Sundén, SOFC modeling considering hydrogen and carbon monoxide as electrochemical reactants, J. Power Sources 232 (15) (2013) 42–54. [9] Y. Jiang, A.V. Virkar, Fuel composition and diluent effect on gas transport and performance of anode-supported SOFCs, J. Electrochem. Soc. 150 (7) (2003) A942–A951. [10] Z. Yang, G.G. Xia, X.H. Li, J.W. Stevenson, (Mn, Co)3O4 spinel coatings on ferritic stainless steels for SOFC interconnect applications, Int. J. Hydrogen Energy 32 (16) (2007) 3648–3654. [11] Z. Yang, G. Xia, J.W. Stevenson, Mn1.5Co1.5 O 4 spinel protection layers on ferritic stainless steels for SOFC interconnect applications, Electrochem. SolidState Lett. 8 (3) (2005) A168–A170. [12] M. Liu, Z. Lü, B. Wei, X. Huang, Y. Zhang, W. Su, Effects of the single chamber SOFC stack configuration on the performance of the single cells, Solid State Ion. 181 (19–20) (2010) 939–942. [13] L. Jin, W. Guan, X. Ma, H. Zhai, W.G. Wang, Quantitative contribution of resistance sources of components to stack performance for planar solid oxide fuel cells, J. Power Sources 253 (2014) 305–314. [14] W.B. Guan, L. Jin, X. Ma, W.G. Wang, Investigation of impactors on cell degradation inside planar SOFC stacks, Fuel Cells 12 (6) (2012) 1085–1094. [15] W.B. Guan, H.J. Zhai, L. Jin, T.S. Li, W.G. Wang, Effect of contact between electrode and interconnect on performance of SOFC stacks, Fuel Cells 11 (3) (2011) 445–450. [16] L. Jin, W. Guan, J. Niu, X. Ma, W.G. Wang, Effect of contact area and depth between cell cathode and interconnect on stack performance for planar solid oxide fuel cells, J. Power Sources 240 (2013) 796–805. [17] S. Liu, C. Song, Z. Lin, The effects of the interconnect rib contact resistance on the performance of planar solid oxide fuel cell stack and the rib design optimization, J. Power Sources 183 (2008) 214–225. [18] L.K. Chiang, H.C. Liu, Y.H. Shiu, C.H. Lee, R.Y. Lee, Thermal stress and thermoelectrochemical analysis of a planar anode-supported solid oxide fuel cell: effects of anode porosity, J. Power Sources 195 (7) (2010) 1895–1904. [19] L.K. Chiang, H.C. Liu, Y.H. Shiu, C.H. Lee, R.-Y. Lee, Thermo-electrochemical and thermal stress analysis for an anode-supported SOFC cell, Renew. Energy 33 (12) (2008) 2580–2588. [20] A. Selimovic, M. Kemm, T. Torisson, M. Assadi, Steady state and transient thermal stress analysis in planar solid oxide fuel cells, J. Power Sources 145 (2) (2005) 463–469. [21] C.K. Lin, T.T. Chen, Y.P. Chyou, L.K. Chiang, Thermal stress analysis of a planar SOFC stack, J. Power Sources 164 (1) (2007) 238–251. [22] S. Zeng, M. Xu, J. Parbey, G. Yu, M. Andersson, Q. Li, B. Li, T. Li, Thermal stress analysis of a planar anode-supported solid oxide fuel cell: Effects of anode porosity, Int. J. Hydrogen Energy 42 (31) (2017) 20239–20248. [23] K.S. Weil, B.J. Koeppel, Thermal stress analysis of the planar SOFC bonded compliant seal design, Int. J. Hydrogen Energy 33 (14) (2008) 3976–3990. [24] M. Xu, T. Li, M. Yang, M. Andersson, Solid oxide fuel cell interconnect design optimization considering the thermal stresses, Sci. Bull. 61 (17) (2016) 1333– 1344.

[25] M. Ni, 2D thermal-fluid modeling and parametric analysis of a planar solid oxide fuel cell, Energy Convers. Manage. 51 (4) (2010) 714–721. [26] M. Andersson, J. Yuan, B. Sundén, SOFC modeling considering electrochemical reactions at the active three phase boundaries, Int. J. Heat Mass Transf. 55 (2012) 773–788. [27] D. Kanno, N. Shikazono, N. Takagi, K. Matsuzaki, Nobuhide Kasagi, Evaluation of SOFC anode polarization simulation using three-dimensional microstructures reconstructed by FIB tomography, Electrochim. Acta 56 (11) (2011) 4015–4021. [28] S.A. Hajimolana, M.A. Hussain, W.M.A.W. Daud, M. Soroush, A. Shamiri, Mathematical modeling of solid oxide fuel cells: a review, Renew. Sustain. Energy Rev. 15 (4) (2011) 1893–1917. [29] J.R. Ferguson, J.M. Fiard, R. Herbin, Three-dimensional numerical simulation for various geometries of solid oxide fuel cells, J. Power Sources 58 (2) (1996) 109–122. [30] Y. Patcharavorachot, A. Arpornwichanop, A. Chuachuensuk, Electrochemical study of a planar solid oxide fuel cell: role of support structures, J. Power Sources 177 (2) (2008) 254–261. [31] M. Ni, Modeling and parametric simulations of solid oxide fuel cells with methane carbon dioxide reforming, Energy Convers. Manage. 70 (2013) 116– 129. [32] O. Razbani, M. Assadi, M. Andersson, Three dimensional CFD modeling and experimental validation of an electrolyte supported solid oxide fuel cell fed with methane-free biogas, Int. J. Hydrogen Energy 38 (2013) 10068–10080. [33] S.H. Chan, K.A. Khor, Z.T. Xia, A complete polarization model of a solid oxide fuel cell and its sensitivity to the change of cell component thickness, J. Power Sources 93 (1–2) (2001) 130–140. [34] M. Ni, M.K.H. Leung, D.Y.C. Leung, Parametric study of solid oxide fuel cell performance, Energy Convers. Manage. 48 (5) (2007) 1525–1535. [35] D.Y. Zhonggang Zhang, Guogang Yang, Jingfeng Chen, Yifeng Zheng, He Miao, Weiguo Wang, Jinliang Yuan, Naibao Huang, Three-dimensional CFD modeling of transport phenomena in multi-channel anode-supported planar SOFCs, Int. J. Heat Mass Transf. 84 (2015) 942–954. [36] B. Todd, J.B. Young, Thermodynamic and transport properties of gases for use in solid oxide fuel cell modelling, J. Power Sources 110 (1) (2002) 186–200. [37] H.K. Sanghyeok Lee, Kyung Joong Yoon, Ji-Won Son, Jong-Ho Lee, Byung-Kook Kim, Wonjoon Choi, Jongsup Hong, The effect of fuel utilization on heat and mass transfer within solid oxide fuel cells examined by three-dimensional numerical simulations, Int. J. Heat Mass Transf. 97 (2016) 77–93. [38] R.C. Reid, J.M. Prausnitz, B.E. Poling, The Properties of Gases and Liquids, fourth ed., R.R. Donnelley & Sons Company, USA, 1986. [39] J. Yuan, Y. Huang, B. Sundén, W.G. Wang, Analysis of parameter effects on chemical coupled transport phenomena in SOFC anodes, Heat Mass Transf. 45 (4) (2009) 471–484. [40] M. Ni, M.K.H. Leung, D.Y.C. Leung, Micro-scale modelling of solid oxide fuel cells with micro-structurally graded electrodes, J. Power Sources 168 (2) (2007) 369–378. [41] M. Andersson, Doctoral Thesis, Lund University, Sweden, 2011. [42] M. Andersson, J. Yuan, B. Sunden, Grading the amount of electrochemical active sites along the main flow direction of an SOFC, J. Electrochem. Soc. 160 (1) (2013) F1–F12. [43] J.Y. Martin Andersson, Bengt Sundén, SOFC modeling considering electrochemical reactions at the active three phase boundaries, Int. J. Heat Mass Transf. 55 (2012) 773–788. [44] W.G. Bessler, D. Jürgen Warnatz, G. Goodwin, The influence of equilibrium potential on the hydrogen oxidation kinetics of SOFC anodes, Solid State Ion. 177 (39–40) (2007) 3371–3383. [45] M. Ni, Modeling of SOFC running on partially pre-reformed gas mixture, Int. J. Hydrogen Energy 37 (2) (2012) 1731–1745. [46] M. Xu, T. Li, M. Yang, M. Andersson, I. Fransson, T. Larsson, B. Sundén, Modeling of an anode supported solid oxide fuel cell focusing on thermal stresses, Int. J. Hydrogen Energy 41 (33) (2016) 14927–14940. [47] L. Liu, G.-Y. Kim, A. Chandra, Modeling of thermal stresses and lifetime prediction of planar solid oxide fuel cell under thermal cycling conditions, J. Power Sources 195 (8) (2010) 2310–2318.