is the specific heat of species  i , which can be calculated as [32]

where a1 , a2   and a3   are  constants。

i i i

In the solid domain, heat conduction is considered, which can be described as

where Q  is the heat source, s is the thermal conductivity of the solid。

Heat source Q in SOFC can be classified as Ohmic heat source Qohm  , activation heat source

Qact

and entropy heat source  Qentr 。

Qohm  is the heat generation from ohmic loss, which can be expressed  as

Qact  is the heat generation from activation loss, which can be expressed as

Qentr is the heat generation from entropy change in the electrochemical reactions, which can    be

expressed as

where  F  is the Faraday constant, act    is the activation loss and S is the entropy。

2。2。Boundary Conditions

As described above, the charge transport equation, momentum transport equation, mass transport  equation and heat  transport are taken  into account  in  this  model。  Boundary conditions are

required for solving these coupled partial differential equations。 For easy description of boundary conditions, some boundaries in our model are tagged as indicated in Fig。 3c。 Many researchers didn’t give each of boundaries setting due to so many boundaries in 3D model, which may lead to some misunderstanding。 So in this work we describe the boundary setting for each boundary of each equation in order to be easy understanding。 Co-flow pattern is adopted in this study, since it has the most uniform temperature distribution than counter-flow pattern or cross-flow pattern [5]。

2。2。1。Boundary Setting for the Charge Transport Equation

The boundaries B1 and B12 are the top surface and the bottom surface of half unit cell, respectively, where electronic potential is set as E at boundary B1 and is specified at boundary B12。

The Nernst potential E is given as

G0 is the Gibbs free energy change at the condition that the partial pressure of H2, H2O

and O2 are all 1 bar。 The boundaries B4 and B8 represent the interfaces between rib and electrode。 Contact resistance is specified at these boundaries, which means the local current densities cross boundary B4 ( iribAn  ) and boundary B8 ( iCarib ) are determined by below expression。

is the contact resistance,  and  are the electric potentials at B4 on

rib side and anode side, respectively。     and  are the electric potentials  at  B8 on  cathode

e,Ca/rib e,rib/Ca

side and rib side, respectively。 The interface between electrolyte and anode is represented by B6,

where the normal electronic current density   is

an trans

due to the reaction at B6 and the normal ionic

current density is  ian     。 The interface between electrolyte and cathode is represented by B7, where   the

normal electronic current density is  ica due to the reaction at B7 and the normal ionic current density

is  ica       。 All others boundaries are set as electric insulation。

itrans and itrans

can be calculated by Butler-Volmer equations as follows [34, 35]:

are the forward and reverse reaction symmetric factor, respectively,  E