Numerical Investigation of Effects of Obstacles in Flow Channels and Depth of Flow Channels for PEMFCs

Abstract

The channel is a crucial component of the polymer electrolyte membrane fuel cell (PEMFC). Since the channel can change the reactant transfer capability, water removal capability, and distribution of the reactant, it affects the performance and durability of PEMFCs. This study investigated the effects of obstacles in the serpentine-type flow channel on the performance of PEMFCs by computational fluid dynamics (CFD). The height of the obstacles was varied to analyze the electrochemical performances of the fuel cells. In addition, the depth of the flow channel was varied to compare the performances of the PEMFCs. To better represent the real-world tendency, the agglomerate model and the Forchheimer inertial effect were used. The results showed that changes in the channel depth caused greater performance improvements compared to the installation of obstacles, due to the enhanced mass transfer and improved water removal. However, the results for the installation of obstacles showed the lower non-uniformity of the current density and a reduced pressure drop compared to the changes in the channel depth, offering advantages in terms of flooding, the fuel cell life, and the operating cost.

1. Introduction

The hydrogen industry has gained attention due to carbon neutrality and the depletion of fossil fuels. The hydrogen industry is the best way to create a sustainable energy economy using alternative energy sources based on renewable energy [1]. In this context, the fuel cell is a key component as the energy conversion device that converts hydrogen energy to electrical energy in the hydrogen society. Since fuel cells produce electrical energy directly, they have excellent energy conversion efficiencies and produce virtually no harmful emissions. There are many types of fuel cells: polymer electrolyte membrane fuel cells (PEMFCs), solid oxide fuel cells, phosphoric acid fuel cells, and molten carbonate fuel cells. Fuel cells are distinguished by their electrolytes, which must conduct ions, not electrons [1]. Among them, the PEMFC electrolyte is a thin polymeric membrane with ionic conductivity based on water-based transport mechanisms. Thus, water management is crucial in the PEMFC performance [1]. PEMFCs have gained attention as the next-generation energy source for transportation and power generators because they have a high performance, good on/off ability, and good durability due to their low operating temperature. However, more research is needed because there are still technical problems, such as the low specific power, high cost of components, and difficult water management. These issues are due to the unoptimized PEMFC components: a bipolar plate (BPP), channel, gas diffusion layer (GDL), microporous layer (MPL), catalyst layer (CL), and membrane [2] The BPP separates the reactants between the anode and cathode, provides electrical paths, and offers structural support. Jang et al. made a polymer bipolar plate to make ultra-light PEMFCs, using sputter to overcome the polymer’s low electron conductivity [3]. Among the PEMFC components, the channel is a key component because it can reduce the concentration loss due to a reactant shortage. The channel facilitates the movement of reactants to the catalyst layer and removes the liquid water, which decreases the performance of PEMFCs because it blocks the reactant and causes irreversible damage [1]. Also, the channel ensures the uniform distribution of reactants at the catalyst layer. This uniform distribution helps to alleviate the degradation of PEMFCs because it helps to make a more uniform temperature and current density distribution [4]. For these reasons, many researchers have continued to improve the mass transfer and water removal capabilities of channels by modifying their structure or adding obstacles. However, experimentally, studying the movement of substances and heat inside the fuel cell has limitations in terms of the cost and technology. CFD can reduce the time spent experimenting to find the optimal value. Therefore, computational fluid dynamics (CFD) is needed to understand and analyze the movement of reactants and heat inside the fuel cell [5].

Resent advances in CFD have led to a lot of studies using it for various purposes. Li et al. developed an optimization model using dynamic mesh technology to investigate the multiphase vortex-induced vibration (MVIV) behavior in tidal power plants. They analyzed the relationship between the displacement responses and flow states and validated the MVIV transition behaviors using a fluid-induced vibration sensing platform [6]. In addition, many studies that modify the channel to improve the performance of fuel cells have been conducted using CFD. Ying et al. investigated various obstacle shapes placed in straight channels through numerical analysis to identify the shape with the best mass transport capability [7]. Yang et al. used CFD to compare the performances of different obstacles and height distributions within triple-channel serpentine flow fields [8]. Son et al. also utilized CFD to compare the fuel cell performances based on different obstacle patterns in parallel flow fields [9]. Yu et al. developed a sinusoidal and intercepted flow field and analyzed its performance through both experimentation and simulation, and they additionally conducted the performance comparison based on the number of inlets [10]. In this study, analyses were conducted using CFD to investigate the effects of various obstacle heights and flow field depths on the mass transfer capability in the single serpentine channel. Additionally, the impacts of these modifications on the distribution of the current density and pressure drop in the flow field were also compared.

2. Model Description

2.1. Analysis Model

The dimensions of the single-cell PEMFC model used in this study were based on the actual shape of the single cell used in experiments. The PEMFC model also consists of a BPP, channel, GDL, MPL, CL, and membrane. The detailed dimensions are listed in

Table 1. Geometrical parameters.

Geometrical Parameters	Value	Unit
Reaction area	5.06	cm2
Channel width	1	mm
Rib width	1.1	mm
BPP thickness	1.5	mm
GDL thickness	160	µm
MPL thickness	90	µm
CL thickness	5	µm
Membrane thickness	25.4	µm
Gap between obstacles	2.625	mm

.

The mass transfer capability of the cathode channel plays a more dominant role in determining the performance of the PEMFC, especially at low voltage, because oxygen has lower diffusivity than hydrogen, water is only produced at the cathode, and the cathode electrochemical reaction is more complex [11,12]. Therefore, in this study, obstacle placement and variation in the channel depth were applied to only the cathode channel. In all cases, the anode channel was the same as the cathode channel used in the reference case (Ref).

Figure 1 shows the PEMFC model schematics and the cross-section view along line A-A. The channel depth and obstacle height dimensions for each case are listed in

Table 2. Channel depth and obstacle height parameters [mm].

	Ref	O0.1	O0.2	O0.3	O0.4	D0.7	D0.6	D0.5	D0.4
d1	0.8	0.8	0.8	0.8	0.8	0.7	0.6	0.5	0.4
d2	-	0.1	0.2	0.3	0.4	-	-	-	-

. d1 is the channel depth and d2 is the obstacle height. The obstacle cases were represented based on their heights as 0.1 mm (O0.1), 0.2 mm (O0.2), 0.3 mm (O0.3), and 0.4 mm (O0.4), and the depth cases were indicated by their depths as 0.7 mm (D0.7), 0.6 mm (D0.6), 0.5 mm (D0.5), and 0.4 mm (D0.4).

2.2. Model Assumptions

The model developed for this study established the following assumptions to improve the analysis speed and reduce the computational load while maintaining the result accuracy:

(1)

The gas is an incompressible ideal gas and the flow is laminar [13];

(2)

The fuel cell operates at the steady state;

(3)

The effects of gravity are neglected;

(4)

Contact resistance is ignored [14];

(5)

The crossover phenomenon does not occur [15].

2.3. Governing Equations

In computational fluid dynamics, various governing equations are used to calculate the physical movement of the reactants and electrochemical reactions. As mentioned earlier, since the system operates at the steady state, the time terms do not exist from the governing equations.

The mass conservation equation is as follows:

$$\nabla ·\left(\rho \overrightarrow{u}\right)=S_{mass}$$

(1)

where $\rho$ is the gas density, and $\overrightarrow{u}$ is the gas velocity. $S_{mass}$ is the source term for mass generation.

The momentum conservation equation is as follows:

$$\nabla ·\left(\rho \overrightarrow{u}\overrightarrow{u}\right)=\nabla ·\left(\mu \nabla \overrightarrow{u}\right)-\nabla P+S_{mom}$$

(2)

where $\mu$ is the dynamic viscosity, and $P$ is the pressure. $S_{mom}$ is the momentum source term. The momentum source is the pressure loss in the porous layers and generally follows Darcy’s law. However, as the velocity of the reactants increases, a nonlinear pressure drop occurs in the porous layer due to inertial effects. In this study, the velocity in the porous layer increased due to the presence of obstacles and changes in the channel depth. Therefore, the Darcy–Forchheimer law, which includes inertial effects in the porous region, was used in the momentum source term to account for this pressure drop. For this reason, the model can better represent the real-world tendency [16].

The energy conservation equation is as follows:

$$\nabla ·\left(\rho c_{p}T\overrightarrow{u}\right)=\nabla ·\left(k_{eff}\nabla \mathrm{T}\right)+S_T$$

(3)

where $c_p$ is the specific heat, $k_{eff}$ is the effective thermal conductivity, and $T$ is the temperature. $S_T$ is the source term for thermal energy and includes the heat generation by electrochemical reactions, ohmic resistance, and phase changes of the substances.

The species conservation equation is as follows:

$$\nabla ·\left(\rho Y_i\overrightarrow{u}\right)=\nabla ·\left(\rho D_{eff,i}\nabla Y_i\right)+S_i$$

(4)

where $Y_i$ is the mass fraction of species (i), and $D_{eff,i}$ is the effective diffusivity of species (i). $S_i$ is the species source term of species (i), which is the species generated or consumed by the electrochemical reaction and phase change.

The charge conservation equation is as follows:

$$\nabla ·\left({\sigma }_{eff,s}\nabla {\varphi }_s\right)+S_s=0$$

(5)

$$\nabla ·\left({\sigma }_{eff,e}\nabla {\varphi }_e\right)+S_e=0$$

(6)

where ${\sigma }_{eff,s}$ is the effective electrical conductivity, and ${\sigma }_{eff,e}$ is the effective ionic conductivity. ${\varphi }_s$ is the solid potential, and ${\varphi }_e$ is the electrolyte potential. $S_s$ and $S_e$ are the potential source terms, and they are current by the electrochemical reaction.

The liquid water transfer equation is as follows:

$$\nabla ·\left({\rho }_{l}f\overrightarrow{u}\right)=\nabla ·\left({\rho }_l{D}_l\nabla s\right)+S_{vl}$$

(7)

where ${\rho }_l,f,D_l$, and $s$ are the liquid density, interfacial drag coefficient, liquid diffusivity, and liquid volume fraction, respectively. $S_{vl}$ is the liquid source term, which is the liquid generated or disappeared by condensation and evaporation.

The dissolved water transfer equation is as follows:

$$-\nabla ·\left({\displaystyle \frac{n_d}{F}}{\sigma }_{eff,e}\nabla {\varphi }_e\right)=\nabla ·\left({\displaystyle \frac{{\rho }_m}{EW}}{D}_{dw}\nabla \lambda \right)+S_{dw}$$

(8)

The above equation is used to calculate the movement of dissolved water in the electrolyte, where $n_d$ is the electro-osmotic drag coefficient, which is the number of water molecules dragged per proton. $F,{\rho }_m$, and $EW$ are the Faraday constant, membrane density, and equivalent weight of the membrane, respectively. $D_{dw}$ is the dissolved water diffusivity, which is the function of the water content ($\lambda$). $S_{dw}$ is the dissolved water source term, which represents the dissolved water generated or consumed by the electrochemical reaction and the water’s absorption/desorption in the ionomer.

The PEMFC model used in this study was the spherical agglomerate model. The spherical agglomerate model adds water and ionomer terms to the cathode current generation equation to reflect the effects of the water and ionomer film. Therefore, the model can represent the electrochemical kinetics due to the slow oxygen reduction reaction in the cathode CL more accurately. The source terms and parameters used in the governing equations are shown in

Table 3. Source terms [9,17].

Location	Description	Unit
Anode CL Cathode CL Other	$S_{mass}=\left\{\begin{array}{c}{S}_{H_2}-S_{vl}-M_{H_{2}O}{S}_{vd}\\ S_{O_2}-S_{vl}-M_{H_{2}O}{S}_{vd}\\ {-S}_{vl}\end{array}\right.$	kg m−3s−1
All	$S_{mom}=-\left({\displaystyle \frac{\mu }{K}}\overrightarrow{u}+\beta {\displaystyle \frac{\rho }{{ϵ}^{1.5}K}}\left|\overrightarrow{u}\right|\overrightarrow{u}\right)$	kg m−2s−2
Anode CL	$S_{H_2}=-{\displaystyle \frac{M_{H_2}{J}_a}{2F}}$	kg m−3s−1
Cathode CL	$S_{O_2}=-{\displaystyle \frac{M_{O_2}{J}_c}{4F}}$	kg m−3s−1
CL Other	$S_{H_{2}O}=\left\{\begin{array}{c}-S_{vl}-M_{H_{2}O}{S}_{vd}\\ -S_{vl}\end{array}\right.$	kg m−3s−1
All except membrane	$S_{vl}=\left\{\begin{array}{c}{\gamma }_{cond}\left(1-s\right)M_{H_{2}O}\left(P_{H_{2}O}-P_{sat}\right)/RT\\ {\gamma }_{evap}s{M}_{H_{2}O}\left(P_{H_{2}O}-P_{sat}\right)/RT\end{array}\begin{array}{c}{P}_{H_{2}O}\ge P_{sat}\\ P_{H_{2}O}<P_{sat}\end{array}\right.$	kg m−3s−1
CL	$S_{vd}\left\{\begin{array}{c}{\gamma }_{ads}{\displaystyle \frac{{\rho }_{mem}}{EW}}\left({\lambda }_{eq}-\lambda \right){\lambda }_{eq}\ge \lambda \\ {\gamma }_{des}{\displaystyle \frac{{\rho }_{mem}}{EW}}\left({\lambda }_{eq}-\lambda \right){\lambda }_{eq}<\lambda \end{array}\right.$	mol m−3
Membrane	$S_T=\left\{\begin{array}{c}{\displaystyle \frac{i_e}{{\sigma }_{eff,e}}}\\ {\displaystyle \genfrac{}{}{0pt}{}{\left({\eta }_a-\frac{T∆S_a}{2F}\right)j_a+\frac{i_s^2}{{\sigma }_{eff,s}}+\frac{i_e^2}{{\sigma }_{eff,e}}+h_{lg}{S}_{vl}}{\left({\eta }_c-\frac{T∆S_c}{4F}\right)j_a+\frac{i_s^2}{{\sigma }_{eff,s}}+\frac{i_e^2}{{\sigma }_{eff,e}}+h_{lg}{S}_{vl}}}\\ {\displaystyle \genfrac{}{}{0pt}{}{\frac{i_s^2}{{\sigma }_{eff,s}}+h_{lg}{S}_{vl}}{\frac{i_s^2}{{\sigma }_{eff,s}}}}\end{array}\right.$	W m−3
Anode CL
Cathode CL
MPL, GDL, chn
BPP
Anode CL	$S_s=\left\{\begin{array}{c}-j_a\\ j_c\end{array}\right.$	A m−3
Cathode CL
Anode CL	$S_e=\left\{\begin{array}{c}{j}_a\\ -j_c\end{array}\right.$	A m−3
Cathode CL
Anode CL	$S_{dw}=\left\{\begin{array}{c}{S}_{vd}\\ S_{vd}+{\displaystyle \frac{j_c}{2F}}\end{array}\right.$	mol m−3
Cathode CL

,

Table 4. Parameters.

Parameter	Value	Unit
Anode transfer coefficient [18], ${\alpha }_a$	0.5	-
Cathode transfer coefficient [18], ${\alpha }_c$	1.0	-
Entropy change by hydrogen oxidation [19], $∆S_a$	0.104	J mol−1 K−1
Entropy change by oxygen reduction [19], $∆S_c$	−326.36	J mol−1 K−1
Electrical conductivity of BPP, GDL, MPL, CL [20,21], ${\sigma }_{BPP},{\sigma }_{GDL},{\sigma }_{MPL},{\sigma }_{CL}$	92,600, 8000, 5000, 5000	S m−1
Thermal conductivity of BPP, MPL, CL, membrane [20,21,22], $k_{BPP},k_{MPL},k_{CL},k_{mem}$	120, 1, 1, 0.16	W m−1K−1
Through-plane thermal conductivity of GDL [20], $k_{GDL}$	1.7	W m−1K−1
In-plane thermal conductivity of GDL [20], $k_{GDL}$	21	W m−1K−1
Porosity of GDL, MPL [23], ${ϵ}_{GDL},{ϵ}_{MPL}$	0.8, 0.7	-
Permeability of GDL, MPL, CL [20], $K_{GDL},K_{MPL},K_{CL}$	1 × 10−11, 1 × 10−12, 1 × 10−13	m2
Contact angle of GDL, MPL, CL [24], ${\theta }_{GDL},{\theta }_{MPL},{\theta }_{CL}$	110, 130, 95	deg
Surface tension [24], $\sigma$	0.0625	N m−1
Latent heat of water phase change [25], $h_{lg}$	2.36 × 106	J kg−1
Membrane density [26], ${\rho }_{mem}$	1970	kg m−3
Membrane equivalent weight [26], $EW$	1.050	Kg mol−1
Condensation rate [17], ${\gamma }_{cond}$	100	s−1
Evaporation rate [17], ${\gamma }_{evap}$	100	s−1
Inertial coefficient [27], $\beta$	2.88 × 10−6
Platinum loading [28], $m_{pt}$	0.4	mg cm−2
Platinum density [29], ${\rho }_{pt}$	21,450	kg m−3
Carbon loading, $m_C$	0.27	mg cm−2
Carbon density [28], ${\rho }_C$	1800	kg m−3
Agglomerate radius, $r_{agg}$	0.2	µm
Ionomer volume fraction in catalyst layer [27], $L_i$	0.4	-
Ionomer volume fraction in agglomerate [27], $L_{i,agg}$	0.5	-

and

Table 5. Equations for electrochemical reaction.

Description	Equation		Unit
Open-circuit voltage [1]	$V_{oc}=1.229-0.0008456\left(T-298.15\right)+{\displaystyle \frac{RT}{2F}}ln{P}_{H_2}{P}_{O_2}^{0.5}$		V
Overpotential [25]	${\eta }_a={\varphi }_s-{\varphi }_e,{\eta }_c={\varphi }_s-{\varphi }_e-V_{oc}$		V
Anode volumetric current density [30]	$j_a=\left(1-s\right)a_{eff}{i}_{o,a}{\left({\displaystyle \frac{P_{H_2}}{C_{H_2}^{ref}{H}_{H_2}}}\right)}^{0.5}\left[\exp \left({\displaystyle \frac{{\alpha }_{a}F{\eta }_a}{RT}}\right)-exp\left({\displaystyle \frac{-{\alpha }_{c}F{\eta }_a}{RT}}\right)\right]$		A m−3
Cathode reaction rate [30]	$k_c={\displaystyle \frac{{\alpha }_{eff}{i}_{0,c}}{4F\left(1-{ϵ}_{CL}\right)C_{O_2}^{ref}}}\left[-\exp \left({\displaystyle \frac{{\alpha }_{a}F{\eta }_c}{RT}}\right)+exp\left({\displaystyle \frac{-{\alpha }_{c}F{\eta }_c}{RT}}\right)\right]$		A m−3
Anode exchange current density [31]	$i_{0,a}=1·{10}^2\exp \left[-{\displaystyle \frac{E_a^{act}}{R}}\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{353}}\right)\right]$		A m−2
Cathode exchange current density [30]	$i_{0,c}=5·{10}^{-4}\exp \left[-{\displaystyle \frac{E_c^{act}}{R}}\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{353}}\right)\right]$		A m−2
Proton conductivity [1]	${\sigma }_e=\left(0.514\lambda -0.326\right)\exp \left[1268\left({\displaystyle \frac{1}{303}}-{\displaystyle \frac{1}{T}}\right)\right]$		S m−1
Effective proton conductivity [18]	${\sigma }_{eff,e}=L_i^{1.5}{\sigma }_e$	Membrane	S m−1
Effective electrical conductivity [18]	${\sigma }_{eff,s}=\left\{\begin{array}{c}{{\left(1-ϵ-L_i\right)}^{1.5}\sigma }_s\\ {\left(1-ϵ\right)}^{1.5}{\sigma }_s\end{array}\right.$	CL MPL, GDL	S m−1
H2 and H2O diffusivity on anode side [32]	$D_{H_2,a}=D_{H_{2}O,a}=1.055\times {10}^{-4}{\left({\displaystyle \frac{T}{333}}\right)}^{1.75}\left({\displaystyle \frac{101325}{P}}\right)$		m2 s−1
O2 diffusivity on cathode side [31]	$D_{O_2,c}=0.2652\times {10}^{-4}{\left({\displaystyle \frac{T}{333}}\right)}^{1.75}\left({\displaystyle \frac{101325}{P}}\right)$		m2 s−1
H2O diffusivity on cathode side [31]	$D_{H_{2}O,c}=0.2982\times {10}^{-4}{\left({\displaystyle \frac{T}{333}}\right)}^{1.75}\left({\displaystyle \frac{101325}{P}}\right)$		m2 s−1
Effective diffusivity [9]	$D_{eff,i}={ϵ}^{1.5}{\left(1-s\right)}^{1.5}{D}_i$		m2 s−1
Dissolved O2 diffusivity in ionomer [31]	$D_{O_{2,i}}=3.1\times {10}^{-7}\exp \left(-{\displaystyle \frac{2768}{T}}\right)$		m2 s−1
Dissolved O2 diffusivity in liquid [31]	$D_{O_2\omega }=1.98\times {10}^{-9}\left({\displaystyle \frac{{\mu }_{l@293K}}{{\mu }_l}}\right)\left({\displaystyle \frac{T}{293}}\right)$		m2 s−1
Saturation pressure [17]	$P_{sat}={10}^{(-2.1794+0.02953(T-273.15)-9.1837\times {10}^{-5}{\left(T-273.15\right)}^2+1.4454\times {10}^{-7}{\left(T-273.15\right)}^3)}$		atm
Water activity [17]	$a={\displaystyle \frac{P_{H_{2}O}}{P_{sat}}}+2s$		-
Interfacial drag coefficient [31]	$f={\displaystyle \frac{K_l{\mu }_g}{K_g{\mu }_l}}$		-
Relative permeability [33]	$K_l=s^3,K_g={\left(1-s\right)}^3$		m2
Capillary pressure [34]	$P_c=\sigma \cos \theta {\left({\displaystyle \frac{ϵ}{K}}\right)}^{0.5}(1.417s-2.12s^2+1.263s^3)$		Pa
Liquid water diffusivity [31]	$D_l=-{\displaystyle \frac{K{K}_l}{{\mu }_l}}{\displaystyle \frac{d{P}_c}{ds}}$		m2 s−1
Equilibrium water content [17]	${\lambda }_{eq}=\left\{\begin{array}{c}0.043+17.81a-39.85a^2+36a^3\\ 14+1.4\left(a-1\right)\\ 16.8\end{array}\right.\begin{array}{c}a\le 1\\ 1<a\le 3\\ a>3\end{array}$		-
Electro-osmotic drag coefficient [28]	$n_d=2.5{\displaystyle \frac{\lambda }{22}}$		-
Dissolved water diffusivity [17]	$D_{dw}=\left\{\begin{array}{c}2.693\times {10}^{-10},(\lambda \le 2)\\ D_{dw}^0\left[0.87\left(3-\lambda \right)+2.95\left(\lambda -2\right)\right],(2<\lambda \le 3)\\ \begin{array}{c}{D}_{dw}^0\left[2.95\left(4-\lambda \right)+1.64\left(\lambda -3\right)\right],(3<\lambda \le 4)\\ D_{dw}^0\left(2.563-0.33\lambda +0.0264{\lambda }^2-0.000671{\lambda }^3\right),(\lambda >4)\end{array}\end{array}\right.$ $where{D}_{dw}^0={10}^{-10}\exp \left[2516\left({\displaystyle \frac{1}{303}}-{\displaystyle \frac{1}{T}}\right)\right]$		m2 s−1
Adsorption rate [28]	${\gamma }_{ads}={\displaystyle \frac{1.14\times {10}^{-5}}{t_{CL}}}{\displaystyle \frac{\lambda V_{\omega }}{V_m+\lambda V_{\omega }}}\exp \left[2416\left({\displaystyle \frac{1}{303}}-{\displaystyle \frac{1}{T}}\right)\right]$		s−1
Desorption rate [28]	${\gamma }_{des}={\displaystyle \frac{4.59\times {10}^{-5}}{t_{CL}}}{\displaystyle \frac{\lambda V_{\omega }}{V_m+\lambda V_{\omega }}}\exp \left[2416\left({\displaystyle \frac{1}{303}}-{\displaystyle \frac{1}{T}}\right)\right]$		s−1
Molar volume [28]	$V_{\omega }={\displaystyle \frac{M_{H_{2}O}}{{\rho }_l}},V_m={\displaystyle \frac{EW}{{\rho }_m}}$		m3 mol−1

. The equations for the spherical agglomerate model are listed in

Table 6. Equations of spherical agglomerate model [18].

Description	Equation	Unit
Cathode volumetric current density	$j_c=4F{\displaystyle \frac{P_{O_2}}{H_{O_2}}}{\left[{\displaystyle \frac{1}{\xi \left(1-{ϵ}_{CL}\right)k_c}}+{\displaystyle \frac{r_{agg}+{\delta }_i+{\delta }_{\omega }}{r_{agg}}}\left({\displaystyle \frac{{\delta }_i}{a_{agg,i}{D}_{O_2,i}}}+{\displaystyle \frac{{\delta }_{\omega }}{a_{agg,\omega }{D}_{O_2,\omega }}}\right)\right]}^{-1}$	A m−3
Pt mass ratio in Pt/C particles	$M{R}_{Pt}={\displaystyle \frac{m_{Pt}}{m_{Pt}+m_c}}$	-
Pt/C volume fraction in catalyst layer	$L_{pt/C}={\displaystyle \frac{m_{Pt}}{t_{CL}}}\left({\displaystyle \frac{1}{{\rho }_{Pt}}}+{\displaystyle \frac{1}{{\rho }_C}}{\displaystyle \frac{1-M{R}_{Pt}}{M{R}_{Pt}}}\right)$	-
Active surface area	$a_{eff}={\displaystyle \frac{{10}^3{m}_{Pt}}{t_{CL}}}\left(222.79M{R}_{Pt}^3-158.57M{R}_{Pt}^2-201.53M{R}_{Pt}+159.5\right)$	m−1
CL porosity	${ϵ}_{CL}=1-L_{Pt/C}-L_i$	-
Number of agglomerate particles per CL volume	$N_{agg}={\displaystyle \frac{3L_{Pt/C}}{4\pi r_{agg}^3\left(1-L_{i,agg}\right)}}$	-
Ionomer film thickness	${\delta }_i=r_{agg}\left[\sqrt[3]{{\displaystyle \frac{L_i\left(1-L_{i,agg}\right)}{L_{Pt/C}}}-L_{i,agg}+1}-1\right]$	m
Water film thickness	${\delta }_{\omega }=\sqrt[3]{{\left(r_{agg}+{\delta }_i\right)}^3+{\displaystyle \frac{3s{ϵ}_{CL}}{4\pi N_{agg}}}}-\left(r_{agg}+{\delta }_i\right)$	m
Effective agglomerate surface area of ionomer film	$a_{agg,i}={\displaystyle \frac{3L_{\frac{Pt}{C^{{ϵ}_{CL}}}}}{r_{agg}^3\left(1-L_{i,agg}\right)}}{\left(r_{agg}+{\delta }_i\right)}^2$	m−1
Effective agglomerate surface area of water film	$a_{agg,\omega }={\displaystyle \frac{3L_{\frac{Pt}{C^{{ϵ}_{CL}}}}}{r_{agg}^3\left(1-L_{i,agg}\right)}}{\left(r_{agg}+{\delta }_i+{\delta }_{\omega }\right)}^2$	m−1
Thiele modulus	$\mathsf{\Phi }={\displaystyle \frac{r_{agg}}{3}}\sqrt{{\displaystyle \frac{k_c}{L_{i,agg}^{1.5}{D}_{O_{2,i}}}}}$	-
Effectiveness factor	$\xi ={\displaystyle \frac{1}{\mathsf{\Phi }}}\left({\displaystyle \frac{1}{\tanh 3\mathsf{\Phi }}}-{\displaystyle \frac{1}{3\mathsf{\Phi }}}\right)$	-

.

2.4. Validation

The PEMFC model used in this study was developed using Ansys Fluent 2023R2. Physical variables not supported by the program were expressed using User-Defined Functions (UDFs). The SIMPLE algorithm was chosen to solve the governing equations, and second-order upwind was used for all physical variables. The convergence criteria were defined with the volumetric current density at the anode and cathode, being within 10⁻⁴ A/cm² and a residual error of 10⁻⁶ for all physical variables. Figure 2 shows the I–V curves obtained from both the experiment and simulation. Since the experimental results match well with the simulation data, the PEMFC model is reliable. The operating conditions of the experiment are summarized in

Table 7. Operating conditions.

Parameters	Value	Unit
Anode, cathode inlet temperature	80/80	°C
Anode, cathode humidity	100/100	%
Anode, cathode stoichiometry	2.0/2.0	-
Reference current density	1.0	A cm−2

.

3. Results and Discussion

In this study, how the obstacle height and channel depth affect the mass transfer and water removal capabilities of the channel was investigated. The simulation was conducted in the 0.7 V to 0.35 V range, where the concentration loss could be well observed, to improve the simulation time efficiency. Specifically, the analysis focused on 0.5 V, which was the maximum power density.

Figure 3 and Figure 4 show the current density and the temperature contours at the interface between the cathode CL and MPL at 0.5 V. The current density has the higher value at the inlet region because there was enough oxygen at the inlet. However, at the outlet region, the current density decreased because oxygen was consumed, and water was produced by the reaction. The temperature distribution shows the tendency to follow the current density contours because a high current density means a lot of reactions and the reaction makes heat. Figure 5 shows the I–V curves for all cases in the 0.7 V–0.35 V range of the cell voltage. It is observed that the current density of the PEMFC increased in all cases, especially in the low-voltage region, where concentration losses dominate. These results show that both the installation of obstacles and changes in the channel depth helped to reduce the concentration losses of the fuel cell. Therefore, both the installation of obstacles and changes in the channel depth improved the performance of the fuel cell.

3.1. Analysis of the Obstacle Height

Figure 6 shows the average oxygen molar concentration in the CL. Figure 7 shows the oxygen molar concentration contour at the interface between the cathode CL and MPL when obstacles were installed in the channel. Due to oxygen consumption by the electrochemical reaction, the oxygen molar concentration was high at the inlet and decreased along the line of the flow channel in all cases. The average oxygen molar concentration increased in all cases. Compared to the reference case, O0.1 shows an increase of about 1.4%, while O0.4 shows an increase of about 24%. This enhancement can be explained by the convection mechanism. The installation of the obstacles changed in the flow direction. Thus, by hitting the obstacles, the flow direction was changed from along the channel to the vertical direction of the channel, which made oxygen move easily to the catalyst layer. In addition, the installation of the obstacles generated an additional pressure drop because changes in the flow direction acted as an increase in the length of the flow path, and the cross-section area of the flow channel also became smaller at the region where the obstacle was installed. As a result, the additional pressure drop made the greater pressure difference between the channels, and it made a greater by-pass flow, which increased the oxygen concentration under the rib [25]. These results show that the installation of obstacles improved the mass transfer capability of the channel by the enhancement of the convection mechanism, and the higher the obstacle, the greater the increase in the oxygen molar concentration, because the higher the obstacle, the greater the change in the direction of the oxygen flow and the greater the pressure drop.

Figure 8 shows the average liquid saturation in the CL. Figure 9 shows the liquid saturation contour at the interface between the cathode CL and MPL when obstacles were installed in the channel. Since water was generated by the electrochemical reaction, the liquid saturation was low at the inlet and increased along the line of the flow channel in all cases. The installation of obstacles reduced the average liquid saturation in the CL. Compared to the reference case, O0.1 shows a decrease of about 1.4% and O0.4 shows a decrease of about 17%. This is because of the velocity increases in the porous region. Figure 10 shows the velocity contour at the interface between the GDL and the MPL. As mentioned above, the velocity in the porous region increases due to the enhancement of the convection flow, not only at the rib part but also at the channel part in the porous region, and the higher the obstacle, the greater the velocity. Since the high velocity makes it difficult to accumulate liquid water, the installation of the obstacles decreases the average liquid saturation in the CL [10].

Consequently, the installation of obstacles enhances the mass transfer and water removal capabilities of the channel by changing the flow direction and increasing the by-pass flow.

3.2. Analysis of Channel Depth

Figure 11 shows the average oxygen molar concentration in the CL. Figure 12 shows the oxygen molar concentration contour at the interface between the cathode CL and MPL when the channel depth was varied. Similar to the effect of the obstacle installation, reducing the channel depth in all cases increased the oxygen molar concentration. Compared to the reference, D0.7 shows an increase of about 5.4% and D0.4 shows an increase of about 30%. This improvement can be explained by the convection and diffusion mechanism. Similar to the obstacle cases, reducing the channel depth generated an additional pressure drop because reducing the channel depth made the cross-section area of the flow channel small. Therefore, reducing the depth enhances the by-pass flow. In addition, reducing the depth improves the effect of the diffusion because the distance that needs to be traveled by diffusion decreases. These results show that reducing the channel depth improves the mass transport capability of the channel, and the shallower the depth, the greater the increase in the oxygen molar concentration, because the shallower the depth, the greater the pressure drop and the shorter the distance of the diffusion.

Figure 13 shows the average liquid saturation in the CL. Figure 14 shows the liquid saturation contour at the interface between the cathode CL and MPL when the channel depth was varied. A shallower channel depth results in lower liquid saturation values. Compared to the reference, D0.7 shows a decrease of about 4.4% and D0.4 shows a decrease of about 23%. This can also be explained by the velocity contour. Figure 15 shows the velocity contours at the interface between the cathode GDL and MPL when the channel depth was varied. The velocity increased only at the rib part in the porous region, and the shallower the depth, the greater the velocity. The velocity at the channel part in the porous region had no difference compared to the reference case because there was no obstacle.

Consequently, reducing the channel depth enhanced the mass transfer and water removal capabilities of the channel by shortening the diffusion path and enhancing the by-pass flow.

3.3. Non-Uniformity

The uneven distribution of the current density decreases the stability and durability of the PEMFC because the current density non-uniformity resulting from the uneven distribution of the reactants leads to the formation of hot spots [4]. Therefore, evaluating the current density non-uniformity is essential to investigate the effect of the modification of the channel. However, evaluating the non-uniformity by contours is difficult to measure quantitatively. Therefore, Equation (9) was used to quantitatively evaluate the current density non-uniformity [8]:

$$E=\sqrt{{\displaystyle \frac{\int {\left(J-\overline{J}\right)}^{2}dS}{{\overline{J}}^2\int dS}}}$$

(9)

$$\overline{J}={\displaystyle \frac{\int JdS}{\int dS}}$$

(10)

where $J$ is the current density, and $\overline{J}$ is the average current density over the interface of the MPL and the CL. The calculated $E$ indicates the non-uniformity, with lower values representing a more uniform distribution of the substance across the plane. The non-uniformity obtained by the equation is shown in the bar graph in Figure 16. In all cases with modification, the non-uniformity decreased compared to the Ref case, and a greater modification degree led to lower non-uniformity. Only the O0.1 case had higher non-uniformity compared to the D0.7 case. As a result, the cases with the installation of the obstacle had higher stability and durability than the other cases, except for the O0.1 case.

3.4. Pressure Drop

The high pressure drop in the PEMFC system creates a high average pressure and a large pressure difference between the anode and cathode, which reduces the lifetime of the PEMFC and damages the membrane [35,36]. In addition, the increased pressure drop results in higher pump power required to supply the reactants [25]. Therefore, a comparison of the pressure drop is required. Figure 17 shows the pressure drop results obtained for each case in the bar graph. The pressure drop is calculated from the difference between the outlet pressure and the inlet pressure. The installation of the obstacle and the reduction in the channel depth increase the pressure drop in the fuel cell system. It was observed that as the height of the obstacles increased and the channel depth decreased, the pressure drop also increased. In addition, changes in the channel depth had a more significant impact on increasing the pressure drop compared to the installation of the obstacle. As a result, obstacle cases have advantages in terms of the pressure drop compared to depth cases.

4. Conclusions

In this study, how the installation of the obstacle and the reduction in the channel depth in the serpentine channel affect the mass transport capability and water removal capacity of the channel was investigated by CFD. The results are as follows:

(1)

The installation of the obstacle and the reduction in the channel depth improved the performance of the fuel cell, with a greater degree of modification leading to higher performance improvements. Especially, the performance was remarkably improved in the low-voltage region, where the concentration loss is dominant;

(2)

The changes in the channel depth were more effective than the installation of obstacles at enhancing the mass transport and water removal capabilities of the channel;

(3)

The installation of obstacles can increase the performance of the fuel cell with lower non-uniformity and lower pressure losses compared to reductions in the channel depth. Therefore, the installation of obstacles has higher stability and durability and a long lifetime compared to the reduction in the channel depth;

(4)

However, the O0.1 case had higher non-uniformity compared to the D0.7 case. Therefore, in this case, the obstacle must be at least 0.2 mm in height to have the benefit of low non-uniformity when installed.

Hopefully, the results of this study will help researchers in optimizing PEMFC flow channels. In addition, the gap between the obstacles, the shape of the obstacles, and the flow channel type are important characteristics that can influence the transport mechanisms of the reactants and pressure drop. The effects of these parameters will be investigated in future work.