Numerical investigations for a solid oxide electrolyte cell stack

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/he

Numerical investigations for a solid oxide electrolyte cell stack Zonglei Xu, Na Ren, Mao Tang, Xiongwen Zhang*, Fengming Wang, Guojun Li Key Laboratory of Thermal-Fluid Science and Engineering of MOE, School of Energy & Power Engineering, Xi'an Jiaotong University, Xi'an, 710049, China

article info

abstract

Article history:

The performance of a solid oxide electrolyzer cell (SOEC) is highly related to the stack

Received 15 November 2018

structure and operating conditions. This work presents numerical investigations on a 10-

Received in revised form

cell SOEC stack under different operating mass flow rates. The distribution characteris-

25 March 2019

tics of the flow field, pressure field and temperature field with different mass flow rates are

Accepted 15 April 2019

compared. The cell performances of SOEC at different layers in the stack are investigated

Available online 15 May 2019

and discussed. The simulation results show that the performances of SOECs in the stack are differentiated by mass flow rate and species mass fraction of fuel channels at different

Keywords:

heights. Under the small mass flow rate condition, the Nernst voltages and operating

Numerical simulation

voltages of cells at the lower positions are slightly higher than that of cells at upper po-

Performance investigation

sitions. Under the large mass flow rate condition, the Nernst voltages and operating volt-

Solid oxide electrolyzer cell stack

ages of the cells in the middle stack is slightly higher than that of the upper and lower cells. © 2019 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved.

Introduction The increasing consumption of fossil fuels has caused serious environmental problems like air pollution and global warming. Growing demands for reducing the dependence on fossil fuels and looking for a replacement fuel [1e4]. Hydrogen is considered to be a highly attractive alternative fuel for its potential to address the environmental problems [5]. However, currently hydrogen is mainly produced from fossil fuelbased methods, which usually generate greenhouse gases as by-products. In a long term, water electrolysis powered by renewable energy or nuclear energy is a sustainable, efficient and environmental-friendly way to produce hydrogen on a large scale [6e10]. Among different types of electrolyzer cells, the high temperature solid oxide electrolyzer cell (SOEC) has

attracted growing interest for its lower electrical energy requirement compared with low temperature electrolyzer cells [11]. In particular, the systems that integrate SOEC with nuclear energy and solar energy have been investigated by several researches [12e14]. An SOEC can be viewed as the reverse operation of a solid oxide fuel cell (SOFC) [15e17]. A number of studies including experiments [18e21] and model simulations [22e30] have been performed on the performance characterization of SOEC. At the modeling aspects, early research focused on models at the cell level. Navasa et al. [22] proposed a three-dimensional computational fluid dynamics (CFD) model for an SOEC and examined the effects of different operating voltages. Ni et al. [23,24] proposed a 2-dimensional (2-D) model for a planar SOEC and conducted extensive parameter studies in terms of

* Corresponding author. E-mail address: [email protected] (X. Zhang). https://doi.org/10.1016/j.ijhydene.2019.04.141 0360-3199/© 2019 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved.

20998

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

electrode microstructure, inlet gas compositions, temperature and pressure. Xu et al. [25] conducted numerical studies on SOECs with variant flow configurations and analyzed the detailed distribution of species, temperature and electrical parameters. Menon et al. [26] established a numerical model coupled kinetics at the micro-scale level with macro-scale phenomena for SOECs to investigate the electrochemical behavior and irreversible losses. Jin et al. [27] developed a 2-D SOEC model to examin the performance sensitivity to delaminations occurred at oxygen electrode/electrolyte interface. Lay et al. [28] proposed a micro model to describe the electrochemical mechanisms in H2 and O2 electrode of an SOEC and analyzed the effects of microstructure parameters on the electrode behaviour. Nerat et al. [29] modeled an electrolyte supported SOEC and investigated the influence of anode delamination on the electrochemical performance. Recently, Zhang et al. [30] modeled the steady state behavior by extreme learning machine algorithm and analyzed the performances of SOEC under various gas compositions. Besides, modeling of SOECs for Syngas and alternative fuels production [31e34] and modeling of regenerative SOFCs [35e38] are also receiving increased attention. However, to achieve practical commercial large-scale hydrogen production, single SOECs need to be assembled in series to form a stack. The performance of SOEC stack have been explored by several researches though experimental methods [39e43]. At the stack level modeling, fewer studies have been published on SOEC stack modeling compared to SOFC stack modeling [44e48] and SOEC cell level modeling. Hawkes et al. [49] performed a 3D simulation on an SOEC stack with the FLUENT and presented the effects of the variation of operating voltages on the temperature, gas composition and Nernst voltage. Laurencin [50] et al. proposed a 2D in-housemodel that can provide the distribution of temperature, current densities, gas compositions and overpotentials to analyse the performances of an SOEC stack. The thermal equilibrium analysis of the stack and parametric studies of the irreversible losses were also provided in their work. Recently, Banerjee et al. [51] presented a hierarchical model to examin the coelectrolysis performance of an SOEC stack. However, few studies focused on the effects of operating mass flow rates on the cell performance variation in the SOEC stack. Experimental results [39] revealed that the performance differences among the cells may be an important factor causing the

degradation of the SOEC stack. But the reasons why the variation occurs are difficult to identify through experimental methods [39]. The cell performance variation and flow uniformity in the SOFC stack have been extensively investigated based on numerical studies [52e56]. Similar to the SOFC stack, the cell performance variation in an SOEC stack may be caused by the differences in the inlet gases entering the channels of each cell, which are closely related to design and structure of the manifolds and the working conditions [56,57]. On the other hand, the endothermic behavior of the electrolysis reaction may impose changes in the temperature distribution thus result in the cell performance variations [49,52]. The purpose of this investigation was to conduct 3D CFD numerical simulations of an SOEC stack to study the effects of operating mass flow rates on cell performance variation. The distribution characteristics of flow field, temperature field and pressure field in the manifolds, the flow distribution in the stack, the inlet gas composition distribution, and the operating characteristics of cells in the stack under two conditions were investigated and analyzed.

Electrochemical model Fig. 1 shows the electrochemical processes of an SOEC with an oxygen ion-conducting electrolyte. The vapor decomposes at the electrolyte/cathode interface of each cell as H2 O þ 2e /H2 þ O2

(R1)

The generated hydrogen is diffused through the cathodes into channels. The generated oxygen ions are conducted through the electrolytes to the anodes and then oxidized to form oxygen: O2  2e /0:5O2

(R2)

The generated oxygen diffuses into the gas channels while the electrons are conducted to the cathodes through the external circuit. The overall electrochemical reaction can be expressed as: H2 O / H2 þ 0:5O2

(R3)

The Nernst voltage of the electrochemical reaction takes the form:

Fig. 1 e Electrochemical processes of an oxide-ion conducting SOEC.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

Ev ¼ En þ hact þ hconc þ hohm

Table 1 e Coefficients of specific heat capacity for the gases [58]. a1 =J mol1 $K1 H 2O O2 H2

a2  103 =Jmol1 $K2

a3  106 =Jmol1 $K3

9.6212 12.9874 0.8373

1.1848 3.8644 2.0138

30.3794 25.8911 29.0856

hact ¼

  2RT i 1 sinh ne F 2i0

where     gi T; pi ¼ hi ðTÞ  Tsi T; pi ¼

ZT cpi ðTÞdT  T

T0

cpi ðTÞ dT T

T0

 RT ln

p0 pi

(2)

cpi in Eq. (2) is given by: cpi ðTÞ ¼ a1 þ a2 T þ a3 T2

(5)

where the exchange current density i0 is modeled as [60]:

(1)

ZT

(4)

The activation polarization can be expressed as [59,60]:

i0 ¼ gc  1
gH2 þ 0:5gO2  gH2 O En ¼ 2F

20999

(3)

The constant parameters a1, a2 and a3 are provided in Table 1. The voltage of each cell between the interconnectors can be calculated as:



pH2 pref 

!

! pH2 O exp pref

 Ec at the cathode=electrolyte interfaces RT

!0:25 pO2 i0 ¼ ga exp pref   Ea at the anode=electrolyte interfaces  RT

(6)

(7)

where Ec ¼ 1:0  105 J=mol;Ea ¼ 1:2  105 J=mol;gc ¼ 7:0  107 A= m2 ; ga ¼ 11:2  108 A=m2 [25], which are evaluated by experiments. The ohmic losses hohm includes losses due to the ionic resistance of the electrolyte, the electronic resistances of both electrodes and the contact resistances between the electrodes and the interconnectors. It can be expressed as:

Fig. 2 e The depiction of the overall SOEC stack model.

Fig. 3 e The superposition of two SOECs with an interconnector.

21000

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

Fig. 4 e Grid quality for the whole model (a) Mesh used at each cross section in Y-direction, (b) Corner view of the mesh for a cell.

21001

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

given in Section Numerical methods. The contact resistance Rcontact is due to the current collection between electrodes and interconnectors. The value of this parameter was empirically specified and adjusted by fitting the numerical results to the experimental data. In this model, it was specified as 0.03 U cm2 [25,62,63]. The concentration polarizationhconc accounts for the difference between the actual gas partial pressures at reaction site and the ones in the electrode surfaces. The electrochemical reactions occur at the three phase boundary. In this study, the gas partial pressures at the electrolyte-electrode interface are used for calculating the Gibbs free energy. Thus, the concentration polarization term in Eq. (4) can be ignored [25].

Table 2 e Structure parameters for the planar SOEC stack [25,49]. Parameter

value

Stack Stack length Stack width Stack height Manifolds Manifold length Manifold width Manifold height Manifold wall thickness Manifold inlet/outlet diameter Cells Number of Cells Number of channels Gas flow configuration Active area Anode thickness Cathode thickness Electrolyte thickness Cathode-support thickness Interconnector thickness Gas channel width Gas channel height Rib width Pore diameter Grain diameter Tortuosity Porosity

60 mm 60 mm 25 mm 40 mm 10 mm 25 mm 0.5 mm 6 mm 10 13  13 cross flow 16 cm2 25 mm 10 mm 10 mm 400 mm 1 mm 1.5 mm 0.5 mm 2 mm 3.0 mm 1.5 mm 5 0.4

!    n hohm ¼ hohm;el þ hohm;a þ hohm;c þ Rcontact  i ,!

Numerical methods The conservation equations of mass, momentum, energy, species, and electric potential are summarized below [25] . (11)

    vp div εr! n nj ¼ div εm grad nj þ Sj  ε vwj

(12)

  divðr! n hÞ ¼ div leff grad T þ Sq

(13)

   ! div εrYi V ¼ div rDi;eff grad Yi þ Si

(14)

divðs grad FÞ ¼ 0

(15)

(8)

The ohmic losses through the electrolytes hohm;el can be calculated as ! .   n s hohm;el ¼ del  i ,!

(9)

wheredel denotes the electrolyte thickness. The conductivity s can be modeled as follows: s¼

 ! div εr V ¼ Sm

  A0 E0 exp  T T

(10)

where A0 and E0 for anodes, cathodes, electrolytes and interconnectors are constant and referred to References [61]. The ohmic losses in the electrodes and interconnectors are computed by solving the electric potential equation, which is

where ε is the porosity, h signifies the gas enthalpy, leff indicates the effective thermal conductivity; Di;eff denotes the effective diffusion coefficient of species i, which can be found in the previous publications [25]. The source terms of these equations are summarized as [25]: m Sj ¼  Vj for the porous electrodes and a Sj ¼ 0 for gas channels

(16)

2

Sq ¼ sðgrad FÞ for the electric zones

(17)

Table 3 e Operating conditions and material properties for the planar SOEC stack [15,16]. Operating conditions Inlet gas temperature (K)

Pressure (bar)

Inlet gas compositions

1073

3

Air: O2 21%, N2 79%

Case A Case B

Operating voltage(v)

Fuel:H220%,H2O 80%

s(S m1) 1

1

l(Wm K )

Anode

Electrolyte

  9:5  107 1150 [16] exp  T T



10300 3:34  104 exp  T

4

2

11

 [16]

Inlet fuel flow rate (Kg/s)

Inlet air flow rate (Kg/s)

1.67  105

3.38  105

1.67  104

3.38  104

Cathode

Interconnector

  9:7  104 2100 [15] exp  T T

7.768eþ5

4

30

21002

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

Fig. 5 e Flow field of H2/H2O inlet manifold and air inlet manifold in two cases (a) case A, (b) case B.

Sq ¼ ACV



!   !  i ,n 2F

!

i¼

T,DSE þ hohm i þ hact i

for grid cells adjacent to the electrode=electrolyte interfaces (18)

SO2 ¼ 0:5MO2 ACV

SH2 ¼MH2 ACV

!   !  i ,n 2F

!   !  i ,n 2F

for the anode

and SH2 O ¼ MH2 O ACV

(19) !   !  i ,n 2F

for the cathode (20)

wherea is the permeability, which can be found in the previous publications [25,64]; DSE denotes the change of entropy in the reaction, ACV signifies the area per unit volume of the grid cells adjoining the electrode/electrolyte interfaces [25]; Mi expresses the molecular weight for species i. The source term Sj of the momentum conservation equation (Eq. (12)) for the porous electrodes is modeled by Darcy's Law as Eq. (16). The source terms Sq for energy equations in the electric zones are due to the ohmic losses. The source terms Si for species are set to be zero in the gas channels. The boundary conditions of electrolytes were treated as wall in this study. The current density of the electrolyte is written as:

ðFa  Fc Þ  En  hact rel

(21)

All external surfaces were regarded as adiabatic. The temperature, the mass flow rates and the gas compositions were assigned at the manifold inlets. The pressure was specified at the manifold outlets. Uniform potential was assumed and defined at the upper and lower terminals of the stack. It was assumed that the SOECs in the stack should produce the same amount of current density as they are connected in series. This requirement was satisfied with an iterative adjustment of the cell voltage that was set as the potential boundary condition at the current inlet and outlet of each cell. The numerical simulations were conducted by using the CFD software FLUENT 14.5, where the SIMPLEC algorithm and the laminar model were used. The electrochemical model and source terms were programmed and provided by user-defined functions (UDFs).

Results and discussions Model geometry and numerical grids The SOEC stack model investigated in this work was developed based on the geometry of Ref. [45], as shown in Fig. 2. The

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

21003

Fig. 6 e Pressure field at the inlet section of H2/H2O channels and air channels in two cases (a) case A, (b) case B.

overall SOEC stack model comprises two major parts: a series of cells and the external manifolds. The H2O/H2 gas mixture and air enter the gas-in manifolds though inlets on the top and then flow into the cells at different heights. After the reaction, the gas flows out from the cell channels and enters the gas-out manifolds, and then eventually exit the manifolds through the outlets on the top. Fig. 3 shows a part of the stack that is the superposition of two SOECs with an interconnector in the middle for connecting different cells. The gas channels are arranged in both sides of the interconnectors with a cross flow design. The computational domain consists of the anodes, the cathodes, the gas flow channels, the inter-connectors, and the manifolds. The electrolytes and the wall surfaces of the manifols are regarded as solid walls and the interconnectors are solid zones. The porous zones comprise the anodes and the cathodes. Approximately 3.6 million mesh elements are included in the whole computational domain. Fig. 4(a) shows the top view of the grid used for the cells and manifolds. This grid was used at each level throughout the entire height of the whole model. There are 108  108 elements for the active area in the X and Z directions. The gas-in manifolds are portrayed by finer grids to accurately describe the flow field. Fig. 4(b) is a corner view of the grid for a cell. There are six layers stacked in the Y-direction: the interconnector (the top layer), the air channel (embeded in the connector), the anode (the third

layer), the cathode (the fourth layer), the fuel channel (embeded in the connector) and the other interconnector (the bottom layer), with five, three, two, four, three and five elements, respectively. Two grid density tests were performed. In the first study, the number of grids in all three directions of a single cell was doubled, and the results obtained such as voltages, temperatures and species compositions were identical to the present model. The second test was performed to verify the grids of the manifolds are sufficiently fine. In this test, the number of the grids of the gas-in manifolds and the gas-out manifolds was doubled, and the results obtained were also identical to the present model.

Materials and operating conditions Table 2 provides the structural parameters of the stack. The materials for electrolyte, cathode and anode are the Y2O3stabilized ZrO2 (YSZ), the Ni-stabilized ZrO2 (YSZ), and the Sr-doped LaMnO3 (LSM), respectively. Stainless steel is the material of interconnectors. The operating conditions and material properties for the SOEC stack are provided in Table 3. It was assumed that the mass flow rate may be the crucial factor influencing the flow field in the manifolds, thus two cases with different mass flow rates were considered in this study. It is suggested that the case A and case B can represent the small mass flow rate condition and the large mass flow

21004

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

Fig. 7 e Temperature field at the inlet section of H2/H2O channels and air channels in two cases (a) case A, (b) case B.

rate condition, respectively. In order to examin the influence of the endothermic behavior on the temperature distribution, an operating voltage 11v that is lower than the thermal neutral voltage was applied.

Model validation A primitive model verification has been conducted in Ref. [25]. The simulation results of a one-cell stack with an active area of 89 mm  89 mm were compared with experimental data. The measured cell consists of the YSZ electrolyte, the Ni/YSZ cathode and the LSM anode. The one-cell stack was operated at 1073 K and supplied with H2O/H2 gas mixture at 0.217 L/min and air at 0.33 L/min. The simulation results based on the present model agree well with the experimental measurements. The maximum absolute error with the model is approximately 309 A/m2 [25]. Because the steam electrolysis reaction is endothermic and accompanied with irreversible heat generation, the net heat generation in an SOEC may be negative, zero, or positive depending on the operating voltage. At the thermal neutral voltage, the net heat generation is zero and the temperature of the cell might be equal to the operating temperature. The thermal neutral voltage can be calculated though the application of energy balance principle. It is 1.2867 V at the

operating temperature of 1073 K [63]. The model was also validated with thermal neutral operations. The model simulations at the thermal neutral voltage 1.2867 V were conducted to show the deviations of the simulated temperature distribution from the theoretical value of 1073 K. The maximum relative error and absolute error between the model simulation and theoretical calculation are 1.52% and 16.3 K [25].

Results of the manifolds Fig. 5 shows the flow field of H2/H2O inlet manifold and air inlet manifold in two cases. The flow field of inlet manifold has a critical influence on the distribution of pressure field and temperature field in the manifold, thus influence the flow distribution uniformity of the stack. Moreover, this may be an important factor relating to the cell performance variation in the stack. In case A, the flow field in the H2/H2O inlet manifold appears to be chaotic and the gas almost spreads throughout the whole space of the manifold, as indicated in Fig. 5(a). As the mass flow rate of air is higher, it can be seen from Fig. 5(a) that the air impinges on the bottom after entering the air inlet manifold; then flows backward to both sides and forms two vortices. In case B, as shown in Fig. 5(b), both H2/H2O gas mixture and air flow hit the bottom at a relatively higher velocity after entering the inlet manifolds, and forms two

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

21005

Fig. 8 e The inlet mass flow distribution at H2/H2O side and air side in two cases (a) case A, (b) case B.

vortices respectively. It is noted that the presence of vortices may increase the gradient of the pressure field. Fig. 6 shows the pressure field at the inlet section of H2/H2O channels and air channels in two cases. In case A, it can be observed from Fig. 6(a) that the inlet pressures of the H2/H2O channels for different cells are very uniform. For the air side, the inlet pressures of the intermediate cells are smaller than that of the cells at both ends due to the vortices. In case B, it can be identified in Fig. 6(b) that the maximum value of pressure is concentrated at the position where the gas flow collides with the gas distributor bottom, because the flow velocity is sharply reduced to almost zero there. The pressure at the center of the vortices is lower. Because the pressure fields in the gas outlet manifolds are all very uniform and not displayed here, the pressure difference between the inlet and the oulet of a channel is nearly positively dependent on the inlet pressure. Therefore, the inlet pressure of the cells closer to the center of the vortices is smaller, leading to the inlet mass flow rates of cells in the middle stack to be smaller than that of cells at both ends in case B. Fig. 7 shows the temperature field at the inlet section of H2/ H2O channels and air channels in two cases. In case A, it can be identified from Fig. 7(a) that the mean temperature is lower than the inlet temperature 1073 K. This is because the electrolysis reaction appears to be endothermic as the operating voltage of each cell is lower than the thermal neutral voltage (1.2867 V at 1073 K). It can be noted that the temperature

decreases significantly along the gas flow direction in case A. A possible explanation for this is that the rate of heat absorption of the reaction is grater than that of heat diffusion in the manifold when the mass flow rate is small. In case B, it can be observed from Fig. 7(b) that the temperature field distribution on both the H2/H2O side and the air side are relatively uniform, because the high mass flow rate might benefit the heat diffusion. It can also be observed from Fig. 7 that the temperature gradient in the gas inlet manifold decreases with an increase in the inlet mass flow rate. The influence of the endothermic behavior of the electrolysis reaction may be more significant under the small mass flow condition.

Results of the stack Fig. 8 shows the inlet mass flow distribution at H2/H2O side and air side in two cases. The inlet mass flow distribution is mainly related to the pressure distribution in the manifold. The inlet mass flow rate of a cell increases as the inlet pressure increases. In case A, as shown in Fig. 8(a), the inlet mass flow distribution at the H2/H2O side and the air side are both very uniform due to the uniform pressure distribution. In case B, a similar distribution of the inlet mass flow at the H2/H2O side and the air side can be found in Fig. 8(b). The cells in the middle position have lower inlet mass flow rates because the pressure in the middle of the inlet manifold is smaller as indicated in Fig. 6.

21006

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

Fig. 9 e The average water mass fraction of inlet H2/H2O gas mixture and the average oxygen mass fraction of inlet air of each cell in two cases (a) Water mass fraction, (b) Oxygen mass fraction.

Fig. 9 shows the average water mass fraction of inlet H2/ H2O gas mixture and the average oxygen mass fraction of inlet air for each cell in two cases. For the H2/H2O side, when the temperature changes, the density change in hydrogen as the ideal gas will be more significant than that in steam which is superheated as the pressure is about 3 bar and the temperature is about 1000 K in the two cases. Therefore, as Fig. 9(a) shows, the average water mass fraction of each cell in case A has a larger difference due to a greater temperature gradient in the manifold. The average water mass fraction for cells at higher positions is larger in case A. Though numerically small, the difference in the average water mass fraction may be sufficient to cause a difference in the Nernst voltage. Compared to case A, the difference in the average water mass fraction for each cell in case B is not obvious. For the air side, the average oxygen mass fraction of each cell is very similar in both cases as Fig. 9(b) shows. Fig. 10 shows the average Nernst voltage of each cell in two cases. The Nernst voltage is mainly related to the mass flow rate and composition of the gas in the flow channels. In case A, the average Nernst voltage is negatively correlated to the

Fig. 10 e The average Nernst voltage of each cell in the two cases (a) case A, (b) case B. water mass fraction of inlet H2/H2O gas mixture. As Fig. 10(a) shows, the Nernst voltages of cells at the lower positions are slightly higher due to the smaller water mass fraction of the inlet flow. In case B, the inlet mass flow rate becomes the dominant factor influencing the average Nernst voltage. The average Nernst voltage decreases when the inlet mass flow rate increases. As Fig. 10 (b) shows, the average Nernst voltages of cells in the middle are slightly higher than that of cells at both ends. It also can be observed that the average Nernst voltage of the cell at the same layer in case B is smaller than that in case A due to a higher inlet mass flow rate. Fig. 11 displays the operating voltage of each cell in two cases. Because the cells at different layers are connected in series, the current density of each cell should be the same (1585.1 A/m2 for case A and 3397.4 A/m2 for case B). It can be observed from Figs. 10 and 11 that the cells with higher average Nernst voltage have higher operating voltages. In case A, as Fig. 11(a) shows, the operating voltages of cells at the lower positions are higher than that of cells at the upper positions. In case B, the operating voltages of the intermediate cells are higher than that of cells at the two sides, as Fig. 11(b) shows. In general, the distribution of operating voltage is similar to the Nernst voltage. However, it is notable that the

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

21007

analyzed. The difference of working performance between the cells at different heights in the stack was compared and discussed. Under the small mass flow rate condition (case A), the pressure gradient in the gas inlet manifold is small, and the tempertature gradient in the gas inlet manifold is larger due to the more significant effect of endothermic behavior. In this case, the performance difference between cells at different heights in the stack is mainly due to the difference in water mass fractions of inlet flows caused by the large temperature gradient. In the case A, the inlet mass flow rate at different heights is fairly uniform. The cells at lower positions has higher Nernst voltages. Under the large mass flow rate condition (case B), the pressure gradient in the gas inlet manifold is large while the tempertature gradient in the gas inlet manifold is small. The performance difference for cells at different heights in the stack is mainly due to the different inlet mass flow rates. In this case, the Nernst voltages of the middle cells is slightly higher than that of the upper and lower cells. The distribution of operating voltage is similar to the Nernst voltage.

Acknowledgement This work was supported by the National Natural Science Foundation of China(Grant No: 51737011, 51676161, and 51776172）,the State of Grid (Grant No: SGSDJN00FZQT1700446) and the National Key Research and Development Program of China (2016YFB0901900).

Nomenclature

Fig. 11 e The average Voltage of each cell in the two cases (a) case A, (b) case B.

difference in the operating voltage of each cell is not numerically significant, because the difference in the inlet mass flow rate or the water mass fraction of H2/H2O gas mixture for each cell is not large and does not cause a very significant gap in the Nernst voltage. And this uniformity is favorable to the longterm durability of the stack. It is possible that if the active area changes, the number of cells stacked changes, or the structure of the gas manifolds changes, the distribution of pressure field and temperature field in the manifolds will be different, so that the variation in mass flow rate or the water mass fraction at inlets for cells at different heights may be more significant. And then a different cell performance variation in cells might be observed .

Conclusions A numerical investigation on a planar SOEC stack has been presented and discussed. Effects of the mass flow rates are shown for this stack. The flow fields, pressure fields and temperature fields of gas inlet manifolds were presented and

A0 a1 a2 a3 ACV cp D E En E0 Ev F g h i i0 M ! n ne p p0 R Rcontact s

coefficient of conductivity (S m1 K) fitting function coefficient of specific heat capacity (J mol1 K1) fitting function coefficient of specific heat capacity (J mol1 K2) fitting function coefficient of specific heat capacity (J mol1 K3) specific surface area of control volume (m1) specific heat capacity (J mol1 K1) diffusion coefficient (m2 s1) activation energy (J mol1) Nernst voltage(V) coefficient of conductivity (K) voltage (V) Faraday constant (¼96485 A s mol1) molar Gibbs free energy (J mol1) molar enthalpy (J mol1) current density (A m2) exchange current density (A m2) molecular mass (kg mol1) unit normal vector of reaction interface electrons transferred per reaction pressure (Pa) standard pressure (Pa) universal gas constant (¼8.314 J mol1 K1) contact resistance(U cm2) molar entropy (J mol1 K1)

21008

S DSE T T0 V wj xi Yi

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

source term entropy change of the electrochemical reaction (J mol1) temperature (K) standard temperature (K) velocity (m s1) coordinate at the jth (x, y, z) direction Eq. (12) mole fraction of species i mass fraction of species i

Greek letters
a permeability of porous material (m2) ~ a kinetic parameters in Eqs. (6) and (7) (A m2) € a thickness z inertial resistance factor (m1) c¸ polarization (V) € e thermal conductivity (W m1 K1) ı` fluid viscosity (Pa s) ~ n fluid density (kg m3)
o conductivity (S m1) ^ o tortuosity € O electric potential (V) Subscripts a anode act activation c cathode conc concentration contact contact resistance eff effective el electrolyte i index of chemical species j direction of coordinate k index of chemical species m mass ohm ohmic ref reference q quantity of heat

references

[1] Schnoor JL. Energy and global warming: the great convergence. Environ Sci Technol 2004;38:119A. [2] Hosseini SE, Wahid MA, Aghili N. The scenario of greenhouse gases reduction in Malaysia. Renew Sustain Energy Rev 2013;28:400e9. [3] Zervas PL, Sarimveis H, Palyvos JA, et al. Model-based optimal control of a hybrid power generation system consisting of photovoltaic arrays and fuel cells. J Power Sources 2008;181(2):327e38. [4] Hosseini SE, Andwari AM, Wahid MA, BagheriG. A review on green energy potentials in Iran. Renew Sustain Energy Rev 2013;27:533e45. [5] Mazloomi K, Gomes C. Hydrogen as an energy carrier: prospects and challenges. Renew Sustain Energy Rev 2012;16(5):3024e33. [6] Jensen SH, Sun X, Ebbesen SD, et al. Hydrogen and synthetic fuel production using pressurized solid oxide electrolysis cells. Int J Hydrogen Energy 2010;35(18):9544e9. [7] Hosseini SE, Wahid MA. Hydrogen production from renewable and sustainable energy resources: promising

[8]

[9]

[10]

[11]

[12]

[13]

[14] [15] [16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

green energy carrier for clean development. Renew Sustain Energy Rev 2016;57:850e66. Dincer I, Acar C. Review and evaluation of hydrogen production methods for better sustainability. Int J Hydrogen Energy 2015;40(34):11094e111. Dincer I. Technical, environmental and exergetic aspects of hydrogen energy systems. Int J Hydrogen Energy 2002;27(3):265e85. Yildiz B, Kazimi M. Efficiency of hydrogen production systems using alternative nuclear energy technologies. Int J Hydrogen Energy 2006;31(1):77e92. Ni M, Leung MKH, Leung DYC. Technological development of hydrogen production by solid oxide electrolyzer cell (SOEC). Int J Hydrogen Energy 2008;33(9):2337e54. Kim JS, Boardman RD, Bragg-Sitton SM. Dynamic performance analysis of a high-temperature steam electrolysis plant integrated within nuclear-renewable hybrid energy systems. Appl Energy 2018;228:2090e110. Balta MT, Kizilkan O, Yılmaz Fatih. Energy and exergy analyses of integrated hydrogen production system using high temperature steam electrolysis. Int J Hydrogen Energy 2016;41(19):8032e41. Dincer I. Green methods for hydrogen production. Int J Hydrogen Energy 2012;37(2):1954e71. EG&G Technical Services, Inc. Fuel cell handbook. 2004. Ferguson JR, Fiard JM, Herbin R. Three-dimensional numerical simulation for various geometries of solid oxide fuel cells. J Power Sources 1996;58(2):109e22. Wang Y, Yu J, Weng S. Numerical investigation of different loads effect on the performance of planar electrode supported SOFC with syngas as fuel. Int J Hydrogen Energy 2011;36(9):5624e31. Zhang Xiaoyu, O'Brien James E, Tao Greg, Zhou Can, Gregory K. Housley. Experimental design, operation, and results of a 4 kW high temperature steam electrolysis experiment. J Power Sources 2015;297:90e7. Stoots Carl M, O'Brien James E, Condie Keith G, Hartvigsen Joseph J. High-temperature electrolysis for largescale hydrogen production from nuclear energyeexperimental investigations. Int J Hydrogen Energy 2010;35(10):4861e70. Schefold J, Brisse A, Poepke H. 23,000h steam electrolysis with an electrolyte supported solid oxide cell. Int J Hydrogen Energy 2017;42(19):13415e26. Wang Z, Mori M, Araki T. Steam electrolysis performance of intermediate-temperature solid oxide electrolysis cell and efficiency of hydrogen production system at 300 Nm3 h1. Int J Hydrogen Energy 2010;35(10):4451e8.
n B. Computational fluid dynamics Navasa M, Yuan J, Sunde approach for performance evaluation of a solid oxide electrolysis cell for hydrogen production. Appl Energy 2015;137(137):867e76. Ni M, Leung MKH, Leung DYC. Parametric study of solid oxide steam electrolyzer for hydrogen production. Int J Hydrogen Energy 2007;32(13):2305e13. Ni M. Computational fluid dynamics modeling of a solid oxide electrolyzer cell for hydrogen production. Int J Hydrogen Energy 2009;34(18):7795e806. Xu Z, Zhang X, Li G, et al. Comparative performance investigation of different gas flow configurations for a planar solid oxide electrolyzer cell. Int J Hydrogen Energy 2017;42(16). Menon V, Janardhanan VM, Deutschmann O. A mathematical model to analyze solid oxide electrolyzer cells (SOECs) for hydrogen production. Chem Eng Sci 2014;110(110):83e93. Jin X, Xue X. Computational fluid dynamics analysis of solid oxide electrolysis cells with delaminations. Int J Hydrogen Energy 2010;35(14):7321e8.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 0 9 9 7 e2 1 0 0 9

[28] Lay-Grindler E, Laurencin J, Delette G, et al. Micro modelling of solid oxide electrolysis cell: from performance to durability. Int J Hydrogen Energy 2013;38(17):6917e29. ic
. Modelling of anode delamination in [29] Nerat M, ÐaniJuric solid oxide electrolysis cell and analysis of its effects on electrochemical performance. Int J Hydrogen Energy 2018. S0360319918306955. [30] Zhang C, Liu Q, Wu Q, et al. Modelling of solid oxide electrolyser cell using extreme learning machine. Electrochim Acta 2017;251:137e44. [31] Kazempoor P, Braun RJ. Hydrogen and synthetic fuel production using high temperature solid oxide electrolysis cells (SOECs). Int J Hydrogen Energy 2015;40(9):3599e612. [32] Stempien JP, Ni M, Sun Q, et al. Production of sustainable methane from renewable energy and captured carbon dioxide with the use of Solid Oxide Electrolyzer: a thermodynamic assessment. Energy 2015;82:714e21. [33] Ebbesen SD, Graves C, Mogensen M. Production of synthetic fuels by Co-electrolysis of steam and carbon dioxide. Int J Green Energy 2009;6(6):646e60. [34] Chen B, Xu H, Ni M. Modelling of SOEC-FT reactor: pressure effects on methanation process. Appl Energy 2017;185:814e24. [35] Kazempoor P, Braun RJ. Model validation and performance analysis of regenerative solid oxide cells for energy storage applications: reversible operation. Int J Hydrogen Energy 2014;39(11):5955e71. [36] Ferrero D, Lanzini A, Leone P, et al. Reversible operation of solid oxide cells under electrolysis and fuel cell modes: experimental study and model validation. Chem Eng J 2015;274:143e55. [37] Wendel CH, Gao Z, Barnett SA, et al. Modeling and experimental performance of an intermediate temperature reversible solid oxide cell for high-efficiency, distributedscale electrical energy storage. J Power Sources 2015;283:329e42. [38] Jin X, Xue X. Mathematical modeling analysis of regenerative solid oxide fuel cells in switching mode conditions. J Power Sources 2010;195(19):6652e8. [39] Zhang Xiaoyu, O'Brien, et al. Improved durability of SOEC stacks for high temperature electrolysis. Int J Hydrogen Energy 2013;38(1):20e8. [40] Zheng Y, Li Q, Chen T, et al. Comparison of performance and degradation of large-scale solid oxide electrolysis cells in stack with different composite air electrodes. Int J Hydrogen Energy 2015;40(6):2460e72. [41] Jin L, Guan W, Ma X, et al. Quantitative contribution of resistance sources of components to stack performance for planar solid oxide fuel cells. J Power Sources 2014;253(253):305e14. [42] Reytier M, Iorio SD, Chatroux A, et al. Stack performances in high temperature steam electrolysis and co-electrolysis. Int J Hydrogen Energy 2015;40(35):11370e7. [43] Barelli L, Bidini G, Cinti G. Air variation in SOE: stack experimental study. Int J Hydrogen Energy 2018;43(26):11655e62. [44] Recknagle KP, Williford RE, Chick LA, et al. Threedimensional thermo-fluid electrochemical modeling of planar SOFC stacks. J Power Sources 2003;113(1):109e14. [45] Goodman SB, Wright TM. A 3D CFD model for predicting the temperature distribution in a full scale APU SOFC short stack

[46]

[47]

[48]

[49]

[50]

[51]

[52]

[53]

[54]

[55]

[56]

[57]

[58] [59]

[60]

[61]

[62] [63]

[64]

21009

under transient operating conditions. Appl Energy 2014;135(12):539e47. Nishida RT, Beale SB, Pharoah JG. Comprehensive computational fluid dynamics model of solid oxide fuel cell stacks. Int J Hydrogen Energy 2016;41(45):20592e605. Li A, Song C, Lin Z. A multiphysics fully coupled modeling tool for the design and operation analysis of planar solid oxide fuel cell stacks. Appl Energy 2017;190:1234e44. Nakajo A, Mueller F, Brouwer J, et al. Mechanical reliability and durability of SOFC stacks. Part I : modelling of the effect of operating conditions and design alternatives on the reliability. Int J Hydrogen Energy 2012;37(11):9249e68. Hawkes G, O'Brien J, Stoots C, et al. 3D CFD model of a multicell high-temperature electrolysis stack. Int J Hydrogen Energy 2009;34(9):4189e97. Laurencin J, Kane D, Delette G, et al. Modelling of solid oxide steam electrolyser: impact of the operating conditions on hydrogen production. J Power Sources 2011;196(4):2080e93. Banerjee A, Wang Y, Diercks J, et al. Hierarchical modeling of solid oxide cells and stacks producing syngas via H2O/CO2 Co-electrolysis for industrial applications. Appl Energy 2018;230:996e1013. Lin B, Shi Y, Cai N. Numerical simulation of cell-to-cell performance variation within a syngas-fuelled planar solid oxide fuel cell stack. Appl Therm Eng 2017;114:653e62. Burt AC, Celik IB, Gemmen RS, et al. Cell-to-cell variations with increasing SOFC stack size[C]//ASME 2004. In: 2nd international conference on fuel cell science, engineering and technology. American Society of Mechanical Engineers; 2004. p. 31e8. Boersma RJ, Sammes NM. Computational analysis of the gasflow distribution in solid oxide fuel cell stacks. J Power Sources 1996;63(2):215e9. Zhao C, Yang J, Zhang T, et al. Numerical simulation of flow distribution for external manifold design in solid oxide fuel cell stack. Int J Hydrogen Energy 2017;42(10):7003e13. Bi W, Chen D, Lin Z. A key geometric parameter for the flow uniformity in planar solid oxide fuel cell stacks. Int J Hydrogen Energy 2009;34(9):3873e84. Chen D, Zeng Q, Su S, et al. Geometric optimization of a 10cell modular planar solid oxide fuel cell stack manifold. Appl Energy 2013;112:1100e7. Liu GY, Liu ZG, Ying JM, He YL. Engineering thermodynamics. Beijing: Higher Education Press; 1998 (In Chinese). Zhang X, Li G, Li J, et al. Numerical study on electric characteristics of solid oxide fuel cells. Energy Convers Manag 2007;48(3):977e89. Chan SH, Khor KA, Xia ZT. A complete polarization model of a solid oxide fuel cell and its sensitivity to the change of cell component thickness. J Power Sources 2001;93(1e2):130e40. Costamagna P, Selimovic A, Del Borghi M, et al. Electrochemical model of the integrated planar solid oxide fuel cell (IP-SOFC). Chem Eng J 2004;102(1):61e9. Fergus JW. Metallic interconnects for solid oxide fuel cells. Mater Sci Eng, A 2005;397(1e2):271e83. Hawkes GL, Brien JEO, Stoots CM, et al. CFD model of a planar solid oxide electrolysis cell for hydrogen production from nuclear energy. Nucl Technol 2005;158(2):132e44. Xu P, Yu B. Developing a new form of permeability and KozenyeCarman constant for homogeneous porous media by means of fractal geometry. Adv Water Resour 2008;31(1):74e81.