Numerical Analysis on Temperature Distribution in a Single Cell of PEFC Operated at Higher Temperature by1D Heat Transfer Model and 3D Multi-Physics Simulation Model

Akira Nishimura^*, Kyohei Toyoda, Daiki Mishima, Eric Hu

Division of Mechanical Engineering, Graduate School of Engineering, Mie University, Tsu, Japan.

DOI: 10.4236/epe.2023.155010                

Abstract

This study is to understand the impact of operating conditions, especially initial operation temperature (T_ini) which is set in a high temperature range, on the temperature profile of the interface between the polymer electrolyte membrane (PEM) and the catalyst layer at the cathode (i.e., the reaction surface) in a single cell of polymer electrolyte fuel cell (PEFC). A 1D multi-plate heat transfer model based on the temperature data of the separator measured using the thermograph in a power generation experiment was developed to evaluate the reaction surface temperature (T_react). In addition, to validate the proposed heat transfer model, T_react obtained from the model was compared with that from the 3D numerical simulation using CFD software COMSOL Multiphysics which solves the continuity equation, Brinkman equation, Maxwell-Stefan equation, Butler-Volmer equation as well as heat transfer equation. As a result, the temperature gap between the results obtained by 1D heat transfer model and those obtained by 3D numerical simulation is below approximately 0.5 K. The simulation results show the change in the molar concentration of O₂ and H₂O from the inlet to the outlet is more even with the increase in T_ini due to the lower performance of O₂ reduction reaction. The change in the current density from the inlet to the outlet is more even with the increase in T_ini and the value of current density is smaller with the increase in T_ini due to the increase in ohmic over-potential and concentration over-potential. It is revealed that the change in T_react from the inlet to the outlet is more even with the increase in T_ini irrespective of heat transfer model. This is because the generated heat from the power generation is lower with the increase in T_ini due to the lower performance of O₂ reduction reaction.

Keywords

PEFC, Heat Transfer Model, Temperature Distribution, Numerical Simulation, High Temperature Operation

1. Introduction

According to the Japanese Energy and Industry Technology Development Organization (NEDO) road map 2017 in Japan, a high-temperature operation such as 363 K and 373 K is requested for the stationary and mobile application use of polymer electrolyte fuel cell (PEFC), respectively, during the duration from 2020 to 2025 [1] . On the other hand, the PEFC system using Nafion membrane as a polymer electrolyte membrane (PEM) is usually operated under 353 K [2] [3] [4] . When PEFC system is operated at a higher temperature such as 363 K and 373 K, we can obtain the following merits: 1) the kinetics improvement of the electrode; 2) for the vehicle usage, the cooling system can be smaller because of the increase in the temperature gap between PEFC stack system and the coolant; and 3) the durability enhancement of CO contained in the H₂ reformed from hydrocarbon [5] . However, we should solve the following demerits: 1) damage of PEM; 2) electrode elution; 3) performance drop due to uneven distribution of gas flow, pressure, temperature, voltage and current in PEFC [6] . It can be believed that the even distribution of H₂, O₂, H₂O, temperature and current density provide not only the higher power generation performance but also the longer lifetime when we operated the PEFC system at higher temperature [6] [7] .

The temperature distribution in a single cell of PEFC is crucial to the performance of PEFC. Uneven temperature distribution could cause degradation of PEM and catalyst layer. Localized temperature rise would cause thermal decomposition of PEM. PEM could also be broken by thermal stress caused by an uneven temperature distribution [8] [9] . Therefore, it is important to understand the temperature distribution in a single cell of PEFC in order to improve the power generation performance and realize the long life span, which is the aim of this study.

According to the literature survey, some studies have been conducted on high temperature of PEFC (HTPEFC), which focuses on the development of components consisting of PEFC.

As to PEM, several studies investigated to develop a new material for HTPEFC. The PEFC using phosphoric acid-developed polybenzimidazole membrane could be operated from 393 K to 433 K [10] . From this report, the power density was 0.254 W/cm², 0.299 W/cm² and 0.389 W/cm² at the current density of 0.7 A/m², 0.8 A/m² and 0.9 A/m², respectively when operated at 393 K, 413 K and 433 K, respectively, resulting from the improvement of proton conductivity at higher temperature operation. The other study has developed the polybenzimidazole/graphene oxide composite membrane [11] . The PEFC stack consisting of 12 individual cells equipped with the membrane performed the power density of 3.6 W/cm² at the current density of 0.67 A/cm² when operated at 433 K, resulting from the improvement of proton conductivity at higher temperature operation.

As to catalyst layer, some studies investigated the structure and preparation process. The catalyst layer having different microstructures and the effect of Pt loading on the performance and degradation of HTPEFC was investigated, reporting that the mass transfer was affected remarkably by the impacts of microstructures and Pt loadings [2] . From this report, the catalyst preparation process is important to obtain the higher power generation characteristics of HTPEFC [12] . The performance of membrane electrode assembly (MEA) having the anode electrode modified by Pt pulse electrodeposition was 437.2 mW/mg-Pt, which was almost 1.36 times higher than that of the pristine MEA.

As to gas diffusion layer (GDL), some studies investigated the structure such as porosity and thickness. The numerical study revealed that the effect of uneven porosity distribution was more considerable when the current densities were higher [13] . From this report, the reaction kinetics were hardly affected by changing the porosity configurations. The other numerical study revealed the thickness and porosity exhibited the opposite impact on diffusion flux, which reduced with the increase in GDL thickness but increase with the increase in porosity [14] . According to this report, the optimum thickness for anode GDL and cathode GDL would be 50 μm - 120 μm and 140 μm - 170 μm respectively, and the optimum value for GDL porosity ranged from 35% to 45%.

As to separator, several studies investigated the structure. The previous study carried out 3D numerical simulation to understand the effect of interdigitated flow field on not only the mass transfer characteristics but also the power generation characteristics [15] . According to this report, the interdigitated flow provided not only the better power generation performance compared with the parallel flow but also the similar characteristics as the serpentine flow. Additionally, there was the optimum ratio of channel to rib to obtain the higher power density. The other numerical study focused on the rib design [16] . The ratio of channel to rib influenced the distributions of gas diffusion, electron conduction and current density in the porous electrodes significantly. Moreover, the optimum ratio of channel to rib was 1 which provided the peak power density of 0.428 W/cm² and the current density of 1.2 A/cm². The widths of top and bottom edges of the anode and cathode flow channels were investigated as an independent variable with a constrained range for the optimization of the performance of HTPEFC [17] . From this report, the trapezoidal structure of cross-sectional area of the flow channel was the best shape to obtain the highest power generation performance. It also revealed that the pressure drop and the outlet power of the optimal model were larger by 1.7% and 6.5% than those of the original model at 0.4 V, respectively.

However, only a few papers [16] [18] investigated the temperature distribution near the interface between PEM and catalyst layer at the cathode, which is defined as a reaction surface in this study, excluding other studies by the authors [19] - [24] . The authors’ studies [19] - [24] investigates the effect of PEM’s thickness, GDL’s thickness and separator’s thickness on the distribution of the temperature at the reaction surface (T_surf), in a single cell of PEFC at a higher temperature such as 363 K and 373 K by 1D heat transfer model using the experimental temperature distribution data obtained by means of a thermograph. However, this model investigated the heat transfer phenomena in a single cell of PEFC only. Therefore, it is necessary to compare the temperature distribution which is obtained considering the mass transfer phenomena and the electrochemical reaction as well as heat transfer phenomenon in order to verify the heat transfer model proposed by the authors.

The aim of this study is to clarify and verify the distribution of T_react at higher temperatures, i.e. 363 K and 373 K calculated by 1D heat transfer model proposed by the authors. This study carries out the numerical simulation using a 3D model by COMSOL Multiphysics composed of multi-physics simulation codes considering the mass transfer phenomenon, the electrochemical reaction and heat transfer phenomenon to verify the distribution of T_react at higher temperatures. If we can verify 1D heat transfer model by 3D model, it can be said that the 1D heat transfer model is effective to predict the distribution of T_react without complex calculation and long calculation time. The operation temperature is changed by 353 K, 363 K and 373 K. As to 353 K, this study has selected it to exhibit the characteristics at a standard operating temperature condition compared with the characteristics at a higher temperature. The relative humidity (RH) of supply gas at anode of 80 %RH and cathode of 80 %RH (A80%RH, C80%RH), that at anode of 80 %RH and cathode of 40 %RH (A80%RH, C40%RH), that at anode of 40 %RH and cathode of 80 %RH (A40%RH, C80%RH) and that at anode of 40 %RH and cathode of 40 %RH (A40%RH, C40%RH) is also investigated. The distributions of O₂, H₂O and current density on the interface between PEM and catalyst later at the cathode, which is obtained by 3D model, are investigated to support the distribution on the temperature distribution.

2. Calculation Procedure

2.1. 1D Multi-Plate Heat Transfer Model

Figure 1 illustrates the multi-plate single cell of PEFC module (1D) used in this study. In the module, the separator’s back is the opposite side of the surface contacting the GDL. The separator’s back surface temperature T_surf,c and T_surf,a were measured using thermograph.

The heat transfer across the module is assumed to be in 1D direction only. In the module, the cell is divided into a gas channel and a rib part. In Figure 1, the upper and the lower parts represents rib part and channel part, respectively. For

Figure 1. 1D multi-plate heat transfer model.

both parts, the heat transfer was assumed to be in the through-plane direction. The reaction heat generated on the reaction surface is transferred to the cathode and anode sides separately. Although the gas flowing through the gas channel from the inlet to the outlet of the cell carries away some heat, the amount of heat taken is less than 1% of the estimated reaction heat of approximately 20 W [25] . Therefore, the heat carried away by the gas flow was neglected in this model. Additionally, the mass flow rate of gas flowing through the gas channel is very small ranging from 10⁻⁸ to 10⁻⁶ kg/s, resulting that the thermal conduction of gas in the gas channel is assumed since the gas is thought to be static.

2.2. Heat Generation Rate by Reaction

The heat generation rate H_react as a reaction product is calculated as follows:

$H_{\text{react}}=E_{\text{i}}-W_{\text{E}}$ (1)

where E_i is the ideal (total) energy generation rate by the water formation from H₂ and O₂ based on higher heating value except the initial temperature of cell (T_ini) = 373 K. The lower heating value is adopted for T_ini = 373 K. W_E is the electric work generated by PEFC. E_i and W_E are expressed as follows:

$E_{\text{i}}=m_{{\text{H}}_{\text{2}}}\times q_{\text{HHV}}\text{or}{q}_{\text{LHV}}$ (2)

$W_{\text{E}}=I\times V$ (3)

where I is the load current obtained by the experiment (=20 A). In this study, the power generation data from a load current of 20 A (=0.80 A/cm²) were used for the heat transfer modeling. $m_{{\text{H}}_{\text{2}}}$ is the molar flow rate of supplied H₂, which is equal to the ideal reaction consumption rate of H₂ required for the generation of 20 A, i.e., the stoichiometric ratio of 1.0. Here, the stoichiometric ratio is the ratio of the feed amount of H₂ or O₂ to that required to generate a current of 20 A. The flow rate of the supply gas (H₂) at the stoichiometric ratio of 1.0 is defined as follows.

$m_{{\text{H}}_{\text{2}}}=I/nF$ (4)

where $m_{{\text{H}}_{\text{2}}}$ is the molar flow rate of the supplied H₂ [mol/s], n is the valence of the ion (=2 for H₂) [−], and F is the Faraday constant (=96500) [C/mol]. $m_{{\text{O}}_{\text{2}}}$ is the molar flow rate of the supplied O₂ [mol/s] and is calculated as follows:

${\text{H}}_{\text{2}}+1/2{\text{O}}_{\text{2}}={\text{H}}_{\text{2}}\text{O}$ (5)

The actual stoichiometric ratio of the supply gas was confirmed, using the mass flow controller installed at the inlet of the single cell and the mass flow meter installed at the outlet of the cell in the power generation experiment [26] .

2.3. Heat Balance Equations for Calculating Reaction Surface Temperature

The reaction heat at rib and channel are expressed by the following equations:

$H_{\text{rib},\text{c}}=K_{\text{rib},\text{c}}A\left(T_{\text{react},\text{rib}}-T_{\text{surf},\text{c}}\right)/\text{2}$ (6)

$H_{\text{chan},\text{c}}=K_{\text{chan},\text{c}}A\left(T_{\text{react},\text{chan}}-T_{\text{surf},\text{c}}\right)/2$ (7)

$H_{\text{rib},\text{a}}=K_{\text{rib},\text{a}}A\left(T_{\text{react},\text{rib}}-T_{\text{surf},\text{a}}\right)/2$ (8)

$H_{\text{chan},\text{a}}=K_{\text{chan},\text{a}}A\left(T_{\text{react},\text{chan}}-T_{\text{surf},\text{a}}\right)/\text{2}$ (9)

$H_{\text{react}}=H_{\text{rib},\text{c}}+H_{\text{chan},\text{c}}+H_{\text{rib},\text{a}}+H_{\text{chan},\text{a}}$ (10)

where A is the heat transfer area, which is the active are of MEA (i.e., power generation area = 0.0025 m²). The overall heat transfer coefficients K_rib,c, K_chan,c, K_rib,a and K_chan,a are defined as follows:

$1/K_{\text{rib},\text{c}}=d_{\text{cat}}/k_{\text{cat}}+d_{\text{GDL}}/k_{\text{GDL}}+d_{\text{rib}}/k_{\text{rib}}+d_{\text{sep}}/k_{\text{sep}}$ (11)

$1/K_{\text{chan},\text{c}}=d_{\text{cat}}/k_{\text{cat}}+d_{\text{GDL}}/k_{\text{GDL}}+d_{\text{chan}}/k_{\text{chan},\text{c}}+d_{\text{sep}}/k_{\text{sep}}$ (12)

$1/K_{\text{rib},\text{a}}=d_{\text{PEM}}/k_{\text{PEM}}+d_{\text{cat}}/k_{\text{cat}}+d_{\text{GDL}}/k_{\text{GDL}}+d_{\text{rib}}/k_{\text{rib}}+d_{\text{sep}}/k_{\text{sep}}$ (13)

$1/K_{\text{chan},\text{a}}=d_{\text{PEM}}/k_{\text{PEM}}+d_{\text{cat}}/k_{\text{cat}}+d_{\text{GDL}}/k_{\text{GDL}}+d_{\text{chan}}/k_{\text{chan},\text{a}}+d_{\text{sep}}/k_{\text{sep}}$ (14)

Table 1 lists the specification of cell components used in the model. In Table 1, the effective thermal conductivity of porous media k, are the values of the cell components used in the present experiment and in references [26] [27] . Since the effective thermal conductivities given in Table 1 are obtained when the cell component pores are filled with the air at room temperature, the corrected effective thermal conductivities are calculated for the cell components pores filled with H₂ or O₂ at 353 K or 363 K or 373 K, which were the T_ini value assumed in this study. In this calculation, the thermal conductivity of each gas is from The Japan Society of Mechanical Engineers [28] .

Table 1. Specifications of PEFC components.

In order to solve Equations (6)-(9), the temperatures measured using the thermograph were substituted into these equations as T_surf,c and T_surf,a. The operation conditions used for power generation in order to measure temperatures with thermograph are given in Table 2.

Regarding a cathode gas, this study selects O₂. It can be expected that H₂, which is produced from a renewable energy via H₂O electrolyzer, will be used as a fuel for PEFC in order to realize a zero-CO₂-emission society in the near future. When H₂ is produced by H₂O electrolysis, O₂ is also produced as a by-product. This study suggests that not only H₂ but also O₂ produced from H₂O electrolysis are used for PEFC. This study also proposes that the total system consisting of renewable energy, H₂O electrolyzer, and PEFC system operated using H₂ and O₂ produced by H₂O elecrtrolyzer. Therefore, in this study, O₂ is adopted as the cathode gas for the numerical simulation. If O₂ was adopted as a cathode gas, a higher current density on the interface between PEM and the catalyst layer could be expected, especially under the rib, compared to the case using an air [29] .

Analysis using 1D model as well as 3D model is carried out by means of the data obtained under the conditions listed in Table 2. The experimental procedure for measuring temperature during the power generation has been explained in the reference [26] . In the heat transfer analysis, it was assumed that T_surf,c on the rib side was equal to T_surf,c on the channel side as well as T_surf,a because the difference between them could not be recognized by the measured data.

By the comparison of temperature distribution between in-plane and through-plane, the difference between T_react,rib and T_react,chan was found to be small, i.e., less than 1 K [30] [31] [32] . Consequently, it is believed that the heat flow in the through-plane direction dominates the heat transfer in the cell.

Table 2. Operating conditions of power generation for temperature measurement by thermograph.

Considering the above described assumptions and Equations (6)-(14), the reaction surface temperature T_react is expressed as follows:

$\begin{array}{l}{T}_{\text{react}}=T_{\text{react},\text{rib}}=T_{\text{react},\text{chan}}\\ =\left\{\text{2}{H}_{\text{react}}/A+\left(K_{\text{rib},\text{c}}+K_{\text{chan},\text{c}}\right)T_{\text{surf},\text{c}}+\left(K_{\text{rib},\text{a}}+K_{\text{chan},\text{a}}\right)T_{\text{surf},\text{a}}\right\}\\ /\left(K_{\text{rib},\text{c}}+K_{\text{chan},\text{c}}+K_{\text{rib},\text{a}}+K_{\text{chan},\text{a}}\right)\end{array}$ (15)

2.4. 3D Numerical Simulation Model

In this study, the 3D numerical simulation has been conducted using a multi-physics software COMSOL Multiphysics. This software has the simulation code for PEFC composed of the continuity equation, the Brinkman equation for a momentum transfer, the Maxwell-Stefan equation for a diffusion transfer and Butler-Volmer equation for an electrochemical reaction. This simulation code for PEFC has been validated well by many previous studies [18] [33] [34] .

The continuity equation which considers the gas species in porous media, e.g. catalyst layer, MPL and GDL as well as the gas channel is expressed as follows:

$\frac{\partial }{\partial t}\left({\epsilon }_{\text{p}}\rho \right)+\nabla \cdot \left(\rho u\right)=Q_{\text{m}}$ (16)

where ε_p indicates the porosity [−], ρ indicates the density [kg/m³], $u$ indicates the gas velocity vector [m/s], Q_m indicates the mass source term [kg/(m³∙s)] and t indicates the time [s].

Brinkman equation considering the relationship between the pressure and gas flow velocity, which is solved in porous media, e.g. catalyst layer, MPL and GDL as well as in gas channel, is expressed as follows:

$\begin{array}{l}\frac{\rho }{{\epsilon }_{\text{p}}}\left\{\frac{\partial u}{\partial t}+\left(u\cdot \nabla \right)\frac{u}{{\epsilon }_{\text{p}}}\right\}\\ =-\nabla p+\nabla \cdot \left[\frac{1}{{\epsilon }_{\text{p}}}\left\{\mu \left(\nabla u+{\left(\nabla u\right)}^{\text{T}}\right)-\frac{2}{3}\mu \left(\nabla \cdot u\right)I\right\}\right]-\left({\kappa }^{-1}\mu +\frac{Q_{\text{m}}}{{\epsilon }_{\text{p}}^2}\right)u+F\end{array}$ (17)

where p indicates the pressure [Pa], μ indicates the viscosity [Pa·s], $I$ indicates the unit vector [−], κ indicates the permeability [m²] and $F$ indicates the force vector [kg/(m²∙s)], e.g., gravity.

Maxwell-Stefan equation which considers the mass transfer such as the diffusion, ion transfer and convection transfer is expressed as follows:

$N_{\text{i}}=-D_{\text{i}}\nabla C_{\text{i}}-z_{\text{i}}{u}_{\text{m},\text{i}}F{C}_{\text{i}}\nabla {\phi }_{\text{l}}+C_{\text{i}}u=J_{\text{i}}+C_{\text{i}}u$ (18)

$\frac{\partial C_{\text{i}}}{\partial t}+\nabla \cdot N_{\text{i}}=R_{\text{i},\text{tot}}$ (19)

where $N_{\text{i}}$ indicates the vector of the molar flow rate on the interface between PEM and electrode [mol/(m²·s)], D_i indicates the diffusion coefficient [m²/s], C_i indicates the concentration of ion i [mol/m³], z_i indicates the valence of ion [−], u_m,i indicates the mobility of ion i [(s∙mol)/kg], F indicates the Faraday constant [C/mol], φ_l indicates the electrical potential of liquid [35] [V], $J_{\text{i}}$ indicates the molar flow rate of the convection transfer [mol/(m²∙s)], and R_i,tot indicates the reaction rate of species [mol/(m³∙s)].

Butler-Volmer equation calculates the electrochemical reaction as follows:

$i=i_0\left\{\exp \left(\frac{{\alpha }_{\text{a}}F\eta }{RT}\right)-\exp \left(\frac{-{\alpha }_{\text{c}}F\eta }{RT}\right)\right\}$ (20)

$\eta ={\phi }_{\text{s}}-{\phi }_{\text{l}}-E_{\text{eq}}$ (21)

where i indicates the current density [A/m²], i₀ indicates the exchange current density [A/m²], α_a indicates the charge transfer coefficient at the anode [−], η indicates the activation over-potential [35] [V], R indicates the gas constant [J/(mol·K)], T indicates the temperature [K], α_c indicates the charge transfer coefficient at the cathode [−], φ_s indicates the electrical potential of solid [35] [V], E_eq indicates the equilibrium electric potential [35] [V].

Heat transfer equation considering electrical reaction is expressed as follows:

$\rho C_{\text{p}}u\cdot \nabla T=\nabla \cdot \left(k\nabla T\right)+Q_{\text{jh}}+{\displaystyle \underset{m}{\sum }{a}_{\text{v}}}{Q}_{\text{e}}$ (22)

$Q_{\text{jh}}=-\left(i_{\text{s}}\cdot \nabla {\phi }_{\text{s}}+i_{\text{l}}\cdot \nabla {\phi }_{\text{l}}\right)$ (23)

$Q_{\text{e}}=\left(\eta +T\frac{\delta E_{\text{eq}}}{\delta T}\right)i$ (24)

where C_p indicates the specific heat [J/(kg∙K)], u is gas velocity [m/s], k indicates the thermal conductivity [W/(m∙K)], a_v indicates the activation specific area [1/m], $i_{\text{s}}$ indicates the current density vector in electrode [A/m²] and $i_{\text{l}}$ indicates the current density in electrolyte [A/m²].

Figure 2 illustrates 3D model of single cell of PEFC for the numerical simulation used in this study [20] [34] . This structure follows the commercial single cell used in the experimental studies carried out by the authors [36] [37] . The roof of the gas separator at anode side and cathode side is omitted in this model. The cell has a gas separator with a serpentine flow channel consisting of five gas channels with the width of 1.0 mm and the depth of 1.0 mm as well as a rib with the width of 1.0 mm. The size of cell components listed in Table 1 is adopted for this numerical simulation. Operation conditions listed in Table 2 are also adopted for this numerical simulation. Table 3 lists physical parameters adopted for this numerical simulation. To investigate and compare the distribution of T_react between 1D heat transfer model and 3D numerical simulation, this study selects the analysis points of A to K as shown in Figure 3. The average value on the cross sectional area of the interface between PEM and cathode catalyst layer at each point, which covers the part under gas channel and that under rib, is calculated.

This study set the following assumptions [20] [34] .

1) The distributions of the inlet gas flow rate at the anode side and the cathode side are uniform, respectively.

2) The pressure of the outlet of the gas channel is the atmospheric pressure.

3) No slip on the gas channel wall excluding the inlet and the outlet of the gas channel is considered.

Figure 2. 3D model for numerical simulation in this study.

Figure 3. Analysis points for the quantitative evaluation along with the gas flow through gas channel.

Table 3. Physical parameters [19] [20] [21] [24] [26] [28] [29] [34] [38] - [44] .

4) The cell voltage obtained by the power generation experiment is set at the cathode electrode and the earth ground in set at the anode electrode. The in-plane distribution of cell voltage at the cathode electrode is uniform.

5) Reactant gases are treated as an ideal gas and incompressible Newton fluid.

6) H₂O is treated as a vapour.

7) The cell temperature is uniform and the outside boundary of the 3D model is set at T_ini.

8) The effective porosity and the permeability of the porous media are isotropic. The conductivity in the porous media is also isotropic.

The impacts of T_ini and RH of supply gas on the distribution of T_react have been investigated by 1D model and 3D model. The impacts of T_ini and RH of supply gas on the distributions of H₂, O₂, H₂O and current density have been also investigated by 3D model. In this paper, we focus on the distributions of O₂, H₂O and current density as well as distribution of T_react exhibited on the reaction surface, which are shown later.

3. Results and Discussion

3.1. Comparison of Distribution of T_react between 1D Heat Transfer Model and 3D Numerical Simulation

Figure 4 and Figure 5 show distributions of T_react calculated by 1D heat transfer

Figure 4. Distributions of T_react calculated by 1D heat transfer model at T_ini = 353 K.

Figure 5. Distributions of T_react calculated by 3D heat transfer model at T_ini = 353 K.

model and 3D numerical simulation at T_ini = 353 K, respectively. Figure 6 and Figure 7 show distributions of T_react calculated by 1D heat transfer model and 3D numerical simulation at T_ini = 363 K, respectively. Figure 8 and Figure 9 show distributions of T_react calculated by 1D heat transfer model and 3D numerical simulation at T_ini = 373 K, respectively. In these figures, RH of supply gas is changed.

Figure 6. Distributions of T_react calculated by 1D heat transfer model at T_ini = 363 K.

Figure 7. Distributions of T_react calculated by 3D heat transfer model at T_ini = 363 K.

Figure 8. Distributions of T_react calculated by 1D heat transfer model at T_ini = 373 K.

Figure 9. Distributions of T_react calculated by 3D heat transfer model at T_ini = 373 K.

According to Figures 4-9, it is seen that the change in T_react from the inlet to the outlet is more even with the increase in T_ini irrespective of the investigated model. Since the saturation pressure of H₂O vapour increase with the temperature exponentially [44] , it is easy to dehydrate PEM at higher temperature. In other words, it is easy to decrease the proton conductivity of PEM at higher temperature. As a result, the power generation performance is dropped at higher temperature due to large ohmic loss, resulting in the lower generated heat. Therefore, the change in T_react from the inlet to the outlet is more even with the increase in T_ini. Since this study has set the excess gas which is over s.r. = 1.0 as the inlet gas flow rate, the generated heat is accumulated along with the excess gas flow through the gas channel. Consequently, it is thought T_react increases from the inlet to the outlet largely, especially at T_ini = 353 K since the power generation performance is better.

Comparing the results obtained by 1D heat transfer model with those obtained by 3D numerical simulation, the temperature gas between them is below approximately 0.5 K. Therefore, it can be claimed that 1D heat transfer model predicts the distribution of T_react well even though we think the heat transfer in single cell of PEFC only. This study has calculated the amount of heat taken by the gas flow along through the gas channel from the inlet to the outlet of the cell from the results obtained by 3D numerical simulation, resulting that it is approximately 0.01% of the heat generated. Therefore, it is thought that 1D heat transfer model can predict the distribution of T_react well. However, the conditions validated by this study are 353 K, 363 K and 373 K only. In the near future, this study will validate under the other operation condition to verify the accuracy of 1D Heat Transfer Model proposed by the authors. In the following section, we discuss the other distributions to clarify the phenomena.

3.2. Distributions of Molar Concentration of O₂ Calculated by 3D Numerical Simulation

Figure 10 shows distributions of molar concentration of O₂ calculated by 3D

Figure 10. Distributions of molar concentration of O₂ calculated by 3D numerical simulation. (a) 353 K; (b) 363 K; (c) 373 K.

numerical simulation at T_ini = 353 K, 363 K and 373 K, respectively. In this figure, RH of supply gas is changed.

According to Figure 10, it is seen that the change in the molar concentration of O₂ from the inlet to the outlet is more even with the increase in T_ini. It is known from the previous studies [45] [46] that the proton conductivity of PEM increases with the increase in temperature as well as the increase in RH. On the other hand, the saturation pressure of H₂O vapour increases with the temperature exponentially [44] , resulting that it is easy to dehydrate PEM at T_ini = 373 K compared with T_ini = 353 K. As a result, the proton conductivity of PEM decreases at T_ini = 373 K. If the proton conductivity of PEM decreases, the performance of the O₂ reduction reaction drops by the lack of proton. Since the hydration of PEM is not enough at T_ini = 373 K, the high O₂ partial pressure is needed to progress the O₂ reduction reaction [44] . Therefore, the change in the molar concentration of O₂ from the inlet to the outlet decreases with the increase in T_ini.

3.3. Distributions of Molar Concentration of H₂O Calculated by 3D Numerical Simulation

Figure 11 shows distributions of molar concentration of H₂O calculated by 3D numerical simulation at T_ini = 353 K, 363 K and 373 K, respectively. In this figure, RH of supply gas is changed.

According to Figure 11, it is seen that the change in the molar concentration of H₂O from the inlet to the outlet is more even with the increase in T_ini. As discussed above, the proton conductivity of PEM increases with the increase in temperature as well as the increase in RH [45] [46] . On the other hand, since the saturation pressure of H₂O vapour increases with the temperature exponentially [44] , it is easy to dehydrate PEM at T_ini = 373 K compared with T_ini = 353 K. Therefore, the proton conductivity of PEM decreases at T_ini = 373 K. If the proton conductivity of PEM decreases, the performance of O₂ reduction reaction drops by the lack of proton. Moreover, the dehydration of PEM is not enough at T_ini = 373 K, resulting that high O₂ partial pressure is needed to progress the O₂ reduction reaction [44] . Consequently, it can be claimed that the change in the molar concentration of H₂O from the inlet to the outlet becomes more even with the increase in T_ini. According to Figures 4-9, it is observed that the change in T_react from the inlet to the outlet is more even with the increase in T_ini. Since the generated heat from the power generation is lower with the increase in T_ini due to the lower performance of O₂ reduction reaction [47] , it can be thought that the change in T_react from the inlet to the outlet is more even with the increase in T_ini.

It can be seen from Figure 11 that the change in the molar concentration of H₂O from the inlet to the outlet for A40%RH, C40%RH is more even compared with the other RH conditions. Since A40%RH, C40%RH is the dry condition, PEM and catalyst layer are dehydrated easily [48] . The proton conductivity of PEM is smaller under a dry condition [49] . In addition, the RH influences the performance of O₂ reduction reaction carrying out on the ionomer in the catalyst layer at the cathode [40] . Therefore, there is the optimum H₂O saturation for ionomer in the catalyst layer at cathode [21] , indicating that the performance of O₂ reduction reaction which produces H₂O is lower for A40% RH, C40% RH. On the other hand, it is seen from Figures 4-9 that the distribution of T_react for A40%RH, C40%RH is relatively even compared to the other RH conditions.

Figure 11. Distributions of molar concentration of H₂O calculated by 3D numerical simulation. (a) 353 K; (b) 363 K; (c) 373 K.

Since the generated heat from the power generation is lower for A40%RH, C40%RH due to the lower performance of O₂ reduction reaction [47] , it can be thought that the distribution of T_react for A40%RH, C40%RH is more even.

3.4. Distributions of Current Density Calculated by 3D Numerical Simulation

Figure 12 shows distributions of current density calculated by 3D numerical

Figure 12. Distributions of current density calculated by 3D numerical simulation. (a) 353 K; (b) 363 K; (c) 373 K.

simulation at T_ini = 353 K, 363 K and 373 K, respectively. In this figure, RH of supply gas is changed.

According to Figure 12, it is seen that the change in the current density from the inlet to the outlet is more even with the increase in T_ini and the value of current density is smaller with the increase in T_ini. As discussed above, the proton conductivity of PEM decreases with the increase in T_ini, resulting that the performance of O₂ reduction reaction drops due to the lack of proton. In addition, the hydration of PEM is not enough at high temperature, high O₂ partial pressure is needed to progress the O₂ reduction reaction [44] . Consequently, it can be claimed that the current density decreases with the increase in T_ini due to the increase in the ohmic over-potential and the concentration over-potential [44] .

According to Figures 4-9, it is observed that the change in T_react from the inlet to the outlet is more even with the increase in T_ini. Since the generated heat from the power generation is lower with the increase in T_ini due to the lower performance of O₂ reduction reaction [47] , it can be thought that the change in T_react from the inlet to the outlet is more even with the increase in T_ini.

From the investigation of this study, we can claim that it is necessary to control the hydration of PEM and catalyst layer in order to obtain the high power generation performance at higher temperature such as 363 K and 373 K. As a procedure to control the hydration of PEM and catalyst layer under higher temperature operation condition, this study suggests recirculating the H₂O which is emitted from the cell and promoting the heat transfer in order to cool the cell. This study would like to investigate these trials in the next step.

4. Conclusions

The temperature distribution on the reaction surface simulated/predicted by the 1D heat transfer model has been validated by the 3D numerical simulation using COMSOL Multiphysics. The effects of T_ini, especially higher temperature than usual operation condition and RH of supply gas on the distribution of T_react, have been investigated. In addition, the impacts of T_ini and RH of supply gas on the distributions of O₂, H₂O and current density have been also investigated by the 3D model. The following conclusions can be drawn from the study:

1) It is revealed that the change in T_react from the inlet to the outlet is more even with the increase in T_ini irrespective of the investigated model. This is because the generated heat from the power generation is lower with the increase in T_ini due to the lower performance of O₂ reduction reaction.

2) It is confirmed that the temperature gap between the results obtained by 1D heat transfer model and those obtained by 3D numerical simulation is below approximately 0.5 K. It can be claimed that 1D heat transfer model predicts the distribution of T_react well.

3) According to the 3D numerical simulation, the change in the molar concentration of O₂ and H₂O from the inlet to the outlet is more even with the increase in T_ini due to the lower performance of O₂ reduction reaction.

4) According to the 3D numerical simulation, the change in the current density from the inlet to the outlet is more even with the increase in T_ini and the value of current density is smaller with the increase in T_ini since the ohmic over-potential and the concentration over-potential increase.

It is necessary to control the hydration of PEM and catalyst layer in order to obtain high power generation performance at higher temperatures such as 363 K and 373 K. As an example of the procedure, this study suggests recirculating the H₂O which is emitted from the cell, and promoting the heat transfer in order to cool the cell.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	NEDO (New Energy and Industry Technology Development Organization) (2022). http://www.nedo.go.jp/content/100871976.pdf
[2]	Zhang, J., Wang, H., Li, W., Zhang, J., Lu, D., Yan, W., Xiang, Y. and Lu, S. (2021) Effect of Catalyst Layer Microstructures on Performance and Stability for High Temperature Polymer Electrolyte Membrane Fuel Cells. Journal of Power Sources, 505, Article ID: 230059. https://doi.org/10.1016/j.jpowsour.2021.230059
[3]	Zhang, G. and Kandlikar, S.G.A. (2012) Critical Review of Cooling Technique in Proton Exchange Membrane Fuel Cell Stacks. International Journal of Hydrogen Energy, 37, 2412-2429. https://doi.org/10.1016/j.ijhydene.2011.11.010
[4]	Agbossou, K., Kolhe, M., Hamelin, J. and Bose, T.K. (2004) Performance of a Stand-Alone Renewable Energy System Based on Energy Storage as Hydrogen. IEEE Transactions on Energy Conversion, 19, 633-640. https://doi.org/10.1109/TEC.2004.827719
[5]	Li, Q., He, R., Jensen, J.O. and Bjerrum, N.J. (2003) Approaches and Recent Development Polymer Electrolyte Membrane for Fuel Cells Operating above 100 °C. Chemical of Materials, 15, 4896-4915. https://doi.org/10.1021/cm0310519
[6]	Lee, C.Y., Weng, F.B., Kuo, Y.W., Cheng, C.H., Cheng, C.K. and Lin, J.T. (2016) In-Situ Measurement of High-Temperature Resistant Integrated Microsensor Embedded in High Temperature Proton Exchange Membrane Fuel Cell Stack. Sensors, 16, Article No. 1731. https://doi.org/10.3390/s16101731
[7]	Lee, C.Y., Weng, F.B., Kuo, Y.W., Cheng, Y.T., Cheng, C.K., Tsai, C.H. and Lee, T.J. (2016) Persistent Effect Test for High Temperature Resistant Integrated Microsensor Embedded in High Temperature Proton Exchange Membrane Fuel Cell Stack. Sensors and Actuators A: Physical, 250, 202-209. https://doi.org/10.1016/j.sna.2016.09.026
[8]	Wang, M., Guo, H. and Ma, C. (2006) Temperature Distribution on the MEA Surface of a PEMFC with Serpentine Channel Flow Bed. Journal of Power Sources, 157, 181-187. https://doi.org/10.1016/j.jpowsour.2005.08.012
[9]	Tsuji, K. (2008) Domestic Fuel Cell Co-Generation System Entering Real Commercial Stage. Hydrogen Energy System, 33, 93-96.
[10]	Ryu, S.K., Vinothkannan, M., Kim, A.R. and Yoo, D.J. (2022) Effect of Type and Stoichiometry of Fuels on Performance of Polybenzimidazole-Based Proton Exchange Membrane Fuel Cells Operating at the Temperature Range of 120-160 °C. Energy, 238, Article ID: 121791. https://doi.org/10.1016/j.energy.2021.121791
[11]	Budak, Y. and Devrim, Y. (2022) Micro-Cogeneration Application of a High-Temperature PEM Fuel Cell Stack Operated with Polybenzimidazole Based Membranes. International Journal of Hydrogen Energy, 45, 35198-35207. https://doi.org/10.1016/j.ijhydene.2019.11.173
[12]	Kim, D.K., Kim, H., Park, H., Oh, S., Ahn, S.H., Kim, H.J. and Kim, S.K. (2019) Performance Enhancement of High-Temperature Polymer Electrolyte Membrane Fuel Cells Using Pt Pulse Electrodeposition. Journal of Power Sources, 438, Article ID: 227022. https://doi.org/10.1016/j.jpowsour.2019.227022
[13]	Kanchan, B.K., Randive, P. and Pati, S. (2021) Implications of Non-Uniform Porosity Distribution in Gas Diffusion Layer on the Performance of a High Temperature PEM Fuel Cell. International Journal of Hydrogen Energy, 46, 18571-18588. https://doi.org/10.1016/j.ijhydene.2021.03.010
[14]	Xia, L., Ni, M., He, Q., Xu, Q. and Cheng, C. (2021) Optimization of Gas Diffusion Layer in High Temperature PEMFC with the Focuses on Thickness and Porosity. Applied Energy, 300, Article ID: 117357. https://doi.org/10.1016/j.apenergy.2021.117357
[15]	Agarwal, H., Thosar, A.U., Bhat, S.D. and Lele, A.K. (2022) Interdigitated Flow Field Impact on Mass Transport and Electrochemical Reaction in High-Temperature Polymer Electrolyte Fuel Cell. Journal of Power Sources, 532, Article ID: 231319. https://doi.org/10.1016/j.jpowsour.2022.231319
[16]	Xia, L., Xu, Q., He, Q., Ni, M. and Seng, M. (2021) Numerical Study of High Temperature Proton Exchange Membrane Fuel Cell (HT-PEFC) with a Focus on Rib Design. International Journal of Hydrogen Energy, 46, 21098-21111. https://doi.org/10.1016/j.ijhydene.2021.03.192
[17]	Huang, T., Wang, W., Yuan, Y., Huang, J., Chen, X., Zhang, J., Kong, X., Zhang, Y. and Wan, Z. (2021) Optimization of High-Temperature Proton Exchange Membrane Fuel Cell Flow Channel Based on Genetic Algorithm. Energy Reports, 7, 1374-1384. https://doi.org/10.1016/j.egyr.2021.02.062
[18]	Zhang, J., Zhang, C., Hao, D., Ni, M., Hung, S., Liu, D. and Zheng, Y. (2021) 3D Non-Isothermal Dynamic Simulation of High Temperature Proton Exchange Membrane Fuel Cell in Start-Up Process. International Journal of Hydrogen Energy, 46, 2577-2593. https://doi.org/10.1016/j.ijhydene.2020.10.116
[19]	Nishimura, A., Kono, N., Toyoda, K., Mishima, D. and Kolhe, M.L. (2022) Impact of Separator Thickness on Temperature Distribution in Single Cell of Polymer Electrolyte Fuel Cell Operated at Higher Temperature of 90°C and 100°C. Energies, 15, Article No. 4203. https://doi.org/10.3390/en15124203
[20]	Nishimura, A., Toyoda, K., Mishima, D., Ito, S. and Hu, E. (2022) Numerical Analysis on Impact of Thickness of PEM and GDL with and without MPL on Coupling Phenomena in PEFC Operated at Higher Temperature Such as 363 K and 373 K. Energies, 15, Article No. 5936. https://doi.org/10.3390/en15165936
[21]	Nishimura, A., Kono, N., Toyoda, K., Kojima, Y. and Kolhe, M.L. (2021) Impact Analysis of MPL on a PEFC Cell’s Temperature Distribution with Thin PEM and GDL for Operating at Higher Temperature than Usual. Journal of Energy and Power Engineering, 15, 39-51. https://doi.org/10.17265/1934-8975/2021.02.001
[22]	Nishimura, A., Yamamoto, K., Okado, Y., Kojima, Y., Hirota, M. and Kolhe, M.L. (2020) Impact of Analysis of MPL and PEM Thickness on Temperature Distribution with PEFC Operating at Relatively Higher Temperature. Energy, 205, Article ID: 117875. https://doi.org/10.1016/j.energy.2020.117875
[23]	Nishimura, A., Sato, Y., Kamiya, S., Okado, T., Yamamoto, K., Hirota, M. and Hu, E. (2019) Impact of Thickness of Polymer Electrolyte Membrane and Gas Diffusion Layer on Temperature Distribution in Polymer Electrolyte Fuel Cell Operated at Temperature around 90 °C. Journal of Energy and Power Engineering, 13, 97-115. https://doi.org/10.17265/1934-8975/2019.03.002
[24]	Nishimura, A., Sato, Y., Yoshimura, M., Kamiya, S. and Hirota, M. (2018) Impact of Thickness of Polymer Electrolyte Membrane on Temperature Distribution in Single Cell of Polymer Electrolyte Fuel Cell Operated at High Temperature. Journal of Energy and Power Engineering, 12, 80-92. https://doi.org/10.17265/1934-8975/2018.02.004
[25]	Nishimura, A., Iio, K., Baba, M., Yamauchi, T., Hirota, M. and Hu, E. (2014) Modeling of Heat Transfer in Single Cell of Polymer Electrolyte Fuel Cell by Means of Temperature Data Measured by Thermograph. Journal of Chemical Engineering of Japan, 47, 521-529. https://doi.org/10.1252/jcej.13we275
[26]	Nishimura, A., Shibuya, K., Morimoto, A., Tanaka, S., Hirota, M., Nakamura, Y., Kojima, M., Narita, M. and Hu, E. (2012) Dominant Factor and Mechanism of Coupling Phenomena in Single Cell of Polymer Electrolyte Fuel Cell. Applied Energy, 90, 73-79. https://doi.org/10.1016/j.apenergy.2011.01.003
[27]	Khandelwah, M. and Mench, M.M. (2006) Direct Measurement of Through-Plane Thermal Conductivity and Contact Resistance in Fuel Cell Materials. Journal of Power Sources, 161, 1106-1115. https://doi.org/10.1016/j.jpowsour.2006.06.092
[28]	The Japan Society of Mechanical Engineers (1993) JSME Heat Transfer Handbook. Maruzen, Tokyo, 387.
[29]	Penga, Z., Tolj, I. and Barbir, F. (2016) Computational Fluid Dynamics Study of PEM Fuel Cell Performance. International Journal of Hydrogen Energy, 41, 17585-17594. https://doi.org/10.1016/j.ijhydene.2016.07.092
[30]	Kawase, M., Inagaki, T., Kawashima, S. and Miura, K. (2009) Effective Thermal Conductivity of Gas Diffusion Layer in Through-Plane Direction. ECS Transactions, 25, 1529-1537. https://doi.org/10.1149/1.3210709
[31]	Oshima, A., Nishimura, A., Morimoto, A., Tanaka, S., Hirota, M. and Narita, M. (2010) Theoretical Investigation on Influence of Inflow Gas Condition and Gas Channel Structure of Separator on Mass and Temperature Distribution in Single Cell of Polymer Electrolyte Fuel Cell. Preprints of Mechanical Engineering Congress, 203-204.
[32]	Jung, C.Y., Shim, H.S., Koo, S.M., Lee, S.F. and Yi, S.C. (2012) Investigations of the Temperature Distribution in Proton Exchange Membrane Fuel Cell. Applied Energy, 93, 733-741. https://doi.org/10.1016/j.apenergy.2011.08.035
[33]	Das, S.K. and Gibson, H.A. (2021) These Dimensional Multi-Physics Modeling and Simulation for Assessment of Mass Transport Impact on the Performance of a High Temperature Polymer Electrolyte Membrane Fuel Cell. Journal of Power Sources, 499, 161-188. https://doi.org/10.1016/j.jpowsour.2021.229844
[34]	Nishimura, A., Toyoda, K., Kojima, Y. and Kolhe, M.L. (2021) Numerical Simulation on Impacts of Thickness of Nafion Series Membranes and Relative Humidity on PEMFC Operated at 363 K and 373 K. Energies, 14, Article No. 8256. https://doi.org/10.3390/en14248256
[35]	Chen, H., Guo, H., Ye, F. and Ma, C.F.A. (2021) Numerical Study on Oriented-Type Flow Channels with Porous-Blocked Baffles of Proton Exchange Membrane Fuel Cells. International Journal of Hydrogen Energy, 46, 29443-29458. https://doi.org/10.1016/j.ijhydene.2020.12.178
[36]	Nishimura, A., Okado, T., Kojima, Y. and Hu, E. (2021) Impact of MPL on Temperature Distribution in Single Polymer Electrolyte Fuel Cell with Various Thickness of Polymer Electrolyte Membrane. Energies, 13, Article No. 2499. https://doi.org/10.3390/en13102499
[37]	Nishimura, A., Kamiya, S., Okado, T., Sato, Y., Hirota, M. and Kolhe, M.L. (2019) Heat and Mass Transfer Analysis in Single Cell of PEFC Using Different PEM and GDL at Higher Temperature. International Journal of Hydrogen Energy, 44, 29631-29640. https://doi.org/10.1016/j.ijhydene.2019.05.192
[38]	Copper, N.J., Santamaria, A.D., Becton, M.K. and Park, J.W. (2017) Neutron Radiography Measurement of In-situ PEMFC Liquid Water Saturation in 2D & 3D Morphology Gas Diffusion Layers. International Journal of Hydrogen Energy, 42, 16269-16678. https://doi.org/10.1016/j.ijhydene.2017.05.105
[39]	Merck (2023). https://www.sigmaaldrich.com/SG/en/applications/materials-science-and-engineering/batteries-supercapacitors-and-fuel-cells
[40]	Rostami, L., Nejad, P.M.G. and Vatani, A.A. (2016) Numerical Investigation of Serpentine Flow Channel with Different Bend Sizes in Polymer Electrolyte Membrane Fuel Cells. Energy, 97, 400-410. https://doi.org/10.1016/j.energy.2015.10.132
[41]	Senn, S.M. and Poulikakos, D. (2004) Polymer Electrolyte Fuel Cells with Porous Materials as Fluid Distributions and Comparisons with Traditional Channelled Systems. Transactions of ASME, 126, 410-418. https://doi.org/10.1115/1.1738424
[42]	Takayama, T. (2018) Numerical Simulation of Transient International States of PEFC Cell and Stack Considering Control of Anode System. Research Report of Mizuho Research Technology, 9, 1-14.
[43]	TORAY (2023). http://www.torayca.com/en/lineup/composites/com_009_01.html
[44]	Xing, L., Das, P.K., Song, X., Mamlouk, M. and Scott, K. (2015) Numerical Analysis of the Optimum Membrane/Ionomer Water Content of PEMFCs: The Interface of Nafion Ionomer Content and Cathode Relative Humidity. Applied Energy, 138, 242-257. https://doi.org/10.1016/j.apenergy.2014.10.011
[45]	Akimoto, K., Sasabe, T., Yoshida, T., Naito, H., Kawamura, K. and Hirai, S. (2019) Investigation of Effects of High Temperature and Pressure on a Polymer Electrolyte Fuel Cell with Polarization Analysis and X-Ray Imaging of Liquid Water. Journal of Power Sources, 431, 205-209. https://doi.org/10.1016/j.jpowsour.2019.04.115
[46]	Jia, T., Shen, S., Zhao, J., Jin, J., Pan, B., Duan, X., Meng, C. and Che, Q. (2020) Ultrathin Membranes Formation via the Layer by Layer Self-Assembly of Carbon Nanotubes-Based Inorganics as High Temperature Proton Exchange Membranes. International Journal of Hydrogen Energy, 45, 14517-14527. https://doi.org/10.1016/j.ijhydene.2020.03.175
[47]	Miao, T., Tongsh, C., Wang, J., Cheng, P., Liang, J., Wang, Z., Chen, W., Zhang, C., Xi, F., Du, Q., Wang, B., Bai, F. and Jiao, K. (2022) Current Density and Temperature Distribution Measurement and Homogeneity Analysis for a Large-Area Proton Exchange Membrane Fuel Cell. Energy, 239, Article ID: 121922. https://doi.org/10.1016/j.energy.2021.121922
[48]	Springer, T.E., Zawodzinski, T.A. and Gottesfeld, D. (1991) Polymer Electrolyte Fuel Cell Model. Journal of the Electrochemical Society, 138, 2334-2341. https://doi.org/10.1149/1.2085971
[49]	Zhou, J., Shukla, S., Putz, A. and Secanell, M. (2018) Analysis of the Role of the Microporous Layer in Improving Polymer Electrolyte Fuel Cell Performance. Electrochimica Acta, 268, 366-382. https://doi.org/10.1016/j.electacta.2018.02.100