Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas

Energy xxx (xxxx) xxx

Contents lists available at ScienceDirect

Energy journal homepage: www.elsevier.com/locate/energy

Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas Chao Yang a, b, Xiuhui Jing b, He Miao a, Yu Wu c, Chen Shu b, Jiatang Wang a, Houcheng Zhang d, Guojun Yu b, Jinliang Yuan a, * a

Faculty of Maritime and Transportation, Ningbo University, Zhejiang, China Merchant Marine College, Shanghai Maritime University, Shanghai, China College of Engineering Science and Technology, Shanghai Ocean University, Shanghai, China d Department of Microelectronic Science and Engineering, Ningbo University, Zhejiang, China b c

a r t i c l e i n f o

a b s t r a c t

Article history: Received 24 April 2019 Received in revised form 10 October 2019 Accepted 15 October 2019 Available online xxx

A two-dimensional mathematical model is developed for a single-cell based on the planar configuration and validated by relevant experimental data, with an aim to describe the coupling phenomena of the multiphysics transport processes and the meso-scale elementary reactions. It is revealed that desorption and adsorption reactions in the electrode mostly take place near the electrolyte and the channel, respectively; the distribution of the surface species depends on the gas diffusion in the porous electrode affected by the thickness and microstructure of the electrode. The electrochemical reactions are centralized in about 100 mm thick electrode from the electrolyte. Nis and COs are the major surface species in both fuel cell (FC) and electrolysis cell (EC) modes. Os is higher in the FC mode, particularly near the electrolyte due to the desorption and charge transfer reactions; The microscopic structure properties, including average porosity, tortuosity and particle size, are also influential on the elementary reactions due to the gas diffusion through the tortuous pathways and the active sites on the catalyst surfaces. It is also found that the performance predicted in the global models is often overestimated, because the limitations of the local elementary reactions are not considered in the global model. © 2019 Elsevier Ltd. All rights reserved.

Keywords: rSOFC Elementary reactions Dual-mode operation Surface species Transport phenomena

1. Introduction Reversible solid oxide fuel cell (rSOFC) is able to work with dualfunctions of electricity generation and energy storage as shown in Fig. 1(a). It can produce electricity as solid oxide fuel cell (SOFC) with flexible fuels such as H2, syngas or other hydrocarbon gases. Furthermore, it can also electrolyze H2O or CO2/H2O mixture by the external electricity as solid oxide electrolysis cell (SOEC), to produce H2 or CO/H2 (syngas), respectively [1e3], as expressed in Eqs. (1)e(3). The syngas can be converted into other fuels with a further approach such as Fischer-Tropch processes [1,4]. rSOFC is appropriated to balance the conflicts between energy supply and demand in both traditional and new energy systems, with advantages of compact, high power density and simple design [5,6]. The reactions coupled transport processes in the dual-functional electrodes are crucial for the reversible performance [7,8]. It is more complicated

when the syngas is employed involving methane reforming, watergas shift (WGS) in Eq. (4) and electrochemical reactions simultaneously [9]. Electrochemical reactions in fuel electrode: FC mode


! H2 þ O2


























H2 O þ 2e EC mode

(1)

FC mode


! CO þ O2


























CO2 þ 2e EC mode

(2)

Electrochemical reaction in air electrode: FC mode


! O2 þ 4e


























2O2 EC mode

(3)

Water-gas shift (WGS) reaction in fuel electrode: FC mode

* Coreponding author. E-mail address: [email protected] (J. Yuan).


! CO þ H2 O


























CO2 þ H2 EC mode

(4)

https://doi.org/10.1016/j.energy.2019.116379 0360-5442/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

2

C. Yang et al. / Energy xxx (xxxx) xxx

Fig. 1. (a) The schematic diagrams of rSOFC in dual modes, (b) the meso-scale elementary reactions in Ni/YSZ electrode.

The experimental approaches are possible but limited to provide the coupling details of the complex reactions and transport processes, consequently the mathematical modeling approaches have been developed to understand the fundamental but important processes in rSOFC [8e15]. Notably, the global kinetics were simply considered as the gas reaction mechanism in the work above. However, the reactions occur mostly in the porous electrode, where the adsorption/desorption and other surface reactions take place heterogeneously rather than uniformly [15]. For example, on the Ni/ YSZ electrode, the charge transfer on the triple-phase boundary (TPB) and the thermo-catalytic reforming on the surface of catalyst Ni take place as shown in Fig. 1(b). Therefore, the meso-scale elementary reaction mechanism should be applied for understanding the microscopic phenomena and effects in the electrodes of rSOFC. The mechanism of the elementary reforming and electrochemical reactions were investigated on the surface of Ni, which is commonly applied for the fuel electrodes. (1) Reforming reactions: A meso-scale elementary reaction mechanism was developed by Hecht et al. for CH4/CO2/H2O in 800  C, in SOFC [16]. Janardhanan et al. adjusted the parameters for the above mechanism model for the thermodynamic consistency in 500e2000  C [17]. Deutschmann et al. developed a multi-step steam reforming mechanism of CH4 in the reforming reactor with Ni, which can be thermodynamically consistent in 546e1546  C [18]. Delgado et al. extended the mechanism model, the effect of H2 and H2O on the dry reforming with CO2 was also evaluated [19]. The fuel reforming reactions were also investigated for the catalysts Rh, Pt and Pd [20e22]. (2) Electrochemical reactions: The elementary mechanisms of the charge transfer reactions were also investigated on the TPB such as Ni/YSZ. Hydrogen and oxygen spillover were commonly recognized as the two paths of the charge transfer reactions [23e25]. Hs was the major adsorbates on the TPB, which was controlled by the rate-limiting steps involving the mid-product OHs [23]. The adsorption/desorption of H2 and the diffusion of O2 on YSZ were considered as the rate-limiting steps in the H2/H2O electrochemistry [26]. CO/CO2 elementary electrochemistry was investigated in Ni/YSZ or Pt/YSZ in two paths: The oxygen spillover

occurs with the charge transfer on the TPB [27], or CO is participated in the charge transfer as adsorbates directly [28,29]. The elementary reaction mechanism of hydrocarbon fuel and carbon were also investigated for different electrode and electrolyte [30e33]. The mechanisms mentioned above were commonly proved and applied in many mathematical modeling work for SOFC [34e37], SOEC [38,39], and rSOFC [40e42]. However, there are still deficiencies according to the open published literature: (1) The meso-scale elementary reactions are reversible theoretically in SOFC and SOEC operations, but still different in terms of the detailed pathways of the surface reactions. It was reported that Hs and COs were the major intermediate surface species for the co-electrolysis reactions in SOEC [10], while COs and Os were considered as the major ones in SOFC [30]. There was almost no study specific for the reversible elementary reactions between SOFC and SOEC, which is necessary for rSOFC. (2) The Butler-Volmer equation based on the global kinetics with H2/H2O/ O2 was commonly applied to describe the electrochemical reactions, which is not appropriated for syngas as fuel. Fueyo et al. also pointed that the exchange current density and charge transfer coefficient were not identical in the reversible operations [43]. (3) The energy transport were neglected in most work [30,41], which may not be appropriated under the high temperature and the mode-switching operating conditions in rSOFC. (4) Some work presented above were based on the button cell along the thickness direction. The gas channel, the convection flux, the pressure gradient in the porous electrode were neglected [10,41]. It is not appropriated for the planar or tubular configurations in the real applications. In this study, a two-dimensional single-channel model is developed for a planar rSOFC single cell based on the cell sample prepared specialized in the experiment. A 20-steps meso-scale reversible reforming and charge exchange reactions mechanism are considered with syngas (H2/H2O/CO/CO2) for both FC and EC modes, the multiphysics transport coupled with the elementary reactions are also considered. The thickness of the electrode and microstructure properties are investigated to understand the interactions between the multiphysics transport and elementary

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

C. Yang et al. / Energy xxx (xxxx) xxx

3

Fig. 2. Geometry of (a) planar rSOFC unit cell, and (b) a 2D single channel model along the main flow direction.

reactions. The predictions by the global reaction model are also presented to compare with those by the elementary reactions model. 2. Mathematical modeling description 2.1. Geometry model and model assumptions An electrolyte-supported planar cell with 4 cm  4 cm active area is prepared for the experiment and modeling. A twodimensional geometrical model (yellow shadow face in Fig. 2(a)) is built in the mid-plane along the main flow direction. Five layers

are considered in Fig. 2(b): the self-supported YSZ electrolyte layer (ELE) with 200 mm thickness, one Ni/YSZ fuel diffusion layer (FDL) and Ni/GDC fuel active layer (FAL) on the fuel side, one LSM air diffusion layer (ADL) and one LSM/GDC air layer (AAL) on the oxygen side. Each of the all four electrode layers is 30 mm thick. 500 mm thick gas channels are built on both sides. Notably, the rib of the interconnector is not considered due to the absence of the cross-sectional plane in this two-dimensional model. The following assumptions are applied: (1) the ideal gases with laminar flow on both electrodes; (2) the local temperature equilibrium conditions inside the cell and heat insulation on all the external boundaries; (3) the steady-state conditions, homogeneous

Table 1 The elementary reforming reaction mechanism applied for syngas on Ni-based catalyst in rSOFC [10,17]. Reaction Adsorption reactions

Desorption reactions

Surface reactions

a b c

1f. H2þNis þ Nis/Hs þ Hs 2f. O2þNis þ Nis/Os þ Os 3f. H2O þ Nis/H2Os 4f. CO2þNis/CO2s 5f. CO þ Nis/COs 1b. Hs þ Hs/H2þNis þ Nis 2b. Os þ Os/O2þNis þ Nis 3b. H2Os/H2O þ Nis 4b. CO2s/CO2þNis 5b. COs/CO þ Nis 6f. Os þ Hs/OHs þ Nis 6b. OHs þ Nis/Os þ Hs 7f. OHs þ Hs/H2Os þ Nis 7b. H2Os þ Nis/OHs þ Hs 8f. OHs þ OHs/Os þ H2Os 8b. Os þ H2Os/OHs þ OHs 9f. Os þ COs/CO2s þ Nis 9b. CO2s þ Nis/Os þ COs 10f. Os þ Cs/COs þ Nis 10b. COs þ Nis/Os þ Cs

Aa b

1.000E-02 1.000E-02b 1.000E-01b 1.000E-05b 5.000E-01b 5.593Eþ19 2.508Eþ23 4.579Eþ15 9.334Eþ07 4.041Eþ11 5.000Eþ22 2.005Eþ21 3.000Eþ20 2.175Eþ21 3.000Eþ21 5.432Eþ23 2.000Eþ19 3.214Eþ23 5.200Eþ23 1.354Eþ22

n

Ea

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 3

0 0 0 0 0 88.12 470.39 62.68 28.8 112.85e50XCOsc 97.9 37.19 42.7 91.36 100 209.37 123.6e50XCOsc 86.5 148.1 115.97e50XCOsc

Arrhenius parameters for the reaction constant rate given by: k¼ATnexp(-E/(RT)). Sticking coefficient. The active energy dependent on the surface coverage X.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

4

C. Yang et al. / Energy xxx (xxxx) xxx

microstructures of electrodes; (4) the elementary reactions occur on the surface of Ni in the porous electrode, while the electrochemical reactions occur only on TPB; (5) the surface species uniformly distributed on the surface of Ni. 2.2. Elementary reforming reaction mechanism and model For the catalyst Ni and mixture H2/H2O/CO/CO2 in 800e850  C specific for rSOFC in this model, a 20-steps elementary kinetics mechanism, validated in similar work for both SOFC and SOEC [10,17,30,41,44], is employed. In which, there are 5 gas species and 8 surface species, and the gas adsorption/desorption, surface elementary reaction and thermo-catalytic chemistry are concerned, as shown in Table 1. The rate of the elementary reaction Rj is given as [16]:

Rj ¼

Kr  X

yij  y’ij



j¼1

kj

KgþKs Y

Ci

y’ij

sffiffiffiffiffiffiffiffiffiffiffiffi RT 2 p Mi Ni

 kj ¼ Aj T bj



Ej RT

FC





! Hs þ OH



H2 OO þ e þ Nis O

k¼1

  εkj Xk m X k kj exp  RT

(7)

where Aj, bj, and Ej are Arrhenius parameters listed in Table 1, mkj and εkj are the correction coefficients to describe the species coverage dependency of the rate constants for the surface species k in the reaction j. Xk is the surface coverage of the surface species k. Moreover, the reaction rate of the gas and surface species are written as [17]:

Rj

(8)

j¼1

2.3. Electrochemical reaction model

FC

Oxo

EC

Vo”

’’  OxO þ Nis



!



O s þ e þ VO EC

(9)

where and are the oxygen interstitial and vacancy, respectively, in the electrolyte-phase. As an assumption, only chargetransfer steps take place on the TPB [10].

(12)

FC

OxO þ Nis



!



Os þ 2e þ V ’’O EC

(13)

In Eq. (13), Nis and Os are the active surface site and the surfacephase on the catalyst Ni, respectively, which also participate into the meso-scale reforming reactions in Table 1. Then the reaction rate Re is given as:  Re ¼ kþ e COxO CNis  ke COs CV ’’

(14)

O

where COxO ¼ 4.45  104 mol/m3 and CV } ¼ 4.65  103 mol/m3 are O the concentrations of the oxygen interstitial and vacancy in YSZ [10], which are several orders of magnitude higher than the concentrations of Nis and Os (<1 mol/m3) [16]; therefore, it is appropriated to keep COxO and CV } as the constants for the simplification of O the calculation [10]. Furthermore, kþ e and ke are Arrhenius reaction rate constants given as [36,39]:

kþ e

¼ kþ e;0

k e

¼ k e;0

exp

afuel F hfuel act b

! (15)

RT

0 exp@ 

afuel F hfuel act f RT

1 A

(16)

where af and ab are the charge-transfer coefficients, hact is the activation overpotential, kþ e,0 and ke,0 are the pre-exponential factors given as [36,39]:

kþ e;0 ¼

Different elementary reaction mechanisms were reported in the literature [23,25,32]. For Ni-YSZ electrode with H2/H2O, it is given as [25]:





!



H2 O þ 2e þ V ’’O

(11)

FC

(6)

Y Ks

(10)

(5)

where S is the initial sticking coefficient; GNi ¼ 2.6  109 mol/cm2 is the total site density of Ni [16,17]; y is the sum of the surface reactants stoichiometric coefficients on catalyst Ni; Mi is the molecular weight of the gas species. For the desorption and the surface reactions [17]:

H2 þ OxO

EC

i¼1

S kj ¼ y G

Kr X

FC

 Hs þ OxO



!



OH O þ e þ Nis

EC

where Kr and j are the total number and index of the reactions, respectively; i is the index of species including gas-phase (Kg) and surface-phase (Ks); yij and y’ij are the stoichiometric coefficients of the reactants and products, respectively; Ci is the molar concentration for the gas-phase (mol/cm3) and surface-phase (mol/cm2); kj is Arrhenius reaction rate with the following equations: For the adsorption reactions [17]:

Ri ¼

Furthermore, a two-step charge-transfer reactions (Eqs. (10) and (11)) involving two spillover reactions of hydrogen from Ni surface to oxide ions and hydroxyl ions on the surfaces were agreed well with the experiment data with H2/H2O [25]. In the contrast, the charge-transfer process with Os spillover (e.g. Eq. (12) and (13)) was dominant than Hs spillover in the case with syngas, which may be attributed to the co-redox reactions with high carbon-based fuel [45]. Considering the syngas with the co-redox reactions referred in this model, the one-step O spillover mechanism in Eq. (13) is applied for the charge-transfer reaction in this work [46].

fuel

i0

nF,C 0Ox ,C 0Nis

(17)

O

k e;0 ¼

ifuel 0 nF,C 0Os ,C 0V }

(18)

O

where C 0Ox , C 0Nis , C 0Os and C 0V } are the molar concentrations of the O O surface-phase species in the equilibrium conditions, F is Faraday constant, n is the number of charges transferred, ifuel is the equi0 librium exchange current density [36,39]:

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

C. Yang et al. / Energy xxx (xxxx) xxx

fuel i0

fuel

¼ gfuel 

X aH2



X bH2 O

E  exp  act RT

! (19)

where gfuel is the pre-exponential coefficient; a and b indicate the concentration dependencies of the reactants and products, and Eact is the activation energy. On the other hand, only O2 participates into the electrochemical reaction in the oxygen electrode. The global kinetics electrochemistry is applied in Eq. (3) for simplification [41]. Then the ButlerVolmer equation is appropriated for the air side electrochemical kinetics [41]:

¼ iair 0

I air act

! !# " aair F hair aair F hair XO2 act f act b  exp   exp XO2 ;ref RT RT

(20)

where XO2 and XO2,ref are the local molar fractions of oxygen. The exchange current density iair 0 on the air electrode is given as [7,13]:

iair 0

¼ gair 

X cO2

Eair  exp  act RT

According to the Faraday’s Law, the electrochemical reaction rate Re on the air electrode is given as:

Re ¼

2.4.1. Mass transport properties and source terms The effective gas diffusion coefficient is given as [13]:

ε Di;m  Di;k gas Di;eff ¼  t Di;m þ Di;k

Iair act nF

(22)

P Di;m ¼

¼ Fe  Fion

isi’

 Xi’ D i;i’

(27)

1  Xi

where Xi is the molar fraction of the gas specie i, Di,i’ is known as the binary diffusion coefficient describing the gas diffusion between gas i and i’, as given by [13]:

Di;i’ ¼

0:0026  T 1:5 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 P  M1i þ M1i’  di þ2di’  UD

UD ¼

1:06036

0:1561 þ

ti;i’ T

þ

0:193 1:03587

þ

exp 0:47635ti;i’ T exp 1:52996ti;i’ T

1:76474

exp 3:89411ti;i’ T (29)

(23) k

hair act

¼ Fe  Fion  OCV

(24)

where Fe and Fion are the potential of electron and ion, respectively; OCV is the open circuit voltage on the cell according to the Nernst equation [7]:

0 1 air RT @ P O2 A ln OCV ¼ fuel nF P

(25)

O2

fuel where Pair O2 and PO2 are the partial pressure of oxygen in the air and fuel side, respectively, under the equilibrium conditions. Notably, Pfuel O2 is determined by the elementary reaction of oxygen desorption occurring in the fuel electrode in Table 1. Other electrochemical parameters are listed in Table 2.

b ti;i’ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi li  li’

(30)

where li is the characteristic Lennard-Jones energy of species i, kb is the Boltzmann constant. Further details are shown in Table 4. Di,k represents the Kundsen diffusion for the gas diffusion in the porous electrode, as determined by [13]:

dp Di;k ¼  3

sffiffiffiffiffiffiffiffiffi 8RT pMi

(31)

where dp is the pore diameter of the porous electrode. The source term Sm and Sg,i represent the gas production and consumption caused by both the meso-scale elementary reforming and electrochemical reactions, as given by:

Sm ¼

2.4. Transport processes model

Kg X

Sg;i

(32)

i¼1

The Multiphysics transport processes model is listed in Table 3.

Table 2 Electrochemical reaction parameters applied in the model [7,13]. Electrochemical parameters

Fuel electrode

Air electrode

Forward charge transfer coefficient, af Backward charge transfer coefficient, ab Pre-exponential coefficient, ga (A/m2) Concentration dependency coefficient Activation energy, Eact (J/mol)

0.1 1

0.1 0.4

gfuel ¼ 8  108

gair ¼ 7  108

a

(28)

where P is the gas pressure, d is the characteristic parameter of the binary diffusion, UD is the dimensionless diffusion collision given as [13]:

In addition, the local activation overpotential hact is expressed for the fuel and air electrode [7]:

hfuel act

(26)

where ε and t are the porosity and tortuosity, respectively. Di,m is the molecular diffusion of gases in the mixture, given as [13]:

! (21)

5

a ¼ 0.11, b ¼ 0.67, c ¼ 0.25 Efuel Eair act ¼ 120000 act ¼ 120000



Sg;i ¼ ± Ri , AV;Ni ± Re , AV;TPB  Mi

(33)

where þ and - represent the production and consumption of the mass species, respectively. AV,Ni and AV,TPB are the specific active areas per volume for Ni and TPB, respectively [7,13]:

AV;TPB ¼ p sin2 q  r 2e  ntot  ne  nion 

Ze  Zion  pe  pion Z (34)

Fitted by the experimental data.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

6

C. Yang et al. / Energy xxx (xxxx) xxx

Table 3 The governing equations applied in the model [7,13]. Mass and gas species transport Continuity equation

.

V,ðrf V Þ ¼ Sm .

Gas species equation

Vðrf Yi Þ þ Sg;i V,ðrf V Yi Þ ¼ V,½Dgas i;eff

Momentum transport Momentum equation

.

V,ðrf V V Þ ¼ V,ðmeff VV Þ þ Sd;V 

Energy transport Energy equation

.

vP vn

"

V,ðrV cp;f TÞ ¼ V, keff Vðrcp;f TÞ 

Electron and ion transport Electron equation

 V,ðse;eff VFe Þ ¼ Se

Ion equation

 V,ðsion;eff VFion Þ ¼ Sion ¼  Se

Surface species transport Diffusion equation

 V,ðDsur i;eff VCi Þ ¼ Si

Table 4 Parameters applied in evaluation of diffusion coefficient [7,13].

# P hi Ji þ ST i

Table 5 Structure parameters of electrode [13].

Parameters

H2

H2O

CO

CO2

O2

N2

Structure parameters

FDL

di li/kb

2.827 106.7

2.614 809.1

3.69 91.7

3.941 195.2

3.467 106.7

3.798 71.4

Thickness of layers (mm) Porosity, ε Tortuosity, t Pore diameter of electrode, dp (mm) Contact angle between ion and electron particles, q Effective radius of particles, r (mm) Volume fraction of electronic or ionic particles, f

30 30 30 0.34 0.22 0.3 3 3 3 4.8 0.7 0.3  q¼15 re ¼ 0.5, rion ¼ 0.5 fe ¼ 0.5  (1-ε), fion ¼ 0.5  (1-ε) 6

AV;Ni ¼ pr 2e ntot ne

sin2 q  nion  Ze  Zion sin2 q  ne  Ze  Ze  4 Z Z

!

Average coordination number, Z

FAL

AAL

ADL 30 0.3 3 0.4

(35) where q is the contact angle between the electronic and ionic conductor particles; re and rion are the effective radius of conductor particles; ntot, ne and nion are the numbers of the total, electronic and ionic particles per unit volume, respectively [7]:

1ε

3 n þn 3 p r e ion  g e 3

ntot ¼ 4

(36)

g ¼ rion=r

(37)

" pe=ion ¼ 1 

Ze=ione=ion ¼

 #0:4  4:236  Ze=ione=ion 2:5 2:472 ne=ion  Z   ne=ion þ 1  ne=ion  g2

(38)

(39)

The structure parameters for the electrodes are listed in Table 5.

ne ¼

e

f e  g3 1  fe þ fe  g3

(38)

nion ¼ 1  ne

(39)

where fe is the volume fraction of the electronic particles. Z in Eqs. (34) and (35) is the average coordinate number for a random packing of spheres (e.g., particle Ni) [13]. Ze and Zion are the coordinate numbers of the electronic and ionic particles, respectively [7]:

Ze ¼ 3 þ

 Sd;Vn ¼ 

b¼

C¼

Z3

(36)

2

ne þ nion  ðrion =re Þ

Zion ¼ 3 þ

2.4.2. Momentum transport properties and source terms The source terms Sd,Vx and Sd,Vy in the momentum equations demonstrate the momentum losses caused by the viscosity and inertial resistances in the porous electrodes, given by the Brinkman-Forchheimer-Darcy Model [7,13]:

Z3 ne þ nion  ðrion =re Þ2

. mVn þ εrCjV jVn b

ð2rÞ2  ε3 72t  ð1  εÞ2 ε

b

0:5

1:8 

0:5 180ε5

 (40)

(41)

(42)

where b is the permeability given by the Kozeny-Carman relationship, C is the inertial resistance factor.

 ðrion =re Þ2

(37)

pe and pion in Eqs. (34) and (35) are the probabilities of the two kinds of conductor particles, which are given as [7,13]:

2.4.3. Energy transport properties and source terms The thermal diffusion coefficient in the electrode is contributed by the combined heat conduction in the gas and solid phase as:

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

C. Yang et al. / Energy xxx (xxxx) xxx Table 6 Electric conductivity of different materials [7,13]. Materials

Electronic conductivity (S/m)

Ni YSZ GDC LSM

3.27  106e1065.3T

2.45  105-97.5T

Ionic conductivity (S/m) 3.34  104exp(-10300/T) 3.5  103exp(-6471/T) 2  1013exp(0.0245T)

keff ¼ εkf þ ð1  εÞks

(43)

The source term ST in Table 3 represents the heat production and consumption caused by the reactions and the overpotential losses:

ST ¼ Qsur þ Qpol

Qpol ¼ Re  h þ

7

jIe j2

se;eff

Qsur ¼ ±Re ,AV;TPB  DHe ±

Kr X

Rj ,AV;Ni  DHj

(45)

j¼1

where DHe and DHj are the reaction enthalpy for the electrochemical and reforming reactions, respectively [16e18]. For the heat production caused by the overpotential losses, the active and ohmic polarizations are considered with the source term Qpol:

jIion j2

(46)

sion;eff

where Ie and Iion are the current densities of electron and ion, respectively; se,eff and sion,eff are given in 2.4.4. Furthermore, the P heat production caused by the gas diffusion is given by hiJi in Table 3, in which Ji is the mass flux of the specie i.

2.4.4. Charge transport properties and source terms The diffusion of the charges are driven by the difference of the potentials between the electrons and ions, i.e., Ue and Uion; se,eff and sion,eff are the effective conductivity of the electronic and ionic conductors [13]:

(44)

For the heat production caused by the reactions, the source term Qsur is given as:

þ

se=ion;eff ¼ se=ion 

fe=ion  0:294 ε  0:294 1 þ 1þε

!2 (47)

where f is the volume fraction of the electronic or ionic particles; se and sion are the conductivity in Table 6. The source term Se and Sion represent the reaction rates of the electron and ion on the TPB, respectively:

Se ¼ AV;TPB  Re

(48)

Sion ¼  Se

(49)

Table 7 Boundary and interface conditions. Boundaries

Energy(K)

Species i

Electron(V)

Ion (V)

Inlet in fuel side Interface between Interface between Interface between Interface between Interface between Inlet in air side

1103 Continuity Continuity Continuity Continuity Continuity 1103

Xi Continuity Continuity Insulation Continuity Continuity Xi

None 0 Continuity Insulation Continuity Uoperating None

None Insulation Continuity Continuity Continuity Insulation None

fuel channel and FDL FDL and FAL FAL/AAL and electrolyte AAL and ADL ADL and air channel

Fig. 3. The illustration of rSOFC experimental testing system.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

8

C. Yang et al. / Energy xxx (xxxx) xxx

2.4.5. Surface-phase species diffusion properties and source terms The diffusion of the surface-phase species on the particles is also modeled. However, the diffusion of surface species takes place in the microscale level on the Ni surface, and the diffusion coefficient of the most surface species ranges between 108 and 1010 m2/s, others are even smaller [25]. These values are several orders of magnitude smaller than that of the gas diffusion, which could be neglected in the cell level model in this work [33,36]. The production and consumption of the surface specie is also given as the surface reaction rate:

Si ¼ AV;Ni  Ri

(50)

where Ri is given by Eq. (8). 2.4.6. Boundary conditions The boundary conditions based on the experiment are listed in Table 7. The rib and interconnector are not considered in this 2-D mid-plane model as shown in Fig. 2. The electronic conductivity in the interconnector is much larger than that in the diffusion layer, then the interface between the gas channel and electrode is treated as the boundary of electron potential. Furthermore, the utilization rate of gaseous reactants can be obtained from the species at the inlet and outlet of channels as below [6,7]:



 Ui ¼ Xi;in  Xi;out Xi;in

(51)

where Xi,in and Xi,out are the molar fraction of the gas specie i at inlet and outlet, respectively. 3. Experimental testing procedure and modeling validation 3.1. Testing procedure A 4  4 cm2 electrolyte-supported planar cells are prepared for the testing and modeling validation, as discussed in our previous work [47,48]. The single-cell testing system is shown in Fig. 3. The testing procedures are illustrated as: (1) Cell installation: The single planar cell is fixed in the electrical resistance furnace, Ni and Ag meshes are provided as the current collectors on the fuel and air electrode, respectively; (2) Seal preparation: The cell is heated up to 1123 K with 1 K/min heating rate and sealed with molten glass; (3) Reduction for fuel electrode: H2 with 3% H2O is provided for 3 h, then the cell is heated up to the required temperature (1103 K for FC mode); the current density is kept as a constant at 0.63 A/cm2 for 24 h to ensure the active and stable performance for all region of the fuel electrode [47]. (4) Running test: Fuel and air are provided with the constant flow rates (200 sccm for H2/H2O, 300 sccm for syngas, 500 sccm for air). After the testing, N2 is sent to the cell for cooling. (5) Data recording: The real-time data of the current and voltage are recorded by the electrochemical workstation (CHI 760e, CH Instrument Inc. USA) with four wire approach on the testing. Moreover, the fuel supply and loads are changed when the working condition is switched between FC and EC modes. 3.2. Modeling validation The mathematical model is implemented into the commercial CFD software Ansys/Fluent 16.0 coupled with User-Defined functions (UDFs). The meso-scale reforming/electrochemical reactions, and the multiphysics transport property/source terms are defined in the UDFs. The non-uniformed structural meshes are employed and the grid independence analysis is also performed. According to Fig. 4(a), the voltage and the current density are obtained for different composition of H2/H2O gases in the IeV

Fig. 4. Comparison of IeV curves from experiment and simulation approaches fueled by: (a) H2/H2O, (b) syngas.

curves. It is found that different voltages are observed under the high current density conditions between the experiment and simulation work. Because the 2-D model in Fig. 2, which is not able to reflect the gas diffusion and concentration differences along the cross-sectional direction, leading to an underestimation of the concentration overpotential losses at the high current density [37]. Considering the length of channel (4 cm) and the ranges of the current density (small than ±0.3 A/cm2), the differences of the concentration overpotential are small enough in Fig. 4(a). In overall, the simulation results agree well with that obtained in the experiments in both FC and EC modes. Similar phenomena are obtained when the syngas is applied as the fuel in Fig. 4(b). In addition, the global and elementary reforming reaction models are both implemented and compared with syngas. The global reaction model includes the overall electrochemical and WGS reactions as referred in Eqs. (1)e(4), as presented in our previous work [7,13]. In Fig. 4(b), it is found that the experimental data agrees better with the simulation results obtained by the elementary reaction model (solid lines) than that predicted by the global reaction model (dot lines), which indicates that the elementary reforming and one-step charge-transfer reaction

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

C. Yang et al. / Energy xxx (xxxx) xxx

9

mechanism are more suitable in this study. 4. Results and discussion 4.1. Distribution of surface-phase species in the dual-mode operation

Fig. 5. The average coverage of surface-phase species in dual-modes.

In this work, the distribution of the surface species including Nis, COs, Hs, Os, H2Os, OHs, CO2s, Cs are investigated to evaluate the effects of the meso-scale elementary reactions. According to the assumption, the elementary reactions take place on the surface of Ni with a constant total sites density (2.6  109 mol/cm2), the diffusion of surface species and the gaseous reactions are both neglected. Correspondingly, the coverage of the surface specie Xi indicates that how many Ni sites have been occupied by these specific species during the elementary reactions. According to the average Xi in

Fig. 6. The distribution of temperature in the (a) FC and (b) EC modes.

Fig. 7. The distribution of the surface species coverage in the midline (x ¼ 0.01 m) of the fuel electrode along the thickness direction in dual-mode: (a) XNis and XCOs, (b) XHs and XOs, (c) XOHs and XH2Os, (d) XCO2s and XCs.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

10

C. Yang et al. / Energy xxx (xxxx) xxx

Fig. 5, Nis, COs, Hs and Os are the major intermediate species which occupies 99% of the total available sites on the surface of catalyst Ni in dual-mode. It is also found that XNis is much higher than that of others, and XNis in the FC mode is higher than that in the EC mode. COs is the second major adsorbate; XCOs is about 7%e8% higher in the EC mode than in the FC mode due to the oxidized co-electrolysis reactions, i.e., shifting CO2 to CO. Consequently, XOs is higher in the FC mode than in the EC mode (1.6% vs. 0.2%). XH2Os, XOHs, XCO2ss, XCs are several orders of magnitude smaller than others. It indicates that H2Os, OHs, CO2s and Cs are favored to shift to COs, Hs and Os in dual-mode, which agrees with the experiment and modeling work in references [10,36]. Furthermore, XH2Os, XOHs, and XCO2s are larger in the FC mode than those in the EC mode, respectively, which is attributed to the difference of XOs appeared in the FC and EC modes. XCs is smaller in the FC mode, which may be caused by both higher XOs (shift Cs to COs) and higher temperature. In Fig. 6, the average temperature in the FC mode is above 4.5 K higher than in the EC mode, which is attributed to the differences between the exothermic FC mode and endothermic EC mode. In Fig. 7, the distribution of the surface species at the middle-line (x ¼ 0.01 m) is presented along the thickness direction. FFL and FDL layers with both 30 mm thickness are considered. In Fig. 7(a), XNis decreases from the electrolyte side (the left-hand side) to the gas channel side (the right-hand side), which indicates that more Ni sites are occupied by the surface species near the gas channel, and

more Ni sites are released near the electrolyte. Due to this fact, the absorption reactions are stronger near the gas channel, in contrast the desorption reactions with the gas production may take place mostly near the electrolyte. XCOs increases from the electrolyte side to the fuel side, because COs is produced in Eqs. (5f) and (9b) which are also Nis reduction reactions. Therefore, as the major adsorbate, COs increases with decreasing Nis along the thickness of the electrode. Moreover, the elementary reactions of COs reduction are endothermic according to Table 1, which are more favored to the high temperature near the electrolyte according to Fig. 6. In Fig. 7(b), XHs decreases from the electrolyte side to the fuel channel side in dual-mode. According to Table 1, Hs is produced from the adsorption of H2 (Eq. (1f)), or from the resolving of H2Os and OHs (Eqs. (6b) and (7b)), due to the fact that more H2 are produced in the meso-scale reforming reaction near the electrolyte, while the resolving of H2Os and OHs also shifts to Hs here mostly, which agrees with the finding in Fig. 5. XOs also decreases obviously from the electrolyte side along the thickness direction in the FC mode, because O2 on the TPB releases Os on the Ni surface in the charge transfer reaction in Eq. (13), leading to a larger XOs near the electrolyte in the FC mode. In contrast, XOs is smaller near the electrolyte side in the EC mode, because Os is continuously consumed in the backward charge transfer reaction in Eq. (13). In Fig. 7(c), XH2Os and XOHs decrease from the inside to outside along the thickness of electrode in the FC mode, as opposed to the finding

Fig. 8. The distribution of (a) XNis in the FC mode, (b) XNis in the EC mode, (c) XCO in the FC mode, (b) XCO2 in the FC mode.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

C. Yang et al. / Energy xxx (xxxx) xxx

in the EC mode, which is attributed to the difference between the reductive FC mode and oxidized EC mode as mentioned in Fig. 5. It is also found that XCs increases from the inside to the outside of the electrode in Fig. 7(d), which indicates that the carbon deposition occurs near the gas channel mostly due to the lower temperature in Fig. 6. The distribution of the surface specie is also presented in 2-D view. The distribution of XNis is compared in Fig. 8(a) and (b) for the FC and EC modes, respectively. The gradient of XNis is observed in both modes along the thickness direction. For the FC mode, the smallest XNis is located in the inlet region near the gas channel, because more Ni sites are occupied duo to the gas adsorptions on the top side of the fuel electrode; while the largest XNis is located at the electrolyte side, which indicates that the reforming and desorption reactions are stronger at the bottom of the electrode, as mentioned in Fig. 7. Similar phenomena are observed in the EC mode, but higher gradient are achieved along the main flow direction. According to Table 1, the distribution of Nis is affected by the gas adsorption and desorption reactions, and more Ni sites are occupied by the adsorption at the top side of the electrode due to the easier diffusion of fuel near the gas channel, while more Ni sites are released due to the stronger desorption reactions at the bottom side of the electrode. In other words, the elementary reactions are interdependent with the diffusion of gases. More details are provided for the distribution of gases. The distributions of the molar fraction CO (XCO) and CO2 (XCO2) in the FC mode are shown in Fig. 8(c) and (d), respectively. As CO being the

11

major gaseous reactant in the FC mode, the gradient of XCO is observed along the thickness direction in Fig. 8(c), and the range of XCO is between 22.6% and 25%. The distribution of XCO in the diffusion and functional layers is opposed to that of XNis in Fig. 8(a), which may be attributed to the CO adsorption (i.e., consuming more CO and occupying more Ni free sites near the gas channel). On the other hand, the distribution of XCO2 is similar with that of XNis in the diffusion and functional layers in the FC mode, because CO2 is produced in the desorption reaction with Ni sites released. Therefore, the distribution of Nis (also other surface species) are affected by the diffusion and distribution of gaseous species involving the adsorption/desorption reactions.

4.2. The effect of the electrode thickness The effect of the thickness of the electrode is investigated in the dual-mode operation. Notably, the density of the total sites in the Ni surface are assumed as 2.6  109 mol/cm2 in all the cases, which indicates that more available active sites are provided in the thick electrode cases. In addition, the global reaction model with the classical gaseous reaction mechanism is also considered to compare the effects of the electrode thickness. In Fig. 9, XNis, XCOs, XHs, XOs are presented on the middle-line (x ¼ 0.01 m) along the thickness direction for different thick FDL (30 mm, 200 mm and 600 mm). In Fig. 9(a), XNis decreases with increasing the thickness of FDL in the FC mode, and XNis in the 600 mm thick FDL is 12% lower than that in the 30 mm thick FDL. In

Fig. 9. The distribution of surface species in the middle-line along the thickness direction of fuel electrode in dual-mode, for: (a) XNis, (b) XCOs, (c) XH2Os, (d) XOs.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

12

C. Yang et al. / Energy xxx (xxxx) xxx

other words, the more Ni sites are occupied in the thick electrode, when the same fuel gas (300 sccm) is provided. On one hand, the elementary reactions include the gas adsorption (in Eqs. (1f)-(5f)) and the surface-phase reforming reaction (in Eqs. (6b), (7b) and (9b), 10b), which are favored with increasing the active surface region in the thick electrode; on the other hand, the gas desorption, which occurs mostly near the electrolyte side, are not favored in the thick electrode, because the diffusion of the gas products is more difficult duo to the longer pathway. Therefore, more Ni sites are occupied by the elementary reactions (e.g. gas adsorption), as opposed to that, the less Ni sites are released because the gas desorption are suppressed in the thick fuel electrode. Consequently, XNis are lower in the 600 mm FDL case than others. Similarly, XNis decreases when the thickness of FDL increases in the EC mode. For example, XNis decreases from 80% to 58% when FDL increases from 30 mm to 200 mm. Furthermore, it needs at least 100 mm thick electrode for XNis to stabilize, i.e., to balance the meso-scale elementary reactions. As opposed to XNis, XCOs increases with increasing thickness of FDL in Fig. 9(b). For example, the largest XCOs is obtained in 600 mm FDL in both FC and EC modes. XCOs in the case of 600 mm thickness FDL is 14% higher than that in the 30 mm thick FDL in the FC mode. In contrast, XHs decreases with increasing the thickness of FDL as found in Fig. 9(c), which is opposite to that of XCOs. It may be attributed to the differences on both reaction and diffusion ability between CO/CO2 and H2/H2O. For example, COs comes from the adsorption (from CO) and reforming reactions (in Eqs. (5f), (9r)), then desorbs as CO2 (in Eqs. (9b), (4r)) in the FC mode. The adsorption of CO is improved with the thickness of electrode because more active sites are provided in thick FDL case. However, CO2 desorption, which occurs at the bottom of the electrode mostly, are not favored for the thick electrode, because the bigger diffusion resistance makes it harder to transport the product CO2 from inside to outside. Then the desorption of CO2 is suppressed leading to a higher XCOs in the thick electrode case. Due to the fact that the suppressing of CO2 desorption dominates the changes of XCOs here. In contrast, Hs is produced in Eqs. (1f), (6b) and (7b), which occurs in the H2 adsorptionand the H2Os resolving reactions, respectively. Indeed, the diffusion of H2/H2O is much better than that of CO/CO2 as mentioned in Section 4.2.1, which is not suppressed obviously by increasing the thickness of the electrode. Due to the fact in Fig. 9(c), the increasing available sites of Ni in thick electrode case is dominant, to decrease the coverage of Hs (XHs) totally. In Fig. 9(d), XOs also decreases with increasing the thickness of FDL in the FC mode. There is a slop of decrease of XOs within about 100 mm thickness from the electrolyte side in the FC mode. Because Os is the product in the forward charge transfer reaction according to Eq. (13), which is then reformed to COs, H2Os and CO2s in the elementary reactions mostly. However, there are not enough active surface sites for the reforming reactions in 30 mm FDL case, leading to more Os left in the FC mode. In contrast, there are enough active sites provided in 200 mm and 600 mm FDL cases for the reforming surface reactions, leading to a higher consumption of Os (e.g. shifting COs and Hs to CO2s and H2Os) and lower XOs. Similarly, XOs also decreases with increasing the thickness of FDL in the EC mode. Furthermore, XOs stabilizes in 100 mm thickness of the fuel electrode for the most cases, which indicates that the charge transferring reactions and other Os elementary reactions are centralized within 100 mm thick fuel electrode near the electrolyte. In respect to the charge transfer reactions in different thick fuel electrodes, the distribution of the active current density Ia is evaluated and presented. According to Eqs. (13) and (14), R¼Ia/(nF) represents the reaction rate of charge transfer reaction between the ionic and electronic particles, thus Ia can be used to evaluate the charge transfer and electrochemical performance in the active

region. In Fig. 10(a), Ia decreases down to 0 A/m3 even in 100e150 mm thickness near the electrolyte for the most cases, which indicates that the charge transfer reactions are centralized in this region. It agrees with the distribution of Nis and Os which are major reactants on the surface of Ni in Eq. (13). Similar phenomena are observed in the EC mode shown in Fig. 10(b). It is also found that Ia in the 200 mm thick FDL case is the largest in the FC mode and the smallest (negative value) in the EC mode according to the enlarged shots near the electrolyte, which indicates that better performance on the charge transfer reactions are achieved in 200 mm thick FDL case. The global reaction models are also employed based on Eqs. (1)e(4) to investigate the effects of the electrode thickness. For example, Ia in different thick fuel electrode for both the global and elementary reaction model are presented for the FC mode in Fig. 10(a). It is found that the largest Ia increases by 9.1%, 12.1% and 14.7% in the 30 mm, 200 mm and 600 mm FDL cases when the global reaction model is applied. In the global model, the electrochemical reactions are assumed to take place between the gaseous reactants in the electrode according to Eqs. (1) and (2); In the elementary reaction model, the charge transfer reactions are assumed to occur on the TPB according to Eq. (13), which is affected by the diffusion

Fig. 10. The distribution of the active current density Ia along the midline (x ¼ 0.01 m) of the fuel electrode in different thickness of the fuel electrodes for the elementary and global models (a) in the FC mode, (b) in the EC mode.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

C. Yang et al. / Energy xxx (xxxx) xxx

and reactions on both gaseous and surface species in the fuel electrode. As mentioned in Fig. 4(b), the IeV curves with the worse electrochemical performance in the elementary reaction model agrees better with the experimental data than that by the global reaction model, which indicates that the electrochemical reaction performance is overestimated by the global reaction model. Similarly, the higher Ia (negative) is observed by the global reaction model than by the elementary model in the EC mode in Fig. 10(b). In Fig. 11, the effect of the electrode thickness is also presented for the gas utilization rate Ui¼(Xi,in-Xi,out)/Xi,in Eq. (51) in both global and elementary reaction models. For the elementary reaction model in Fig. 11(a), UCO increases with the thickness of FDL in the FC mode. It is clear that CO adsorption and COs reforming (to CO2) reactions take place mostly in the FC mode, which are also favored to the larger active area in the thick FDL case. On the other hand, UH2 decreases with increasing the thickness of FDL in the FC mode. The decrease of UH2 is unlikely caused by suppressing the adsorption of H2, because the available sites of Ni for the surface reactions increases with the thickness of FDL, which may be attributed that the reforming reactions to shift H2O to H2 on the surface of Ni, which are enhanced when the active area is increased. Similar phenomena are observed on UCO2 and UH2O in the EC mode in Fig. 11(b). Furthermore, the global reforming and electrochemical reaction models are also implemented to evaluate and compare their differences. In the FC mode, higher UCO and lower UH2 are observed by

13

the global reaction model than in the elementary reaction model. The changes of UCO and UH2 are also more obvious with the thick FDL predicted by the global reaction model. UCO increases by 12.2%, 10.5%, and 7.8% in the global reaction model than in the elementary reaction model. Similar phenomena are observed in the EC mode. It is revealed that the gaseous reactants are consumed more in the global reaction model than in the meso-scale elementary reaction model. As mentioned in Fig. 10(b), only gaseous diffusion and reactions are considered for the whole co-redox reactions in the global reaction model. In other words, the performance of the coredox reactions in the global reaction model is overestimated than that in the elementary reaction model. 4.3. The effect of the fuel electrode microstructure properties The microstructure properties of the fuel electrode including the porosity ε, tortuosity t and effective radius of particles re are studied, respectively, with an aim to evaluate the interdependency between the microstructure properties and the performance of the meso-scale elementary reactions. The global and elementary reaction models are both implemented to investigate the differences of these models. In Fig. 12(a), Ia increases in both FC and EC modes when ε increases from 0.21 to 0.35. According to Eq. (26), the diffusion of the gaseous species are improved with increasing the porosity, leading to a lower resistance and better gas diffusion performance. Therefore, higher Ia is obtained when the

Fig. 11. The utilization of fuel gases in different thickness of FDL (a) in the FC mode, (b) in the EC modes.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

Fig. 12. The distribution of the active current density Ia along the middle-line (x ¼ 0.01 m) of the fuel electrode in different (a) porosity ε, (b) tortuosity t, (c) effective radius of particles re in dual-mode.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

C. Yang et al. / Energy xxx (xxxx) xxx

electrochemical reactions are improved with the better gas diffusion. Furthermore, the electrochemical reaction performance is also compared between the elementary and global reaction models. It is found that the largest Ia increases from 1.2  109 A/m3 to 1.58  109 A/m3 by the elementary reaction model, and from 1.38  109 A/m3 to 1.96  109 A/m3 by the global reaction model, when ε increases from 0.21 to 0.35 in the FC mode. Consequently, it is revealed that the electrochemical reaction performance is better by the global reaction model, which is similar with the effects on the electrode thickness in Fig. 10. Similar phenomena are observed when t decreases from 5 to 3 in Fig. 12(b). Because the change of the tortuosity only affects the gas diffusivity and permeability according to Eqs. (26) and (41), then the diffusivity of the gaseous species are increased with decreasing the tortuosity, which is favored to the improvement on the diffusion of gases through the porous electrode. Consequently, the concentration losses decreases, which is beneficial for both gaseous elementary reactions and charge transfer reactions. In addition, a better performance with higher Ia is still obtained in the global reaction cases in both FC and EC modes in Fig. 12(b). In Fig. 12(c), Ia increases when the effective radius of the particles (including electronic and ionic conductor particles) decreases from 0.5 mm to 0.3 mm in both FC and EC modes. It is found that the highest Ia increases from 1.2  109 A/m3 to 2.9  109 A/m3 in the FC

15

mode, and from 4.1  109 A/m3 to 7.2  109 A/m3 in the EC mode, which is much higher than the increase of Ia in Fig. 12(a) and (b). Because the specific area AV,Ni and AV,TPB increase dramatically with the radius of the particles according to Eqs. (34) and (35), while the permeability also increases according to Eq. (41), leading to a much improvement on both reaction and diffusion processes. Furthermore, larger Ia is obtained in the elementary model than in the global model when re increases, e.g., Ia,elementary ¼ 2.9  109 A/m3 vs. Ia,global ¼ 2.55  109 A/m3 in the FC model, while Ia,ele9 3 9 3 mentary ¼ 7.2  10 A/m vs. Ia,elementary ¼ 6.83  10 A/m in the EC mode. It is revealed that the effects on the improvement of the electrochemical reactions are more obvious by the elementary model. Because the whole elementary reactions such as adsorption/desorption, surface reforming and charge transfer reactions are highly correlated to the specific active surface AV,Ni and AV,TPB according to Eqs. (33) and (50), which is not considered in the global reaction model (only gaseous reactions considered). Therefore, the effects of the meso-scale reactions with increasing the particle radius are omitted in the global reaction model, which is not reasonable as discussed above. The effects of the microstructure properties are also investigated for the gas utilization rate Ui. In Fig. 13(a), UCO increases from 0.204 to 0.212 when ε increases from 0.21 to 0.35, which is caused by the improvement of the gas diffusion with higher ε, as discussed in

Fig. 13. The utilization of fuel reactants in different microstructure properties of the fuel electrode (a) in the FC mode, (b) in the EC mode.

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

16

C. Yang et al. / Energy xxx (xxxx) xxx

Fig. 14. The coverage of surface species Nis and COs varies with (a) the porosity ε, (b) the tortuosity t, (c) the particle radius of the fuel electrode re in dual-mode.

Fig. 12(a). Similar phenomena are observed for UCO when t reduces from 5 to 3, due to the better gas diffusion in the lower tortuous pathway case. While UCO increases from 0.204 to 0.255 when re decreases from 0.5 to 0.3 mm, which is more obvious than that by the porosity and tortuosity variation. As discussed in Fig. 11, CO is adsorbed mostly shifting to COs which is the major adsorbate in the FC mode, the active area increases dramatically with decreasing re according to Eqs. (34) and (35), which improves the adsorption of CO mostly. In contrast to that, UH2 increases with higher porosity and lower tortuosity due to the improvement of the gas diffusion. However, UH2 decreases from 0.0695 to 0.0537 when re decreases from 0.5 to 0.3 mm, as opposed to that of UCO. Because H2 is the one of the major gaseous reactants, which is not only provided by the inlet fuel, but also provided by the reforming reactions from H2O according to the elementary reactions mechanism in Table 1 (Eqs. (1b), (3f) and (7b), etc.). Consequently, more H2 are provided in smaller particle radius case with larger specific active surface (AV), leading to a decrease of UH2. According to Fig. 13(b), similar phenomena on the utilization of the gaseous reactants CO2 and H2O are observed in the EC mode. Furthermore, the gas utilization by the global model is also evaluated. In the FC mode, it is found that both UCO and UH2 increase in higher porosity or lower tortuosity cases, because of the improvement of the gas diffusion, as mentioned in the elementary reaction model. However, UCO increases from 0.219 to 0.246 when re decreases from 0.5 to 0.3 mm, which increases by 14.6%; in the contrast, UCO increases by 25% in the elementary reaction model, which is much higher than that in the global model. On the other hand, UCO decreases by 11.9% and 22.8% in the global and elementary models, respectively. It is revealed that the effect of decreasing re is more apparent in the elementary reaction model. As discussed in Fig. 11, CO and H2 are considered to participate into the

elementary reactions in the elementary model, and all the elementary reactions involving the gaseous and surface-phase species occur on the active surfaces of catalyst, which is directly correlated to the decrease of re according to Eqs.(34) and (35). Similar phenomena are observed in the EC mode in Fig. 13(b). Therefore, the meso-scale elementary reactions have been improved in the small particles case with higher active area. In contrast, only gaseous species are considered in the global reaction model. The interdependency of the surface coverage and the microstructure properties is evaluated and presented also for the dualmode operation. In Fig. 14(a), the XNi in the fuel electrode decreases when ε increases from 0.15 to 0.35, then increases when ε > 0.35. The performance of the gases diffusion is improved with higher ε, while the active area decreases simultaneously, which is not favored to the surface reactions. It is clear that the improvement on the gases diffusion is dominant when ε varies between 0.15 and 0.35, to improve the surface reactions inside the electrode with more Ni sites occupied. On the other hand, when ε > 0.35, the active sites are not enough for the surface reactions like adsorption, which is dominant to suppress the surface reactions (e.g., gas adsorption), to decrease the Ni sites occupied by other adsorbates and to lead to an increase of XNis. In other words, there is an optimized porosity (about 0.35 in this work) based on the trade-off between the gas diffusion and active area. In the contrast, as the major intermediate surface species, XCOs changes conversely to that of XNis. XCOs increases from 0.262 to 0.273 in the FC mode, and from 0.401 to 0.413 in the EC mode, respectively, when ε increases from 0.15 to 0.35. Then XCOs decreases to 0.259 in the FC mode, and to 0.398 in the EC mode, respectively, when ε increases from 0.35 to 0.6, because the elementary reactions are suppressed due to the decrease of the active area as mentioned above. In Fig. 14(b), XNis in the fuel

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

C. Yang et al. / Energy xxx (xxxx) xxx

electrode increases with t from 2 to 7. Because the diffusivity of gases and permeability of the electrode increase with decreasing the tortuosity according to Eqs. (26) and (41), then the gas reactants and products are easier to diffuse through the tortuous pathways into or out from the electrode, which leads to a reduction of the concentration overpotential as mentioned in Fig. 12(b). The surface reactions are also improved with more Ni active sites occupied, results in a smaller XNis and larger XCOs with a smaller tortuosity. Due to the same fact, XNis increases with re, as opposed to the decrease of XCOs in Fig. 14(c). 5. Conclusion In this work, a two dimensional single channel model of rSOFC was developed. The syngas with the meso-scale elementary reactions were considered and coupled with the multiphysics transport processes. The detailed information of the surface species and elementary reactions were evaluated and presented. Parameters including the thickness and the micro-structure properties of the fuel electrode were investigated, with an aim to provide inspiration of the optimization on the design and operating in rSOFC. Some conclusions are summarized as below: (1) The distribution of surface phases: Nis, COs, Hs and Os counts about 99% of the total available sites on the surface of Ni, and other adsorbates are favored to shift to the former surface species. Nis decreases from the inner side of the electrode, because adsorption and desorption reactions take place mostly near the gas channel and electrolyte, respectively. Os is much higher in the reductive FC mode than in the oxidized EC mode; Cs is higher in the electrode near the fuel channel in the EC mode, and the carbon deposition may be improved in the dual-mode operation. (2) The effect of thickness: The meso-scale elementary reactions are affected by the thickness of the fuel electrode, with combined contributions of the active area and gas diffusion performance. 100 mm thick electrode is required at least, to provide enough region for the stabilization of the elementary reactions (charger exchanging and reforming) in dual-mode. Higher active current density and utilization rate of gas species are obtained in the global model, due to the overestimation on the electrochemical performance and gas utilization. (3) The effect of microstructure: Higher current density and gas utilization are obtained with larger porosity, lower tortuosity and smaller particle size, because of the improvement on the gaseous diffusion and specific active areas. While the effect of the particle size on the reaction performance is underestimated by the global model. The largest XNis and lowest XCOs are achieved when ε ¼ 0.35 based on the trade-off between the diffusion and active surface. Acknowledgements This work was supported by the National Natural Science Foundation of China (NSFC, 91634102), Ningbo Innovation team grant (Ningbo [2017]74), National Natural Science Foundation of China (NSFC, 51871126) and Natural Science Foundation of Ningbo (2018A610018). References [1] Wang Y, Leung D, Xuan J, Wang H. A review on unitized regenerative fuel cell technologies, part B: unitized regenerative alkaline fuel cell, solid oxide fuel cell, and microfluidic fuel cell. Renew Sustain Energy Rev 2017;75:775e95. [2] Santhanam S, Heddrich MP, Riedel M, Friedrich KA. Theoretical and experimental study of Reversible Solid Oxide Cell (r-SOC) systems for energy storage. Energy 2017;141:202e14. [3] Butera G, Jensen S, Clausen L. A novel system for large-scale storage of electricity as synthetic natural gas using reversible pressurized solid oxide cells. Energy 2019;166:738e54.

17

[4] Becker WL, Braun RJ, Penev M, Melaina M. Production of FischereTropsch liquid fuels from high temperature solid oxide co-electrolysis units. Energy 2012;47:99e115. [5] Zhang H, Kong W, Dong F, Xu H, Chen B, Ni M. Application of cascading thermoelectric generator and cooler for waste heat recovery from solid oxide fuel cells. Energy Convers Manag 2017;148:1382e90. [6] Zhang H, Wang J, Wang F, Zhao J, Miao H, Yuan J. Performance assessment of an advanced triple-cycle system based upon solid oxide fuel cells, vacuum thermionic generators and absorption refrigerators. Energy Convers Manag 2019;193:64e73. [7] Yang C, Wang J, Zhao J, Wu Y, Shu C, Miao H, Wang F, Ye W, Yuan J. CFD modeling and performance comparison of solid oxide fuel cell and electrolysis cell fueled with syngas. Int J Energy Res 2019;43:2656e77. [8] Gomez S, Hotza D. Current developments in reversible solid oxide fuel cells. Renew Sustain Energy Rev 2016;61:155e74. [9] Er-rbib H, Bouallou C. Modeling and simulation of CO methanation process for renewable electricity storage. Energy 2014;75:81e8. [10] Li W, Shi Y, Luo Y, Cai N. Elementary reaction modeling of CO2/H2O coelectrolysis cell considering effects of cathode thickness. J Power Sources 2013;243:118e30. [11] Mohammadi A, Mehrpooya M. A comprehensive review on coupling different types of electrolyzer to renewable energy sources. Energy 2018;158:632e55. [12] Khazaee I, Rava A. Numerical simulation of the performance of solid oxide fuel cell with different flow channel geometries. Energy 2017;119:235e44. [13] Yang C, Shu C, Miao H, Wang Z, Wu Y, Wang J, Zhao J, Wang F, Ye W, Yuan J. Dynamic modelling and performance analysis of reversible solid oxide fuel cell with syngas. Int J Hydrogen Energy 2019;44:6192e211. [14] Barelli L, Bidini G, Ottaviano A. Solid oxide fuel cell modelling: electrochemical performance and thermal management during load-following operation. Energy 2016;115:107e19. [15] Luo Y, Shi Y, Li W, Cai N. Dynamic electro-thermal modeling of co-electrolysis of steam and carbon dioxide in a tubular solid oxide electrolysis cell. Energy 2015;89:637e47. [16] Hecht E, Gupta GK, Zhu H, Dean AM, Kee RJ, Maier L, Deutschmann O. Methane reforming kinetics within a NieYSZ SOFC anode support. Appl Catal Gen 2005;295:40e51. [17] Janardhanan VM, Deutschmann O. CFD analysis of a solid oxide fuel cell with internal reforming: coupled interactions of transport, heterogeneous catalysis and electrochemical processes. J Power Sources 2006;162:1192e202. [18] Maier L, Schaodel B, Herrera Delgado K, Tischer S, Deutschmann O. Stream reforming of methane over nickel: development of a multi-step Surface reaction mechanism. Top Catal 2011;54:845e58. [19] Herrera Delgado K, Stotz H, Maier L, Tischere S, Zellner A, Deutschmann O. Surface reaction kinetics of steam- and CO2-reforming as well as oxidation of methane over Nickel-based catalysts, vol. 5; 2015. p. 871e904. [20] Karakaya C, OtterstAotter R, Maier L, Deutschmann O. Kinetics of the watergas shift reaction over Rh/Al2O3 catalyst. Appl Catal Gen 2014;470:31e44. [21] Karakaya C, Maier L, Deutschmann O. Kinetics for oxidation and reforming of CH4 over Rh/Al2O3 catalyst. Int J Chem Kinet 2016;48:144e60. [22] Hofmann P, Panopoulos KD, Fryda LE, Kakaras E. Comparison between two methane reforming models applied to a quasi-two-dimensional planar solid oxide fuel cell model. Energy 2009;34:2151e7. [23] Mizusaki J. Model for solid electrolyte gas electrode reaction kinetics; key concepts, basic model construction, extension of models, new experimental techniques for model confirmation, and future prospects. Electrochemistry 2014;82:819e29. [24] Utz A, Stormer H, Leonide A. Degradation and relaxation effects of Ni patterned anodes in H2-H2O atmosphere. J Electrochem Soc 2010;157: B920e30. [25] Vogler M, Hutter A, Gauckler L, Warnatz J, Bessler W. Modelling study of surface reactions, diffusion, and spillover at a Ni/YSZ patterned anode. J Electrochem Soc 2009;156:B663e72. [26] Bieberle A, Meier LP, Gauckler LJ. The electrochemistry of Ni pattern anodes used as solid oxide fuel cell model electrodes. J Electrochem Soc 2001;148: A646e56. [27] Yurkiv V, Starukhin D, Volpp HR, Bessler WG. Elementary reaction kinetics ofthe CO/CO2/Ni/YSZ electrode. J Electrochem Soc 2011;158:B5e10. [28] Hanna J, Lee WY, Ghoniem AF. Kinetics of carbon monoxide electrooxidation in solid-oxide fuel cells from Ni-YSZ patterned-anode measurements. J Electrochem Soc 2013;160:F698e708. [29] Sukeshini AM, Habibzadeh B, Becker BP, Stoltz CA, Eichhorn BW, Jackson GS. Electrochemical oxidation of H2, CO, and CO/H2 mixtures on patterned Ni anodes on YSZ electrolytes. J Electrochem Soc 2006;153:A705e15. [30] Peksen M. Numerical thermomechanical modelling of solid oxide fuel cells. Prog Energy Combust Sci 2015;48:1e20. [31] Yang L, Choi Y, Qin W, Chen H, Blinn K, Liu M. Promotion of water mediated carbon removal by nanostructured barium oxide/nickel interfaces in solid oxide fuel cells. Nat Commun 2011;2:1e9. [32] Hanna J, Lee W, Shi Y, Ghoniem A. Fundamentals of electo- and thermochemistry in the anodes of solid-oxide fuel cells with hydrocarbon and syngas fuels. Prog Energy Combust Sci 2014;40:74e111. [33] Luo Y, Shi Y, Li W, Cai N. Comprehensive modeling of tubular solid oxide electrolysis cell for co-electrolysis of steam and carbon dioxide. Energy 2014;70:420e34. [34] Zhu H, Kee RJ, Janardhanan VM, Deutschmann O, Goodwin DG. Modeling

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379

18

[35]

[36]

[37]

[38] [39]

[40]

[41]

[42]

[43]

[44]

[45]

[46]

[47]

[48]

C. Yang et al. / Energy xxx (xxxx) xxx elementary heterogeneous chemistry and electrochemistry in solid-oxide fuel cells. J Electrochem Soc 2005;152:A2427e40. Lee WY, Wee D, Ghoniem AF. An improved one-dimensional membraneelectrode assembly model to predict the performance of solid oxide fuel cell including the limiting current density. J Power Sources 2009;186:417e27. Li W, Shi Y, Luo Y, Wang Y, Cai N. Carbon deposition on patterned nickel/yttria stabilized zirconia electrodes for solid oxide fuel cell/solid oxide electrolysis cell modes. J Power Sources 2015;276:26e31. Menon V, Banerjee A, Dailly J, Deutschmann O. Numerical analysis of mass and heat transport in proton-conducting SOFCs with direct internal reforming. Appl Energy 2015;149:161e75. Narasimhaiah G, Janardhanan VM. Modeling CO2 electrolysis in solid oxide electrolysis cell. J Solid State Electrochem 2013;17:2361e70. Yu X, Shi Y, Wang H, Cai N, Li C, Tomov R, Hanna J, Glowacki B, Ghoniem A. Experimental characterization and elementary reaction modeling of solid oxide electrolyte direct carbon fuel cell. J Power Sources 2013;243:159e71. Luo Y, Li W, Shi Y, Wang Y, Cai N. Reversible H2/H2O electrochemical conversion mechanisms on the patterned nickel electrodes. Int J Hydrogen Energy 2017;42:25130e42. Menon V, Fu Q, Janardhanan VM, Deutschmann O. A model-based understanding of solid-oxide electrolysis cells (SOECs) for syngas production by H2O/CO2 co-electrolysis. J Power Sources 2015;274:768e81. Schluckner C, Suboti V, Lawlor V, Hochenauer C. CFD-simulation of effective carbon gasification strategies from high temperature SOFC Ni-YSZ cermet anodes. Int J Hydrogen Energy 2017;42:4434e48. Garcia-Camprubi M, Izquierdo S, Fueyo N. Challenges in the electrochemical modelling of solid oxide fuel and electrolyser cells. Renew Sustain Energy Rev 2014;33:701e18. Maier L, Schadel B, Herrera Delgado K, Tischer S, Deutschmann O. Steam reforming of methane over nickel: development of a multi-step surface reaction mechanism. Top Catal 2011;54:845e58. Yrukiv V, Utz A, Weber A. Elementary kinetic modeling and experimental validation of electrochemical CO oxidation on Ni/YSZ pattern anodes. Electrochemical Acta 2012;59:573e80. Laurencin J, Huberrt M, Couturier K, Bihan T, Cloetens P, Lefebvre-Joud F, Siebert E. Reactive mechanisms of LSCF single-phase and LSCF-CGO composite electrodes operated in anodic and cathodic polarisations. Electrochim Acta 2015;174:1299e316. Xue J, Miao H, He C. Electrolyte supported solid oxide fuel cells with the super large size and thin ytterbia stabilized zirconia substrate. J Power Sources 2015;279:610e9. Miao H, Li S, Wang Z, Sun S, Kuang M, Liu Z, Yuan J. Enhancing the pyridinic N content of Nitrogen-doped graphene and improving its catalytic activity for oxygen reduction reaction. Int J Hydrogen Energy 2017;42:28298e308.

Nomenclature

English letters a: concentration dependency coefficient Aj: Arrhenius parameter Av: specific area per volume (m2m3) b: concentration dependency coefficient C: inertial resistance factor Ci: molar concentration (mol/cm3 or mol/cm2) Ccharge: specific double-layer capacitance (Fm2) Cp: specific heat capacity (Jkg1K1) c: concentration dependency coefficient D: gas diffusion coefficient (m2s1) dp: pore diameter (mm) Ej: activation energy (Jmol1) F: Faraday constant (Cmol1) DG: Gibbs free energy (Jmol1) DH: enthalpy (Jmol1) I: current density (Am3) i0: exchange current density (Am2) Ji: mass flux of species K: equilibrium constant, specie index, reaction number k: Arrhenius reaction rate (mol, cm, s), heat conductivity (Wm1K1) kb: Boltzmann constant (JK1) M: molecular weight (kgmol1) n: number of charge transferred, number of particles per volume P: gas pressure (Pa) p: probabilities of conductor particle Q: heat source (Wm3) R: reaction rate (molm3s1 or molcm2s1), gas constant (Jmol1K1)

r: radius (mm) S: source term, initial sticking coefficient DS: entropy change (JK1) T: temperature (K) U: gas utilization rate Uoperation: operating cell voltage (V) V: velocity (ms1) X: molar fraction, surface coverage Y: mass fraction Z: average coordination number Greek letters

a: charge transfer coefficient b: permeability (m2), Arrhenius r: pre-exponential coefficient d: collision diameter or characteristic length (m) ε: porosity εkj: correction coefficients for coverage dependency h: overpotential (V) q: contact angle between electron/electrolyte particles ( ) l: characteristic Lennard-Jones energy of gas species m: viscosity (kg  m-1s-1) mkj: correction coefficients for coverage dependency y: stoichiometric coefficients p: PI s: charge conductivity (U  m-1) r: density (kg  m-3) t: tortuosity f: volume fraction of conductor particles G: total surface site density (mol/cm2) F: charge potential (V) UD: dimensionless diffusion collision Subscripts a: exponential constant act: active layer air: air electrode b: backward reaction, exponential constant c: exponential constant e: electron charge, electrochemical reaction eff: effective f: forward reaction, fluid fuel: fuel electrode g: gas i: index of species ideal: ideal voltage ion: ionic charge j: index of elementary reactions k: Kundsen diffusion, index of surface species m: mass pol: polarization r: total number of the reactions s: surface, solid Abbreviation AAL: air functional layer ADL: air diffusion layer EC: electrolysis cell ELE: electrolyte layer FC: fuel cell FAL: fuel functional layer FDL: fuel diffusion layer GDC: gadolinia-doped ceria LSM: Sr-doped LaMnO3 OCV: open circuit voltage or Nernst voltage (V) rSOFC: reversible solid oxide fuel cell SOFC: solid oxide fuel cell SOEC: solid oxide electrolytic cell TPB: triple phase boundary UDFs: user defined functions WGS: water gas shift YSZ: yttrium-stabilized zirconia

Please cite this article as: Yang C et al., Analysis of effects of meso-scale reactions on multiphysics transport processes in rSOFC fueled with syngas, Energy, https://doi.org/10.1016/j.energy.2019.116379