## Abstract

In this research, the impact of structural parameters of bipolar plates on the proton exchange membrane (PEM) fuel cell performance has been investigated using numerical method, and this model incorporates all the essential fundamental physical and electrochemical processes occurring in the membrane electrolyte, cathode catalyst layer, electrode backing, and flow channel, with some assumptions in each part. In formulation of this model, the cell is assumed to work under steady state conditions. Also, since the thickness of the cell is negligible compared to other dimensions, one-dimensional and isothermal approximations are used. The structural parameters considered in this paper are: the width of channels (W_{c}), the width of support (W_{s}), the number of gas channels (n_{g}), the height of channels (h_{c}), and the height of supports (h_{p}). The results show that structural parameters of bipolar plates have a great impact on outlet voltage in high current densities. Also, the number of gas channels, their surface area, the contacting area of bipolar plates, and electrodes have a great effect on the rate of reaction and consequently on outlet voltage. The model predictions have been compared with the existing experimental results available in the literature, and excellent agreement has been demonstrated between the model results and the experimental data for the cell polarization curve.

## Introduction

The proton exchange membrane fuel cell (PEMFC) is an electrochemical energy converter that converts chemical energy of fuel directly into DC electricity. A fuel cell using a very thin polymer membrane as an electrolyte has been considered as a promising candidate for future power sources, especially for transportation applications and residential power, due to their high efficiency, high power density, quick startup, quiet operations, and, most importantly, zero emissions, thereby reducing air pollution and greenhouse gas emissions.

The first demonstration of a fuel cell was by lawyer and scientist William Grove in 1839 [1], although it appears that a Swiss scientist Christian F. Shoenbein independently discovered the very same effect at about the same time (or even a year before) [2]. Optimizing of the geometrical parameters is one of the ways to minimize ohmic losses. Some studies have been developed to investigate the effect of using different channel configurations on cell performance. Sun et al. [3] developed a numerical study, which suggested using a trapezoidal channel cross-section rather than a rectangular or square cross-section to improve cell performance. Chiang and Chu [4] conducted numerical simulations of a 3D isothermal model of a PEMFC for various channel configurations (channel length, height, and shoulder width), with particular emphasis on the effect of channel height on cell performance. They carried out simulations with different combinations of channel width and height, while maintaining the same channel cross-sectional area and constant flow rate. Their results indicated that flat channels (i.e., smaller height) gave better performance.

Yi and Nguyen [5] developed an along-the-channel model for evaluating the effects of various designs and operating parameters on the performance of a PEM fuel cell. The results showed that the humidification of the anode gas is required to improve the conductivity of the membrane and the liquid injection and higher humidification temperature can enhance the cell performance. Ge and Yi [6] developed a two-dimensional model to investigate the effects of operation conditions and membrane thickness on the water transport. The results revealed that the cell performance can be enhanced by increasing the cell temperature. Yi and Nguyen [7] employed a two-dimensional isothermal model of a porous electrode to simulate the hydrodynamics of gas flow through the pore volume of the electrode of a PEMFC. It is concluded that, from the predictions of the PEM fuel cells with interdigitated flow field, the higher gas flow rate through the electrode improves the electrode performance. He et al. [8] and Kumar and Reddy [9,10] studied the influences of electrode and flow field design on the performance of a PEM fuel cell with a half-cell model. It was found that higher differential pressure between inlet and outlet channels would enhance the electrode performance. Kim et al. [11] developed a curve-fitting scheme based on experimental data of the first polarization curve of a PEM fuel cell. Ahmed and Sung [12] performed the simulations of PEMFCs with a new design for the channel shoulder geometry, in which the membrane electrode assembly is deflected from shoulder to shoulder. Bernardi and Verbrugge [13,14] used an analytical approach. They developed a mathematical model of a PEMFC from fundamental transport properties. Where the losses incurred by the activation over potential of the anode and cathode reaction, the ohmic losses incurred by the membrane and the ohmic losses due to the electrodes are subtracted from the reversible cell voltage. Marr and Li [15,16] developed a simplified engineering model of a PEMFC based on the catalyst layer of Weisbrod et al. [17] and the membrane model of Bernardi and Verbrugge [13,14]. This is the model that will be used to incorporate the effect of bipolar plate structural parameters on the PEMFC performance. Marr and Li [15] investigated the composition and performance optimization of cathode catalyst platinum and catalyst layer structure in a proton exchange membrane fuel cell by including both electrochemical reaction and mass transport process. They found that electrochemical reactions occur in a thin layer within a few micrometers thick, indicating ineffective catalyst utilization for the present catalyst layer design.

Also Baschuk and Li [18] investigated a polymer electrolyte membrane fuel cell with variable degrees of water flooding of the membrane.

Therefore, the issue of bipolar structural parameters in a PEMFC has not been fully examined in the previous work. In the following sections, the formulation of the model will be presented. The result is then compared to experimental data available in the literature for validation.

## Model Formulation

The operation of a PEM fuel cell is based on converting the chemical energy of the fuel, such as hydrogen, and oxidant, such as oxygen, into electrical energy, as shown in Fig. 1. The following electrochemical reactions happening on the anode and the cathode are the base of fuel cell function [19].

More precisely, the reactions above happen on a border between the ionically conductive electrolyte and electrically conductive electrode. The electrodes must be porous, allowing the gases to arrive to, as well as product water to leave, the reaction sites, because there are gases involved in fuel cell electrochemical reactions.

In the formulation of this model, the cell is assumed to work under steady state conditions. Also, since the thickness of the cell is negligible compared to its other dimensional, one-dimensional and isothermal approximation is used.

where *E _{r}* is the reversible cell voltage and η

_{act}is the loss due to resistance to the mass transfer limitations and electrochemical reactions in the cathode catalyst layer.

_{ohmic,e}, η

_{ohmic,p}, and η

_{ohmic,m}are the voltage losses due to the ohmic resistance of the electrodes, flow channel plate, and membrane layer. E

_{r}is calculated from a modified version of the Nernst equation, which is modified with an extra term to account for changes in temperature from the standard reference temperature [20]. This is given by Eq. (4).

*G*is the change in Gibbs free energy,

*F*is the Faraday constant, Δ

*S*is the change in the entropy, and

*R*is the universal gas constant, while $PH2$ and $PO2$ are the partial pressures of hydrogen and oxygen in the atmosphere. Taking into account the quantity of parameters, Eq. (4) is simplified to the form below. Table 1 shows the operating and structural parameters.

### Electrode.

where *h _{p}* is the thickness of the solid portion of the plate and

*h*is the height of the flow channels and supports.

_{c}*W*and

*L*are the width and length of the cell, respectively, and

*w*is the width of the channels supports.

_{s}*w*(width of support) and

_{s}*w*(width of channel) can be determined by the equation below,

_{c}where *w*_{c} refers to the width of the flow channels and *w _{s}* is the width of support, while

*n*is the number of the flow channels in the cell.

_{g}Parameters | Notation | Value |
---|---|---|

Temperature, °C | T | 60 |

Gas mixture pressure, atm | p | 5 |

Oxygen mole fraction | Xo_{2} | 0.21 |

CL thickness, μm | L_{C} | 50 |

GDL porosity, vol. % | ɛ_{g} | 40 |

Volume fraction of membrane in the CL, vol. % | L_{m, c} | 0.4 |

Volume fraction of GDL in the CL, vol. % | L_{g, c} | 0.1 |

Mass loading of platinum per unit area of the cathode, K gm^{−2} | m_{pt} | 0.01 |

Mass loading of carbon per unit area of the cathode, K gm^{−2} | m_{c} | 0.003 |

Platinum density, kg m^{−3} | P_{pt} | 21,400 |

Carbon density, kg m^{−3} | P_{c} | 1800 |

Reference oxygen concentration, mol m^{−3} | Co_{2} | 1.2 |

Cathodic transfer coefficient | d_{c} | 1 |

Anodic transfer coefficient | d_{a} | 0.5 |

Bulk protonic conductivity, (Ω m)^{−1} | k | 17 |

Bulk electronic conductivity, (Ω m)^{−1} | k_{s} | 7.27 × 10^{4} |

Oxygen diffusion coefficient within the liquid water, m^{2} s^{−1} | D_{02 –w} | 9.19 × 10^{−9} |

Parameters | Notation | Value |
---|---|---|

Temperature, °C | T | 60 |

Gas mixture pressure, atm | p | 5 |

Oxygen mole fraction | Xo_{2} | 0.21 |

CL thickness, μm | L_{C} | 50 |

GDL porosity, vol. % | ɛ_{g} | 40 |

Volume fraction of membrane in the CL, vol. % | L_{m, c} | 0.4 |

Volume fraction of GDL in the CL, vol. % | L_{g, c} | 0.1 |

Mass loading of platinum per unit area of the cathode, K gm^{−2} | m_{pt} | 0.01 |

Mass loading of carbon per unit area of the cathode, K gm^{−2} | m_{c} | 0.003 |

Platinum density, kg m^{−3} | P_{pt} | 21,400 |

Carbon density, kg m^{−3} | P_{c} | 1800 |

Reference oxygen concentration, mol m^{−3} | Co_{2} | 1.2 |

Cathodic transfer coefficient | d_{c} | 1 |

Anodic transfer coefficient | d_{a} | 0.5 |

Bulk protonic conductivity, (Ω m)^{−1} | k | 17 |

Bulk electronic conductivity, (Ω m)^{−1} | k_{s} | 7.27 × 10^{4} |

Oxygen diffusion coefficient within the liquid water, m^{2} s^{−1} | D_{02 –w} | 9.19 × 10^{−9} |

where *I*_{δ} is the cell current density.

_{e}is the thickness of the electrode. Effective resistivity of the electrode, $\rho R,eeff$, can be derived from the bulk resistivity of the electrode, ρ, and the void fraction of the electrode, ϕ, by the Bruggeman’s correction [20].

### Catalyst.

The catalyst layer is modeled with a set of coupled differential equations. In modeling the cathode catalyst layer, it is assumed to be isothermal, one-dimensional, and uniformly distributed. In this model, the other important assumption is that the negated space within the catalyst layer is large enough so that Knudsen diffusion is unimportant [18].

where *I* is the protonic current density, *C* is the concentration of oxygen, η_{act} is the over potential caused mostly by the resistance to the electrochemical reactions and due to the finite rate of mass diffusion, and *z* is the distance into the catalyst layer measured from the electrode/catalyst layer interface. *i*_{o,ref} denotes the reference current density that is a function of the cell temperature. Also, it is an experimentally derived parameter and is given as a function of temperature for Nafion on platinum in Parthasarathy et al. [21]. Also, γ, reaction order, can be found analytically from the procedure in Newman [22], with a value of 0.5 resulting.

*C*

_{ref}is the reference oxygen concentration, and it is related to

*i*

_{o,ref}and, in this model, has a value of 12 mol/m

^{3}. α

_{c}and α

_{a}denote the cathodic and anodic transfer coefficient, which, in this model, have the values of 1 and 0.5, respectively.

*A*

_{v}is the specific reaction surface and can be derived from the below equation [15],

*m*_{pt} is the catalyst mass loading per unit area of the cathode, δ_{c} is the thickness of the catalyst layer, and *A*_{s} is the catalyst surface area per unit mass of the catalyst.

*K*

_{m}and

*K*

_{s}are the bulk conductivities of the membrane and solid catalyst,

*l*

_{m}is the fraction of membrane in the catalyst layer, and ϕ

_{c}is the void fraction of the catalyst layer, which can be calculated as Eq. (19) depending upon the type and thickness of the catalyst layer.

ρ_{pt} and ρ_{c} are the density of the platinum and its carbon support, and *f _{pt}* represents the amount of platinum catalyst on its carbon support.

*I*is the protonic current density and $DO2eff$ denotes the effective diffusion coefficient of oxygen in the catalyst layer, which is determined as the following equation [23]:

In this model, the mass transfer of oxygen from the flow channel to the reaction site of the cathode catalyst layer will determine the shape of the polarization curve in the concentration over the potential region [18].

### Membrane.

In this equation, δ_{m} is the thickness of the membrane, *K*_{m} is the electrical conductivity of the membrane, which is a function of cell temperature, *K*_{p} is the hydraulic permeability, *K*_{E} is the electrokinetic permeability, ΔP_{a-c} is the pressure differential across the membrane, *C*_{H}^{+} is the fixed charge concentration, and $\mu H2O$ denotes the viscosity of liquid water. The values for the *K*_{m}, *K*_{p}, *C*_{H}^{+}, and *K _{E}* for Nafion can be found in Bernardi and Verbrugge [13,14].

### Bipolar Plate.

In a typical frame and PEM cell, the bipolar plate is a plate of graphic with a serpentine flow channel cut into it. These channels are used to deliver reactants to the electrode and to collect current from the electrode. Two different phenomena are modeled in the plate:

- (1)
Mass transport in the flow channel

- (2)
Electric resistance in the plate

#### Channel Flow.

_{2}O and N

_{2}by using a bulk diffusion correction. It is assumed that there is a constant concentration of O

_{2}at the electrode surface. Additionally, inlet and secondary flow patterns at the corners of the channel were ignored. These assumptions will give more conservative predictions, because adding in correlations for inlet and flow around the corners would increase the mass transfer rate. The mass transport has the form

*N*is the molar flow rate and Δ

*C*is the logarithmic average of concentration change between bulk gas concentration and electrode surface. The exposed area of the electrode on the channel is symbolized with

_{ln}*A*, the diffusion coefficient,

_{c}*D*

_{ibulk}, is the bulk-corrected diffusion, and

*d*is the hydraulic diameter of the channel. However, the Sherwood number,

_{h}*Sh*, is a dimensionless parameter that is defined as

*h*is the mass transfer coefficient (mol/m

_{m}^{2}.s) and

*D*is the diffusion coefficient (m

_{ij}^{2}/s) of the species

*i*in

*j*. For the flow channel, the Sherwood number was derived by the analogy between heat and mass transfer. The Nusselt number correlation [26] for a three-sided adiabatic square duct with one constant heat flux wall was converted into a Sherwood number by using the following conversion formula [25]:

In this equation, Δ*C*_{1} = *C*_{s} – *C*_{inlet} represents the inlet concentration difference and Δ*C*_{2} = *C*_{s} – *C*_{outlet} represents the outlet concentration difference.

*W*equal to the total width of the end supports and the width of the channel. The length is measured from the center of the channel. For the distance in the

*W*direction, we have

*W*is the width of the cell,

*l*is the length of the channel in the

_{w}*W*direction, and

*w*and

_{s}*w*are the widths of the support and channel, respectively. The channel length in the L direction is more complicated, because the inlets and outlets are a different length than the inner channels. The length of the inlet or outlet is

_{c}*L*direction is

For the second group of terms, it is assumed that *w _{s}* has the same distance with respect to both length and width. Also, ohmic resistance, which is calculated in Sec. 2.1, is assumed constant within the temperature range 300–370 K, and diffusion is assumed to spread O

_{2}evenly throughout the electrode surface such that there is a constant O

_{2}concentration at the electrode catalyst layer interface.

### Boundary Conditions.

*P*is the gas mixture pressure at the GDL/CL interface. The Henry constant Ho

_{2}(in atm m

^{3}mol

^{−1}) is typically assumed to be a function of temperature [13] according to

*T*is measured in degrees Kelvin. Because the GDL is very thin, the pressure conditions at the GDL/CL interface are set to the cathode channel values, which, in this model, have been taken to be Xo

_{2}= 0.21 and P = 5 atm. The boundary condition for the protonic current density comes from the assumption that all protons are consumed before they reach the GDL/CL boundary, so that

*z*=

*l*, the protonic current density at the membrane approaches its ultimate value (i.e., the cell current density,

_{c}*I*

_{δ}), so that

The procedure for solving the formulated model to yield a cell voltage for given structural and operational parameters is as follows: first, the reversible voltage of the cell is calculated for the reference temperature and pressure (T = 60 °C, P = 5 atm) using Eq. (5). Then, the losses are calculated in all parts of the cell for calculating overall cell voltage in Eq. (3). Also, for calculating the catalyst influence on overall cell voltage, it is needed to solve ordinary differential equations consisting of Eqs. (14), (16), and (21) along with boundary conditions, which are stated for the catalyst layer in Eqs. (38), (41), and (43) that represent a two-point boundary value problem. The built-in shooting method algorithm within the boundary value problems function of the Matlab software package has been used for the present computation [27].

## Validation

The polarization curve is obtained from the present calculations in base case condition (T = 353.15 K, P = 5 atm), except with *m _{Pt}* = 0.35 mg/cm

^{2}(f = 10 Pt/C Mass %) and

*R*

_{Ohmic}= 0.225 (Ω/cm

^{2}), and is compared with the experimental Ticianelli et al. [28]. As shown in Fig. 4, the results predicted by the numerical are in good agreement with the experimental data.

## Results

The study on the structural parameters of the bipolar plates, such as *W _{c}*,

*W*, ϕ,

_{s}*n*,

_{g}*h*, and

_{p}*h*, which represent the channel width, support width, void fraction, height of support, and channel, respectively, has been done by using the Marr [15] model (Fig. 2). The amounts of other parameters are listed in Table 1.

_{c}Considering any structural changes in channels makes change in the amounts of pressure losses values, so it was assumed that its losses will be relieved by external tools and the pressure will be constant(5 bar). By comparing the results obtained from this research with other models and other sizes of the PEMs, it is revealed that these results are extendable for other fuel cells with any dimensions and structural and operational parameters.

### Width of Channels (*W*_{c}).

_{c}

According to obtained results, when the other parameters are constant and *W _{c}* is changed, the outlet voltage will be changed according to Fig. 5. As shown in this diagram, it is obvious that, when the number of channels or grooves in the plate is constant, the voltage increases with the decrease of

*W*.

_{c}Of course, this change is very insignificant and, according to Fig. 6, the changes in W_{c} are more sensible in high current densities. It can be a good manner to increase total voltage in high current densities.

Figure 5 shows the change of voltage in variable sizes of channel width in *n _{g}* = 12 and

*I*

_{δ}= 10 000 (A/m

^{2}). Figure 6 shows the voltage and power difference when W

_{c}changes from 0.003 (m) to 0.00001 (m) and other parameters are constant.

More surface area is necessary in high current densities; however, when the width of channels decreases, the support area will increase.

The increment of the support area will change the void fraction of the electrode as compression decreases and consequently change the ability of reactant to diffusion to the catalyst surface directly underneath the support changes. So, with decreasing W_{c}, the ohmic resistance will decrease and also changes of void fraction of the electrode will have influence on cell voltage. The electrode void fraction influence on cell voltage will be discussed in Sec. 4.2.

### Void Fraction.

Electrode void fraction is one of the important structural parameters that has influence on cell voltage. According to Marr [20], when void fraction in the electrode increases, mass transfer increases as well and, similarly, electrode resistance increases. The corresponding increase in resistivity from the increase of the electrode void fraction is swamped by increasing performance from better gas diffusion. As shown in Fig. 7, the increasing of the void fraction will lead to enhancement of cell voltage.

### Width of Support (*Ws*).

According to Eq. (8), it can be seen that, for a constant number of channels in the bipolar plate with defined dimensions, W_{s} is increased as W_{c} decreases. So the outlet voltage is increased by increase of W_{s}, and this increase in voltage is due to the decrease of ohmic resistance in addition to changes in void fraction, like in Sec. 4.2. Also, this increase of outlet voltage, due to increase of W_{s}, is more obvious in high current densities.

### Number of Channels (*n*_{g}).

_{g}

It was seen that the outlet voltage is increased because of decrease in *W _{c}*. Also, it is clear that, in high current densities, the effects of structural parameters are more sensible and more area of reaction is needed in these high current densities. Therefore, by increase of

*n*, the surface of contact can be increased, but on the other side, the width of supports will decrease by the increase of the number of channels. This causes changes in the void fraction and, accordingly, the ohmic resistance of the plate. Therefore, the number of channels must have an optimized amount. It has been clearly shown in Fig. 10. According to the voltage-

_{g}*n*diagram, it is obvious that, when n

_{g}_{g}increases in constant amounts of

*W*and

_{c}*I*

_{δ}, the voltage increases and has a maximum amount in a definite amount for each size of

*W*.

_{c}### Height of Channels (h_{c}).

Through obtained data by numerical solution method and according to Fig. 11, it is clear that, as h_{c} (the height of channels shown in Fig. 2) is decreasing, the voltage is increased. This is due to the decrease of ohmic resistance. In this diagram, the number of channels used for the transfer of gas, the width of channel, and mass flow rate of channel are constant and *h _{c}* is variable.

Also, in Fig. 12, the polarization curve for two different h_{c} demonstrates the increase of voltage in higher current densities because of decreasing height of channels. That change in h_{c} has more impact on outlet voltage in higher current densities.

### Height of Support (*h*_{p}).

_{p}

By study effect of the height of supports, which is shown in Fig. 2, it is evident that, by reducing h_{p}, the outlet voltage increases due to reduction of ohmic resistance. Figure 13 shows the influence of *h _{p}* on the outlet voltage.

Other structural and operational parameters remain constant in the investigation of h_{p} parameter.

Figure 14 shows the voltage and power difference when h_{p} changes from 0.009 (m) to 0.00001 (m) and other parameters are constant.

Furthermore, it can be seen that the effect of the channel height is more sensible in high current densities, and according to Fig. 15, the effect of channel height is more sensible than the height of support on cell voltage.

## Conclusion

By investigation of the isothermal, steady, and one-dimensional PEM fuel cell, it is shown that the present model is in excellent agreement with experimental data. Also, it was obtained that structural parameters of bipolar plates and the method of their designing and construction has impact on outlet voltage in high current densities. Also, the result indicates decreasing width of the gas channels, increasing outlet voltage. Similarly, it is evident that the height of channel, the height of support, and the number of gas channels has influence on outlet voltage. Also, the result shows that the contacting area of bipolar plates and the electrode and surface area of gas channels have great effect on the rate of reaction and, consequently, on cell voltage. Changes in the number of gas channels show that there are an optimum number of them for each size of channel width. Also, the increase in number of bipolar plate channels is not a benefit to reach maximum power density and depends on other parameters like contact area or width of channels and width of supports.