international journal of hydrogen energy 34 (2009) 7795–7806

Available at [Link]

journal homepage: [Link]/locate/he

Computational fluid dynamics modeling of a solid oxide

electrolyzer cell for hydrogen production

Meng Ni*
Department of Building and Real Estate, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong, China

article info abstract

Article history: A 2D computational fluid dynamics (CFD) model was developed to study the performance
Received 20 May 2009 of a planar solid oxide electrolyzer cell (SOEC) for hydrogen production. The governing
Received in revised form equations for mass continuity, momentum conservation, energy conservation and species
16 July 2009 conservation were discretized with the finite volume method (FVM). The coupling of
Accepted 23 July 2009 velocity and pressure was treated with the SIMPLEC (Semi-Implicit Method for Pressure
Available online 8 August 2009 Linked Equations – Consistent) algorithm. Simulations were performed to investigate the
effects of operating/structural parameters on heat/mass transfer and the electric charac-
Keywords: teristics of a planar SOEC. It is found that the gas velocity at the cathode increases
Heat and mass transfer significantly along the main flow channel, as the increase in H2 molar fraction decreases
Solid oxide fuel cell (SOFC) the density and viscosity of the gas mixture at the cathode. It is also found that increasing
High temperature steam the inlet gas velocity can enhance the SOEC performance. Another important finding is
electrolysis that the electrode porosity has small effect on SOEC performance. The results of this paper
Electrochemistry provide better understanding on the coupled heat/mass transfer and electrochemical
Porous media reaction phenomena in an SOEC. The model developed can serve as a useful tool for SOEC
Transport phenomena design optimization.
ª 2009 Professor T. Nejat Veziroglu. Published by Elsevier Ltd. All rights reserved.

1. Introduction The success of hydrogen energy relies on the development

of hydrogen production technologies. Presently, hydrogen is
The use of fossil fuels as energy sources has caused serious mainly produced from hydrocarbon reforming, such as
energy crisis and environmental pollution problems at a global natural gas steam reforming [1,2]. However, reforming of
level. Renewable energy sources, such as solar and wind, are hydrocarbon fuels is neither renewable nor clean. Hydrogen
promising to solve these problems. However, renewable can be produced in a cleaner manner by thermochemical
energy sources are usually site-specific and intermittent. cycles for water splitting, photocatalytic water splitting, or
Therefore, hydrogen has been identified as an ideal energy water electrolysis [3,4]. So far, both thermochemical cycles
carrier and a promising fuel to better utilize the renewable and photocatalysis have efficiency too low to be economically
energy sources. Hydrogen can be converted into electricity via sound. Water electrolytic hydrogen is presently the most
fuel cells in an efficient and environmentally friendly manner. viable technology for large-scale renewable hydrogen
Due to the great promise of hydrogen fuel, there is increasingly production [3,4]. Electrolyzer cells can be driven by solar cells,
more interest from government, industry, and academy to or wind turbines together with waste heat from industry to
develop hydrogen-related technologies. produce hydrogen in an environmentally friendly manner.

* Tel.: þ86 852 2766 4152; fax: þ86 852 2764 5131.
E-mail addresses: memni@[Link], bsmengni@[Link]
0360-3199/$ – see front matter ª 2009 Professor T. Nejat Veziroglu. Published by Elsevier Ltd. All rights reserved.
doi:10.1016/[Link].2009.07.080
7796 international journal of hydrogen energy 34 (2009) 7795–7806

Nomenclature rp Mean pore radius of electrode, mm

R Universal gas constant, 8.3145 J/(mol K)
Bg Permeability of electrode
Sm Source term in continuity equation
cp Heat capacity
Sx, Sy Source terms in momentum equations
dH2 Thickness of anode, mm
T Operating temperature, K
dO2 Thickness of cathode, mm
U Velocity in x direction, m/s
Deff
i Effective diffusion coefficient of species i, cm2/s
Uin Gas velocity at the SOFC inlet, m/s
DH2 O;k Knudsen diffusion coefficient of H2O, cm2/s
V SOFC operating potential (V); Velocity in y
DH2 H2 O Binary diffusion coefficient of H2 and H2O, cm2/s
direction, m/s
E Equilibrium potential, V
3 Electrode porosity
E0 Reversible potential at standard condition, V
hact,a Activation overpotential at anode, V
F Faraday constant, 9.6485  104 C/mol
hact,c Activation overpotential at cathode, V
J Current density, A/m2
hohmic Ohmic overpotential of the electrolyte, V
k Thermal conductivity
m Viscosity of gas mixture
L Thickness of electrolyte, mm
x Electrode tortuosity
Mi Molecular weight of species i, kg/mol
r Density of the gas mixture, kg/m3
P Operating pressure, bar
si,j Mean characteristic length of species i and j
PIi Partial pressure (bar) of species i (i ¼ H2O, H2, O2) at
UD Dimensionless diffusion collision integral
electrode–electrolyte interface

Among all the electrolyzer types, high temperature solid oxide modeling, despite of complicated design for the gas supply
electrolyzer cells (SOEC) are more advantageous than lower sub-system. The effect of flow mode on SOEC system perfor-
temperature proton exchange membrane (PEM) electrolyzers mance will be investigated in future studies. Different from
and alkaline electrolyzers because SOEC can produce the Hawkes work, the model is developed and solved with an
hydrogen at a high electrochemical reaction rate with a low in-house code, which allows easy modification and contrib-
electrical energy requirement [5–10]. In recent years, the US utes to the continuing development of high-level software. In
Department of Energy (DOE) has initiated several projects to addition, the focus of the present study is to investigate the
develop the SOEC hydrogen production technology. In addi- effects of operating/structural parameters on the distributions
tion, there are growing interests in integrating SOEC with of current density, gas composition, velocity, and temperature
nuclear energy and geothermal energy technologies to in a planar SOEC.
improve the hydrogen production efficiency [11,12].
Most research works done on SOEC hydrogen production
were experimental in nature, with an emphasis on the
development of new materials to cope with the high operating 2. Model development
temperature and strong reducing/oxidizing environments
required for stable and efficient operation [13–32]. Compa- 2.1. Geometry and operation
rably, only a few mathematical modeling studies on SOEC
hydrogen production are available in literature. Ni et al. The 2D view of a planar SOEC and the computation domain
[33–35] developed 1D electrochemical models for SOEC based are shown in Fig. 1, which includes two flow streams and solid
on oxygen ion conducting electrolyte and proton conducting area that locates the SOEC (electrodes and electrolyte) and the
electrolyte, respectively. A 1D dynamic model of a cathode- interconnects. U and V are the velocity components in the x
supported SOEC was developed by Udagawa et al. [36,37]. The and y directions, respectively. During operation, an electric
1D models are capable of predicting the current density– potential sufficient for electrolysis is applied to the SOEC.
voltage characteristics of an SOEC and are useful tools to Steam/hydrogen gas mixture and air are supplied to the
understand the working mechanisms and to identify the key cathode and anode flow channels, respectively. At the cathode
overpotential losses of an SOEC. In order to better understand side, steam molecules are transported through the porous
the heat/mass transfer and electrochemical characteristics of cathode layer to the triple phase boundary (TPB) at the
an SOEC, computational fluid dynamics (CFD) model is cathode–electrolyte interface, where they receive electrons
needed. Hawkes et al. [7] at Idaho National Laboratory have and decompose to hydrogen gas and oxygen ions. The
done some preliminary CFD study of SOEC using the produced hydrogen molecules are transported out of the
commercial software – FLUENT. The effects of operating porous cathode to the cathode fluid channel and oxygen ions
voltage on the distribution of gas composition, current are subsequently transported through the dense electrolyte to
density, and temperature in a SOEC stack were simulated and the anode side. At the anode-electrolyte, oxygen ions lose
discussed. However, the effects of other operating parameter electrons and become oxygen molecules. The generated
(gas velocity) and electrode microstructure (porosity, pore size oxygen molecules are transported out of the porous anode
et al.) on SOEC performance remain to be investigated. In this layer to the anode flow channel while the electrons produced
paper, a 2D CFD model is developed for a co-flow planar SOEC. are transported to the cathode side via the external circuit to
The co-flow mode is adopted because it is widely used in SOFC complete the cycle. The net reaction of SOEC can be written as
 international journal of hydrogen energy 34 (2009) 7795–7806 7797

x=0 x = xL
y = yM
Interconnect 0.5 mm
H2O; H2 H2; H2O
H2; H2O 1.0 mm
H2O H2

Porous cathode 0.5 mm

Electrolyte 0.1 mm
Porous anode 0.1 mm
O2
Air (O2; N2) 1.0 mm
Air (O2; N2)

y, V Interconnect 0.5 mm
y=0

x, U

Fig. 1 – Two-dimensional view of a planar SOEC.

H2 O/H2 þ 0:5O2 (1) considered. The ohmic overpotential of the SOEC electrolyte
can be easily calculated according to Ohm’s law [39]. In the
2.2. Electrochemistry present study, the concentration overpotentials are implicitly
included in the reversible potential (E ), since the partial pres-
In an SOEC, a current collector is placed along the main flow sures of gaseous species in Eq. (3) were evaluated at the elec-
stream (not shown in the above figure). Due to the high elec- trode–electrolyte interface. In order to obtain these interfacial
trical conductivity of the current collector, the potential of the partial pressures at different locations along the main flow
SOEC along the main flow stream remains constant but the stream, the characteristics of heat/mass transfer in the porous
local current density varies [8,14]. In this study, the current electrodes and the gas flow channels must be obtained. Details
density–potential relationship will be studied with the elec- of the heat/mass transfer are presented in the next section.
trochemical model and CFD model presented in the subse-
quent section. 2.3. Heat and mass transfer in an SOEC
During operation, the electrochemical potential (V) of an
SOEC can be expressed as [33–35,38]: The heat and mass transfer phenomena in an SOEC are gov-
V ¼ E þ hact;c þ hact;a þ hohmic (2) erned by the conservation laws of momentum, energy, and
species. Steady-state laminar flow assumption is used for this
 0:5 ! 2D CFD model. This assumption is valid for an SOEC because
RT PIH2 PIO2 of small Reynolds number in an SOEC. In addition, radiation
E ¼ E0 þ ln (3)
2F PIH2 O heat transfer is not included due to its small effect on heat
transfer processes in SOECs.
where E is the equilibrium potential; hohmic is the ohmic Continuity equation can be written as:
overpotential; hact,a and hact,c are the activation overpotentials
at the anode and cathode, respectively; PI is the partial pres- vðrUÞ vðrVÞ
þ ¼ Sm (4)
sure at the electrode–electrolyte interface, and subscripts H2O, vx vy
H2 and O2 represent the reactant (steam) and hydrogen and
where r is the density of the gas mixture.
oxygen products, respectively; R is the universal gas constant
The volumetric source term Sm in the above equation
(8.3145 J/(mol K)); F is the Faraday constant (96,485 C/mol); and
accounts for the mass balance due to the electrochemical
T is the local operating temperature.
reactions at the active surface (electrode–electrolyte interface)
The equilibrium potential at standard pressure (E0) can be
[40,41]. Therefore, this source term is zero in most of the
calculated as Gibbs free energy change for Eq. (1) over 2F [33].
computation zone and non-zero only in the control volumes
The activation overpotentials are related to the electrode
(in the porous electrode) neighboring the electrode–electrolyte
kinetics at the reaction sites. They represent the overpotential
interface.
incurred due to the activation necessary for charge transfer.
Momentum equations are written as:
The electrode activation overpotential–current density rela-
tionship can be expressed by the well-known Butler–Volmer    
vðrUUÞ vðrVUÞ vP v vU v vU
equation. The ohmic overpotential represents the electric þ ¼ þ m þ m þ Sx (5)
vx vy vx vx vx vy vy
voltage loss due to the resistance to the transport of electrons or
ions. As the electrodes and interconnects have electrical
   
conductivity much higher than the electrolyte, therefore the vðrUVÞ vðrVVÞ vP v vV v vV
þ ¼ þ m þ m þ Sy (6)
ohmic losses at the electrodes and interconnects were not vx vy vy vx vx vy vy
7798 international journal of hydrogen energy 34 (2009) 7795–7806

The inclusion of the source terms Sx and Sy allows the For Knudsen diffusion, as gas molecules frequently collide
momentum equations to be valid for both the porous elec- with the walls of the pores, the transport of molecules can be
trodes and gas flow channels. For example, the source term Sx modeled using the kinetic theory,
can be written as [42]: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
4 8RT
mU DH2 O;k ¼ rp (11)
Sx ¼  (7) 3 pMH2 O
Bg
where rp is the mean pore radius and MH2 O is the molecular
where Bg is the permeability of the porous electrodes and m is
weight of steam. The molecular binary diffusion coefficient
the viscosity of the gas mixture in the electrodes.
DH2 OH2 can be obtained by the Chapman–Enskog theory of
In the gas flow channel, the source terms Sx and Sy are zero
ideal gas [50],
by setting the permeability (Bg) to be infinite. As there are no
reliable permeability data for SOEC, the data for SOFC is  1=2
1 1 T3=2
adopted for SOEC study. Since the permeability value of about DH2 OH2 ¼ 0:00133 þ 2
(12)
MH2 MH2 O PsH2 O;H2 UD
2.0  1010 m2 is widely used in SOFC modeling studies [43–47],
it is used in the present study. where sH2 O;H2 is the mean characteristic length of species H2O
Based on thermal equilibrium assumption for the porous and H2; and UD is a dimensionless diffusion collision integral.
electrodes, only one energy equation is solved for the gas The values of UD and sH2 O;H2 can be obtained from reference [50].
species and the solid structure: Eqs. (10)–(12) describe how to calculate the effective diffu-
        sion coefficient for species diffusion in porous media. For
v rcp UT v rcp VT v vT v vT diffusion in gas channels, the Knudsen diffusion coefficient
þ ¼ k þ k þ ST (8)
vx vy vx vx vy vy should be removed from Eq. (10) since the Knudsen diffusion
where cp and k are the heat capacity and thermal conductivity, effect is negligible. For the effective diffusion coefficient of O2
respectively. at the anode, its calculation procedure is similar to the H2O at
Equation (8) balances the convected thermal energy, the the cathode (Eqs. (10)–(12)) and is not presented here.
heat conduction through the gas species and the solid struc-
ture and a source term ST. The heat source term mainly results 2.4. Boundary conditions
from the entropy change of the electrochemical reaction [48].
In addition, Joule heating due to overpotentials can also At x ¼ 0, a constant flow velocity U ¼ Uin is specified at the
contribute to the heat generation in an SOEC [49]. inlets of the gas flow channels, while U ¼ 0 is specified at the
The species equation can be written as: inlet for the solid part. Similarly, constant gas compositions
    are specified at the inlets of the gas flow channels and
vðrUYi Þ vðrVYi Þ v vYi v vYi constant temperature is specified at the inlet for both the gas
þ ¼ rDi;m þ rDi;m þ Ssp (9)
vx vy vx vx vy vy flow channels and the solid part. At the bottom and the top of
where Yi is the mass fraction of species i and Ssp is the source the computation zone ( y ¼ 0 and y ¼ yM), the interconnect (see
term accounting for the generation or consumption of species Fig. 1), thermally adiabatic conditions are assumed:
i due to electrochemical reaction, which can be determined by ðvT=vyÞ ¼ 0. At the outlet of the computation domain (x ¼ xL),
a similar way with Sm. Di,m is the effective diffusion coefficient the gas flow channels have ðvU=vxÞ ¼ 0, ðvT=vxÞ ¼ 0 and
of species i in the gas mixture. ðvYi =vxÞ ¼ 0, while for the solid part, U ¼ 0, V ¼ 0, ðvT=vxÞ ¼ 0
The diffusion in the porous SOEC electrodes is mainly and ðvYi =vxÞ ¼ 0 are employed.
based on two mechanisms, namely, molecular diffusion and
Knudsen diffusion [33–35]. Molecular diffusion is the domi- 2.5. Numerical methodology
nant mechanism if the pore size is much larger than the mean
free path of the molecular species. In this case, the molecule– The governing equations presented above are discretized into
molecule interaction governs the diffusion process. On the algebraic equations by using the Finite Volume Method (FVM)
other hand, if the pore size is much smaller than the mean [51], which is widely used in computational fluid dynamics
free path of the species, the molecule-pore wall interaction (CFD).
dominates over the molecule–molecule interaction. Thus, In the 2D Cartesian coordinate system, the momentum
Knudsen diffusion becomes an important mechanism. In equations are solved for the velocity components on a stag-
most porous structures, both mechanisms are significant. gered grid arrangement. The convection term and diffusion
At the cathode side, the effective diffusion coefficient of term are treated with upwind and central difference schemes,
steam can be expressed by combining these two diffusion respectively. The pressure–velocity coupling is treated by the
mechanisms using the Bosanquet formula [50], SIMPLEC algorithm (Semi-Implicit Method for Pressure Linked
  Equations – Consistent) with the incompressible form of
1 x 1 1
¼ þ (10) pressure-correction equation. The discretized equations are
Deff
H2 O
3 DH2 H2 O DH2 O;k
solved with the ADI (Alternating Direction Implicit) scheme
where x/3 is the ratio of cathode tortuosity to porosity; based on TDMA (Tridiagonal matrix method) algorithm. After
x=ð3DH2 H2 O Þ is the reciprocal of the effective molecular diffu- the velocity field is obtained, the discretized energy equation
sion coefficient for H2–H2O binary system; and x=ð3DH2 O;k Þ is the and species equations are solved respectively. The properties
reciprocal of the effective Knudsen diffusion coefficient for of the gas mixtures are updated according to the computed
steam. gas composition, such as viscosity, thermal conductivity and
 international journal of hydrogen energy 34 (2009) 7795–7806 7799

effective diffusion coefficients. It should be mentioned that

Table 1 – Parameters used in 2D CFD modeling analyses
the source terms in the governing equations are related to the of a planar SOEC [33,52].
local current density. The local current density varies with x
Parameter Value
while the cell potential is constant. The electrochemical
model developed in the previous section is used to determine Operating temperature, T (K) 1173
the local current density for a given potential. In the iteration, Operating pressure, P (bar) 1.0
Electrode porosity, 3 0.4
the local current density is computed based on the updated
Electrode tortuosity, x 3.0
gas composition and temperature for a given terminal cell
Average pore radius, rp (mm) 0.5
potential. Computation is repeated till convergence is Cathode-supported SOEC
obtained for the momentum equations, energy equation and Anode thickness da (mm) 100
species equations. The in-house 2D CFD code was written and Electrolyte thickness, L (mm) 100
developed in Fortran. Cathode thickness, dc (mm) 500
Height of gas flow channel (mm) 1.0
Length of the planar SOEC (mm) 20
Thickness of interconnect (mm) 0.5
3. Code evaluation Inlet velocity: Uin (m/s) 1.0
Anode inlet gas molar ratio: O2/N2 0.21/0.79
To validate the present CFD code, simulations of 2D fluid flow Cathode inlet gas molar ratio: H2O/H2 0.9/0.1
and heat transfer in a square channel have been performed SOEC operating potential (V) 1.4
for comparison with the data obtained from the commercial Thermal conductivity of SOEC
component (W/(mK))
CFD software: FLUENT. The velocity field and temperature
Cathode 11.0
field were examined and excellent agreement between the
Electrolyte 2.7
present code and FLUENT was found, thus validated Anode 6.0
the present code. To check the independence of the results on Interconnect 1.1
the number of grid points, numerical experiments have been
carried out to ensure that the results are grid independent. It
was found that the results do not change significantly when molar fraction, because: (1) the generation rate of H2 is twice of
the number of control volumes is increased beyond 100  148 O2 (Eq. (1)) and (2) H2O is consumed at the cathode while the
(100 for x direction and 148 for y direction). In the subsequent amount of N2 at the anode does not decrease, diluting the
study, the control volumes of 100  148 were used. concentration of O2. As a result, the molar fraction of H2
increases from 0.1 at the inlet to be higher than 0.28 at the
outlet, while the molar fraction of O2 only increases from 0.21
4. Results and discussion at the inlet to be about 0.25 at the outlet. At the same x, the
molar fractions of H2 and O2 in the porous electrodes are
The dimensions and input parameters used in the present higher than in the gas flow channel. The molar fraction
study were obtained from literature and summarized in Table gradient in the porous electrodes represents the resistance of
1 [33,52]. As the operating potential, inlet velocity and elec- the porous structure to the transport of gas species and is the
trode porosity are important operating/structural parameters source of concentration overpotentials. It should be noted that
for SOEC operation, their effects on the SOEC performance are this difference is more significant for the anode than the
investigated in the subsequent sections. cathode, despite of thin anode used in the simulation. This is
because the large molecular weight of O2 and N2 results in
4.1. Distributions of local current density, gas
composition, temperature and velocity in an SOEC at 12000
different operating potentials
10000
Current density, J (A/m )

Fig. 2(a) and (b) show the distribution of local current density
2

in the planar SOEC at operating potentials of 1.2 V and 1.4 V,

8000
respectively. As expected, the SOEC current density increases
with increasing operating potential. In addition, the local
6000
current density decreases with an increase in x. This is
because the consumption of the reactants (H2O) and genera-
tion of product (H2 and O2) along the channel increases the 4000
Applied potential, (V)
Nernst potential, which in turn decreases the total over- 1.2
potential for a given SOEC operating potential (Eq. (2)). 2000 1.4
Accordingly, the local current density is decreased, based on
the electrochemical model developed in the previous section. 0
The molar fractions of H2 (cathode) and O2 (anode) at an 0.000 0.005 0.010 0.015 0.020
operating potential of 1.2 V are shown in Fig. 3(a). As expected, Local along the main flow stream, X (m)
both H2 molar fraction and O2 molar fraction increase with
increasing x due to the electrochemical reaction (Fig. 2). The Fig. 2 – Distribution of local current density in a planar
change in H2 molar fraction is much greater than that in O2 SOEC at different operating potentials.
7800 international journal of hydrogen energy 34 (2009) 7795–7806

lower effective diffusion coefficients (Eq. (10)–(12)). At a higher seen from Table 2. At the same x, the temperature gradient of
SOEC operating potential (1.4 V), the molar fractions of both H2 the gas mixture is higher than the solid part, as the solid part
and O2 are found higher along the main flow stream than at with high thermal conductivity can conduct heat at a fast rate.
a potential of 1.2 V (Fig. 3(b)). As shown in Fig. 2, higher Axial velocity contours in an SOEC are presented in Fig. 5.
operating potential leads to higher current density and higher The interesting finding is that the axial velocity contours for the
rate of electrochemical reaction, which in turn higher gener- anode are quite different from the cathode. At the anode side,
ation rates of H2 and O2. the air speed accelerates near the inlet due to the development
Fig. 4 shows the distribution of temperature in the SOEC. At of boundary layer and then the velocity profile remains almost
an operating potential of 1.2 V, the temperature of the gas unchanged at the downstream, which is similar to common
mixtures and the solid part decreases slightly with increasing channel flows without electrochemical reactions. Due to small
x (Fig. 4(a)). For comparison, the SOEC temperature increases variation in the gas composition at the anode side (Fig. 3), the
along the main flow stream at an operating potential of 1.4 V transport properties of the gas mixture do not vary significantly
(Fig. 4(b)). As mentioned in model development, the heat along the main flow stream, thus the electrochemical reaction
source term of the energy equation is composed of two parts: does not significantly influence the velocity field. At an oper-
(1) heat sink due to the endothermic nature of the electro- ating potential of 1.2 V, the core gas velocity at the outlet of the
chemical reaction [48] and (2) heat generation due to irre- anode is about 40% higher than at the inlet (Fig. 5(a)). For
versible overpotential losses. Therefore, in an SOEC, the net comparison, the core velocity of H2/H2O mixture at the cathode
heat generation can be negative, zero, or positive, depending side increases considerably along the gas flow channel. It
on the current density (or operating potential) [7]. The net heat
flux is negative at a low operating potential, zero at the
thermal neutral voltage (TNV), and positive at an operating a 0.004
potential higher than the TNV. As the TNV of SOEC is only
a weak function of operating temperature and its value is 1172.87
between 1.2 and 1.3 V (i.e. 1.287 V at 1073 K) [7], the net heat 1172.73
flux is slightly negative at an operating potential of 1.2 V while 0.003 1172.58
1172.44 1172.15
it is positive at a potential of 1.4 V. As a result, the SOEC
1172.29
temperature decreases at 1.2 V but increases at 1.4 V, as can be 1172.0
y (m)

0.002
a 0.004

H2 molar fraction
0.003 0.126 0.152 0.001
0.178
0.204
0.23 0.256 0.282
y (m)

0.002

0.000
0.252
0.000 0.005 0.010 0.015 0.020
0.216 0.222 0.228 0.234 0.24 0.246
0.001
Location along the main flow stream, x (m)

O2 molar fraction
0.000 b 0.004
0.000 0.005 0.010 0.015 0.020
Location along the main flow stream, x (m) 1174.06
1175.13
1176.21
b 0.004 0.003
1177.29
1178.37 1179.44
H2 molar fraction 1180.52
0.003 0.139 0.178
y (m)

0.217 0.002
0.256
0.295 0.334 0.373
y (m)

0.002

0.273 0.001
0.219 0.228 0.237 0.246 0.255 0.264
0.001

O2 molar fraction
0.000 0.000
0.000 0.005 0.010 0.015 0.020 0.000 0.005 0.010 0.015 0.020

Location along the main flow stream, x (m) Location along the main flow stream, x (m)

Fig. 3 – Molar fractions of H2 and O2 in a planar SOEC at Fig. 4 – Temperature (K) distribution in a planar SOEC at
different operating potentials – (a) 1.2 V and (b) 1.4 V. different operating potentials – (a) 1.2 V and (b) 1.4 V.
 international journal of hydrogen energy 34 (2009) 7795–7806 7801

12000
Table 2 – Comparison of the SOEC outlet temperature at
different operating potentials.
10000

Current density, J (A/m )

Operating potential (V) Average SOEC temperature

2
at the outlet (K)
8000
1.2 1171.5
TNV 1173
1.4 1181.5 6000

Note: TNV is defined as the operating potential at which the SOEC

4000
temperature at the outlet and the inlet are the same. Inlet gas velocity, Uin (m/s)
0.5
2000
should be noted that the core velocity influences the thickness 5.0
of the boundary layer and hence influences the heat transfer in
0
the SOEC. As shown in Fig. 3, the hydrogen molar fraction
0.000 0.005 0.010 0.015 0.020
increases significantly along the main flow stream. The varia-
Location along the main flow stream, X (m)
tion in gas composition considerably changes the transport
properties of the gas mixture, as the density and viscosity of H2 Fig. 6 – Distribution of local current density in an SOEC at
are much lower than H2O. As a result, the core gas velocity at the different inlet velocities.
outlet of the cathode is about 70% higher than at the inlet
(Fig. 5(a)). With an increase in operating potential (Fig. 5(b)), the
gas velocity along the cathode channel is increased due to 4.2. Effect of inlet velocity
enhanced electrochemical reaction rate. For comparison, the
gas velocity at the anode shows little dependence on the The local current density distributions at different inlet
operating potential. velocities (0.5 m/s and 5.0 m/s) and an operating potential of
1.4 V are shown in Fig. 6. The local current densities at the
inlet are the same for different inlet velocities. At the
a 0.004

a 0.004
0.0
0.003 H2 molar fraction
1.26 1.47 1.68
1.0 0.003 0.15
0.2
0.25 0.4 0.47
y (m)

0.21
0.002 0.3 0.35
y (m)

0.0 0.002 0.54

1.0 1.26
0.001 1.47
0.225 0.24 0.255 0.285
0.001 0.27 0.3 0.315
0.0
O2 molar fraction
0.000
0.000 0.005 0.010 0.015 0.020 0.000
0.000 0.005 0.010 0.015 0.020
Location along the main flow stream, x (m) Location along the main flow stream, x (m)

b 0.004

b 0.004
0.0
H2 molar fraction
0.003
1.204 1.447
1.0 1.69 0.003
0.113 0.126 0.139 0.152 0.165
y (m)

0.233
0.002
y (m)

0.178
0.002
0.0 0.191

1.0 1.204 0.225

0.001 1.447
0.213 0.216 0.219 0.222
0.001

0.0
O2 molar fraction
0.000 0.000
0.000 0.005 0.010 0.015 0.020 0.000 0.005 0.010 0.015 0.020
Location along the main flow stream, x (m) Location along the main flow stream, x (m)

Fig. 5 – Velocity (m/s) contours in an SOEC at different Fig. 7 – Molar fractions of H2 and O2 in a planar SOEC at
operating potentials – (a) 1.2 V and (b) 1.4 V. different inlet velocities – (a) 0.5 m/s and (b) 5.0 m/s.
7802 international journal of hydrogen energy 34 (2009) 7795–7806

a 0.004 a 0.004

1173.7
1174.5
0.0
0.003 1175.3
1176 0.003
0.6 0.9
1177.5 0.75 1.05
1176.8 1178.3
y (m)

y (m)
0.002 0.15
0.002
0.0
0.45
0.001 0.6
0.001

0.0
0.000
0.000 0.005 0.010 0.015 0.020 0.000
0.000 0.005 0.010 0.015 0.020
Location along the main flow stream, x (m)
Location along the main flow stream, x (m)
b 0.004
b 0.004

1173.8

0.003 1174.6 0.0

1175.4 0.003
1176.2 1177 6.3
7.2
1177.8
1178.7
y (m)

0.002
y (m)

0.9
0.002
0.0

0.001
6.3
0.001 7.2

0.000 0.0
0.000 0.005 0.010 0.015 0.020
0.000
Location along the main flow stream, x (m) 0.000 0.005 0.010 0.015 0.020
Location along the main flow stream, x (m)
Fig. 8 – Temperature (K) distribution in a planar SOEC at
different inlet velocities – (a) 0.5 m/s and (b) 5.0 m/s. Fig. 9 – Velocity (m/s) contours in a planar SOEC at different
inlet velocities – (a) 0.5 m/s and (b) 5.0 m/s.

downstream, the local current density at an inlet velocity of

0.5 m/s decreases much more significantly with increasing x 12000
than that at an inlet velocity of 5.0 m/s. At a higher inlet
velocity, the variation in gas composition along the main flow 10000
Current density, J (A/m )
2

channel is less significant than at a lower velocity (Fig. 7(a) and

(b)). For example, the hydrogen molar fraction at the SOEC 8000
outlet is higher than 0.5 at an inlet velocity of 0.5 m/s, while it
is less than 0.2 at an inlet velocity of 5.0 m/s. As mentioned in 6000
the previous section, the lower hydrogen molar fraction
decreases the Nernst potential, leading to higher current 4000 Electrode porosity
density at a given operating potential. The higher current
0.3
density will increase the ohmic heating due to overpotential 0.6
2000
losses, which in turn tends to increase the SOEC temperature.
However, the SOEC temperature only slightly increases with
0
increasing inlet gas velocity (Fig. 8(a) and (b)). This is because 0.000 0.005 0.010 0.015 0.020
the cooler gas mixture can bring away more heat from the
Location along the main flow stream, X (m)
SOEC at the downstream at a higher inlet velocity (5.0 m/s),
despite of higher heat generation rate. At a low inlet gas Fig. 10 – Distribution of local current density in a planar
velocity (0.5 m/s), the anode gas velocity reaches its maximum SOEC at different electrode porosities.
 international journal of hydrogen energy 34 (2009) 7795–7806 7803

a 0.45 a 0.004

Molar fraction of H2 at the

0.40 H2 molar fraction
electrode-electrolyte interface
Current density, J (A/m )
2

0.003 0.14 0.18 0.22 0.26 0.3

0.35 0.34 0.38

0.30

y (m)
0.002
0.25
0.22 0.23 0.24 0.25 0.26 0.27
0.20 0.001

0.15 Electrode porosity

O2 molar fraction
0.3
0.10 0.000
0.6 0.000 0.005 0.010 0.015 0.020

0.05 Location along the main flow stream, x (m)

0.000 0.005 0.010 0.015 0.020
Location along the main flow stream, X (m) b 0.004
b 0.29 H2 molar fraction
0.003 0.14 0.18 0.22 0.26 0.3
0.28 Molar fraction of O2 at the 0.34 0.38
Current density, J (A/m )
2

electrode-electrolyte interface
0.27
y (m)
Electrode porosity 0.002
0.3
0.26
0.6
0.22 0.23 0.24 0.25 0.26 0.27
0.25 0.001

0.24
O2 molar fraction
0.23 0.000
0.000 0.005 0.010 0.015 0.020
0.22 Location along the main flow stream, x (m)

0.21 Fig. 12 – Molar fractions of H2 and O2 in a planar SOEC at

0.000 0.005 0.010 0.015 0.020 different electrode porosities – (a) 0.3 and (b) 0.6.
Location along the main flow stream, X (m)

Fig. 11 – Gas molar fraction at the electrode–electrolyte

interface – (a) H2 molar fraction and (b) O2 molar fraction.
An important finding of this study is that the SOEC
performance is insensitive to the change of electrode porosity,
as can be seen from the local current density in the SOEC
value near the inlet and then become almost unchanged at the
(Fig. 10). This is different from previous 1D modeling analyses,
downstream; while the cathode gas velocity increases
where significant effect of electrode porosity on SOEC
considerably with a further increase in x (Fig. 9(a)). At a high
performance is observed [33–35]. In 1D modeling analyses,
inlet gas velocity (5.0 m/s), the gas velocity at the cathode is
higher electrode porosity favors gas transport in the porous
similar to the anode, as the variation of gas mixture properties
electrodes, leading to reduced concentration overpotential
does not vary significantly along the main flow stream
and higher SOEC performance. In the present 2D modeling
(Fig. 9(b)). From the present study, it can be seen that
analyses, an increase in electrode porosity decreases the
increasing the inlet gas velocity can enhance the SOEC
gradient of gas molar fraction in the porous electrode. This
performance. However, in practice, the inlet gas velocity
can decrease the molar fractions of H2 and O2 at the electrode–
cannot be infinitely increased as this can leads to decreased
electrolyte interface (Fig. 11 (a) and (b)), which is consistent
H2O utilization and may reduce the system efficiency at a too
with the 1D modeling results. Obviously, the lower interfacial
large velocity. Therefore, a thermodynamics analysis of the
molar fractions of H2 and O2 at a higher electrode porosity
SOEC plant is needed to optimize the operating parameters for
tends to decrease the Nernst potential and thus increase the
overall performance enhancement.
local current density at a given operating potential. However,
as can be seen from the previous analysis, the increase in
4.3. Effect of electrode porosity current density tends to increase the molar fractions of H2 and
O2 at the downstream because of higher rate of electro-
Electrode porosity is an important structural parameter deter- chemical reactions. The combined effects result in small
mining the SOEC performance. In this section, the effects of effect of electrode porosity on SOEC performance. This finding
electrode porosity on the distributions of local current density, is consistent with another study on SOFC in literature, which
gas composition, temperature and velocity are studied at also reported small effect of electrode porosity on SOFC
a given SOEC potential of 1.4 V and an inlet velocity of 1.0 m/s. performance [53]. It can be seen that the gas compositions are
7804 international journal of hydrogen energy 34 (2009) 7795–7806

a 0.004
a 0.004

1173.9
1174.8
0.003 1175.7 0.003
1176.6 1.38
1178.4 1.61 1.84
1177.5 1179.3
1180.2
y (m)

y (m)
0.23
0.002 0.002

0.0

1.15
1.38
0.001 0.001

0.0

0.000 0.000
0.000 0.005 0.010 0.015 0.020 0.000 0.005 0.010 0.015 0.020

Location along the main flow stream, x (m) Location along the main flow stream, x (m)

b 0.004 b 0.004

1174.3 1175.6 0.0

0.003
0.003 1176.9
1.38 1.61 1.84
1178.2
1179.5 1180.8
1182.1
y (m)

0.23
0.002
y (m)

0.002
0.0

1.15
1.38
0.001
0.001

0.0

0.000
0.000 0.000 0.005 0.010 0.015 0.020
0.000 0.005 0.010 0.015 0.020
Location along the main flow stream, x (m)
Location along the main flow stream, x (m)
Fig. 14 – Velocity (m/s) contours in a planar SOEC at
Fig. 13 – Temperature (K) distribution in a planar SOEC at
different electrode porosities – (a) 0.3 and (b) 0.6.
different electrode porosities – (a) 0.3 and (b) 0.6.

different in the porous electrodes but are identical in the flow

agreement was found, thus the present code is suitable for
channels at different electrode porosities (Fig. 12(a) and (b)). It
studying the heat/mass transfer phenomena in an SOEC.
should be mentioned that 1D electrochemical models did not
Parametric analyses were conducted to understand the
consider the changes in gas properties along the main flow
coupled electrochemical reaction and heat/mass transfer
stream and therefore failed to include the adverse effect of
phenomena in a planar SOEC.
current density on gas composition. Due to small difference in
It is found that the velocity contours at the anode and
current density, the SOEC working at 1.4 V shows similar
cathode are quite different. Due to electrochemical reaction,
temperature and velocity contours for different electrode
H2 molar fraction and O2 molar fraction increase along the
porosities (Figs. 13 and 14). It should be mentioned that the
main flow channel. The large increase in H2/H2O ratio along
electrode porosity not only influences the transport of gas
the cathode gas flow channel decreases the density and
species in the porous electrodes, but also influences the rate of
viscosity of the gas mixture, leading to higher velocity at the
electrochemical reactions [54]. The effect of electrode porosity
downstream of the cathode gas channel. For comparison, the
on activation overpotentials will be included in the CFD model
velocity distribution at the anode shows little dependence on
in future study to obtain optimized electrode microstructure.
electrochemical reaction, due to very small change in gas
composition along the anode flow channel. Increasing the
5. Conclusions SOEC operating potential increases the local current density
and the molar fractions of H2 and O2, which in turn increases
A 2D CFD model of a planar SOEC was developed. Simulation the gas velocity at the cathode downstream. At an operating
results were compared with data from FLUENT and good potential lower than the thermal neutral voltage (TNV), the
 international journal of hydrogen energy 34 (2009) 7795–7806 7805

SOEC temperature decreases along the main flow channel. For [10] Shin Y, Park W, Chang J, Park J. Evaluation of the high
comparison, the SOEC temperature increases at a potential temperature electrolysis of steam to produce hydrogen.
higher than TNV. International Journal of Hydrogen Energy 2007;32(10–11):
1486–91.
It is also found that the local current density can be enhanced
[11] Sigurvinsson J, Mansilla C, Lovera P, Werkoff F. Can high
by increasing the inlet gas velocity. Increasing the inlet gas temperature steam electrolysis function with geothermal
velocity decreases the molar fraction ratio of H2/H2O at the heat? International Journal of Hydrogen Energy 2007;32(9):
downstream of the cathode, leading to lower Nernst potential 1174–82.
and higher current density. However, the SOEC temperature is [12] Balta MT, Dincer I, Hepbasli A. Thermodynamic assessment
only slightly higher at a higher inlet gas velocity, as the cooler of geothermal energy use in hydrogen production.
International Journal of Hydrogen Energy 2009;34(7):2925–39.
gas can convect more heat from the SOEC at the downstream
[13] Iwahara H, Esaka T, Uchida H, Maeda N. Proton conduction
despite of higher heat generation from the ohmic heating.
in sintered oxides and its application to steam electrolysis
Another important finding is that the electrode porosity for hydrogen production. Solid State Ionics 1981;3–4:359–63.
shows small effect on SOEC performance. A high electrode [14] Donitz W, Erdle E. High-temperature electrolysis of water
porosity facilitates the transport of gas species in the porous vapor-status of development and perspectives for
electrodes, tending to decrease the Nernst potential and thus application. International Journal of Hydrogen Energy 1985;
increase current density. However, an increase in current 10(5):291–5.
[15] Iwahara H, Uchida H, Yamasaki I. High-temperature steam
density tends to increase the molar fractions of H2 and O2,
electrolysis using SrCeO3-based proton conductive solid
which in turn limits the increase in local current density. The
electrolyte. International Journal of Hydrogen Energy 1987;
combined effects result in small dependence of current 12(2):73–7.
density on electrode porosity. Correspondingly, the tempera- [16] Quandt KH, Streicher R. Concept and design of a 3.5 MW pilot
ture field and velocity field of the SOEC do not vary signifi- plant for high temperature electrolysis of water vapor.
cantly with a change in electrode porosity. International Journal of Hydrogen Energy 1986;11(5):309–15.
The results presented in this paper provide better under- [17] Donitz W, Dietrich G, Erdle E, Streicher R. Electrochemical
high temperature technology for hydrogen production or
standing on the working mechanisms of SOECs. The model
direct electricity generation. International Journal of
developed can serve as a useful tool for design optimization. Hydrogen Energy 1988;13(5):283–7.
[18] Arashi H, Naito H, Miura H. Hydrogen production from high-
temperature steam electrolysis using solar energy.
references International Journal of Hydrogen Energy 1991;16(9):603–8.
[19] Perfiliev MV. Problems of high-temperature electrolysis of
water vapor. International Journal of Hydrogen Energy 1994;
[1] Sherif SA, Barbir F, Veziroglu TN. Wind energy and the 19(3):227–30.
hydrogen economy – review of the technology. Solar Energy [20] Hino R, Miyamoto Y. Hydrogen production by high
2005;78(5):647–60. temperature electrolysis of water vapor steam – test results
[2] Ni M, Leung DYC, Leung MKH, Sumathy K. An overview of obtained with an electrolysis tube. Journal of the Atomic
hydrogen production from biomass. Fuel Processing Energy Society of Japan 1995;37(11):1042–9.
Technology 2006;87(5):461–72. [21] Momma A, Kato T, Kaga Y, Nagata S. Polarization behavior of
[3] Ni M, Leung MKH, Leung DYC, Sumathy K. A review and high temperature electrolysis cells (SOEC). Journal of the
recent developments in photocatalytic water-splitting using Ceramic Society of Japan 1997;105(5):369–73.
TiO2 for hydrogen production. Renewable and Sustainable [22] Hino R, Haga K, Aita H, Sekita K. R&D on hydrogen
Energy Reviews 2007;11(3):401–25. production by high temperature electrolysis of steam.
[4] Lee MS, Goswami DY, Stefanakos EK. Immobilization of Nuclear Engineering and Design 2004;233(1–3):363–75.
calcium oxide solid reactant on a yttria fabric and [23] Hong HS, Chae US, Choo ST, Lee KS. Microstructure and
thermodynamic analysis of UT-3 thermochemical hydrogen electrical conductivity of Ni/YSZ and NiO/YSZ composites for
production cycle. International Journal of Hydrogen Energy high temperature electrolysis prepared by mechanical
2009;34(2):745–52. alloying. Journal of Power Sources 2005;149:84–9.
[5] O’Brien JE, McKellar MG, Stoots CM, Herring JS, Hawkes GL. [24] Yildiz B, Kazimi MS. Efficiency of hydrogen production
Parametric study of large-scale production of syngas via systems using alternative nuclear energy technologies.
high-temperature co-electrolysis. International Journal of International Journal of Hydrogen Energy 2006;31(1):77–92.
Hydrogen Energy 2009;34(9):4216–26. [25] Jensen SH, Larsen PH, Mogensen M. Hydrogen and synthetic
[6] Stoots C, O’Brien J, Hartvigsen J. Results of recent high fuel production from renewable energy sources.
temperature coelectrolysis studies at the Idaho National International Journal of Hydrogen Energy 2007;32(15):3253–7.
Laboratory. International Journal of Hydrogen Energy 2009; [26] Hong HS, Chae US, Choo ST. The effect of ball milling
34(9):4208–15. parameters and Ni content on a YSZ-coated Ni composite for
[7] Hawkes G, O’Brien J, Stoots C, Hawkes B. 3D CFD model of a high temperature electrolysis cathode. Journal of Alloys
a multi-cell high temperature electrolysis stack. and Compounds 2008;449(1–2):331–4.
International Journal of Hydrogen Energy 2009;34(9): [27] Brisse A, Schefold J, Zahid M. High temperature water
4189–97. electrolysis in solid oxide cells. International Journal of
[8] Ni M, Leung MKH, Leung DYC. Technological development of Hydrogen Energy 2008;33(20):5375–82.
hydrogen production by solid oxide electrolyzer cell (SOEC). [28] Yu B, Zhang WQ, Xu JM, Chen J. Microstructural
International Journal of Hydrogen Energy 2008;33(9):2337–54. characterization and electrochemical properties of Ba0.5Sr0.
[9] Herring JS, O’Brien JE, Stoots CM, Hawkes GL, Hartvigsen JJ, 5Co0.8Fe0.2O3d and its application for anode of SOEC.
Shahnam M. Progress in high-temperature electrolysis for International Journal of Hydrogen Energy 2008;33(23):6873–7.
hydrogen production using planar SOFC technology. [29] Liang M, Yu B, Wen M, Chen J, Xu J, Zhai Y. Preparation of
International Journal of Hydrogen Energy 2007;32(4):440–50. LSM–YSZ composite powder for anode of solid oxide
7806 international journal of hydrogen energy 34 (2009) 7795–7806

electrolysis cell and its activation mechanism. Journal of [42] Yuan JL, Lv XR, Sunden B, Yue D. Analysis of parameter
Power Sources 2009;190(2):341–5. effects on transport phenomena in conjunction with
[30] Laguna-Bercero MA, Skinner SJ, Kilner JA. Performance of chemical reactions in ducts relevant for methane reformers.
solid oxide electrolysis cells based on Scandia stabilized International Journal of Hydrogen Energy 2007;32(16):
zirconia. Journal of Power Sources 2009;192(1):126–31. 3887–98.
[31] Hauch A, Jensen SH, Ramousse S, Mogensen M. Performance [43] Tseronis K, Kookos IK, Theodoropoulos C. Modeling mass
and durability of solid oxide electrolysis cells. Journal of the transport in solid oxide fuel cell anodes: a case for
Electrochemical Society 2006;153(9):A1741–7. a multidimensional dusty gas-based model. Chemical
[32] Hauch A, Ebbesen SD, Jensen SH, Mogensen M. Solid oxide Engineering Science 2008;63(23):5626–38.
electrolysis cells: microstructure and degradation of the Ni/ [44] Yuan JL, Rokni M, Sunden B. Three-dimensional
yttria-stabilized zirconia electrode. Journal of the computational analysis of gas and heat transport
Electrochemical Society 2008;155(11):B1184–93. phenomena in ducts relevant for anode-supported solid
[33] Ni M, Leung MKH, Leung DYC. Parametric study of solid oxide oxide fuel cells. International Journal of Heat and Mass
steam electrolyzer for hydrogen production. International Transfer 2003;46(5):809–21.
Journal of Hydrogen Energy 2007;32(13):2305–13. [45] Ji Y, Yuan K, Chung JN, Chen YC. Effects of transport scale on
[34] Ni M, Leung MKH, Leung DYC. Mathematical modeling of the heat/mass transfer and performance optimization for solid
coupled transport and electrochemical reactions in solid oxide fuel cells. Journal of Power Sources 2006;161(1):380–91.
oxide steam electrolyzer for hydrogen production. [46] Xia WS, Yang YZ, Wang QS, Effects of operations and
Electrochimica Acta 2007;52(24):6707–18. structural parameters on the one-cell stack performance of
[35] Ni M, Leung MKH, Leung DYC. Electrochemical modeling of planar solid oxide fuel cell. Journal of Power Source, in press.
hydrogen production by proton-conducting solid oxide DOI: 10.1016/[Link].2009.06.009.
steam electrolyzer. International Journal of Hydrogen Energy [47] Arpino F, Carotenuto A, Massarotti N, Nithiarasu P. A robust
2008;33(15):4040–7. model and numerical approach for solving solid oxide fuel
[36] Udagawa J, Aguiar P, Brandon NP. Hydrogen production cell (SOFC) problems. International Journal of Numerical
through steam electrolysis: model-based steady state Methods for Heat & Fluid Flow 2008;18(7–8):811–34.
performance of a cathode-supported intermediate [48] Ni M, Leung MKH, Leung DYC. Energy and exergy analysis of
temperature solid oxide electrolysis cell. Journal of Power hydrogen production by solid oxide steam electrolyzer plant.
Sources 2007;166(1):127–36. International Journal of Hydrogen Energy 2007;32(18):
[37] Udagawa J, Aguiar P, Brandon NP. Hydrogen production 4648–60.
through steam electrolysis: control strategies for a cathode- [49] Li PW, Chyu MK. Simulation of the chemical/electrochemical
supported intermediate temperature solid oxide electrolysis reactions and heat/mass transfer for a tubular SOFC in
cell. Journal of Power Sources 2008;180(1):354–64. a stack. Journal of Power Sources 2003;124(2):487–98.
[38] Chan SH, Khor KA, Xia ZT. A complete polarization model of [50] Reid RC, Prausnitz JM, Poling BE. The properties of gases &
a solid oxide fuel cell and its sensitivity to the change of cell liquids. 4th ed. New York: McGraw-Hill Book Company; 1987.
component thickness. Journal of Power Sources 2001;93(1–2): [51] Patankar SV. Numerical heat transfer and fluid flow. New
130–40. York: McGraw-Hill Book Company; 1980.
[39] Ferguson JR, Fiard JM, Herbin R. Three-dimensional [52] Yuan JL, Sunden B. Analysis of intermediate temperature
numerical simulation for various geometries of solid oxide solid oxide fuel cell transport processes and performance.
fuel cells. Journal of Power Sources 1996;58(2):109–22. ASME Journal of Heat Transfer 2005;127:1380–90.
[40] Zhu HY, Kee RJ. A general mathematical model for analyzing [53] Suzuki M, Shikazono N, Fukagata K, Kasagi N. Numerical
the performance of fuel cell membrane-electrode analysis of coupled transport and reaction phenomena in an
assemblies. Journal of Power Sources 2003;117(1–2):61–74. anode-supported flat-tube solid oxide fuel cell. Journal of
[41] Chan SH, Xia ZT. Polarization effects in electrolyte/electrode- Power Sources 2008;180(1):29–40.
supported solid oxide fuel cells. Journal of Applied [54] Chan SH, Xia ZT. Anode micro model of solid oxide fuel cell.
Electrochemistry 2002;32(3):339–47. Journal of the Electrochemical Society 2001;148(4):A388–94.