Analysis of heat and mass transport characteristics in anode-supported solid oxide fuel cells at various operating conditions

Abstract

In this study, the heat and mass transport characteristics of anode-supported solid oxide fuel cells (SOFCs) are numerically investigated by using a two-dimensional model. The mathematical model is validated by comparing the numerical results with experimental data found in the open literature. The species and temperature distributions of SOFCs at different cell voltages are presented and compared. Effects of operating temperature, flow direction arrangement, and flow velocity on the overall cell performance and local temperature distribution are also analyzed. It is concluded that the local temperature is increased with decreasing operating cell voltage, increasing operating temperature, and decreasing cathode flow velocity. The temperature distribution is significantly changed when counter-flow arrangement is used instead of coflow arrangement. In addition, the effect of anode flow velocity on temperature distribution is negligible.

1. Introduction

With the rapid development of global economy, the energy demand has been gradually increasing during the past decades. Fossil fuels including coal, petroleum, and natural gas are still the major energy sources around the world. Due to the environment and energy issues, most societies are trying to use clean and sustainable energy instead of fossil fuels. Fuel cells are regarded as one of the most promising alternative energy technologies [1]. According to the electrolyte, fuel cells can be classified into different types [2]. Proton exchange membrane fuel cells (PEMFCs) are typical low-temperature fuel cells [3–9]. And solid oxide fuel cells (SOFCs) are typical high-temperature fuel cells. SOFCs have attracted wide attention due to their advantages compared to other types of fuel cells. Tubular and planar SOFCs are two typical configurations of such cells. For the planar SOFCs, there are three types of cells including anode-supported cells, electrolyte-supported cells, and cathode-supported cells.

Cell performance of anode-supported and electrolyte-supported SOFCs under varying operating temperature and pressure was experimentally investigated by Hsieh et al. [10]. The cell performance was significantly improved with increasing operating temperature and pressure. A spiral-like interconnector was proposed and applied to improve the mass transfer inside SOFCs [11]. Experiments were carried out to evaluate the mass transfer performance. Compared to experimental measurements, numerical simulations can provide detailed information of transport phenomena inside fuel cells. Numerous mathematical models have been developed and applied to improve the understanding of fuel cells. A mathematical model was developed to predict the cell performance and corresponding potential losses [12]. It was reported that the potential losses were mainly attributed to the Ohmic loss and this needed to be minimized to improve the cell performance. The effects of operating and design conditions on cell performance were systematically investigated. Ni et al. [13] performed numerical analyses on the cell performance of SOFCs with functionally graded electrodes, where nonuniform porosity and particle size were linearly distributed. Similarly, nonlinear porosity distribution has been applied in the electrodes [14]. The cell performance of SOFCs with parabolic, linear, and inverse parabolic porosity distributions was compared with that of SOFCs with uniform porosity distribution. It was found that the cell performance was strongly influenced by the porous microstructure distributions. The local transport characteristics of SOFCs with square and rectangular channels were numerically studied and compared by Andersson et al. [15]. Recently, the effect of flow channel geometry on the cell performance of SOFCs was evaluated by Khazaee et al. [16]. Numerical results showed that the best cell performance was provided by the rectangular channel, followed by the triangular channel and the trapezoid channel. Different flow field designs can also be found in the papers presented by Danilov and Tade [17], Canavar and Timurkutluk [18], Saied et al. [19], and Shen et al. [20]. Due to the fuel flexibility, there are also many investigations related to internal reforming reaction occurring in SOFCs [21–23].

The lifetime of SOFCs is significantly affected by the thermal stress, which is caused by the prominent temperature gradient. The thermal stress distribution within SOFCs has already been extensively investigated by different researchers [24–26]. In addition, there are also some investigations focusing on the temperature distribution within SOFCs. The species and temperature distributions of SOFCs were predicted by a mathematical model developed by Mahcene et al. [27]. The effect of several heat sources was systematically investigated. The temperature distributions of SOFCs were numerically simulated by Yahya et al. [28] using the lattice Boltzmann method (LBM). The effects of thermal conductivity, cell configuration, operating temperature, porosity, and current density on temperature distribution were comprehensively discussed and analyzed. The heat and mass transfer characteristics of SOFCs with varying electrolyte thicknesses and operating temperatures were investigated and discussed by Park et al. [29]. Similarly, the effect of fuel utilization on heat and mass transfer within SOFCs was investigated by Lee et al. [30].

Although there are various numerical investigations on the cell performance and transport processes inside SOFCs, a comprehensive study on the temperature distribution is still limited. The existing papers mainly focus on the chemical reactions occurring in the electrodes. The aim of this study was to study the heat and mass transport characteristics in anode-supported SOFCs at various operating conditions. In the present study, a two-dimensional, nonisothermal, and steady-state model based on the finite volume method (FVM) was developed and used to investigate the effects of operating voltage, flow direction arrangement, operating temperature, and flow velocity on the heat transfer characteristics of anode-supported planar SOFCs.

2. Model description

2.1. Physical model and assumptions

In the present study, a two-dimensional mathematical model was developed for anode-supported planar solid oxide fuel cells. As shown in Figure 1, the computational domain consists of several layers, i.e., the anode/cathode interconnect, the anode/cathode flow channel, the anode/cathode diffusion layer, the anode/cathode functional layer, and the electrolyte layer between the anode and cathode electrodes. The geometric parameters of the physical model and the operating conditions are summarized in Table 1.

Figure 1. Schematic of anode-supported solid oxide fuel cells with coflow arrangement.

Table 1. Fuel cell geometric parameters and operating conditions [31].

Parameter	Value	Unit
Fuel cell length	40	mm
Interconnect thickness	0.3	mm
Flow channel thickness	1	mm
Anode diffusion layer thickness	0.38	mm
Anode functional layer thickness	0.02	mm
Electrolyte thickness	0.01	mm
Cathode functional layer thickness	0.02	mm
Cathode diffusion layer thickness	0.05	mm
Operating pressure, Pa/Pc	1.0	atm
Operating temperature, Ta/Tc	1073	K

The assumptions made in the mathematical model are as follows: steady state; laminar flow; incompressible gas; isotropic and homogeneous electrodes; neglected radiation heat transfer. The conservation equations of mass, momentum, species, energy, and charge are solved to describe the transport processes in fuel cells. The transport equations are connected and coupled with each other through the source/sink terms.

2.2. Governing equations

Mass conservation equation: (1) $$\nabla \cdot (\rho \overrightarrow{u})=S_{\text{mass }}$$(1) where ρ and $\overrightarrow{u}\mathrm{ }$are the mixture fluid density and velocity, respectively. S_mass is determined by the mass consumption or generation caused by the electrochemical reaction. (2) $$S_{\mathrm{mass}}=S_{H_2}+S_{H_{2}O}\mathrm{ }\left(\text{anode}\right)$$(2) (3) $$S_{\mathrm{mass}}=S_{O_2}\left(\text{cathode}\right)$$(3)

Momentum conservation equation: (4) $$\nabla \cdot \left(\rho \overrightarrow{u}\overrightarrow{u}\right)=\nabla ·\left(\mu \nabla \overrightarrow{u}\right)-\nabla P+S_{\mathrm{mom}}$$(4) where P and μ are the mixture pressure and dynamic viscosity, respectively. S_mom is the source term of the momentum equation. In the porous regions, the Darcy’s law is adopted. (5) $$S_{\mathrm{mom}}=-\frac{\mu }{K}\overrightarrow{u}$$(5) where K is the permeability.

Species conservation equation: (6) $$\nabla \cdot \left(\rho \overrightarrow{u}{Y}_i\right)=\nabla ·\left(\rho D_{\mathit{eff},i}\nabla Y_i\right)+S_i$$(6) where Y_i and D_eff,i are the mass fraction and effective diffusivity for the ith species, respectively. S_i is the consumption or generation amount of the ith (H₂, O₂, H₂O) species. (7) $$S_{H_2}\mathrm{ }=-\frac{j_a}{2F}{M}_{H_2}$$(7) (8) $$S_{O_2}\mathrm{ }=-\frac{j_c}{4F}{M}_{O_2}$$(8) (9) $$S_{H_{2}O}\mathrm{ }=\frac{j_a}{2F}{M}_{H_{2}O}$$(9)

Energy conservation equation: (10) $$\nabla \cdot \left(\rho c_p\overrightarrow{u}T\right)=\nabla ·\left(k_{\mathit{eff}}\nabla \mathrm{T}\right)+S_T$$(10) where c_p and k_eff are the specific heat and effective thermal conductivity, respectively. The irreversible, reversible, and Ohmic heat generation terms are included in the source term of the energy equation. (11) $$S_T=j_{a,c}\left|{\eta }_{a,c}\right|+j_{a,c}\frac{T\mathrm{\Delta }S}{\mathit{nF}}+{\sigma }_{\mathit{eff},i}{‖\nabla {\mathrm{\varphi }}_i‖}^2+{\sigma }_{\mathit{eff},s}{‖\nabla {\mathrm{\varphi }}_s‖}^2$$(11)

Charge conservation equation: (12) $$\nabla \cdot \left({\sigma }_{\mathit{eff},e}\nabla {\varphi }_e\right)+S_e=0$$(12) (13) $$\nabla \cdot \left({\sigma }_{\mathit{eff},i}\nabla {\varphi }_i\right)+S_i=0$$(13) where σ_eff,e is the effective electronic conductivity, and σ_eff,i is the effective ionic conductivity. The source terms are described by the Butler–Volmer equation. (14) $$S_e\mathrm{ }=-j_a\mathrm{ }{S}_i\mathrm{ }=+j_a\mathrm{ }\left(\mathrm{anode}\right)$$(14) (15) $$S_e\mathrm{ }=+j_c\mathrm{ }{S}_i\mathrm{ }=-j_c\mathrm{ }\left(\mathrm{cathode}\right)$$(15)

The corresponding expressions and parameters used in the mathematical model are summarized and presented in Tables 2 and 3, respectively. More detailed information can be found in [29].

Table 2. Complementary expressions.

Description
Effective mass diffusivity	$D_{eff,i}\mathrm{ }=\frac{\epsilon }{\tau }\times \frac{D_{i,m}\times D_{i,m}}{D_{i,m}+D_{i,k}}$
Tortuosity	$\tau \mathrm{ }={\epsilon }^{-0.5}$
Bulk diffusivity	$D_{i,m}=D_0\left(\frac{P_0}{P}\right){\left(\frac{T}{T_0}\right)}^{1.5}$
Knudsen diffusivity	$D_{i,k}=\frac{d_p}{3}\times \sqrt{\frac{8RT}{\pi M_i}}$
Effective thermal conductivity	$k_{eff}=\left(1-\epsilon \right)k_s+\epsilon k_f$
Effective electronic conductivity	${{\sigma }_{eff,e}=\sigma }_e(1-\epsilon )V_e$
Effective ionic conductivity	${{\sigma }_{eff,i}=\sigma }_i(1-\epsilon )V_i$
Anode electronic conductivity	${\sigma }_{a,e}=\frac{9.5\times {10}^7}{T}exp(-\frac{1150}{T})$
Cathode electronic conductivity	${\sigma }_{c,e}=\frac{4.2\times {10}^7}{T}exp(-\frac{1200}{T})$
Ionic conductivity	${\sigma }_i=3.34\times {10}^{4}exp(-\frac{10300}{T})$
Anode volumetric current density	$j_a\mathrm{ }={i_a^{ref}a}_{eff}{\left(\frac{C_{H_2}}{C_{H_2}^{\mathit{ref}}}\right)}^{{\gamma }_a}\left[e^{{\alpha }_{a}F{\eta }_a/RT}-e^{{-\alpha }_{c}F{\eta }_a/RT}\right]$
Cathode volumetric current density	$j_c\mathrm{ }={i_c^{ref}a}_{eff}{\left(\frac{C_{O_2}}{C_{O_2}^{\mathit{ref}}}\right)}^{{\gamma }_c}\left[-e^{{\alpha }_{a}F{\eta }_c/RT}+e^{{-\alpha }_{c}F{\eta }_c/RT}\right]$
Anode over-potential	${\eta }_a={\mathrm{\varphi }}_e-{\mathrm{\varphi }}_i$
Cathode over-potential	${\eta }_c={\mathrm{\varphi }}_e-{\mathrm{\varphi }}_i-U_o$

Table 3. Parameters used in the mathematical model.

Parameter	Value	Units
Porosity of diffusion layer	0.3	–
Porosity of functional layer	0.3	–
Volume fraction of Ni/LSM	0.4	–
Volume fraction of YSZ	0.6	–
Anode concentration dependence	0.5	–
Cathode concentration dependence	0.5	–
Reference hydrogen concentration	10.78	mol m^–3
Reference oxygen concentration	2.38	mol m^–3
Anode permeability	1.76 × 10^−11	m^2
Cathode permeability	1.76 × 10^−11	m^2
Anode thermal conductivity	11	W m^–1 K^–1
Cathode thermal conductivity	6	W m^–1 K^-1
Electrolyte thermal conductivity	2.7	W m^–1 K^–1
Interconnect thermal conductivity	20	W m^–1 K^–1

2.3. Numerical implementation and boundary conditions

The mathematical model is implemented using the commercial software ANSYS FLUENT. The governing equations of charge transport are implemented by using the user-defined scalar (UDS) equations. In addition, the user-defined functions (UDF) are used to define and calculate the source terms and the parameters in the iterative calculation process.

At the inlet of the flow channels, the velocity, temperature, and species mass fractions are prescribed. At the outlet of flow channels, a pressure-outlet boundary condition is assigned. Constant electric potentials (ϕ_e = 0 and ϕ_e = V_cell) are specified at the anode and cathode terminal walls, respectively. And the adiabatic boundary condition is specified at all the surrounding walls. And detailed boundary conditions are summarized in Table 4.

Table 4. Boundary conditions.

Description	Conditions	Value	Units
Top wall	Electrical potential	0	V
Bottom wall	Electrical potential	0.7	V
	Velocity	0.3	m/s
Fuel inlet	Mass fraction	H2:H2O = 0.95:0.05	–
	Temperature	1073	K
	Velocity	3	m/s
Air inlet	Mass fraction	O2:N2 = 0.233:0.767	–
	Temperature	1073	K
Fuel outlet	Pressure	0	Pa
Air outlet	Pressure	0	Pa

The Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm is used to link the pressure and velocity fields. The second-order upwind discretization scheme is adopted for momentum, species, energy, and UDS equations to ensure numerical stability and accuracy. Four mesh systems (x×y), namely, mesh I (40 × 34), mesh II (80 × 46), mesh III (120 × 58), and mesh IV (160 × 70), are used to investigate the mesh independence. In order to balance the accuracy and computational resources, the mesh system III has been chosen for the following simulations.

3. Results and discussion

Prior to the numerical simulations, a model validation was carried out and it is found that the numerical results agree well with the experimental data reported by Zhang et al. [31], as shown in Figure 2. The experimental operating temperature is 1073 K, and the operating pressure is 1.0 atm. Air and humidified hydrogen were introduced into the cathode and anode sides with the flow rate of 500 and 200 sccm, respectively.

Figure 2. Comparison between the numerical results and experimental data.

Figures 3 and 4 show the hydrogen and water mass fraction distributions along the cell length direction at cell voltages 0.9, 0.7, and 0.5 V, respectively. For all three cases, the hydrogen mass fraction is gradually decreased, while the water mass fraction is gradually increased along the cell length direction. This is because hydrogen is consumed and water is simultaneously produced due to the electrochemical reactions occurring in the functional layers. In addition, the oxygen mass fraction distributions along the cell length direction at three cell voltages are illustrated in Figure 5. The mass fraction is gradually decreased due to the consumption of oxygen. It is also observed that the magnitude of the mass fraction is significantly affected by the cell voltage. The corresponding current densities for three cell voltages are 0.31, 0.69, and 1.11 A/cm², respectively. Larger amounts of hydrogen and oxygen are consumed at a lower cell voltage. Accordingly, the amount of water is increased with decreasing cell voltage.

Figure 3. Hydrogen mass fraction distribution at the anode diffusion layer and functional layer interface under different operating voltages.

Figure 4. Water mass fraction distribution at the anode diffusion layer and functional layer interface under different operating voltages.

Figure 5. Oxygen mass fraction distribution at the cathode diffusion layer and functional layer interface under different operating voltages.

The temperature distributions at cell voltages 0.9, 0.7, and 0.5 V are shown in Figure 6. For all three cases, the minimum temperature appears at the cell inlet region and the maximum temperature occurs at the cell outlet region. In addition, the local temperature distribution at the anode diffusion layer and functional layer interface is presented in Figure 7. Temperature is nonlinearly increased along the cell length direction for all three cases. The temperature is increased from 1079.8 to 1110.8 K, from 1094.9 to 1204.8 K, and from 1120.9 to 1344.2 K at the cell voltage 0.9, 0.7, and 0.5 V, respectively. The temperature difference between the inlet temperature (x = 0 m) and the outlet temperature (x = 0.4 m) is greatly increased with decreasing cell voltage.

Figure 6. Temperature distribution of SOFCs: (a) operating cell voltage 0.9 V; (b) operating cell voltage 0.7 V; (c) operating cell voltage 0.5 V.

Figure 7. Temperature distribution at the anode diffusion layer and functional layer interface under different operating voltages.

In the following section, the effects of operating temperature, flow direction, anode velocity, and cathode velocity on the temperature distribution are systematically investigated. The cell operating voltage is fixed at 0.7 V in the corresponding numerical simulations. The temperature distributions of SOFCs at an operating temperature varying from 1023 to 1123 K are shown in Figure 8. It is clearly seen that both the local temperature magnitude and the temperature variation along the cell length direction are gradually increased with increasing operating temperature. The temperature is increased from 1034.0 to 1085.8 K, from 1094.9 to 1204.8 K, and from 1163.4.9 to 1360.4 K at the cell operating temperatures 1023, 1073, and 1123 K, respectively. The corresponding current densities for the three cases at the cell voltage 0.7 V are 0.37, 0.69, and 1.17 A/cm², respectively. It is concluded that the cell performance is significantly improved by the increasing operating temperature. In addition, the corresponding temperature inside the fuel cells is greatly increased.

Figure 8. Temperature distribution at the anode diffusion layer and functional layer interface under different operating temperature.

Generally, three different flow arrangements including coflow, counter-flow, and crossflow are used in SOFCs. In this study, a comparison between the coflow and counter-flow arrangements is performed to study the heat transport characteristics in SOFCs. The temperature distributions of SOFCs with coflow and counter-flow arrangements are shown in Figure 9. The anode flow channel and cathode flow channel have the same flow direction in the coflow arrangement. On the other hand, the flow direction in the cathode flow channel is opposite to the anode flow channel in the counter-flow arrangement. It is observed that a lower maximum temperature and a higher minimum temperature are obtained by the SOFC with counter-flow arrangement compared to the SOFC with coflow arrangement. The maximum temperature of the two cases is 1204.8 and 1189.9 K, respectively. The minimum temperature of coflow arrangement is 1094.9 K, while that of counter-flow arrangement is 1131.6 K. For the coflow arrangement, the temperature is monotonically increased along the cell length direction. For the counter-flow arrangement, the temperature is firstly increased and then gradually decreased along the cell length direction. The maximum temperature appears around the x = 0.15 m position, which is at the upstream of the anode electrode. The corresponding current densities for the two cases at the cell voltage 0.7 V are 0.69 and 0.70 A/cm², respectively. The overall cell performance is not affected by the flow arrangements, while the local temperature distribution is significantly changed.

Figure 9. Temperature distribution at the anode diffusion layer and functional layer interface under different flow arrangements.

Figure 10 shows the temperature distributions of SOFCs with anode flow velocity varying from 0.2 to 0.4 m/s. It can be seen that the temperature variation between the different cases is very small. This indicates that the effect of the anode flow velocity on temperature is negligible. The corresponding current densities for the three cases at the cell voltage 0.7 V are 0.68, 0.69, and 0.70 A/cm², respectively. The effect of the cathode flow velocity on the temperature distributions of SOFCs is also examined and shown in Figure 11. The velocity varies from 2 to 4 m/s with an interval of 1 m/s. It is observed that the local temperature is considerably decreased with increasing velocity, especially the maximum temperature. The corresponding current densities for the three cases at the cell voltage 0.7 V are 0.69, 0.69, and 0.70 A/cm², respectively. The corresponding maximum temperature for the three cases at the cell voltage 0.7 V are 1242.7, 1204.8, and 1180.5 K, respectively. This indicates that the local temperature can be greatly decreased with increasing air flow velocity.

Figure 10. Temperature distribution at the anode diffusion layer and functional layer interface under different anode flow velocities.

Figure 11. Temperature distribution at the anode diffusion layer and functional layer interface under different cathode flow velocities.

4. Conclusions

In this study, a two-dimensional, nonisothermal, steady-state model was applied to investigate the heat and mass transport characteristics of SOFCs. The mass, momentum, species, temperature, and charge conservation equations were solved, and the results were analyzed and discussed. The conclusions are as follows:

The temperature is gradually increased along the flow direction, and the temperature is significantly increased with decreasing cell operating voltage. The overall cell performance can be greatly improved by increasing the operating temperature, and the local temperature is also significantly increased. The flow arrangements have strong impact on the local temperature distribution, and a lower maximum temperature is provided by the counter-flow arrangement. The effect of anode flow velocity on temperature distribution can be ignored, while the temperature is gradually decreased with increasing cathode flow velocity.

Nomenclature
a	=	effective surface area, m–1
c	=	mole concentration, mol m–3
cp	=	specific heat, J kg–1 K–1
d	=	diameter, m
D	=	diffusivity, m2 s–1
F	=	Faraday constant, 96485 C mol–1
j	=	volumetric current density, A m–3
k	=	thermal conductivity, W m–1 K–1
K	=	permeability, m2
M	=	molecular weight, kg mol–1
P	=	pressure, Pa
R	=	universal gas constant, 8.314 J mol–1 K–1
S	=	source term
T	=	temperature, K
$\overrightarrow{u}$	=	velocity vector, m/s
U	=	thermodynamic equilibrium potential, V
V	=	voltage, V
Y	=	mass fraction

Greek Symbols
α	=	transfer coefficient
γ	=	concentration dependence
ε	=	porosity
η	=	over-potential, V
μ	=	dynamic viscosity, Pa s
ρ	=	density, kg m–3
σ	=	electronic/ionic conductivity, S m–1
τ	=	tortuosity
ϕ	=	potential, V

Superscripts and Subscripts
a	=	anode
c	=	cathode
e	=	electron
eff	=	effective
f	=	fluid
g	=	gas
i	=	ith species/ion
mass	=	mass equation
mom	=	momentum equation
p	=	pore
ref	=	reference
s	=	solid
T	=	temperature equation

Additional information

Funding

The authors acknowledge the financial supports from National Natural Science Foundation of China (No. 51606013, No. 51779025 and No. 21676040). This work is also supported by the Fundamental Research Funds for the Central Universities of China (No. 3132019191 and No. 3132019187).