A Three-Dimensional Simulation Model for Proton Exchange Membrane Fuel Cells with Conventional and Bimetallic Catalyst Layers

Abstract

A three-dimensional steady-state model has been developed to study the phenomena that occurs during Proton Exchange Membrane Fuel Cell’s (PEMFC) operation. Electrochemical and transport phenomena on both the anode and cathode sides were investigated. Particular emphasis has been given to the composition and structure of the catalyst layers (CLs), considering parameters such as the metal loading, the most effective specific metal surface, the agglomeration, and the particle size. In this context, two types of CLs were investigated. The first type concerns conventional CLs consisting of Pt/C, while the second type refers to bimetallic CLs consisting of Pt-Ru/C. In both cases, the CLs were examined for various loadings of Pt, Ru, and C to define the optimum atomic ratio for an enhanced PEMFC performance, while, in parallel, possible challenges are intedified. The mathematical model for simulating the entire phenomena and the method for modeling the bimetallic catalyst layers are presented. The results show a good agreement between the model and the experimental data reported in the literature. Additionally, the scenario of bimetallic CLs consisting of Pt-Ru/C with a ratio of 50-50 (Pt-Ru) significantly improved the overall PEMFC electrochemical performance.

1. Introduction

Fuel cells have been recognized globally as ideal candidates for future power generation devices in several sectors, such as in the automotive industry, in distributed power generation, and in portable devices [1]. However, their cost remains high, hampering their market roll-out. Proton exchange membrane fuel cells (PEMFCs) are composed of a membrane electrode assembly (MEA) covered by two gas diffusion layers (GDL) that are placed between current collector plates. The bipolar plates usually incorporate the gas flow channels in the anode and cathode sides of the PEMFC. The electrochemical half-cell reactions (oxidation at the anode and reduction at the cathode) take place on the corresponding catalyst layers (CLs), as depicted in Figure 1. The hydrogen fuel in the anode side is oxidized through an electrochemical half-cell reaction called Hydrogen Oxidation Reaction (HOR) to protons (H⁺) and electrons. The electrons flow through an external electric circuit to the cathode side, while protons are conducted through the polymer electrolyte (Nafion) and transferred to the cathode side. In the air-exposed cathode, oxygen, protons, and electrons recombine to form H₂O. The reaction occurring at the cathode side of PEMFC is called Oxygen Reduction Reaction (ORR). Equations (1)–(3) describe the electrochemical reactions that take place on both electrodes, as well as the overall cell reaction:

Anode electrode (HOR): H₂ → 2H⁺ + 2e⁻

(1)

$$\mathrm{Cathode}\mathrm{electrode}\left(\mathrm{ORR}\right):\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.O_2+2H^{+}+2e^{-}\to H_{2}O$$

(2)

$$\mathrm{Overall}\mathrm{reaction}:{\mathrm{H}}_2+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.O_2\to H_{2}O+electricity+heat$$

(3)

From the first years of studying PEMFCs, noble metals have been incorporated into CLs, with platinum (Pt) being the most applied solution for achieving high electrochemical performance. However, the application of Pt implies some drawbacks. Some of the issues are the anode catalyst being poisoned by carbon impurities contained in the fuel mixture, which can gradually decrease PEMFC’s performance. Additionally, using noble metals such as Pt in both anode and cathode electrode layers increases the overall cost of PEMFC. Thus, researchers have focused on lowering the costs of PEMFCs in an effort to accelerate their marketization prospects. This was attempted by minimizing the platinum loading on both the anode and cathode CLs by applying more cost-effective noble metals such as Ruthenium (Ru), while at the same time maintaining, or even improving, the performance and efficiency of PEMFCs. Several studies have highlighted the ongoing research toward this direction.

Albers et al. [3] examined a series of PtRu/C electrode composites by combining several techniques to study and comprehend their morphological structure, such as transmission electron microscopy, scanning transmission electron microscopy, energy dispersive X-ray microanalysis, X-ray diffraction, and inelastic incoherent neutron scattering. Their study showed that the local composition of the precious metal particles with different sizes and compositions indicated a proportional increase in both particles and the Pt/Ru ratio. At the same time, Ruthenium prevented the agglomeration of the platinum particles from retaining smaller particle sizes. Henry et al. [4] studied carbon-supported platinum-ruthenium (Pt-Ru/C) nanoparticles in anode catalysts in proton-exchange membrane fuel cells to examine the enhanced tolerance in carbon monoxide. The anode catalyst was investigated after a 1000 h aging test under 26 ppm CO, and the structure of the catalyst was analyzed via high-resolution transmission electron microscopy (TEM) images and energy-dispersive X-ray Spectroscopy (EDS). Their analyses showed the dissolution of Ru from the PtRu/C within the micro-porous layer and the membrane. The results indicated that Pt and hydrogen catalyze the Ru reduction within the membrane. In addition, the localization of the precipitation band near the cathode showed that Pt came from the dissolution of cathodic Pt/C and that the hydrogen crossover reduced both Pt and Ru ions. Yu Chen et al. [5] examined a self-made 40-cell PEMFC stack with an active area of 112.85 cm² for each membrane electrode assembly. The anode catalyst was composed of Pt-Ru, and tested under different reformate gases of different CO and hydrogen contents. The study showed that a PEMFC stack could tolerate a CO concentration above 50 ppm under non-diluted H₂. However, it tolerated 10 ppm CO under diluted H₂. Saeed et al. [6] evaluated the performance of a PEMFC with an active area of 25 cm², considering several variables such as the flow pattern, the flow rate, and the degradation of Pt-Ru/C catalyst. In their study, scanning electron microscopy (SEM) images were used for the topography of the electrode. The results showed degradation in the elements of CL due to the corrosion phenomenon as a direct result of the electrochemical reaction among Pt-Ru. The authors concluded that the metal alloys presented an increased current at less negative values, while others showed a decrease in the current to less positive values.

In this regard, more recent studies have examined the operation of PEMFCs, emphasizing the compositions of CLs to identify how the operation of PEMFCs could be affected. Brouzgou et al. [7] investigated the CO tolerance of bimetallic CL consisting of PtMe (Me = Iridium or Palladium) for their application in hydrogen PEMFCs. The durability of the catalyst layers was studied via various techniques. The study concluded that incorporating a second metal on the electrode, such as Palladium, improves Pt’s tolerance. In addition, both electrocatalytic activity and durability towards ORR were enhanced by alloying Pt with another metal. Also, both Palladium and Iridium increased the electrocatalytic activity of pure Pt. It is worth mentioning that accelerated durability tests were also performed for both Pt-Ir/C and Pt-Pd/C catalysts, and the results were promising, indicating stable behavior after 5000 accelerated durability test cycles. In the study of Berova et al. [8], the stability of Ru on Pt catalysts was investigated at fuel-cell operating conditions. The scope of this work was to examine how different shell thicknesses of the CLs affect the degradation behavior. To this end, the catalyst layer with the lower thickness initially exhibited better performance. However, it strongly degraded during the stress tests, ending up with lower performance compared with the catalyst layer with higher thickness. Zhao et al. [9] highlighted the catalyst utilization to minimize the Pt loading. For this purpose, a model was developed to determine how the electrode preparation method, porosity, the dispersion degree of carbon agglomerates, the carbon support size, etc., would affect the performance of the catalyst layer. The study concluded that the parameters mentioned above had an optimum value, and that deviating from these values can lead to electron, proton, and mass-transport issues. Seidfar et al. [10] developed a multiphase non-isothermal pseudo-three-dimensional model to comprehend how four deterministic parameters interact between CL characteristics that provide constant CL porosity. The study concluded that low Pt-loadings increase the oxygen-transport resistance, leading to reduced PEMFC performance. Low ionomer content escalates the ionomer potential losses. Thus, the cathodic overpotential increases due to higher transport resistances.

Besides experimental and modeling studies on PEMFCs, durability stress tests have a vital role in assisting researchers to understand the long-term behavior of the examined CL compositions. Hengge et al. [11] investigated the stability, chemical composition, and structure of a Pt/Ru catalyst at room temperature to simulate existing conditions while ramping up PEMFCs. The results showed that lower maximum-potential values increased stability, while dissolution and dealloying were found to be the main degradation mechanisms, with Ru being dissolved usually. Brkovic et al. [12] highlighted the two crucial factors affecting the commercialization of PEMFCs. These are the significant high prices of catalysts and the degradation mainly by CO. For this purpose, the authors developed a tungsten-carbide-oxide as a new non-carbon-based catalyst support for Pt-Ru (Pt-Ru/WxCyOz), which is based on PEMFC anode CLs. The performance of the developed CL was investigated via cyclic voltammetry, linear scan voltammetry, and rotating disk electrode voltammetry. The developed CL was tested as an anode in a PEMFC. A synthetic reformate gas mixture was employed as fuel feedstock in a PEMFC, which presented a significant power drop of 35.3% for a commercial Pt-Ru/C CL, while for the developed CL anode catalyst, this drop was about 16%. Fan et al. [13] investigated the degradation mechanisms of PEMFCs consisting of Pt black and Pt/C CLs after 100 h of operation. The degradation of Pt black CL was more intense than Pt/C due to the different decay mechanisms. More specifically, the degradation of Pt black is caused by Pt agglomeration and oxidation, leading to increased ohmic and mass-transport resistances. On the other hand, Pt/C CL degradation was mainly attributed to the reduction of the electrochemical surface and carbon corrosion.

The present study aims to develop a three-dimensional model to simulate the electrochemical and transport phenomena taking place on both the anode and cathode sides of a PEMFC, comprising monometallic and bimetallic catalyst layers. Furthermore, crucial parameters will be considered regarding the composition and structure of the catalyst layers, such as the metal loading and percentage, the most effective specific metal surface, the agglomeration, and the particle size, to achieve a more detailed description of the catalyst layers. In this sense, two types of CLs will be investigated. In the case of monometallic catalyst layers, the composition of Pt/C was taken into consideration, while for the bimetallic catalysts, Pt-Ru/C was examined. In both cases, the catalyst layers were investigated for the various loading of Pt, Ru, and C to identify the optimum ratio for achieving an enhanced performance of PEMFC, and simultaneously identify possible weaknesses. The structure of this paper is as follows. Firstly, the geometry of the computational domain will be introduced, followed by a description of the governing equations. Then, the bimetallic catalyst-layer modelling approach will be explained. Finally, the numerical results and conclusions will be addressed.

2. Model Geometry and Description

During the last few years, numerous efforts have been employed to comprehend the phenomena occurring during the operation of a PEMFC and to identify possible solutions for optimizing the performance of the technology. As mentioned previously, the scope of this work is to examine different types of CLs and to investigate if it is possible to minimize the platinum loading in CLs with less expensive noble metals such as Ru. The model aims to study the phenomena (electrochemical and transport) in the entire geometry of the PEMFC consisting of gas flow channels (GFC), gas diffusion layers (GDLs), catalyst layers (CLs), and the electrolyte (membrane) on both the anode and cathode sides. The examined geometry is depicted in Figure 2, indicating the three-dimensional geometry of the PEMFC, and Figure 3 presents the zy geometry and the PEMFC’s components. The three-dimensional domain considered in this study is a single channel unit of a PEM fuel cell that is followed by the description of the governing equations.

The reactant gases of hydrogen and oxygen (supplied as air) are fed to the anode and cathode GFCs, respectively. Detailed information regarding the PEMFC geometry and the operating conditions are summarized in

Table 1. PEMFC geometric dimensions.

Parameter	Value	Unit
Fuel cell length	20	mm
GFC height	1	mm
GFC width	0.78	mm
Rib width	0.90	mm
GDL width	0.38	mm
CL thickness	0.23	mm
Membrane thickness	0.1	mm

and

Table 2. Fuel cell operating conditions.

Parameter	Anode	Cathode
Reactant gas	Hydrogen	Air
Operating pressure, Pa/Pc	1 atm	1 atm
Relative humidity, RHa/RHc	100%	100%
Operating temperature, Ta/Tc	343 K	343 K

. The fuel cell operating temperature and pressure are 343 K and 1 atm. Additionally, 100% relative humidity is applied for both the anode and cathode reactant gases. The conservation equations of species, flow, energy, and charge conservation are solved using the commercial CFD package COMSOL MULTIPHYSICS 5.3 using the finite element method (FEM).

The computational mesh is shown in Figure 4. The number of computational cells was limited to 22,356 in all regions of the PEM fuel cell. More specifically, the meshing in CLs, GDLs, and the electrolyte regions is denser for obtaining a more detailed view of the results in both the electrochemical and transport phenomena. In the literature, there are studies in which the number of computational cells is higher than this study; however, the number of cells in the current work can provide results with significant good convergence at a relatively quick time.

3. Mathematical Model

The assumptions in which the model is developed rely on PEMFC operation under steady-state and non-isothermal conditions. The gases are considered to be compressible ideal gases. The flow on the gas channels is considered laminar. Transport mechanisms are described via the Maxwell–Stefan diffusion model. In porous media, additional transport mechanisms are described via diffusion and mass transfer. There is no crossover of gases and water through the electrolyte. The porous media of the PEMFC are isotropic and homogeneous. The anode and cathode catalyst layers are porous electrodes of small agglomerates dispersed homogeneously [14]. Heat transfer coupled with electrochemical reactions accounting: (1) Charge transport in the electrolyte, (2) charge transport in the solid conductor materials, and (3) activation overpotentials in the electrode reactions.

3.1. Brinkman Equations on Anode/Cathode Gas Channels

The combination of continuity and momentum equations form the Brinkman equations. Brinkman equations are applied to solve the free flow on both porous media and non-porous domains [14,15,16].

$$\mathsf{\rho }\left(\mathrm{u}·\text{∇}\right)\mathrm{u}=\nabla [-\mathrm{pI}+\mathsf{\mu }(\nabla \mathrm{u}+{(\text{∇}\mathrm{u})}^{\mathrm{T}})-\frac{2}{3}\mathsf{\mu }(\nabla ·\mathrm{u})\mathsf{{\rm I}}]+\mathrm{F}$$

(4)

$$\nabla ·\left(\mathsf{\rho }\mathrm{u}\right)=\mathrm{0}$$

(5)

where μ is the dynamic viscosity of the fluid (kg/(m·s)), u refers to the velocity vector (m/s), ρ is the density of the fluid (kg/m³), p is the pressure (Pa), F concerns the influence of gravity and other volume forces (kg/m²·s²), and T is the absolute temperature (K).

3.2. Brinkman Equations for the Anode/Cathode Porous Domains

The equations applied for the anode/cathode porous domain are the following.

$$\frac{1}{{\mathsf{\epsilon }}_{\mathrm{p}}}\mathsf{\rho }\left(\mathrm{u}·\nabla \right)\mathrm{u}\frac{1}{{\mathsf{\epsilon }}_{\mathrm{p}}}=\nabla \left[-\mathrm{pI}+\mathsf{\mu }\frac{1}{{\mathsf{\epsilon }}_{\mathrm{p}}}\left(\nabla \mathrm{u}+{\left(\text{∇}\mathrm{u}\right)}^{\mathrm{T}}\right)-\frac{2}{3}\mathsf{\mu }\frac{1}{{\mathsf{\epsilon }}_{\mathrm{p}}}\left(\nabla ·\mathrm{u}\right)\mathsf{{\rm I}}\right]-\left(\mathsf{\mu }{\mathsf{\kappa }}^{-1}+{\mathsf{\beta }}_{\mathrm{F}}\left|\mathrm{u}\right|+\frac{{\mathrm{Q}}_{\mathrm{m}}}{{\mathsf{\epsilon }}_{\mathsf{\rho }}^2}\right)\mathrm{u}+\mathrm{F}$$

(6)

$$\nabla ·\left(\mathsf{\rho }\mathrm{u}\right)={\mathrm{Q}}_{\mathrm{m}}$$

(7)

where ${\mathsf{\epsilon }}_{\mathrm{p}}$ is the porosity, k is the permeability tensor of the porous medium, Q_m is the mass source, and ${\mathsf{\beta }}_{\mathrm{F}}$ is the Forchheimer drag (kg/m⁴).

The mass source term denotes the mass deposit and creation within the domains and is calculated according to Equation (8).

$${\mathrm{Q}}_{\mathrm{m}}={\displaystyle \sum _{\mathrm{i}}^{}}{\mathrm{R}}_{\mathrm{i},\mathrm{mass}}$$

(8)

where,

$${\mathrm{R}}_{\mathrm{i},\mathrm{mass}}={\mathrm{M}}_{\mathrm{i}}·{\mathrm{R}}_{\mathrm{i},\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{a}\mathrm{r}}$$

(9)

where M_i is the molar mass (kg/mol) and, for the porous electrode, the electrochemical reactions result in species source term and are calculated based on the following equation:

$${\mathrm{R}}_{\mathrm{i},\mathrm{molar}}=-{\sum }_{\mathrm{m}}^{}{\mathrm{a}}_{\mathrm{v},\mathrm{m}}\frac{{\mathrm{v}}_{\mathrm{i},\mathrm{m}}{\mathrm{i}}_{\mathrm{m}}}{{\mathrm{n}}_{\mathrm{m}}\mathrm{F}}$$

(10)

where ${\mathrm{a}}_{\mathrm{v},\mathrm{m}}$ is the specific surface area (m²/m³), ${\mathrm{v}}_{\mathrm{i},\mathrm{m}}$ represents the stoichiometric coefficients of products and reactants, ${\mathrm{i}}_{\mathrm{m}}$ is the current density (A/m²), ${\mathrm{n}}_{\mathrm{m}}$ is the number of participating electrons, and F is the Faraday’s constant.

3.3. Boundary Conditions

The boundary condition applied on the inlet of the anode and cathode gas channels is defined according to Equation (11).

u = u₀

(11)

The following boundary condition is presented at the outlet of the channels.

$$\widehat{{\mathrm{p}}_0} \text{≤} {\mathrm{p}}_0;\mathrm{u}·\mathrm{t}=0;{\mathrm{n}}^{\mathrm{T}}\left[-\mathrm{p}\mathrm{I}+\mathsf{\mu }\frac{1}{{\mathsf{\epsilon }}_{\mathrm{p}}}\left(\nabla \mathrm{u}+{\left(\nabla \mathrm{u}\right)}^{\mathrm{T}}\right)-\frac{2}{3}\mathsf{\mu }\frac{1}{{\mathsf{\epsilon }}_{\mathrm{p}}}\left(\nabla ·\mathrm{u}\right)\mathsf{{\rm I}}\right]\mathrm{n}=-\widehat{p_0}$$

(12)

The symmetry condition describes the no-penetration and shear stress according to the following Equation (13).

$$\mathrm{u}·\mathrm{n}=0;\mathrm{K}-(\mathrm{K}·\mathrm{n})\mathrm{n}=0;\mathrm{K}=\mathsf{\mu }\left(\text{∇}\mathrm{u}+{\left(\text{∇}\mathrm{u}\right)}^{\mathrm{T}}\right)\mathrm{n}$$

(13)

A non-slip boundary condition is assumed for the walls (u = 0).

3.4. Conservation of Species

The governing equation for an individual species is:

$$\nabla ·{\mathrm{j}}_{\mathrm{i}} + \mathsf{\rho }\left(\mathrm{u}·\text{∇}\right){\mathsf{\omega }}_{\mathrm{i}}={\mathrm{R}}_{\mathrm{i}}$$

(14)

where ρ refers to the mixture density (kg/m³), u is the mass averaged velocity of the mixture (m/s), ${\mathsf{\omega }}_{\mathrm{i}}$ is the mass fraction, ${\mathrm{j}}_{\mathrm{i}}$ is the mass flux relative to the mass averaged velocity and R_i is the consumption or production rate of the reactant or product species, respectively.

The multicomponent gas diffusion is described via the Maxwell–Stefan diffusivities. The mass flux relative to the mass average velocity ${\mathrm{j}}_{\mathrm{i}}$ is defined by the generalized Fick equation:

$${\mathrm{j}}_{\mathrm{i}} = -\mathsf{\rho }·{\mathsf{\omega }}_{\mathrm{i}}{\sum }_{\mathrm{k}}^{}{\mathrm{D}}_{\mathrm{i},\mathrm{k}}{\mathrm{d}}_{\mathrm{k}}-{\mathrm{D}}_{\mathrm{i}}^{\mathrm{T}}\nabla \mathrm{l}\mathrm{n}\mathrm{T}$$

(15)

where ${\mathrm{D}}_{\mathrm{i},\mathrm{k}}$ refers to the multicomponent Fick diffusivities (m²/s), d_k is the diffusional driving force on species k, and ${\mathrm{D}}_{\mathrm{i}}^{\mathrm{T}}$ concerns the thermal diffusion coefficients.

For ideal gas mixtures, the diffusional driving force is:

$${\mathrm{d}}_{\mathrm{k}} = \nabla {\mathrm{x}}_{\mathrm{k}}+\frac{1}{\mathrm{p}}\left[({\mathrm{x}}_{\mathrm{k}}-{\mathsf{\omega }}_{\mathrm{k}})\nabla \mathrm{p}-\mathsf{\rho }{\mathsf{\omega }}_{\mathrm{k}}{\mathrm{g}}_{\mathrm{k}}+{\mathsf{\omega }}_{\mathrm{k}}{\sum }_{\mathrm{l}=1}^{\mathrm{Q}}{\mathsf{\rho }\mathsf{\omega }}_{\mathrm{l}}{\mathrm{g}}_{\mathrm{l}}\right]$$

(16)

Considering the ideal gas law and the partial pressures definition for ideal gas mixtures, the diffusional driving force is written as:

p = c·R_g·T

(17)

p_k = x_k·p

(18)

where c is the total molar concentration (mol/m³), R_g is the universal ideal gas constant 8.314 J/(mol·K), p is the total pressure (Pa), p_k is the partial pressure (Pa), ${\mathsf{\rho }}_{\mathrm{k}}$ is the density of species k (kg/m³), and ${\mathrm{g}}_{\mathrm{k}}$ refers to the external force acting on species k. In the case of ionic species, the external force arises due to the electric field, and X_k is the mole fraction.

The mole fraction x_k is given by:

$$\mathrm{x}{\mathrm{k}}_{}=\frac{{\mathsf{\omega }}_{\mathrm{k}}}{{\mathrm{M}}_{\mathrm{k}}}\mathrm{M}$$

(19)

The mean molar mass of the mixture M (kg/mol) is calculated as:

$$\frac{1}{\mathrm{M}}={\sum }_{\mathrm{i}=1}^{\mathrm{Q}}\frac{{\mathsf{\omega }}_{\mathrm{i}}}{{\mathrm{M}}_{\mathrm{i}}}$$

(20)

To account for the effect of the porosity of the porous domains in calculating the diffusivities, the Bruggeman model is applied. The effective transport factor f_e is defined by the porosity ${\mathsf{\epsilon }}_{\mathrm{p}}$ and tortuosity factor ${\mathsf{\tau }}_{\mathrm{F}}$ as is presented below:

$$\mathrm{f}{\mathrm{e}}_{}=\frac{{\mathsf{\epsilon }}_{\mathrm{p}}}{{\mathsf{\tau }}_{\mathrm{F}}}$$

(21)

where the tortuosity factor is calculated as:

$${\mathsf{\tau }}_{\mathrm{F}}={\mathsf{\epsilon }}_{\mathrm{p}}^{-1/2}$$

(22)

The rate expression describing its production and consumption R_i is calculated as follows:

$$\mathrm{R}{\mathrm{i}}_{}={\mathrm{M}}_{\mathrm{i}}{\sum }_{\mathrm{m}}^{}{\mathrm{R}}_{\mathrm{i},\mathrm{m}}-{\mathsf{\omega }}_{\mathrm{i}}{\sum }_{\mathrm{i}}^{}{\mathrm{M}}_{\mathrm{i}}{\sum }_{\mathrm{m}}^{}{\mathrm{R}}_{\mathrm{i},\mathrm{m}}$$

(23)

where i refers to species and m is the type of reaction.

The boundary condition is applied on the inlet of the anode/cathode gas channels, and it is defined considering the mass fractions.

$${\mathsf{\omega }}_{\mathrm{i}}={\mathsf{\omega }}_{0,\mathrm{j}}$$

(24)

where ω_0,j refers to the mass fraction in the entrance of the channels.

The following boundary condition is applied to the domains where convection is assumed to be the predominant driving force of mass through the outflow boundaries.

$$-\mathrm{n}·\mathsf{\rho }{\mathsf{\omega }}_{\mathrm{i}}{\sum }_{\mathrm{k}}^{}{{\stackrel{\mathrm{-}}{\mathrm{D}}}_{\mathrm{ik}}\mathrm{d}}_{\mathrm{k}}=\mathrm{0}$$

(25)

A symmetry boundary condition represents the boundaries where the species concentration is symmetric.

−n·N_i = 0

(26)

Also, a no-flux boundary condition has been added, presenting the boundaries where no mass flows in or out, resulting in a zero total flux.

−n·N_i = 0

(27)

3.5. Conservation of Charge

A charge balance has been applied to describe the electronic and ionic charge transport in the porous domains. Applying Ohm’s law [16] and the effective conductivities using Bruggeman correction, the equations for the anode/cathode catalyst layers are expressed as follows:

$$\nabla ·{\mathrm{i}}_{\mathrm{l}}={\mathrm{Q}}_{\mathrm{l}}+{\mathrm{i}}_{\mathrm{v},\mathrm{total}},{\mathrm{i}}_{\mathrm{l}}=-{\mathsf{\sigma }}_{\mathrm{l},\mathrm{eff}}\nabla {\mathsf{\Phi }}_{\mathrm{l}}\mathrm{and}{\mathsf{\sigma }}_{\mathrm{l},\mathrm{eff}}={\mathsf{\epsilon }}_{\mathrm{l}}^{1.5}{\mathsf{\sigma }}_{\mathrm{l}}$$

(28)

$$\nabla ·{\mathrm{i}}_{\mathrm{s}}={\mathrm{Q}}_{\mathrm{s}}-{\mathrm{i}}_{\mathrm{v},\mathrm{total}},{\mathrm{i}}_{\mathrm{s}}=-{\mathsf{\sigma }}_{\mathrm{s}}\nabla {\mathsf{\Phi }}_{\mathrm{s}}$$

(29)

where i_l is the current density vector for the electrolyte, i_s denotes the current density vector for the solid phase, i_v,total is the local current source, σ_l is the electrolyte conductivity (S/m), Φ_l refers to the electrolyte potential, Q_l source term refers to the electrolyte, σ_s refers to the electrode conductivity (S/m), Φ_s is the electrode potential, and Q_s source terms for the electrode.

$${\mathrm{i}}_{\mathrm{v},\mathrm{total}}={\sum }_{\mathrm{m}}^{}{\mathrm{i}}_{\mathrm{v},\mathrm{m}}$$

(30)

The electrode domains only conduct current and the equation is stated below.

$$\nabla {\mathrm{i}}_{\mathrm{s}}={\mathrm{Q}}_{\mathrm{s}},{\mathrm{i}}_{\mathrm{s}}=-{\mathsf{\sigma }}_{\mathrm{s}}\nabla {\mathsf{\Phi }}_{\mathrm{s}}$$

(31)

The ion current in the membrane is described via the following equation:

$$\nabla {\mathrm{i}}_{\mathrm{l}}={\mathrm{Q}}_{\mathrm{l}},{\mathrm{i}}_{\mathrm{l}}=-{\mathsf{\sigma }}_{\mathrm{l}}\nabla {\mathsf{\Phi }}_{\mathrm{l}}$$

(32)

The Butler–Volmer equation is used to describe the relation between the current and activation overpotential, where the anodic and cathodic terms of the current density expression depend on the local concentrations of the electroactive species at the electrode surface [17,18].

$${\mathrm{i}}_{loc,\mathrm{m}}={\mathrm{i}}_0·\left({\mathrm{C}}_{\mathrm{R}}\exp \left(\frac{{\mathrm{a}}_{\mathrm{a}}\mathrm{Fn}}{\mathrm{RT}}\right)-{\mathrm{C}}_0\exp \left(\frac{{-\mathrm{a}}_{\mathrm{c}}\mathrm{Fn}}{\mathrm{RT}}\right)\right)$$

(33)

where a_a and a_c are the anodic and cathodic charge transfer coefficients, C_R and C₀ are dimensionless expressions, describing the dependence on the reduced and oxidized species in the reaction, n is the overpotential, and F is the Faraday’s constant.

n = Φ_s − Φ_l − E_eq

(34)

where E_eq denotes the equilibrium potential.

The boundary conditions applied for the charge conservation equation are the following. The anode electrode is considered as the electric ground Φ_s = 0. In addition, the electric potential of the cathode electrode is assumed as the cell voltage Φ_s = V_cell. Also, the insulation boundary condition has been added for all the walls. Thus, current vectors for the electrolyte and electrodes are described as follows.

−n·i_l = 0, −n·i_s = 0

(35)

3.6. Conservation of Energy

Heat is generated due to the electrochemical reactions and resistance to the transport of protons, electrons, and species. Heat generation can have a significant impact on the performance of PEMFC. The energy equation governing the heat transport in solids is described as follows [19]:

$$\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}}\mathrm{u}·\text{∇}\mathrm{T}+\text{∇}\mathrm{q}=\mathrm{Q}$$

(36)

$$\mathrm{q}=-\mathrm{k}\text{∇}\mathrm{T}$$

(37)

where ${\mathrm{C}}_{\mathrm{p}}$ is the specific heat (J/kg·K), u is the velocity field (m/sec), q is the heat flux vector, k is the solid thermal conductivity (W/(m·K)), and Q is the heat source (W/m³).

For describing the heat generation in the CLs, a general heat source has been added to calculate the heat generated according to the total power dissipation density.

Q = Q₀

(38)

A heat flux boundary has been added across the boundaries of the membrane electrode assembly (MEA). Convective heat flux has been added and is described by the following equation:

−n·q₀ = q₀

(39)

q₀ = h(T_ext − T)

(40)

where h is the heat transfer coefficient, estimated at 190 (W/(m²·K)), while the external temperature T_ext = 293.15 K and q₀ is the total flux across the selected boundaries.

A thermal insulation boundary was added, indicating no heat flux across the boundary.

−n · q = 0

(41)

More specifically, the heat transfer phenomena are coupled to electrochemical reactions. The irreversible voltage losses can occur due to the following phenomena, (1) charge transport in the electrolyte (joule heating), (2) charge transport in the solid conductor materials (joule heating), and (3) activation overpotentials in the electrode reactions. The charge transport in the electrodes and electrolyte creates a joule-heating source Q_JH.

$${\mathrm{Q}}_{\mathrm{JH}}=-({\mathrm{i}}_{\mathrm{s}}·\nabla {\mathsf{\Phi }}_{\mathrm{s}}+{\mathrm{i}}_{\mathrm{l}}·\nabla {\mathsf{\Phi }}_{\mathrm{l}})$$

(42)

Heating due to the electrochemical reactions can be described via Faraday’s law for the electrodes as follows:

$${\mathrm{Q}}_{\mathrm{m}}=\left(\frac{{\mathsf{\Delta }\mathrm{H}}_{\mathrm{m}}}{{\mathrm{n}}_{\mathrm{m}}\mathrm{F}}-\left(\frac{{\mathsf{\Delta }\mathrm{G}}_{\mathrm{m}}}{{\mathrm{n}}_{\mathrm{m}}\mathrm{F}}-{\mathrm{n}}_{\mathrm{m},\mathrm{tot}}\right)\right){\mathrm{i}}_{\mathrm{m}}$$

(43)

where ${\mathsf{\Delta }\mathrm{H}}_{\mathrm{m}}$ is the enthalpy change of the reaction, and ${\mathsf{\Delta }\mathrm{G}}_{\mathrm{m}}$ is the Gibbs free energy of the reaction. ${\mathsf{\Delta }\mathrm{G}}_{\mathrm{m}}$ is defined as:

$${\mathsf{\Delta }\mathrm{G}}_{\mathrm{m}}={\mathsf{\Delta }\mathrm{H}}_{\mathrm{m}}-T\mathsf{\Delta }{S}_m$$

(44)

where ${\mathsf{\Delta }\mathrm{s}}_{\mathrm{m}}$ is the net entropy change. In addition, the equilibrium potential is related to the Gibbs free energy via the following equation.

$${\mathrm{E}}_{\mathrm{eq},\mathrm{m}}=-\frac{{\mathsf{\Delta }\mathrm{G}}_{\mathrm{m}}}{{\mathrm{n}}_{\mathrm{m}}\mathrm{F}}$$

(45)

The total overpotential is defined as:

n_m,toto = Φ_s − Φ_l − E_eq,m

(46)

Considering the relation:

$$\frac{\mathsf{\partial }{\mathrm{E}}_{\mathrm{eq},\mathrm{m}}}{\mathsf{\partial }\mathrm{T}}=\frac{{\mathsf{\Delta }\mathrm{S}}_{\mathrm{m}}}{{\mathrm{n}}_{\mathrm{m}}\mathrm{F}}$$

(47)

by substituting Equations (44) and (45) into Equation (43) results in the following expression:

$${\mathrm{Q}}_{\mathrm{m}}=\left({\mathsf{\Phi }}_{\mathrm{s}}-{\mathsf{\Phi }}_{\mathrm{l}}-{\mathrm{E}}_{\mathrm{eq},\mathrm{m}}+\mathrm{T}\frac{\mathsf{\partial }{\mathrm{E}}_{\mathrm{eq},\mathrm{m}}}{\mathsf{\partial }\mathrm{T}}\right){\mathrm{i}}_{\mathrm{m}}$$

(48)

The overpotential expression presents irreversible activation losses, and the last term is the reversible heat change due to the net change of entropy in the conversion process.

For an electrode surface, the sum of all individual heat sources of the electrode is calculated as follows:

$${\mathrm{Q}}_{\mathrm{EC}}={\sum }_{\mathrm{m}}^{}{\mathrm{Q}}_{\mathrm{m}}$$

(49)

For a porous electrode, joule heating and electrochemical sources are summed up for a total heat source in the domain.

$${\mathrm{Q}}_{\mathrm{TOT},\mathrm{p}}={\sum }_{\mathrm{m}}^{}{\mathrm{a}}_{\mathrm{v},\mathrm{m}}{\mathrm{Q}}_{\mathrm{m}}+{\mathrm{Q}}_{\mathrm{JH}}$$

(50)

4. Bimetallic Catalyst Layer Model

Efficient catalyst layers are always at the forefront of interest in PEMFCs since producing highly active and durable electro-catalysts is essential, while minimizing the use of noble metals, especially platinum (Pt). On the anode side of PEMFCs, the poisoning of the catalyst layers can have a significant impact that causes a detrimental effect on PEMFC’s overall electrochemical performance. This could happen when the H₂ fuel is not 100% pure and contains carbon monoxide (CO), originating from less costly H₂ production routes_, such as the steam reforming of natural gas. On the cathode side, where O₂ is reduced (ORR), Pt is also the most common solution since it provides fast electrochemical reaction rates. However, the kinetics of ORR is sluggish compared with this HOR. To tackle this drawback, a common solution is the formulation of catalyst layers consisting of high Pt loadings, which increases the cost even more. For this purpose, bimetallic catalysts have a crucial role in turning conventional catalyst layers into a more efficient and economically viable solution. The materials that are typically incorporated in bimetallic catalyst layers are platinum-group metals (PGMs) and iron-group metals. The goal is to enhance the catalytic activity and replace at least some of the Pt amount [20].

In this chapter, a method for simulating bimetallic CLs is proposed. The rationale is to identify how PEMFC’s performance is affected when a second metal, such as Ruthenium, is added to the catalyst layer. Ruthenium is a noble metal that belongs to the platinum group of the periodic table. In this sense, the scope is to examine bimetallic CLs consisting of Pt-Ru that is supported on carbon at different atomic ratios between platinum and ruthenium, and highlight possible discrepancies in the overall performance of PEMFC.

Before starting the simulation procedure, some crucial parameters regarding the technical specification of the bimetallic CL should be calculated. These parameters concern the type of solid solution formed (interstitial or substitutional), the solid mixture’s density (Pt-Ru), the bimetallic CL’s heat capacity, and the electric conductivity. Initially, for identifying the solid solution, formation criteria are considered. One of the main criteria is that a substitutional solid solution is formed when the atomic radius difference between the two metal entities is less than 15%. On the contrary, an interstitial solid solution is formed if the difference is greater than 15%. The atomic radius of the materials examined in the current work is r_Pt = 177 pm and r_Ru = 178 pm. Thus, the difference in the atomic radius between the two materials is less than 15%, forming a substitutional solid solution. In addition, a crucial parameter that is essential to be determined is the crystal structure of the Pt-Ru/C catalyst. According to Angelucci et al. [21], the Pt-Ru/C CL presented a similar crystal structure to Pt/C. More specifically, the Pt/C diffractogram displayed peaks at approximately 39.0°, 46.0°, 67.5°, 81.0°, and 86.0° attributed to the diffraction from (111), (200), (220), (311), and (222) planes of the Pt fcc crystal structure. At the same time, the diffractogram of Pt-Ru/C CL does not present features that correspond to the ruthenium hexagonal structure (hcp). The Ru features displayed peaks at 38.4°, 44.0°, 69.4°, 78.4°, and 85.6° ascribed to the diffraction from (100), (101), (110), (103), and (201) planes of the Ru crystal. Considering the above, Pt-Ru/C appears to have an fcc crystal structure.

The second parameter concerns the density of the Pt-Ru mixture. In the model, three different Pt/Ru atomic ratios will be examined, Pt₅₀Ru₅₀, Pt₇₅Ru₂₅, and Pt₈₀Ru₂₀. In this sense, the values of molecular weight, the atomic ratio between Pt and Ru in the solid mixture, and the particle size should be considered to calculate the density of the solid mixture for each examined case.

Table 3. List of parameters used for the calculation of density [21].

Mixture	EDX (Relative Atomic Composition)		Particle Size (nm)	Lattice Parameter (afcc, nm)
	Pt	Ru
Pt50Ru50	48.65	51.35	3.6	0.387
Pt75Ru25	74.42	25.58	5.9	0.389
Pt80Ru20	78.22	21.78	6.3	0.390

and

Table 4. Molecular weight of Pt and Ru.

Material	Molecular Weight (g/mol)
Pt	195.08
Ru	101.07

list the parameters’ values for all the examined cases. In addition, an example concerning the density calculation for the case of Pt₅₀Ru₅₀ will follow.

Number of moles calculation:

The number of moles is calculated according to Equation (51).

$${\mathrm{X}}_{\mathrm{i}}=\frac{100\mathrm{\%}\mathrm{i}}{{\mathrm{M}}_{\mathrm{i}}}$$

(51)

where ${\mathrm{X}}_{\mathrm{i}}$ refers to the number of moles, %i is the atomic ratio of each material in the solid solution, and M_i is the molecular weight in (gr/mol), while i refers to the metal entity (Pt or Ru). An example for the calculation of Pt moles for the case of Pt₅₀Ru₅₀ is following:

$${\mathrm{X}}_{\mathrm{Pt}}=\frac{100·0.4865}{195.08}=0.249$$

Following the same procedure, the number of Ru moles are calculated equal to ${\mathrm{X}}_{\mathrm{Ru}}$= 0.508.

4.1. Fraction

Mole fraction is calculated according to Equation (52).

$${\mathrm{x}}_{\mathrm{i}}=\frac{{\mathrm{x}}_{\mathrm{i}}}{{\mathrm{x}}_{\mathrm{Pt}}+{\mathrm{x}}_{\mathrm{Ru}}}$$

(52)

where x_i refers to the mole fraction for the metal i (Pt or Ru).

Example of mole fraction calculation for ${\mathrm{x}}_{\mathrm{Pt}}.$

$${\mathrm{x}}_{\mathrm{Pt}}=\frac{0.249}{0.249+0.508}=0.249$$

Following the same process, the mole fraction for Ru is ${\mathrm{x}}_{\mathrm{Ru}}$= 0.671.

4.2. Unit Cell Mass

As mentioned previously, the difference in atomic radius between platinum and Ruthenium is below 15%, indicating that Pt-Ru forms a substitutional solid solution. In this respect, the mass of the unit cell is calculated according to the following equation.

$${\mathsf{{\rm M}}}_{\mathrm{cell}}=\frac{4}{{\mathrm{N}}_{\mathrm{A}\mathrm{v}}}·({\mathrm{x}}_{\mathrm{Pt}}·{\mathrm{M}}_{\mathrm{Pt}}+{\mathrm{x}}_{\mathrm{Ru}}·{\mathrm{M}}_{\mathrm{Ru}})$$

(53)

where N_Av is Avogadro’s number indicating the number of units in one mole of any substance and equals 6.023·10²³. An example of unit cell mass calculation for Pt₅₀Ru₅₀ according to Equation (53) is presented below:

$${\mathsf{{\rm M}}}_{\mathrm{cell}}=\frac{4}{6.023·{10}^{23}}·(0.329·195.08+0.671·101.07)=8.7678·{10}^{-22} \mathrm{g}$$

4.3. Unit Cell Volume

Considering the lattice parameter (a_fcc) from Table 3, the unit cell volume is calculated according to the following equation.

$${\mathrm{V}}_{\mathrm{cell}}={\mathrm{a}}^3$$

(54)

An example of unit cell volume calculation is presented:

$${\mathrm{V}}_{\mathrm{cell}}={(3.87·{10}^{-8})}^3{cm}^3=5.8·{10}^{-23}{cm}^3$$

Considering the previously calculated parameters, the density of the solid mixture is calculated according to the following equation:

$$\mathsf{\rho }=\frac{{\mathsf{{\rm M}}}_{\mathrm{cell}}}{{\mathrm{V}}_{\mathrm{cell}}}$$

(55)

By substituting the values resulting from the calculations of Equations (53) and (54), the density of the solid mixture is ρ = 15.12 g/cm³

The calculated densities for each examined case are presented in

Table 5. Densities for each composition of Pt-Ru/C CL.

Compositions	Density (g/cm3)
Pt50Ru50	15.127
Pt75Ru25	17.779
Pt80Ru20	18.161

.

4.4. Heat Capacity of the Bimetallic CL

Heat capacity is a measurable physical quantity equal to the ratio of the heat added to (or removed from) an object to the resulting temperature change. Specific heat is the amount of heat per unit mass that is required to raise one degree Celsius. The heat capacity of a mixture can be calculated using the rule of mixtures. The new heat capacity depends on the proportion of each component, which can be calculated from mass or volume. For estimating the exact amounts of each material in the CL, crucial information was used from the study of Angelucci et al. [21]. These are the geometric surface area of the electrodes (1.27 cm²) and the metal loading, which in all cases were 0.4 mg of metal per cm⁻² and 20 wt% of metal on carbon [21]. The heat capacity of Pt-Ru/C CL will be calculated by considering Εquation (56).

$${\mathrm{Cp}}_{\mathrm{mixture}}=\left(\frac{{\mathrm{m}}_1}{{\mathrm{m}}_{\mathrm{mixture}}}\right){\mathrm{Cp}}_1+\left(\frac{{\mathrm{m}}_2}{{\mathrm{m}}_{\mathrm{mixture}}}\right){\mathrm{Cp}}_2+\left(\frac{{\mathrm{m}}_3}{{\mathrm{m}}_{\mathrm{mixture}}}\right){\mathrm{Cp}}_3$$

(56)

where m_i is the mass of each material in (mg), ${\mathrm{m}}_{\mathrm{mixture}}$ is the total mass of the mixture in (mg), and cp_i is the heat capacity of each material.

4.5. Electric Conductivity

Electric conductivity is also a crucial parameter for the characterization of CL specifications. It reflects the ease at which an electric charge can pass through a material and is calculated according to Equation (57). In this study, the electric conductivity is calculated considering commercial Pt-Ru/C CLs with the specifications presented in

Table 6. Properties of the commercial CL.

Catalyst Properties	Values
Loading (mg/cm2)	2
Thickness (μm)	215
Electrical resistivity (mΩ·cm2)	10

.

The Electrical Resistivity Is the Inverse of Electric Conductivity:

$$\mathsf{\sigma }=\frac{1}{\mathsf{\rho }}$$

(57)

where ρ is the specific conductivity in (ohm·m), calculated considering the electrical resistivity and the thickness of the CL from Table 4.

$$\mathsf{\rho }=\frac{10\mathrm{m}\mathsf{\Omega }·{\mathrm{cm}}^2}{0.0215\mathrm{cm}}=0.00465ohm·m$$

Considering Equation (57) the electric conductivity is about 215$\frac{\mathrm{S}}{\mathrm{m}}$ for the Pt-Ru/C CL, while the electric conductivity for conventional CLs (Pt/C) is around 222 $\frac{\mathrm{S}}{\mathrm{m}}$. Both CLs present a slight difference in the electric conductivity values, and this is because both materials have similar characteristics since they are categorized in the Platinum Group Metals (PGMs). Additionally, several cases of conventional CLs were simulated in the developed model to identify how the performance and the overall operation of the PEMFC vary when applying different platinum and carbon loadings on both the anode and cathode CLs. The input data were taken from the study of Xing et al. [22]. Each examined case is presented in

Table 7. Physical properties of the examined Pt/C layers.

Parameters	Cases	lCL (m)	rPt (m)	mPt (mg/cm2)	mC (mg/cm2)	%Pt	εPt	εI	εCL	${\mathsf{\alpha }}_{\mathbf{Pt}}^{\mathbf{eff}}\left(\frac{1}{\mathbf{m}}\right)$
Anode	Base Case	5·10−5	2.5·10−9	0.2	0.6	25	0.059	0.3	0.64	1.32·105
	1st Case	5·10−5	2.5·10−9	0.2	0.5	28.5	0.049	0.3	0.65	1.10·105
	2nd Case	5·10−5	2.5·10−9	0.2	1.0	16.6	0.097	0.3	0.60	2.17·105
	3rd Case	5·10−5	2.5·10−9	0.6	1.0	37.5	0.100	0.3	0.59	6.76·105
	4th Case	5·10−5	2.5·10−9	1.0	1.0	50	0.104	0.3	0.59	1.16·106
Cathode	Base Case	5·10−5	2.5·10−9	0.4	0.6	40	0.060	0.3	0.63	2.72·105
	1st Case	5·10−5	2.5·10−9	0.4	0.5	44.4	0.051	0.3	0.64	2.29·105
	2nd Case	5·10−5	2.5·10−9	0.4	1.0	28.5	0.098	0.3	0.6	4.42·105
	3rd Case	5·10−5	2.5·10−9	0.8	1.0	44.4	0.102	0.3	0.59	9.19·105
	4th Case	5·10−5	2.5·10−9	0.4	1.3	23.5	0.127	0.3	0.57	5.70·105

. Furthermore, for each case of Pt-Ru/C CLs, the parameters that were recalculated were (1) the percentage of platinum, (2) the agglomeration formed by the platinum and carbon, (3) the porosity of the CL, and (4) the effective surface area in which the electrochemical reactions of HOR and ORR are taking place, in other words, the active electrochemical zone.

In Table 7, the Platinum and Carbon loadings in anode and cathode catalyst layers are presented for all the examined cases of the monometallic catalysts. More specifically, in the base case, the anode side is composed of Platinum at 0.2 mg/cm² while the platinum loading on the cathode side is at 0.4 mg/cm², and the Carbon loading for both the anode and cathode sides are at 0.6 mg/cm². The 1st case is composed of platinum at 0.2 mg/cm² on the anode side and 0.4 mg/cm² at the cathode side, and carbon loading at 0.5 mg/cm² at both the anode and cathode sides. In the 2nd case, the platinum loading remains the same for both the anode and cathode catalysts as in the 1st case, while the carbon loading in the anode and cathode is at 1.0 mg/cm². In the 3rd case, platinum loading in the anode is at 0.6 mg/cm² and 0.8 mg/cm² in the cathode, while the carbon loading is at 1.0 mg/cm² for both catalyst sides. Finally, in the 4th case, the platinum and carbon loadings are at 1.0 mg/cm² at the anode side^,, while on the cathode side, platinum is at 0.4 mg/cm², and carbon loading is at 1.3 mg/cm².

5. Results and Discussion

5.1. Species Concentration

Figure 5, Figure 6 and Figure 7 present the mass fraction of the species (hydrogen, oxygen, and water) on both the anode and cathode sides along the flow direction. Figure 5 presents the distribution of hydrogen and oxygen mass fractions at a cell Voltage equal to 0.6 V on the anode and cathode compartments (GFC, GDL, and CL) of the PEMFC. The mass fractions of hydrogen and oxygen decrease as the reactant gases pass through the GDL to the CL, indicating the consumption of gases due to the electrochemical reactions. In the cathode compartment, oxygen concentration decreases along with the GFC, indicating oxygen consumption due to ORR in the cathode catalyst layer (Figure 5b). The oxygen concentration in the GDL is much lower than in the GFC. Thus, the oxygen consumed and the amount that diffuses towards the CL that is driven via the concentration gradient balance the oxygen concentration in the CL. On the other hand, hydrogen concentration decreases along with the GFC as it is consumed via HOR. However, the decrease in the GDL is lower than that of oxygen on the cathode side since the diffusivity of hydrogen is higher.

Figure 7 depicts the water molar concentrations at the same cell Voltage = 0.6 V for both anode and cathode compartments. The anode’s water concentration is lower than the cathode side’s concentration. That was expected, considering that water is produced on the cathode side by consuming oxygen due to ORR (Figure 5). Water concentration slightly increases from the gas channel to the catalyst layer on the anode side. This is explained since, on the anode side, the electro-osmotic drag (water molecules that are dragged by the proton traveling from the anode to the cathode side) and the water back diffusion (water is transferred into the membrane due to the water concentration gradient from the cathode to the anode side) is occurred and can cause an increase in water’s concentration.

5.2. Model Validation

In this section, the predicted I-V curve will be compared with data from the study of Atyabi et al. [23], and possible reasons for the observed discrepancies will be identified. Furthermore, a statistical evaluation will be conducted to determine if the discrepancies are within acceptable statistical limits. The statistical parameters considered for the evaluation are the coefficient of determination (R²) and the mean bias error (MBE). The coefficient of determination is used to measure the variability within the observed values from the study and the calculated values that are derived from the developed model, while the mean bias error (MBE) captures the average bias within the observed and calculated values.

Model validation was carried out by comparing the polarization characteristics resulting from the simulation model (blue curve) and the experimental study of Atyabi et al. [23] (green curve), as shown in Figure 8. The numerical model and the data provided by the literature presented comparable results. However, the main difference is observed at the high current density values (at the region of concentration losses), where the developed model provides overestimated results. More precisely, the liquid water is vital at high current densities and can cause liquid water flooding in GDL, disrupting the transmission of reactant gases to the CL. Leading to the deterioration of the overall performance of the PEMFC due to water-flooding phenomena and, at the same time, blocking GDL. Thus, the discrepancy between the developed model and the study of Atyabi et al. [23] relies on the fact that in the present modelling work, it is assumed that the PEMFC operates at single-phase conditions and water remains as steam throughout the entire process. The main reason for this assumption has to do with the fact that water-phase-change phenomena are significantly complex [24,25]. Also, the scope of this work has to do with modeling in detail the CL and the overall performance of PEMFC considering various metal compositions and loadings. Furthermore, the part of activation losses in Figure 8 does not present a steep voltage drop as can be found in the experimental study of Mu et al. [26]. The main reason relies on the fact that in this model, the CL has been modeled as an electrode consisting of small agglomerates. However, in the simulation, the structure-dependent local transport of oxygen inside the agglomerate has not been considered [26].

In Figure 9, the comparison between the values derived from the model and study for estimating the deviation between the two polarisation curves is presented. Considering Figure 8 and Figure 9, it is concluded that the calculated polarization curve is in good agreement with the data provided in the literature. More specifically, R² has been calculated at 0.9989, while the MBE for the current density values is approximately 0.034.

5.3. Polarization Curves Comparison (Pt/C & Pt-Ru/C)

This section will investigate the polarization curves for Pt/C and Pt-Ru/C CLs at various platinum and carbon loadings. The examined loadings for the CL consisting of Pt/C are presented in Table 7. On the other hand, the case of Pt-Ru/C CL has been studied under constant metal loading and metal-on-carbon loading. The parameter that varies is the atomic ratio between platinum and Ruthenium. The data emerged from the study of Angelucci et al. [21], considering a geometric surface area of 1.27 cm² in all cases, a metal loading at 0.4 mg/cm², and 20 wt% of metal in carbon. In Figure 10, the polarization curve for all examined cases, including the base case, is compared to identify how the performance of PEMFC varies. In cases 1, 2, and 3, the platinum loading increases for both anode and cathode CLs from 0.2 to 0.6 mg/cm² for the anode side, while the cathode side increased from 0.4 to 0.8 mg/cm². On the other hand, the carbon loading increases from 0.5 to 1.0 mg/cm² for both the anode and cathode sides. The I-V curve has been improved by increasing platinum and carbon loadings in the first three scenarios. This reflects that higher platinum loadings are causing an increase in the electro-catalytic surface area, which can also boost the electro-activity of both catalyst layers. Additionally, adding more carbon loading can also enhance the electrochemical performance of PEMFC due to the higher amounts of electrons transported from the electrode to the external electrical circuit. In the 4th case, 1 mg/cm² of platinum and carbon were applied on the anode side while the loading was 0.4 mg/cm² for platinum and 1.3 mg/cm² for the cathode side. This scenario reveals that the electro-catalytic surface area increases by constantly increasing the platinum and carbon loading but it does not enhance PEMFC performance. This behaviour indicates an upper limit that improves the overall performance of the PEMFC and, by exceeding these limits, can cause a drop in the performance of the PEMFC. Additionally, in the 4th scenario, the platinum loading at the cathode side was decreased to 0.4 mg/cm². This indicates that higher amounts of platinum loading are more suitable for the cathode CL since ORR is approximately five to six orders of magnitude slower than HOR, which is reflecting in the higher cathodic activation overpotential compared with the activation overpotential of the anode [20,27].

In Figure 11, the polarization curves for the different Pt_X-Ru_100-X/C CLs are presented. In this case, the parameter that was varied was the atomic ratio between platinum and Ruthenium. Three cases were examined: Pt-Ru at a 50-50 ratio, Pt-Ru at 75-25, and Pt-Ru at an 80-20 ratio. From all three scenarios, the case of Pt₅₀-Ru₅₀ presents the best performance compared to the other two scenarios. The case of Pt₇₅-Ru₂₅ and Pt₈₀-Ru₂₀ present almost similar I-V curves, with the Pt₈₀-Ru₂₀ showing a slightly decreased performance compared with the Pt₇₅-Ru₂₅ case. This indicates that by decreasing Ru content, this decreases the overall performance of PEMFC. In all three cases, the parameter that indicates whether the performance of PEMFC improved is the electro-catalytic surface area. In the 1st case, the electro-catalytic surface area is at 6.95·10⁵ 1/m, the 2nd case is at 3.59·10⁵ 1/m, and the 3rd case is at approximately 3.29·10⁵ 1/m. Figure 12 confirms that alloying platinum with a second metal, such as Ruthenium, enhances the catalytic activity and the overall performance of the PEMFC. As indicated in Figure 12, Pt₅₀-Ru₅₀ presents significantly enhanced performance compared with all the other monometallic cases. The only case that had better performance was the 3rd case, in which the platinum and carbon loadings were significantly high, especially the Platinum loading on the cathode side. However, this scenario is not a financially viable solution since these amounts of platinum loading will increase the cost of the PEMFC. To this extent, Pt₅₀-Ru₅₀ could be a promising solution since a high performance is reached using 50% less Platinum when the latter is combined with Ruthenium. More precisely, enhanced performance is feasible in three different ways. The first concerns the bifunctional effects, where the second metal delivers one of the necessary reactants. The second relies on improving the electronic effects in which the second metal operates as a promoter and modifies the electronic properties of the catalytically active metals by affecting the adsorption/desorption capabilities of the catalyst (metal–metal interactions). The third relates to the morphological effects in which the catalytically inert metal changes the distribution of the active sites, opening new reaction pathways [9]. In the current study, the enhancement of PEMFC’s performance is due to changes in the morphological characteristics of the CL, as mentioned previously, in terms of the electro-catalytic surface area.

5.4. Temperature Analysis

Figure 13 depicts the temperature distribution inside the PEMFC’s compartments for all the examined cases of Pt/C and Pt-Ru/C catalyst layers. The temperature distribution was studied under the same operational conditions (cell potential, V = 0.6 V) in all examined cases. The temperature limits presented minor variances. For the Pt/C CLs, as the platinum and carbon loading increases, the operational temperature also increases, attributed to the exothermic character of the cell reaction as well as to the joule-heating that is caused by the increased current densities. However, the temperature limits for the 4th case are slightly decreased, which can be ascribed to the lower platinum loadings (Table 5), causing a lower electrochemical performance and thus a lower production of heat. On the contrary, for Pt-Ru/C CLs, the operational temperature limits are slightly reduced for the cases of 75-25 and 80-20, compared with 50-50. This can be attributed to the enhanced morphological properties (increased electrocatalytic surface area of Pt-Ru/50-50 case), leading to the formation of more active sites, which can result in higher reaction rates and, subsequently, higher temperatures. Nonetheless, it is obvious that higher temperatures are resulted in the anode and cathode catalyst layers compared with other compartments of the PEMFC, since these are the parts of PEMFC in which both the HOR and ORR take place. Thus it is concluded that the rate at which the reactant gases are consumed in the CLs affects the overall temperature in the PEMFC; as the electro-catalytic activity increases the exothermicity of the overall cell reaction, it releases more heat while the increased developed current densities increase the joule-heating and thus the cell temperature. Additionally, the temperature at the cathode side is slightly increased than the temperature on the anode side due to the reversible and irreversible entropy production [28].

5.5. Transport Phenomena

In Figure 14, the convective and diffusive flux in both the anode and cathode compartments of PEMFC via a line graph plot is presented. The green line indicates the diffusive flux magnitude, while the blue one concerns the convective flux magnitude. For this purpose, a cross-section along the z-axis (Figure 15) was applied to examine both fluxes in the various layers of PEMFC. In both cases, it is indicated that the reactant gases in gas channels are transported through convection. This can be obtained by considering the increased values regarding the convection flow. Moreover, it is observed that diffusive flux also coexists to transport gases to the active layers. Nonetheless, as the reactant gases flow through the GLDs and CLs, the convective flux is significantly reduced until it drops to zero, where the mass transport of the reactant gases via diffusion prevails. Furthermore, as the reactant gases flow through the CLs, diffusive flux tends to be zero. This implies that the gases are reaching CL, where the electrochemical half-cell reactions of HOR and ORR occur. The significantly low convective and diffusive flux values indicate the laminar flow of reactant gases. Moreover, the flow on the cathode side is faster than on the anode in order to establish the stoichiometric balance of the reactants on the CL.

6. Conclusions

A three-dimensional computational model was developed to examine the electrochemical and transport phenomena in PEMFC using Pt/C and Pt_XRu_100-X/C electrocatalysts. The model is based on steady-state, single-phase, and non-isothermal conditions. Nonetheless, a significant effort has been made to simulate bimetallic CLs for investigating the operation and performance of the PEMFC when the electrodes consist of Pt-Ru bimetallic composites that are supported on carbon paper. For this purpose, user-defined equations have been employed to calculate several parameters regarding the specifications of CLs, which were presented in detail in Section 3. Initially, the operation of PEMFC was studied considering conventional Pt/C CL with various platinum and carbon loadings. The results showed that as the platinum loading increases, the overall performance of the PEMFC also increases. Furthermore, after exceeding the limit of 1.0 mg/cm², the performance of PEMFC was reduced. Moreover, the operation of the PEMFC was examined for various Pt-Ru atomic ratios, and the results were compared with that of CLs consisting of convectional Pt/C electrodes. The results showed an improved behaviour of the PEMFC when CLs consist of Pt-Ru with a ratio of 50-50. However, the cases of 75-25 and 80-20 did not enhance the performance of the PEMFC and presented a slight decrease in performance. Last but not least, the results were compared with studies in the literature and presented a good agreement with the data that were provided by the developed model.