A new one dimensional steady state model for PEM fuel cell

In consideration of the thermodynamic parameters of reactants changed along the channel, a new calculation method for the output voltage of PEM (proton exchange membrane) fuel cell was introduced in this paper based on the empirical equations. Then a one dimension steady state model for PEM fuel was established based on this new calculation method. Simulation studies on the distributions of the partial pressure of gases and current density along the channel were implemented with the model. The simulation results showed that, increasing the pressure of inlet gases can increase the partial pressure differences among the sections of the channel, which would deteriorate the uniformity of current density along the channel; A proper increment of the stoichiometric ratio of inlet gases can alleviate the spatial variations of the parameters of reactants, and the uniformity of current density can be improved by such alleviation. In order to validate the effectiveness of the model, the comparison between simulation results and known experimental data was conducted, which indicates that the presented model has proper accuracy.


Introduction
Electric Proton exchange membrane (PEM) fuel cells are electrochemical energy conversion devices that directly convert the chemical energy of supplied reactants (hydrogen and oxygen) into electricity. Of the many types of fuel cells, the PEM fuel cell is considered a commercially implementable frontrunner for its low temperature operation and high efficiency. Establishing accurate numerical model for fuel cells is significant for the optimization of their design and operation. The models of fuel cell generally include mechanism models and semi-empirical models. The mechanism models based on transport and electrochemical equations can apply to any kind of fuel cell, but the computation process is time-consuming and some parameters in such model are difficult to be measured. The semi-empirical model based on the combination of mechanism model and experimental data can only use for some specific fuel cells, but the calculation is convenient and the numerical results agree with experimental results well [1][2][3][4] . According to the Butler-Volmer equation, the partial pressures of reactants at active layer directly affect the activation overpotential of fuel cell. The effect of variations of reactants along channel on the output voltage of fuel cell can be ignored when the pressure difference between inlet and outlet of the channel is small. In the case of complicated channels which have big flow resistance, it's essential to consider the effect of the reactants' variations on the fuel cell's performance. With the consideration of the change of the thermodynamic parameters of reactant gases along the channel, a new calculation method for the output voltage of fuel cell has been presented in this paper based on empirical equations. And a one dimension steady state model for PEM fuel cell was established based on this new method. The model can predict the performances of fuel cell under various operating conditions and the distributions of current density, the partial pressure of oxygen and water vapor along the channel. Using the inlet of the channel and flow direction as the origin and positive direction of x axis respectively, the fuel cell was divided into N control volumes along the x axis. Numbers on the right side of Fig.1 is the number of control volumes; numbers on the left side is the number of the section of control volumes. For example x i is the number of the outlet section of control volume i, and also is the inlet section of control volume (i+1). That is to say the outlet parameters of former control volume are the inlet parameters of latter control volume. The parameters on section x 0 and x N are boundary conditions which can be measured. Through diffusion layer, hydrogen and oxygen in control volume transfer from channel to active layer where the chemical reaction occurred. Hydrogen was decomposed into protons and electrons at anodic side, while oxygen combining with the protons and electrons delivered from anode generated water at cathodic side. Water transporting through the membrane usually caused by proton dragging, concentration gradient and pressure gradient across the membrane. As the pressure difference between anodic side and cathodic side of fuel cell is small, water transport caused by pressure difference would be ignored in this paper. Protons traveling through membrane from anode to cathode drag water molecules along with them, while water diffuses from cathode to anode due to the concentration gradient. The net water across membrane from anode to cathode can be described by the net water transfer coefficient α net , which is defined as the net number of water molecules crossing the membrane of fuel cell per proton. The value of α net is determined by the water content of membrane and current density. In order to focus on the analysis of the effect of the spatial variations of reactants on the performance of fuel cell, and without lose generality, the following assumptions are made in developing a one dimension steady state model for PEM fuel cell: a) The operation of fuel cell in the model is steady, and all parameters in cell are independent of time. b) The operating temperature of fuel cell is constant. c) All gases in the channels are viewed as ideal gas. d) No liquid water exists in the channels, single phase. e) According to Ref. [5], the net water transfer coefficient keeps constant, i.e. α net =0.1. f) The pressure gradient along channel is assumed to be constant, i.e. (dP/dx)=constant [6] . g) The spatial variations of reactants in each single control volume could be ignored.

Control equations for control volume
The amounts of hydrogen and oxygen consumed in control volume i are in direct proportion to the corresponding output current I i . When α net keeps constant, water flux transferred from anode to cathode is also in proportion to the local current density. Applying the mass conservation to control volume i can get equations as the following: where n is the molar flow rate of species in the channel, in mole/s; the superscript i is the number of control volume, the subscripts h, o, w, a, c are hydrogen, oxygen, water vapor, anode and cathode, respectively. I i is the output current of control volume i, in A; F is the Faraday constant, 96485 C/mol. According to the ideal gas law, the partial pressure of gas is in proportion to the corresponding molar fraction. Therefore the partial pressure of the species in control volume i can be shown as: where P i y , n i y are the partial pressure and molar flow rate of species y in control volume i. y stands for oxygen, nitrogen and water vapor at cathodic side, respectively. P i c , n i t are the pressure and molar flow rate of the gas mixture in control volume i. As the pressure gradient along the channel is assumed to be constant, the total pressure of gases in control volume i can be written as: where P c , P e are the inlet and outlet pressure of cathodic channel.
where T cell is the operationg temperature of fuel cell, in K; P * h , P * o are the effective partial pressure of hydrogen and oxygen at active layer, respectively. As shown in Fig.1, the control volumes are connected in parallel. Therefore the output voltage (U i ) of each control volume must be same. For control volume i, the calculation of U i can be written as: where E i r , η i act , η i ohmic are thermodynamic potential, activation overvoltage and ohmic overvoltage of control volume i respectively, in V. The values of activation loss and ohmic loss are depended on many factors. In order to simplify the computing process, empirical equations are employed to calculate the losses of fuel cell. The equations were given in Ref. [1] as: where c i o is the concentration of oxygen at active layer in control volume i, in mol/cm 3 . ξ 1~ξ7 are empirical coefficients, obtained from experimental data. Xue et al. calculated the output voltage of fuel cell by substituting the partial pressure of hydrogen/oxygen within corresponding channel into Eq. (7), and the error in thermodynamic potential E r can be compensated by modifying the empirical coefficients [2] . Henk et al. [3] used the same way to compute the output voltage of fuel cell, and the modified empirical coefficients are ξ 1 =-0.944,ξ 2 =3.54×10 -3 ,ξ 3 =8.3×10 -5 ,ξ 4 =-1.96×10 -4 ,ξ 5 =3.30×10 -3 ,ξ 6 =-7.55×10 -6 ,ξ 7 =1.10×10 -6 .Both the simplified approach and empirical coefficients in Ref. [3] were employed to deal with the output voltage in this paper, combining Eq. (8), (9) and (10), the output voltage for control volume i can given by: where a=ξ 7 / (ξ 5 +ξ 6 T cell ), b=-ξ 4 T cell / (ξ 5 +ξ 6 T cell ), c= -(ξ 1 +ξ 2 T cell +ξ 3 T cell ln(c i o ))/ (ξ 5 +ξ 6 T cell ). As mentioned above, control volumes are connected in parallel. Therefore the output voltage of each control volume is equivalent and equal to the output voltage of the whole cell. Duo to the chemical reaction and flow resistance, the partial pressures of reactants must be changed along the channel. According to Eq. (7), the thermodynamic potentials of different control volumes must be different when the operating temperature of fuel cell keeps constant. However, for a specific PEM fuel cell, the output voltage and current must be certain when the operating conditions and load were given. Assuming the output voltage and current of fuel cell are U and I, the value of the output current I i of control volume i should make the corresponding output voltage U i be equal to U according to Eq. (11). Therefore, for each control volume, the point is to find out proper output currents to ensure all corresponding output voltages are equal to U, and the sum of output currents of all control volumes is equal to I. Fig. 2 shows the detail computation procedure for the model. In order to calculate the output current of control volume 1, the output voltage for control volumes was predefined (open circuit voltage U) at first, and then the corresponding output current of control volume 1 can be obtained by solving Eq. (11). It should be noted that Eq. (11) is a nonlinear equation with one variable I i , the Newton iteration method was adopted to solve the equation in this paper. Applying the same way to other (N-1) control volumes to compute the output currents, add them together, and if the sum of N output currents is equal to the given total output current of Fig.2: Procedure of the proposed calculation method fuel cell, the predefined voltage can be accepted as a reasonable output voltage of the cell. If not, the predefined voltage should reduce a step of h, and keep doing this until the condition is satisfied. As the boundary conditions for the model, the parameters on the inlet and outlet sections of channel include inlet pressure of anode (P a ), inlet pressure of cathode (P c ), outlet pressure of channel (P e ), stoichiometric ratio of hydrogen (S a ), stoichiometric ratio of air (S c ), operating temperature (T cell ), humidification temperature (T hum ), relative humidity of humidified reactant gases (φ) and output current of cell (I).

Results and discussion
In the case of supplying air to cathode, the cathode parameters changed more obviously along the channel than the ones of anodic side. Besides, as the activation loss of fuel cell mainly caused by the oxidation of oxygen, we would focus on analyzing the effect of the spatial variations of cathode reactants on the output voltage of fuel cell. The following parameters employed constant value in this paper: P e =1.0bar, S a =1.5, T cell =353K, T hum = 323K, φ=100%. pressure of oxygen always decreased and the partial pressure of water vapor always increased along the flow path. But it should be noted that at the end of the channel the partial pressure of water vapor declined when the inlet pressure was 2.0bar. The reason is that the lower local current density at the end of channel leads to lower generating velocity of water vapor, and then leads to the declining of the partial pressure of water vapor. As shown in Fig.3, with increasing the pressure of inlet gas, the partial pressure at every channel section was increased, so as the pressure differences among them. Increasing the flow rate of inlet gas can increase the flow rate of oxygen entering cell as well as the flow rate of water vapor expelling from cell. Therefore, both the declining rate of the partial pressure of oxygen and the increasing rate of the partial pressure of water vapor along the channel would slow down due to the increment of the stoichiometric ratio of cathode flow. Fig.4 and Fig.5 show the distributions of current density along the channel when the average output current density of cell are 0.1A/cm 2 and 0.9A/cm 2 , respectively. As the figures show, the current densities are monotone decreasing along the channel under all the operating conditions. The best uniformity of current density was obtained when the pressure and the stoichiometric ratio of cathode inlet gas were 1.2bar and 3.0(△j 0.1 =0.0322A/cm 2 ,△j 0.9 =0.229A/cm 2 ), while the worst one occurred when the pressure and the stoichiometric ratio of cathode inlet gas were 2.0bar and 2.0(△j 0.1 =0.0688A/cm 2 , △j 0.9 =0.506 A/cm 2 ). The uniformity of current density mainly depended on the uniformity of the distribution of reactant gases. As mentioned above, the increment of the pressure of inlet gas could increase the pressure differences among the control volumes, and such kind of differences could be decreased by increasing the flow rate of reactant gas. Therefore, it can be concluded that the uniformity of current density would be improved by increasing the stoichiometric ratio of inlet gas, and would be deteriorated by increasing the pressure of inlet gas. The deterioration of the uniformity of current density would lead to an obviously lower current density at the end of channel, that's the reason why the partial pressure of water vapor declined in Fig.3. Excessively nonuniform distribution of chemical reaction in fuel cell could cause the emerging of local hot spot, which would lead to irreversibly damage on the membrane. Consequently, in order to ensure the safety operation of fuel cell, a proper increment of the flow rate of inlet gas should be employed to improve the uniformity of current density when the operating pressure and output current of cell are high.

Polarization of curves
According to the Butler-Volmer equation, declining the effective partial pressure of oxygen at active layer could increase the activation losses of fuel cell when other operating parameters are kept constant. Since the partial pressure of reactant gases in each control volume could be increased with the increment of inlet pressure and stoichiometric ratio of cathode flow, as Fig.6 shows, the performance of cell was improved by increasing the pressure and flow rate of air. referred in Ref. [3] were used as the model's boundary conditions. And the results obtained from the model had been compared with the experimental data, see Fig.7. It indicates that the presented model has proper accuracy.

Conclusions
Duo to the chemical reaction and flow resistance, the thermodynamic parameters of reactants must be changed along the channel. In order to study the effect of the spatial variations of reactant gases on the performance of PEM fuel cell, a new calculation method for the output voltage of fuel cell was presented in this paper, and a one dimensional steady state model was developed based on this new method. The simulation results of the model indicate that, increasing the pressure and stoichiometric ratio of cathode flow could increase the partial pressure of reactants along the channel, which could improve the performance of fuel cell; the variations of the partial pressure of reactants along the channel could be increased by increasing the pressure of cathode flow, and the uniformity of current density would deteriorated by the increased variations; increasing the flow rate of cathode flow could lessen the partial pressure differences of oxygen in channel, and the uniformity of current density would be improved by the lessened differences. The simulation results of the model agree with known experimental data well, which indicates that the presented model not only saves computing time, but also has proper accuracy.