Nonisothermal two‐phase modeling of the effect of linear nonuniform catalyst layer on polymer electrolyte membrane fuel cell performance

In this research, it is investigated to numerically evaluate the performance of a polymer electrolyte membrane fuel cell (PEMFC). The performance is investigated through the nonuniformity gradient loading at the catalyst layer (CL) of the considered PEMFC. Computational fluid dynamics is used to simulate a 2D domain in which a steady‐state laminar compressible flow in two‐phase for the PEMFC has been considered. In this case, a particular nonuniform variation inside the CL along the channel is assumed. The nonuniform gradient is created using a nonisothermal domain to predict the flooding effects on the performance of the PEMFC. The computational domain is considered as the cathode of PEMFC, which is divided into three regions: a gas channel, a gas diffusion layer, and a CL. The loading variation inside the catalyst is defined as a constant slope along the channel. In order to find the optimum slope, different slope angles are analyzed. The results point out that the nonuniform loading distribution of the catalyst (platinum) along the channel could improve the fuel cell performance up to 1.6% and 5% for power density and voltage generation, respectively. It is inferred that it is better to use more catalyst in the final section of the channel if the performance is the main concern.


| INTRODUCTION
Nowadays, heat and power production from various primary energy authorities demands clean and efficient sources such as fuel cells. 1 For fuel cells, recent years numerous improvements have been made over the cell geometry, design, and material selection. Many changes have been occurred, and several technical barriers and limitations are resolved. However, some obstacles still exist in the development of polymer electrolyte membrane fuel cells (PEMFCs) to be commercialized. As researchers advance in technology of fuel cells, some obstacles get really crucial like low cell performance and high cost due to expensive materials. Hence, considerable research topics mostly experimental have been introduced to improve the performance and eliminate huge costs by changing the fuel cell designs. 2 Performance enhancement for a fuel cell is an achievable goal, only by understanding the related chemical reactions in detail. On the other hand, achieving high-efficiency fuel cell is a function of ultimate catalyst utilization and optimal catalyst distribution. Recently, many models to elaborate the efficiency and performance have been proposed. 3,4 In earlier models, 5,6 the majority of experimental data was available mainly about the diffusion data and the drag coefficient of the electro-osmotic of Nafion to attain the water content. Such investigations, for instance, proved the ionic conductivity of the membrane. Moreover, in previous researches, multicomponent diffusion equation (Stephan-Maxwell) was mainly used and the temperature of fuel cell was considered to be uniform and involvement of different phases has been ignored.
In another study, 3 the transport model for the gas of Bernardi 7 was coupled with the membrane model of Verbrugge and Hill. 8 Additionally, the liquid phase was considered, while the electrochemical reactions were modeled by the Butler-Volmer equation. Literature survey shows that Gurau et al 9 were prominent researchers to apply the 2D model of polymer electrolyte membrane (PEM) fuel cell. Later, Um et al 10 by using computational fluid dynamics (CFD) simulations were prosperous to implement a 2D simulation of a fuel cell. In their model, in order to solve the equations simultaneously, the gas-phase transport and electrochemical reactions were coupled. After that, the implementation of CFD tools in simulation of fuel cells was introduced as a comprehensive topic to solve the chemical and fluid dynamics interactions.
Although extensive focus has been made on the modeling of fuel cells, the two-phase flow modeling has not been carefully studied yet. However, there are several 2D studies in this field, among which some of the recent findings are listed in the following. Wang et al 11 investigated the cathode of the PEMFC in a 2D two-phase model and carried out an extensive analytical research on the dependency of the twophase model to the inlet humidity. In their model, there is the lack of modeling of the catalyst layer (CL), and as an input, the current density is given to the model. Moreover, they applied a multiphase mixture model to consider the maximum number of phenomena in the fuel cell. In another work, You et al 12 used the M2 together with one-dimensional models for the two-phase flow and CL to calculate the current density. Besides, Natargan and Nguyen 13 investigated a 2D two-phase model of a polymer-based fuel cell. They used multicomponent diffusion model while the porosity was modeled with the Darcy equation. The generated water in their model was in liquid phase. In addition to the investigation of thermophysical properties of fuel cells, Dupont and Legendre 14 and Zhang et al 15 tried to investigate the flow pattern in PEMFC channels. Wang and Zhang 16 also investigated a 3D approach, where the transport model is simulated using a two-phase model for the PEMFCs. In such model, the liquid water transport in the flow passages is modeled by implementing the generalized Richards equation. In addition, in their model the effects of geometry of the channel were explored on the performance of the fuel cells. Then, they concluded that the aspect ratio of the channel could affect the efficiency of the fuel cell dramatically. Quan and Lai 17 in a numerical approach investigated the flow behavior inside an interdigitated PEMFC. The wettability of surface of the channel and the working pressure effects on the thermochemical behavior of water were examined. You and Liu 18 established a multicomponent mathematical model in which the polarization curves of fuel cell for the cathode pressure changes have been predicted. They also validated their results with a reasonable agreement with the available experimental data. In another research, Hua and Meng 19 did a parametric study on the liquid water transport. It means that, the effect of exterior temperature of the cell boundaries on the liquid behavior of water and the fuel cell performance was observed.
The idea of application of nonuniform loading of the CL was proposed for the first time by Kulikovsky, 20 a generic CL presented in his work with adjustable catalyst loading inside the layer interfaces. By the aforementioned assumptions, he also achieved the optimal loading distribution at the catalyst level. Recently, catalyst structural properties and their effects on the performance of CL have been studied by several research groups. Kamarajugadda et al 21 studied the effect of composition and structure of the cathode layers on the overall performance of the CL by applying the concept of agglomeration at the layers that the results implied the importance of properties of agglomeration. In another approach, the effects of properties of the agglomeration and the CL thickness on the overall performance of fuel cells have been studied by Khajeh-Hosseini et al. 22 Their simulation results proved that the performance of the cathode of the CL could be impressed by the agglomeration properties (polymer volume fraction, thickness of CL, and carbon loading). As another example, Srinivasarao et al 23  the CL as a network of spherical agglomerates instead of a single bulk layer. In their work, the catalyst loading gradient across the layer was in focus. Therefore, the platinum loading was changed from the gas diffusion layer (GDL) to other interfaces, while its total (accumulative) amount was similar to that of an individual layer of CL. They confirmed that the network of multilayer is superior to a single-layer CL. Roshandel and Ahmadi 24 evaluate a particular variation of catalyst loading along two directions. This variable loading was defined from the membrane/CL interface to the GDL and "in the catalyst plane" under the channels and land areas. In another research, Ebrahimi et al 25 employed a CFD approach to perform investigation on the agglomerate model of the cathode catalyst layer (CCL) to enhance the power density of the PEMFC. Such model was implemented by the optimal distribution of the catalyst in the CCL. In their CFD work, a network of flooded spherical cluster model was created. This development of network was done for the CCL to simulate the transport of both chemical species and charges. It was observed that in the optimum case, the density of the PEMFC reaches to a maximum value, which has about 7% increment in comparison with the base case. Ken-Ming Yin et al 26 mainly concentrated on the effect of the CCL physical properties and its structure to enhance the performance of the PEMFC. They studied the assembly strategy and the distribution of material properties in the CL. Three sublayers were defined for the CL with variable material properties. By considering the water content effects, high porous regions near the GDL show superior performance and a higher fraction of electrolyte in the clustered/agglomerated area near the membrane. Havaej et al 27 investigated the loading variations at the catalyst level by changing it along both vertical and horizontal directions of a PEMFC. Accordingly, by proper loading distributions of catalyst only along one direction in the cell, 3.1% improvement was achieved, while 8% enhancement was obtained by distribution in both directions. It was inferred that the distribution of catalyst is more operative at the high current densities rather than lower ones.
Molaeimanesh et al 28 investigated morphology impacts on the CL of PEM fuel cell. In their approach, the cathode modeling was done by a 3D lattice Boltzmann agglomerate network. Five different CLs of cathode with different morphologies were investigated. When a mole fraction comparison was done, their results showed that oxygen and water vapor are distributed in an uneven manner. In addition, by growing of the ellipsoid stretching, this diverse distribution becomes more severe. Particularly about implementing different designs of CL, Therdthianwong et al 29 investigated the effect of polytetrafluoroethylene (PTFE) addition in the catalyst ink as well as the CL loading pattern. It was found that the presence of 10% PTFE in the CL leads to 34% increase in maximum power density. Also, the optimum Pt loading for PTFE-added CL was found to be 0.25 mg/cm 2 . In another research, Yin et al 30 used catalyst sprayed membrane method to prepare a novel Pt-C/Fe-N-S-C cathode dual catalyst layer (CDCL). The Fe-N-S-C catalyst led to a remarkable water vapor sorption capacity. The optimum Fe-N-S-C content was 1.5-2 mg/cm 2 , and the best operating temperature was 70°C. Moreover, Zheng et al 31 mathematically employed gradient CCL to improve the performance and durability of the call. They investigated several gradient CCL structures with Pt/C catalysts and various electrochemical surface area and Pt mass. It was illustrated that the CCL with an optimum performance is gained by choosing the best values of particle size gradient and Pt loading gradient simultaneously, which leads to a more uniform distribution of Pt active surface for improving the end-of-test performance and durability.
Specifically about nonuniformly distributed components of fuel cells, Prasanna et al 32 evaluated the methods of reduction in the costs of PEMFC by the decrease in the loading of the platinum catalyst. In this analysis, they used the catalyst-gradient electrode method and the distributions of water and reactant gases were considered nonuniform. Their findings illustrated that the losses of the PEMFC are reduced effectively with this technique. In addition, Xing et al 33 studied the performance of the inhomogeneous platinum loading within CL and GDL porosity at the cathode of a PEM fuel cells numerically. Their obtained results indicated that when the platinum loading at the outlet of the cathode increases, the performance of PEM fuel cells and uniformity of the density did not improve. Shi et al 34 scrutinized the design and optimization of platinum as a CL and the nonuniform distribution of current density inside a PEM fuel cells with the aspects of loads, mass transport rates, and the improving the cell performance. Moreover, Wang et al 35 developed a numerical model to investigate the effects of the cathode flow channel on the local transport and PEM fuel cells for the interdigitated and parallel flow fields. Also, they evaluated the effect of operating voltages on the cell performance for the flow channel cross-sectional area and aspect ratio. A few geometry optimization aspects of a PEM fuel cell have been analyzed by Huang et al. 36 They also investigated the reactant transport rates to the removal rates of liquid water and CL in a nonisothermal and two-phase condition.
Despite all of the above-mentioned studies, the nonuniform catalyst attracted less attention and even in a few cases, the thermal effects have been rarely a matter of concern. In this paper, the feasibility of using a nonuniform model for the distribution of catalyst along the fuel cell channel is investigated. A multiphase mixture model is established to iteratively simulate the transport phenomenon in a PEMFC. Moreover, in the two-phase model, nonisothermal approach is adopted to track the thermal effects. In the simulated CL, the loading of catalyst follows a linear distribution along the x direction according to Figure 1. In order to achieve the optimum case of catalyst loading, different patterns are analyzed.

| MODEL DESCRIPTION
A PEMFC is depicted schematically in Figure 1. It consists of an electrolyte membrane and an anode attached together from both sides. In each section, it includes gas flow channels, GDL, and CL. The surrounded rectangular area denoted by "A-B-C-D" is the flow channel, and the porous medium is denoted by "D-C-F-E" shows the GDL. In this case, the oxygen is humidified and entered through the A-D line with a uniform velocity profile.
In this research, the two-phase transport phenomena are solved with a multiphase model for the mixture. This model is adopted according to those presented by Wang and Cheng. 37

| Multiphase model for the mixture
The multiphase mixture model for a steady-state compressible laminar flow with homogeneous GDL is presented as follows. 38 The continuity for the multiphase mixture is written in Equation (1).
The applied Navier-Stokes equations in the channel and in the GDL as a porous medium 12,39 are presented in the following.
In these equations, u is the intrinsic average velocity in Equations (2) and (3). Mass conservation of species in phase k is expressed as.
where C k and j k are used for the concentration and the diffusive flux of species in the phase k, respectively. The last term J k in Equation (4) refers the species transfer rate in interphases. In such mathematical model, due to the molecular diffusion and/or hydrodynamic dispersion, by using the Fickian form, j k can be expressed as the Equation (5). In this equation, D k is the effective mass diffusivity.
The production rate and the destruction rate of species in one phase must be coupled with other phases. This can be seen in the following equation.
For the multiphase mixture, in order to calculate the mass-averaged concentration in all phases, the summation proposed in Equation (7) can be used.
By substituting the production and destruction rate of species in Equation (4), this equation can be converted to: where the species advection correction factor, , is introduced in Equation (9). In the above equation, the index k represents phases (liquid and gas) and denotes species (which includes water and oxygen in this investigation). s k is the saturation of phase k, which denotes the volume of occupied space by the pores inside the corresponding phase. Accordingly, the summation of saturation of all phases is equal to unity. The term on the left-hand side of Equation (8) indicates that in mixture level, species are advected by a modified velocity field u . It means that the original velocity field is replaced by the modified field. 38 As introduced in Equations (10-12), mass concentration of species is equal to the sum of its mass concentration in each phase.
The mass fractionation of oxygen in the liquid phase is assumed to be zero (C O 2 l = 0). So, Equation (8) for oxygen is rearranged as written in Equation (3) 12 : According to mixture law, the mass flux diffusion of each phase relative to the entire multiphase mixture can be expressed as Equation (14), 38 in which j l is defined in Equation (15). 12 Equation (16) represents the steady-state energy conservation in phase k, where the thermal equilibrium locality among the phases has been considered and k k represents the thermal conductivity of phase k. Also, q k is the heat flux associated with phase k, 37 which is neglected in this study. The summation of the phase energy equations represented by Equation (16) yields the energy conservation equation for in the mixture, 37 defining the mixture enthalpy h as can be seen in Equation (17).
In addition, the energy advection correction factor h is introduced as Equation (18).
Then, Equation (19) represents the enthalpy equation according to the above-mentioned definition for the mixture. 37 Accordingly, the mixture velocity, mass diffusion coefficient, and kinematic viscosity are described as Equations (20)(21)(22), respectively.
In the porous media, the total volume of the pores region is balanced and divided between the two involved phases. It means that both phases affect each other. This phenomenon results in a relative permeability definition. Such permeability for each phase is considered as the Equation (23) through this research. 40 Moreover, the capillary pressure is defined as per Equation (24). 12 In the equation above, according to the Leverett function, 41 J(s) can be described as follow.

| Catalyst layer
The previous models have been considered the layer of catalyst to be nondirectional (homogeneous). 22,42 In the present research, the mathematical model of the CL (platinum) consists of the equations that comprise oxygen mass-current balance (Equation (26)), the electrochemical reduction rate of oxygen (Equation (27)), and the law of Ohm (Equation (28)). 42 The flux of oxygen is also linked to the concentration gradient relation presented by Fick's law (Equation (29)).
In Equations (27) and (28), is the loss of activation, R is the universal gas constant, F is the constant of Faraday, i ref 0 is the reference for the exchange of current density, and c and a are the coefficients of cathodic and anodic transfer, respectively. 22 Besides, a s is a symbol of the total surface area of catalyst per unit volume of cathode in the considered CL. This value is corresponding to the catalyst loading which is denoted by m pt . The dependency of a s on m pt is given by Equation (30).
Here, a 0 is the total area per unit mass of the catalyst. The catalyst loading gradient is exerted by h(x) as an added function of loading distribution of catalyst to Equation (27). It should be noted that the total quantity of catalyst in all test cases remains constant. In this paper, h(x) is defined by Equation (31).
As presented, the catalyst loading gradient varies linearly with t. Figure 2 shows the variations of catalyst loading distribution along x direction for l = 4 cm and t = −0.2, 0 and 0.35 as representative slopes. The effects of different nonuniform patterns of catalyst distribution are examined by choosing t = −0.2, 0.1, 0.15, 0.25 and 0.35 for the cases 1 to 5, respectively. In addition, CL with uniform catalyst distribution is also considered as the base case.
Equations (26) to (29) contain four variables I y , , N O 2 , and C O 2 which are current density, over potential voltage, oxygen molar flux, and oxygen concentration, respectively. These equations are solved using a 4th-order Runge-Kutta method with the following boundary conditions. Equations (32) and (33) are given at GDL/CL interface, and Equations (34) and (35) are to be applied at CL/membrane interface.
In the above relations, act is the overall potential voltage at the membrane and is taken as the input of the model. I y is the current density at the CL/membrane interface and is equal to the total current density produced at the CL.

| Output voltage
The voltage, as the output of the cell, is given in Equation (36).
In this equation, U oc is the open-circuit voltage that is calculated according to the proposed model by Nernst equation 39 as presented in Equation (37).

| Boundary conditions
At the GDL/CL interface, the boundary conditions are set as water and oxygen fluxes, which are given by Equation (44).
According to Equation (45), half of the generated heat is assumed to be created from the irreversibility of the electrochemical reaction that is exchanged as a heat flux to the cathode at the corresponding boundary.
The uniform inlet flow is assumed to enter the channel while the following boundary conditions are dominant at this region.
At the outlet, the following conditions are used as presented by Equation (47).
Finally, wall boundary with a specified temperature is considered at ED, FC, and AB walls in Figure 1.

| NUMERICAL PROCEDURE
Two-dimensional steady-state simulation of the compressible laminar flow in the PEMFC is conducted. The main governing equations and the related boundary conditions are listed in Table 2. The overpotential voltage at the CL/membrane interface (Equation (35)) is adopted as the input of the model. To initiate the calculations, the oxygen concentration at the CL/membrane (Equation (33)) should be assumed. Then, the governing equations at the CL (Equations 26-29)) are solved by a fourth-order Runge-Kutta method to find the local current densities, water, and oxygen mass fluxes (Equation (44)) at the GDL/CL interface. These fluxes are fed back again into the GDL, and then, the channel has the governing equations according to the Equation (1), (2), (3), (13), and (19) with the boundary conditions that are solved to obtain velocity, temperature, and concentration fields. The finite volume method and the hybrid scheme are implemented simultaneously for the discretization of the convective terms according to references. 44 The SIMPLE approach is implemented to obtain the results of the discretized linear algebraic equations. A special grid network (staggered) is also utilized to eliminate any unwanted physical fluctuations in the pressure and velocity fields. The initial guess for the concentration of oxygen at the GDL/CL interface is then corrected iteratively, and the procedure continues to a converged solution. The convergence criteria are defined such that when the average current density difference between the two consecutive iterations is less than 10 −4 A/cm 2 . After a fully converged solution, the cell averaged current density is calculated as the mean value of all the local current densities along the CL. By changing act , the polarization curves are obtained. A nonuniform mesh size with 270 × 64 cells in x-y directions has been also implemented.

| RESULTS AND DISCUSSIONS
The main involved operational parameters and the geometrical dimensions used are listed in Table 3. In the base case, a uniform distribution has been assumed for the catalyst along the CL. The effect of oxygen mass fraction at the inlet on Logarithmic saturation pressure of water 43 log 10P sat = −2.1794 + 0.02953T − 9.1837 × 10 −5 T 2 Saturated density of water vapor Mass fractionation of the gasphase water Effective diffusion Conductivity T A B L E 1 Physical property equations the performance of the fuel cell is depicted in Figure 3. The predicted polarization curves are also compared with those presented by You and Liu 12 in this figure. As presented, a reasonable agreement has been achieved between the experimental data and the predicted curves extracted from the simulations.

| Reactant distribution and current density
The oxygen mass concentration at the overpotential voltage of 0.7 V is shown in Figure 4 and Figure 5 at t = −0.2 and 0.35, respectively. Oxygen mass fraction profiles for t = −0.2, 0, and 0.35 in the CL are also shown in Figure 6. It is perfectly depicted that the oxygen consumption in the CL is due to the corresponding electrochemical reactions. Therefore, such decrement in the concentration of oxygen happens along the channel and through the GDL. After that, it falls down to the smallest amount at the most downstream section of the CL. According to Equation (27), the oxygen concentration and the generated current are directly correlated. Therefore, a decrement in the fuel cell efficiency due to reduction in oxygen concentration could be expected. At low current densities, such loss has negligible impacts. On the other hand, at high current densities, due to the higher reduction in concentration of oxygen, the loss is considerable. In general, moving downstream the channel, O 2 mass fraction decreases, while mass fraction of water increases by water production in cathode CL and its transportation from anode side. In other words, as we move toward the channel outlet, partial pressure of water approaches to saturation value at local temperature that makes water condensation occurs. As shown in Figures 6 and 7, by comparing water mass fraction in different CL slope conditions, it can be seen that for t = −0.2, CL thickness near channel outlet is larger, and hence, constant water mass fraction lines are located closer to outlet compared with those of the same values for t = 0. On the other hand, for t = 0.35, constant water mass fraction lines are shifted toward inlet of the channel in comparison with those of the same values for t=−0.2 and t = 0 due to the decreasing thickness of CL while moving downstream the channel.
As it is shown in Figure 8, at the beginning of the CL, the oxygen concentration is higher for t = −0.2, which has less amount of catalyst at this region. However, this concentration decreases more rapidly along the CL than other cases. It can be inferred that such falling rate can be explained by the increasing trend of the catalyst loading along the CL. Also, the oxygen concentration profile is smoother slope for t = 0.35.
Current density distribution along the CL is depicted in Figure 9. For t = −0.2, the averaged current density is higher, because more catalyst is placed at the downstream of the CL, which compensates the reduction in oxygen concentration. Although oxygen concentration at the end of the CL is higher for t = 0.35, the current density is lower since this case has less catalyst at the end of CL compared with the two other cases.

| Temperature profiles
Thermodynamics describes that the saturated water pressure is dependent on temperature. Thus, the fluctuations in temperature could automatically affect the saturation distribution. Besides, by increasing the temperature, some decrement can be observed in the activation loss and this improves the fuel cell efficiency. Nevertheless, a thermal limitation must be defined by implementing heat management system for the fuel cell, because the lifetime of the membrane could be reduced at high temperatures for which exothermic reactions are the potential cause inside the fuel cells.   Figure 10 shows the temperature profile at the CL. As it is illustrated in this figure, for t = 0.35 the peak of temperature is a bit higher than the two other cases. Because at the beginning of the CL it has a higher local current density, the reaction rate and the released heat are higher in this region, as well. Figure 11 shows the polarization curves for t = −0.2, 0, and 0.35. At low current densities, all cases have nearly the same performance. On the other side, at high current densities, the difference between the cases gets more noticeable. At high current densities, t = −0.2 results in the best performance. By comparing the working states on polarization curves in Figure 11 with those of the base case shown in Figure 3, up to 5% improvement in the voltage generation of the fuel cell has been observed. The effect of the nonuniformity of the CL on the fuel cell output power is also investigated. Figure 12 shows the power density curves for t = −0.2, 0, and 0.35. As can be seen, the maximum power density can be enhanced up to 1.6% for t = −0.2, while it can be decreased by up to 1.8% for t = 0.35 compared with the base case. For t = −0.2, more catalyst is applied downstream the CL where reactant concentration is low. In this case, at the beginning, the current density of CL is lower than the two other cases. However, increment in loading along the channel direction causes the current density to be decreased with a slower pace. Therefore, the case t = −0.2 has a higher averaged current density.

| CONCLUSION
A 2D laminar steady-state nonisothermal catalyst loading model is implemented for a PEMFC to study the effect of nonuniform distribution of the CL. The model is demonstrated to be able to calculate the concentrations in the GDL and CL. It can also estimate the current density variations along the CL direction. The following bullet points present the foremost outcomes of the present research.
• By increasing the loading of catalyst along the channel to compensate the reduction in oxygen concentration, an improvement could be expected in the overall performance of fuel cell. • By reduction in the oxygen content, a decreasing trend in the fuel cell efficiency is expectable, which is called as the mass transport loss. In high level of current densities, the reduction in oxygen content is much higher, and the mass loss is more evident. • At the beginning of the CL, the oxygen concentration is higher for t = −0.2 which has less amount of catalyst in this region, while it is more uniform for t = 0.35. • In this study, catalyst loading varies linearly along the CL.
When t = 0.35, the fuel cell has lower performance compared with the uniform CL. It means that the idea of using more catalyst at high concentration locations is not necessarily working. • The maximum temperature for t = 0.35 is a bit higher in comparison with the other two cases, and hence, the rate of electrochemical reaction and therefore the released heat are of greater value in this case. • It was found that when the catalyst loading increases with a slope of 0.2 (t = −0.2) along the channel, the fuel cell has the maximum efficiency.