Parametric Analysis of a Serpentine Flow Pattern Proton Exchange Membrane Fuel Cell for Optimized Performance

A three-dimensional model of a PEM fuel cell consists of a number of equations that are used to simulate fluid flow, gas diffusion, electrochemical process and transport of water and other gases.

The simplified model used here for the simulation assumes that the flow is steady and laminar with all the gases and vapors following ideal gas equation. All the water generated as the result of the reaction of hydrogen with oxygen is in gaseous form. The model equations have been developed and used by many researchers and are easily available in the fuel cell module manual of ANSYS fluent 18.1 [7] and also in literature [8-14]. The governing equations and PEM fuel cell sub-models are given in the following sections.

2.1 Continuity equation

The conservation of mass or continuity equation for a three-dimensional compressible flow is given by Eq. (4).

$\frac{\partial(\rho u)}{\partial x}+\frac{\partial(\rho v)}{\partial y}+\frac{\partial(\rho w)}{\partial z}=-\frac{\partial \rho}{\partial t}$     (4)

where, u, v and w are velocities in the x, y and z direction respectively, $\rho$ is density of reactant gases. Taking into account the porosity of electrodes and membrane this equation transforms to Eq. (5).

$\frac{\partial(\rho \epsilon u)}{\partial x}+\frac{\partial(\rho \epsilon v)}{\partial y}+\frac{\partial(\rho \epsilon w)}{\partial z}=S_{m}$     (5)

$S_{m}$ is the mass sink term and is zero except in the catalyst layer zone and are calculated for hydrogen, oxygen and water using Eqns. (6)-(8).

$S_{H_{2}}\left(k g s^{-1} m^{-3}\right)=\frac{M_{H_{2}}}{2 F} R_{a n}$     (6)

$S_{O_{2}}\left(k g s^{-1} m^{-3}\right)=\frac{M_{O_{2}}}{4 F} R_{c a t}$     (7)

$S_{H_{2 O}}\left(k g s^{-1} m^{-3}\right)=\frac{M_{H_{2} O}}{2 F} R_{c a t}$     (8)

Here, $F$ is the Faraday constant and $\mathrm{M}$ is the molecular weight of the species in $\mathrm{kg} / \mathrm{mol} . R_{a n}$ and $R_{c a t}$ is the exchange current density at the anode and cathode respectively.

2.2 Momentum equation

For porous electrodes, the general form of the momentum equation is given by Eq. (9).

$\nabla \cdot(\epsilon \rho \vec{v} \vec{v})=-\epsilon \nabla p+\nabla \epsilon \tau+S_{m o m}$     (9)

In x-direction, this equation transforms into Eq. (10).

$\begin{aligned} u \frac{\partial(\rho u)}{\partial x} &+v \frac{\partial(\rho u)}{\partial y}+w \frac{\partial(\rho u)}{\partial z} \\ &=-\frac{\partial p}{\partial x}+\frac{\partial}{\partial x}\left(\mu \frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\left(\mu \frac{\partial u}{\partial y}\right) \\ &+\frac{\partial}{\partial z}\left(\mu \frac{\partial u}{\partial z}\right)+S_{\operatorname{mom} x} \end{aligned}$     (10)

Here, $\vec{v}, p, \tau$ and $S_{m o m}$ are velocity vector, pressure, shear stress tensor, viscosity and momentum sink term respectively. Momentum term is zero in flow channels and for the porous zones of the model is given by Eq. (11).

$S_{\operatorname{mom}, x}=-\frac{\mu u}{\beta_{x}}$     (11)

where, $\beta$ is the permeability of the electrodes and is same in all having a value of $10^{-12} \mathrm{m}-12 .$ The momentum equations for y and z directions can be easily formulated by replacing appropriate terms in Eqns. (10) and (11).

2.3 Mass transfer equations

The continuity equation for a certain species i is given by Eq. (12).

$\nabla \cdot\left(\rho \vec{v} y_{i}\right)=-\nabla \cdot \vec{\jmath}_{i}+S_{i}$     (12)

The diffusion mass flux j can be calculated using Fick’s law as given by Eq. (13).

$\overrightarrow{J_{i}}=-\sum_{j=1}^{N-1} \rho D_{i j} \nabla \overrightarrow{y_{i}}$     (13)

Within the porous electrodes mass transfer equation for any species can be written as given by Eq. (14).

$\nabla \cdot\left(\rho \epsilon \vec{v} y_{i}\right)=\nabla \cdot\left(\rho \epsilon D_{i j}^{e f f} \nabla y_{i}\right)+S_{i}$     (14)

where, $y_{i}$ is the mass fraction of species i (H₂, O₂, and H₂O.) diffusivity can be estimated using modified Brueggemann equation (Eq. (15)).

$D_{i j}^{e f f}=D_{i j} \times \epsilon^{1.5}$     (15)

The mass transfer equations for all species has the following general form as given by Eq. (16).

$\begin{aligned} u \frac{\partial\left(\rho y_{i}\right)}{\partial x} &+v \frac{\partial\left(\rho y_{i}\right)}{\partial y}+w \frac{\partial\left(\rho y_{i}\right)}{\partial z} \\ &=\frac{\partial\left(j_{x, i}\right)}{\partial x}+\frac{\partial\left(\rho j_{y, i}\right)}{\partial y}+\frac{\partial\left(j_{z, i}\right)}{\partial z} \\ &+S_{i} \end{aligned}$     (16)

2.4 Mass transfer equations

The rate of change of energy can be modeled using Eq. (17)

$\nabla \cdot(\vec{v}(\rho E+P))=\nabla \cdot\left(k_{e f f} \nabla T-\sum_{j} h_{j} \vec{\jmath}_{j}+\right.$

$\left.\left(\tau_{e f f} \cdot \vec{v}\right)\right)+S_{h}$     (17)

Here $E, h, \boldsymbol{\tau}_{e f f}$ are the total energy, enthalpy and effective shear tensor respectively. $k_{e f}$ is the effective thermal conductivity and can be calculated using Eq. (18) from the thermal conductivity of the fluid, $\left(k_{f}\right)$ and solid, $\left(k_{s}\right)$.phases.

$k_{e f f}=\epsilon k_{f}+(1-\epsilon) k_{f}$     (18)

S_h is the sink term and can be calculated using Eq. (19).

$S_{h}=h_{p h a s e}+h_{r e a c t i o n}+R_{o h m} I^{2}$     (19)

where, h_phase is neglected as there is no phase change modeled and h_reaction is zero in bipolar plate as there is no reaction. Hence Eq. (20) can be used to estimate the temperature change.

$\nabla \cdot(k \nabla T)=-R_{o h m} I^{2}$     (20)

Here k is the thermal conductivity of the bipolar plate.

2.5 Charge conservation equations

Electrochemical reactions occur in the fuel cell that generate current and the driving force behind these reactions is the surface activation over-potential, $\phi$ The conservation of electron and ions is modeled in the solid parts and in membrane using Eqns. (21) and (22).

$\nabla \cdot\left(\sigma_{s o l} \nabla \phi_{s o l}\right)+R_{s o l}=0$     (21)

$\nabla \cdot\left(\sigma_{m e m} \nabla \phi_{m e m}\right)+R_{m e m}=0$     (22)

$\sigma$ is Electrical conductivity of membrane or solid phase of the fuel cell model and R_sol and R_mem are the volume sink terms or current density (A/m³) and are given by set of Eqns. (23)-(26).

Anode side: $R_{s o l}=-R_{a n}(<0)$     (23)

and

$R_{m e m}=R_{a n}(>0)$     (24)

Cathode side: $R_{s a l}=R_{c a t}(>0)$     (25)

and

$R_{m e m}=-R_{c a t}(<0)$     (26)

These terms can be calculated using Butler-Volmer equations (Eqns. (27) and (28)).

$R_{a n}=R_{a n}^{e f f}\left(\frac{C_{H_{2}}}{C_{H_{2}}^{r e f}}\right)^{\gamma a n}\left(e^{\left(\frac{\alpha_{a n} F}{R T}\right) \eta_{a n}}-e^{\left(\frac{\alpha_{c a t} F}{R T}\right) \eta_{c a t}}\right)$     (27)

$R_{c a t}=R_{c a t}^{e f f}\left(\frac{C_{O_{2}}}{C_{O_{2}}^{r e f}}\right)^{\gamma c a t}\left(e^{\left(\frac{\alpha_{c a t} F}{R T}\right) \eta_{c a t}}-e^{\left(\frac{\alpha_{a n} F}{R T}\right) \eta_{a n}}\right)$     (28)

Average current density is calculated from Eq. (29).

$i_{a v e}=\frac{1}{A} \int_{V_{a n}} R_{a n} d V=\frac{1}{A} \int_{V_{c a t}} R_{c a t} d V$     (29)

The water content of the membrane can be estimate using Eq. (30).

$\lambda=\left\{\begin{array}{l}0.043+17.81 a+39.85 a^{2}+36 a^{3}, 0 \leq a \leq 1 \\ 14+1.4(a-1)\end{array}, 0 \leq a \leq 1\right.$     (30)

2.6 Water transport through membrane

The water generated as a result of the reaction of hydrogen with the oxygen at the cathode of a PEM fuel cell diffuse towards the anode due to electro-osmosis and back diffusion. To reduce the computation time and load, single phase was assumed and all the water formed was considered to be in gaseous form.

The results of the numerical simulations were interpreted and then analyzed by plotting various output parameters of the fuel cell against the operating pressure, mass flow rates and temperatures. A detail analysis is provided here based on the numerical study of the varying operating conditions on the performance, current output, hydrogen consumption etc. of a PEM fuel cell. The reliability of the numerical simulations is ensured by comparing the results with experimental results obtained from the literature as shown in Figure 3.

3.png

Figure 3. The comparison of polarization curve of the simulated fuel cell with an experimentally measured data

The issues of a PEM fuel cell such as heat and water management, losses due to pressure and temperature gradient have been studied using numerical simulations. The results of PEMFC simulations can help solve challenges and require accurate prediction of the effect of operating conditions on fuel cell operation. In a PEM fuel cell, all three fundamental phenomena such as heat transfer, fluid flow and electrochemical reaction are directly or indirectly dependent on mass flow rates. An increase in mass flow rate at cathode side may have positive impact on heat transfer while on the other hand a reduced mass flow rate may result in less unused fuel at the fuel channel outlet. It may also affect the water content inside the cathode flow channel.

The polarization curve is the basic performance indicator of a PEM fuel cell being simulated numerically in this study. Figure 3 shows the polarization curve based on simulation values showing cell voltage as a function of current density and agrees well experimental results [15]. This validation clearly shows the capability of the PEM fuel cell model used for current simulations. The single value of current density was achieved at each iteration performed and shows a difference of less than 10% from the experimentally determined values.

4.png

Figure 4. The impact of mass flow rate on current density with increasing operating pressure

The effect of increasing mass flow rate and pressure on current density is shown in Figure 4. The results of the simulation show that the combine effect of increasing mass flow rate and pressure has a direct impact on fuel cell performance. The current shows a linear increase with pressure from 1 atm to 3 atm.  The optimum values of mass flow rate and pressure were found to be 0.0008 g/s and 2.5 atm respectively. Any further increase in mass flow rate results in increased velocity and the heat transfer at bipolar plate also increases. This results in fewer hot areas on the anode side as compared to the cathode side that could be due to lower mass flow rate on the cathode side. On the other hand, it decreases efficiency and lowers fuel consumption rate that ultimately impacts the overall performance and also increased pressure drop across the channel. The inflection points on the curves of current density and change in mass fraction of hydrogen at the outlet, thus corresponds to optimum parameters for a fuel cell as analyzed in detailed here.

Another reason may be the recombination reaction occurring at the cathode side releasing more heat as the cathode current density increases. Any further increase in mass flow rate or operating pressure beyond the optimum values reported above show an increase in the current density to a small extent but it also resulted in more water content at cathode flow field that might cause flooding problem on the cathode side channel.

Figure 5 shows the effect of temperature on current density while keeping the mass flow rates constant. As the temperature of the fuel cell increases the current density decreases due to increased resistance to the flow of electrons. The increasing temperature may increase reaction kinetics at catalyst layer but it affects negatively when electrons travel through bipolar plate. The suitable operating temperature for PEM fuel cell when flow rate and operating pressure were kept constant was observed to be 343 k. The current density drops suddenly if temperature is further increased because benefit of higher kinetics is lost as resistances due to hot spots.

The mass fraction of hydrogen and oxygen at the fuel cell outlet is a good indicator of the performance of a PEM fuel cell. The mass fraction of hydrogen at the outlet is shown in Figure 6 at different operating pressures. The result of this numerical study shows that the mass fraction of hydrogen at the outlet decreases at a high rate with increasing pressure until 2.5 atm and then the curve starts to level out if the pressure is increased further. This shows increased consumption of hydrogen in the fuel cell that enhances the performance of the fuel cell. With the increasing pressure there is higher convection of gas through gas diffusion layer resulting in an increase of reaction rate at the catalyst layer. The lower change in mass fraction was observed at 1 atm compared to 2 atm inside anode flow channel as more fuel was consumed at higher operating pressure and decrease in mass fraction was higher. The same trend was observed inside cathode flow channel where at 1 atm less mass fraction change was observed compared to 2 atm as shown in Figure 6. The decrease in mass fraction of reactants improved from 13% to 23% as the operating pressure was changed from 1 to 2 atm. However, pressure drop inside flow channels also increased when operating pressure was raised as shown in Figure 7.

Figure 8 depicts the pressure drop inside flow field hence overall reduction during different operating conditions can be estimated from above distribution. Furthermore, it is evident from the image that during initial turns of serpentine pressure drop is less compared to end of geometry. The reason is that friction becomes dominant hence pressure drop is maximum at end of channels.

5.png

Figure 5. The impact of temperature on current density when flow rate is constant

6.png

Figure 6. The impact of operating pressure on mass fraction drop

7.png

Figure 7. The drop in oxygen mass fraction

8.jpg

Figure 8. The pressure distribution inside multiple serpentine flow channels. The values of pressure are relative to atmospheric pressure and not absolute pressure

The authors would like to thank Sukkur IBA University for providing licensed software and support for this project.