Experimental Validation of a 3D-CFD Model of a PEM Fuel Cell

The growing energy demand is inevitably accompanied by a strong increase in greenhouse gas emissions, primarily carbon dioxide. The adoption of new energy vectors is therefore seen as the most promising countermeasure. In this context, hydrogen is an extremely interesting energy carrier, since it can be used as a fuel in both conventional energy systems (internal combustion engines, turbines) and in Fuel Cells (FC). In particular, PEM (Polymeric Electrolyte Membrane) FC are given growing attention in the transportation sector as a Life-Cycle viable solution to sustainable mobility.
The use of 3D CFD analysis of for the development of efficient FC architectures is extremely interesting since it can provide a fast development tool for design exploration and optimization. The designer can therefore take advantage of a robust and accurate modelling in order to define and develop fuel cell systems in a more time-efficient and cost-efficient way, to optimize their performance and to lower their production costs.
So far, studies available in the scientific literature lack of quantitative validation of the CFD simulations of complete PEM fuel cells against experimental evidence. The proposed study presents a quantitative validation of a multi-physics model of a Clearpak PEM cell. The chemistry and physics implemented in the methodology allow the authors to obtain both thermal and electrical results, characterizing the performance of each component of the PEM.
The results obtained, compared with the experimental polarization curve, show that the model is not only numerically stable and robust in terms of boundary conditions, but also capable to accurately characterize the performance of the PEM cell over almost its entire polarization range.


Introduction
The wide use of internal combustion engines and the analysis of their combustion management strategies as reported in Fontanesi, S analysis [1] [2] [3] , have led the way in the search for effective technologies to reduce emissions or their complete abatement. The most effective and concretely feasible is that of the use of hydrogen in Fuel Cells.
The main applications envisaged for hydrogen concern the use as a fuel for the generation of electricity and for transport. Plants for the centralized production of electricity and internal combustion engines fuelled with hydrogen are already feasible based on existing technologies. The main technology whose diffusion will largely influence the affirmation of hydrogen as a clean energy carrier is fuel cell. The environmental impact due to the use of this gas is almost zero, in fact whether it is combusted, or it is exploited electro-chemically, the only reaction products are water and heat, following: 2 + 1 2 2 → 2 + ℎ . However, when it used as a fuel for internal combustion engines, the high-temperature thermal cycle lead to formation of harmful intermediate products such as NOx, while the use in low-temperature fuel cells avoids also this issue. However, hydrogen is an extremely volatile and light gas: compared to other fuels it has a reduced energy content per unit of volume, posing non trivial issues for its onboard storage (realized using highly pressurized tanks), while it shows the highest energy content per unit of mass (approx. 3 times that of gasoline).
Polymer Electrolyte Membrane (PEM) fuel cells are receiving the greatest efforts in terms of investment, research and development, especially due to their potential use in the automotive industry. They present a polymeric membrane with high proton conductivity as electrolyte, while platinum is typically used as a catalyst and operate at temperatures between 60°C and 100°C. These cells are not sensitive to corrosion, thus making it possible to use atmospheric air for the cathodic supply, although a control and thermal management systems are required as well as a tight filtration system, due to the sensitivity of the catalyst loading to the presence of CO (even small ppm).
For correct operation, the membrane must be kept at a certain degree of humidification to guarantee the desired proton conductivity.
In order to achieve this, the gas stream at anode is humidified, and it is essential that the cell temperature does not exceed 100°C. In any case, given the low operating temperature and the absence of a liquid and/or corrosive electrolyte, they are particularly suitable for electric propulsion and for the small size (1 -250 kW) generation/cogeneration. Current research efforts on PEM technology development are oriented in three main directions: innovative materials, dynamic cell modelling and control systems. The most important aspect affecting the current development of this technology is the dynamic modelling of fuel cells, which is also the central theme of this project. It consists in simulating the operation of the cell considering the interplay between fluid dynamic, thermal and electrochemical aspects, allowing the authors to understand in more detail the phenomena involved both at the level of the membrane and of the auxiliary subsystems (compressor, cooling systems, humidifiers, etc.). This brings advantages in terms of system optimization. It should also be noted that an accurate modelling of the gas distribution and of the fluid dynamics inside the cell greatly improves the ability to predict electrochemical processes, with a direct impact on the planning and design of the cell.
Another pivotal aspect in fuel cells modelling are grid requirements, especially regarding the transverse direction (i.e. along the cell thickness) where the highest gradients (e.g. species concentration, thermal and electric fields) develop. Kamarajugadda and Mazumder [4] carried out a comparison on the accuracy and efficiency of different membrane modelling strategies. A grid independence analysis was conducted, showing that at least 5 cells were needed along the electrode thickness and 40 along the membrane thickness. This is confirmed by the study carried out by Choopanya and Yang [5], reinforcing the concept that in a PEM fuel cell the main phenomena (diffusion of gases in the diffusion layers, ion transport through the membrane and electronic transport in the components) occur mainly along the cell thickness, which is the most representative characteristic direction of the phenomena occurring during cell operation and therefore it represents the calculation direction along which the highest resolution must be guaranteed. On the contrary, less stringent grid requirements are needed for the longitudinal direction (i.e. along the channels), which only affects the gas flow in terms of complete development, as well as for the transversal direction of the channels, whose relevance is mostly limited to boundary layer development. Such strategies prove to be convenient in limiting the overall node number (hence the total CPU cost of simulations), with the only exceptions being the high-curvature sections, for which high-resolution grids are needed to limit false flow directions.
There is a fair amount of literature on CFD modelling of PEM fuel cells although, as pointed out in Wang's analysis [6], there is still a need for considerable development in this field in order to obtain models with a high degree of accuracy, reliability and applicability to advanced fuel cell design. The aspects requiring the largest improvements are 1) the development of specific physico-chemical models, 2) advanced numerical algorithms, 3) an accurate materials characterization and 4) detailed model validation against well measured canonical test cases. In this paper, a comprehensive numerical model is presented and validated, and in the authors' opinion it constitutes a solid basis for future model development.

Model Equations
The fluid dynamic modelling of all phenomena developing in a PEM fuel cell is a complex task, due to the simultaneous evolution of both chemico-physical and electrochemical processes. These are generally described by solving and transport equations for mass, momentum, energy, chemical species and electric current. In addition, multi-phase treatment is required for the possible presence of condensed liquid water at cathode, as well as a dedicated treatment is needed for flows in porous regions (Gas Diffusion Layers, GDLs). Material characterization plays a crucial role and a relevant use of empirical correlations is usually adopted in this context (e.g. the water transport within the solid polymer membrane and the related electro-osmotic water drag). In the next paragraphs, the fuel cell peculiar areas of modelling are described [7].

Flows in Gas Diffusion Layers
The conservation equation solved for the mass fraction of each i-th chemical species present in the simulation is expressed in Eq. 1 in its integral form and it is generally valid for both fluid regions (e.g. anode/cathode gas channels) and porous media (e.g. GDLs), thanks to the presence of macro-scale porosity ( , defined as the ratio between the volume of the empty areas over the total volume and equal to 1 in single-phase fluid zones): In Eq. 1 is the volumetric source/sink term due to the reacting flow and the species diffusion, for which a gradientoriented Fick's approach is used. The modelling of flows in GDLs requires a dedicated treatment for porous regions, characterized via macroscopic (or large-scale) properties as common practice in engineering-grade simulations. The diffusivity term is modified for flows in porous regions as in Eq. 2, considering both porosity and tortuosity ( , defined as the ratio between the actual convoluted path length over the straight-line distance). In the present analysis both properties are considered isotropic, as no information on preferential diffusion was available for the cell materials.
The continuity equation in GDLs is modified as in Eq. 3, taking into account the medium porosity and the possible presence of volumetric source mass ( ): Modifications are required also for the momentum equation in porous media, given the additional flow resistance and the action of capillary forces in i-direction (Eq.4). Porous regions are characterized via macroscopic (or large-scale) properties as common practice in engineering-grade simulations, e.g. the macro-scale porosity (χ, defined as the ratio between the volume of the empty areas over the total volume and equal to 1 in single-phase fluid zones). The resistance contribution is represented as a tensor ̿ composed by a viscous ( ̿ ) and an inertial ( ̿ ) component (Eq. 5): The inertial resistance ̿ concerns the kinetic energy that must be supplied to the flow in order to move through the porous material and which represents the energy loss related to the turbulence itself. The order of magnitude of flow velocities in fuel cells GDLs well motivates discarding this contribution. Conversely, the viscous resistance ̿ represents the energy loss that the flow undergoes through the interconnected voids of the porous medium, due to forces of a viscous nature. The isotropic tensor for the -th phase (vapor/liquid) is represented as a function of the phase molecular viscosity ( ), of the phase volume fraction ( ) and of the material permeability ( ) as in Eq. 6: As for the capillary effect, it is accounted for as a momentum body source for each -th phase. The Leverett et al. [8] equation is used for it, although a considerable effort is present on the development of more appropriate formulations [9].

Fuel Cell Electrochemistry
As for electric charge modelling, the transport of current within a fuel cell is described by the charge conservation equation and can be applied to both electric and ionic current, in the general form of Eq. 7: In Eq. 6 Φ is the generic term source, acting at the interface between the electrode and the membrane where the current transfer between the two components is present, while is the electrical conductivity (expressed in / ) of solid components (i.e. polymer membrane, bi-polar plates, solid-phase of the porous GDLs). A modelling difficulty is known for an accurate representation of hydrated membrane conductivity, being dependent on water content ( ). Among the many empirical formulations proposed in literature, the Springer [10] correlation for (Eq. 8) is implemented in the model, as appropriate for PEMFC membranes at 30°C operation and considering water presence by means of its activity ( , i.e. the ratio of water partial pressure to its saturation value): Based on Eq. 8 and on the membrane temperature, the solid electrical conductivity for the polymer membrane from [10] is adopted (Eq. 9): Current conservation leads to the equality between the total current generated in the anodic catalytic layer ( ) and that consumed in the respective cathode side ( ), as expressed by Eq. 10: The electrode current densities and are the result of the electrochemical reactions taking place in the catalytic layers of the respective electrodes. Local current densities are modelled through Butler-Volmer formulations (Eq. 11 and 12), functions of temperature, activation overpotential and active catalyst area per unit volume * : The exchange current densities 0, and 0, vary according to pressure and temperature, and they can be expressed by the generic expression: Finally, the average current density can be modelled as the total electric current value generated by a fuel cell divided by the area of the active cell surface (Eq. 14): The electrodynamic model adopted in this study allows the authors to model the generated electric field and the current density as a function of the electric potential acting. The dynamic model (unlike the static one) allows the authors to consider also magnetic effects present within the system. The model equation is obtained by integrating the charge conservation on the cell domain and equating the flow through a closed surface to the amount of electric charge present inside the volume delimited by the same, being the charge density (Eq. 15): Conservation equations are used to model the transport of electrochemical species involved in the reaction and to consider the effects of ionic concentration on the process kinetics. In this model protons + and 2 + are considered as electrochemical species, and dedicated transport equations in the form of Eq. 1 are solved for their mass fraction .
The volumetric source term is not null only at electrode anode and cathode catalyst layers (ACL and CCL), where species are consumed or generated by the electrochemical reactions. The relative terms of well/source are expressed through Faraday's law and are equal to Eq. 16-18, where and represent the anode and cathode current densities of the electrode, respectively, while 2 and 2 are the molar mass of hydrogen and oxygen. The negative sign of 2 and 2 are due to their consumption at electrodes, whereas water is produced (positive sign for 2 ).

Membrane Model and Fuel Cell Thermal Balance
As for the solid membrane, a solid ion approach is adopted in conjunction with the mentioned electrochemical surface reaction model, and it permits to include reactions taking place in the fluid but involving ions migrating into or out of the solid region (i.e. reactions at CLs). In fact, it is constituted by a porous structure, through which complex ionic and water transport mechanisms take place. The membrane isolates the two fluid streams at anode and cathode (which are also differently pressurized in the analysed fuel cell); hence, it is impermeable to gases.
The generation of heat in the fuel cell is considered in both its contributions: Ohmic and reaction heating. As for the Ohmic heating ( ℎ ), it is generated via Joule effect in a conductor due to the electric current. This is included considering an additional source term in the energy conservation equation, expressed by Eq. 19: The second contribution is the heat generation associated to electrochemical reactions ( ). This is modelled as in Eq. 20, where the overpotential represents the irreversible heat generation, whereas the second term indicates the equilibrium potential temperature derivative: In the presented comprehensive model, both cell fluid-dynamics and electrochemistry are described using all these submodels. As for the physical models related to the fluid dynamics, the conservation of mass, momentum and energy equations are applied by finite-volume integration on the computational grid in a time-independent (steady-state) approach. Regarding the cell electrochemistry, electrodynamic potential is used to model the electric field generated by the difference of potential agent in the cell as a boundary condition; electrochemical surface reaction is used to model the electrode oxidation and reduction reactions and in particular to calculate the generated current; governing equations for the electrochemical species are necessary to define the various chemical species involved in the reactions and to predict their consumption/generation; the solid ion model permits to include the ionic transfer through the membrane; ohmic heating considers the heat generated by the circulating electric current flow; finally, electrochemical reaction heating is able to determine the heat generation attributable to electrochemical electrode processes.

Numerical Model
A Clearpak cell from Pragma Industries is chosen for testing and modelling. This particular fuel cell is fully functional and it is intended for academic purposes primarily; the classic bipolar plates containing the gas distribution channels are in fact replaced by transparent polycarbonate plastic plates that completely enclose the cell, permitting to observe the production and transport of water to the cathode. A three-dimensional grid of the fuel cell is created in the STAR-CCM+ commercial software, consisting in approximately 726.400 finite volume cells. A representation of the computational grid is presented in Figure 1. Two different calculation grids were defined, both structured (directed mesh with hexahedral cells) but characterized by a different degree of resolution, in order to be able to express themselves regarding the influence of mesh resolution on calculation. Describing the procedure in more detail, during the process of defining a structured mesh, a two-dimensional sketch (called a patch) is created on the so-called source surface of each component. This consists in dividing these surfaces into quadrangular elements within which it is possible to control the size of the cells in the plane. Once the patch is defined, it is possible to manage the volumetric distribution of the surface mesh along a guide surface (in this case the lateral one), defining the number of cell layers and the parameters of a possible non-homogeneous distribution; this is done by means of parameters called Stretching function, such as the Thickness ratio or the ratio between the thicknesses of two adjacent cells.
In this way a structured three-dimensional mesh is obtained, such as the one shown as an example in Figure 1 (a). The choice of the characteristic patch dimensions in the various components is dictated by the desire to obtain a mesh that is as compliant as possible. In this way it is possible to guarantee a correct transfer of information between the cells, avoiding problems of false diffusion of the results and increasing the accuracy of the simulation.
The characterization of the gas flow at anodic and cathodic channels reproduced the conditions of the laboratory tests, where the cell was supplied with hydrogen and air flows, both humidified with a relative humidity of 100%. In Table 1 the physical properties foreseen in the whole cell model are collected, whereas the material properties are summarized in Table 2.  The fuel cell polarization curve was measured in the experiments at 100% relative humidity and pressurized flows: 0.5 bar at anode channel, 2 bar at cathode channel. The operating temperature was 27°C and the circulating current density has been cyclically varied from 0 to 2 A/cm 2 (via incremental steps of 0.25 A/cm 2 ) to obtain the corresponding voltage and thus defining the characteristic curve of the cell.
CFD simulations replicated these tests with the aim to reproduce the measured polarization curve. Four potential values are applied, and the resulting current density at electrode collectors is monitored, as reported in Table 3. The initial and boundary conditions in the four tests provided are summarized in the following tables:

Results
The simulation, in the four variants studied, took an average of 700 iterations to achieve convergence, obtaining residues associated with the solution which in any case are smaller in order of magnitude of 1e-3. For all variables, a 2 nd order differencing scheme in space has been adopted.
Three dimensional fields are analysed for each voltage. As visible in Figure 2, hydrogen is consumed from inlet to outlet on the anode side, whereas the same is verified for oxygen on the cathode side. Images are relative to the 0.7 V configuration and fields of mole fractions are extracted at the ACL (anode catalyst layer) and CCL (cathode catalyst layer), respectively, hence at the contact interface between porous diffusion layers and polymer membrane. Reactants diffusion through GDLs allows electrochemical surface reactions to develop on the entire membrane surface, as shown by the almost monotonic concentration gradient from inlet to outlet. Given the importance of water at both anode and cathode side, the same views are reported for vapor water in Figure 3. For both cell sides an increase in water content is observed: however, at anode side this is due to the drop in hydrogen content moving through the cell active surface, leading to a water-enriched outlet gas mixture. On the contrary, at cathode side this is due to water formation because of the oxygen reduction reaction (ORR), thus constituting a proper reaction product. Given the crucial relevance of reactants distribution on the active membrane surface, velocity fields are analysed at both gas channel and diffusion layers. In Figure 4 the velocity magnitude is reported for the anode side on two transverse sections cutting midway of the gas channel and of the porous GDL, respectively. Despite the generally low values of velocity (justifying the model hypothesis of laminar flows), a clear difference is present between the pressure-drive flow in the winding gas channel and the diffusive-driven flow in the GDL. The local rate of electrochemical reactions is monitored. In Figure 5 the reaction flux (expressed in 2 ) for hydrogen and oxygen is reported: as expected, regions where a high reaction rate is present can be identified from both reactants. Relevant reaction flux non-homogeneity is observed: this is attributed to a non-perfect surface irroration with reactants, essentially linked to their diffusive-only transport from gas channel to the active surface. In order to improve this aspect, a re-design of the GDL is under investigation (e.g. different materials, increase in porosity, etc.) for improved reactants delivery. Finally, the results of the four simulated conditions are reported against the measured polarization curve in Figure 6. As visible, a good agreement is found for the voltage-current slope, expressing an encouraging model sensitivity to the variation of the working operating conditions. Still, a current overestimation is present for all points: this is attributed to the lack of the bi-polar plates in the adopted model. Despite their generally high electrical conductivity, they are known to contribute to Ohmic losses and to decrease the simulated net current. An improved model of the same cell including bi-polar plates will be the object of future studies.

Conclusions
In this paper a 3-D CFD model of a PEM fuel cell is created using STAR-CCM+. The comprehensive model permits to include both physico-chemical aspects and electrochemical processes, permitting further development of the model in order to be able to use it as an effective cell design analysis and optimization tool. The results obtained are satisfactory with respect to the polarization curve slope, although an overestimation of current density is present due to the lack of simulated bi-polar plates (hence underestimating the current ohmic losses). The model can simulate the cell operation in different conditions, giving a correct trend of the results in almost the entire polarization range of the device. Transport equations are solved for species and ions, and a careful characterization of material properties is carried out using provided values as well as from literature survey, e.g. for membrane electrical conductivity and its dependence on temperature and water content.
The obtained results allow the authors to gain a model-based understanding of the complex interplay of all processes as well as to have a visual 3D representation of design limitations (e.g. regions of impeded reactant diffusion showing a reaction rate minimum). The model will be object of further evolution in future studies, and it constitutes the basis for design exploration and optimization of PEM fuel cells for automotive applications.