Prediction of Performance Variation Caused by Manufacturing Tolerances and Defects in Gas Di ﬀ usion Electrodes of Phosphoric Acid (PA)–Doped Polybenzimidazole (PBI)-Based High-Temperature Proton Exchange Membrane Fuel Cells

: The automated process of coating catalyst layers on gas di ﬀ usion electrodes (GDEs) for high-temperature proton exchange membrane fuel cells results inherently into a number of defects. These defects consist of agglomerates in which the platinum sites cannot be accessed by phosphoric acid and which are the consequence of an inconsistent coating, uncoated regions, scratches, knots, blemishes, folds, or attached ﬁne particles—all ranging from µ m to mm size. These electrochemically inactive spots cause a reduction of the e ﬀ ective catalyst area per unit volume (cm 2 / cm 3 ) and determine a drop in fuel cell performance. A computational ﬂuid dynamics (CFD) model is presented that predicts performance variation caused by manufacturing tolerances and defects of the GDE and which enables the creation of a six-sigma product speciﬁcation for Advent phosphoric acid (PA)-doped polybenzimidazole (PBI)-based membrane electrode assemblies (MEAs). The model was used to predict the total volume of defects that would cause a 10% drop in performance. It was found that a 10% performance drop at the nominal operating regime would be caused by uniformly distributed defects totaling 39% of the catalyst layer volume (~0.5 defects / µ m 2 ). The study provides an upper bound for the estimation of the impact of the defect location on performance drop. It was found that the impact on the local current density is higher when the defect is located closer to the interface with the membrane. The local current density decays less than 2% in the presence of an isolated defect, regardless of its location along the active area of the catalyst layer.


Introduction
Advent PBI (type APM) are phosphoric acid (PA)-doped polybenzimidazole (PBI)-based high-temperature membrane electrode assemblies (MEAs) produced and commercialized by Advent Technologies Inc. under license from BASF. They are capable of operating in proton exchange membrane fuel cells (PEMFCs) between 120 • C and 180 • C without external humidification, which renders significant benefits over the low-temperature perfluorosulfonic acid (PFSA) membranes such as Nafion. These benefits include simplified water and thermal management, faster electrode kinetics for both electrode reactions, and considerably improved anode tolerance to CO concentrations up to 3% [1]  The Computational Fluid Dynamics (CFD) methodology is recognized as an essential tool for the complete engineering analysis of PEMFCs systems. Previous numerical simulations of hightemperature PEMFCs containing PA-doped PBI membranes include the models of Cheddie and Munroe [16], Siegel et al. [17,18], and Hu et al. [19].
An objective of this study is to create a CFD model that predicts performance variation caused by manufacturing tolerances and defects of the GDE and which enables the creation of a six-sigma product specification for Advent PBI MEAs. Another objective is to perform numerical simulations that provide quantitative estimations of the impact that various GDE manufacturing defects have on the fuel cell performance. A third objective of the study is to perform a sensitivity analysis of the defect location along the GDE and on the GDE porosity.

Mathematical Model
The 3D computational domain used in the study is shown in Figure 2. It comprises a cathode and an anode gas channel (1 mm  1 mm cross-section, 3 mm distance between adjacent channels), a cathode and an anode GDE each consisting of a macro-porous gas diffusion layer (0.4 mm thick), a micro-porous catalyst layer (30 μm thick), and a membrane (70 μm thick). The domain extends from the channels' inlets to the outlets and from the symmetry plane along the channels to the parallel symmetry plane running between two adjacent channels. This domain-size is sufficient to capture the fundamental issues related to fuel cell operation and the computational results may be extended to describe the operation of an entire HT-PEMFC. The CFD model is presented in Appendixes A and B. It consists of transport equations for mass, momentum, and chemical species (hydrogen in the anode channel and GDE, oxygen, nitrogen, and water vapor in the cathode channel and GDE) and for the conservation of the electrical charge in the ionomer (PA)-phase of the GDE and membrane. At normal operating temperatures and pressures of HT-PEMFCs, water produced at cathode can only The Computational Fluid Dynamics (CFD) methodology is recognized as an essential tool for the complete engineering analysis of PEMFCs systems. Previous numerical simulations of high-temperature PEMFCs containing PA-doped PBI membranes include the models of Cheddie and Munroe [16], Siegel et al. [17,18], and Hu et al. [19].
An objective of this study is to create a CFD model that predicts performance variation caused by manufacturing tolerances and defects of the GDE and which enables the creation of a six-sigma product specification for Advent PBI MEAs. Another objective is to perform numerical simulations that provide quantitative estimations of the impact that various GDE manufacturing defects have on the fuel cell performance. A third objective of the study is to perform a sensitivity analysis of the defect location along the GDE and on the GDE porosity.

Mathematical Model
The 3D computational domain used in the study is shown in Figure 2. It comprises a cathode and an anode gas channel (1 mm × 1 mm cross-section, 3 mm distance between adjacent channels), a cathode and an anode GDE each consisting of a macro-porous gas diffusion layer (0.4 mm thick), a micro-porous catalyst layer (30 µm thick), and a membrane (70 µm thick). The domain extends from the channels' inlets to the outlets and from the symmetry plane along the channels to the parallel symmetry plane running between two adjacent channels. This domain-size is sufficient to capture the fundamental issues related to fuel cell operation and the computational results may be extended to describe the operation of an entire HT-PEMFC. The CFD model is presented in Appendices A and B. It consists of transport equations for mass, momentum, and chemical species (hydrogen in the anode channel and GDE, oxygen, nitrogen, and water vapor in the cathode channel and GDE) and for the conservation of the electrical charge in the ionomer (PA)-phase of the GDE and membrane. At normal operating temperatures and pressures of HT-PEMFCs, water produced at cathode can only be as vapor, therefore the model is single-phase. Unlike low-temperature perfluorosulfonic acid (PFSA)-based membranes such as Nafion ® , PA-doped PBI membranes do not need to be rehydrated during operation to improve their proton conductivity, therefore the hydrogen and air entering the domain are considered dry gasses. Water in the cathode flow-field is considered at equilibrium with the PA aqueous solution in the membrane, and since the electro-osmotic drag coefficient of water in PA-PBI membranes is nearly zero [6], sorption/desorption of water and electro-osmotic discharge of water [20][21][22][23][24][25] are neglected in the model. be as vapor, therefore the model is single-phase. Unlike low-temperature perfluorosulfonic acid (PFSA)-based membranes such as Nafion ® , PA-doped PBI membranes do not need to be rehydrated during operation to improve their proton conductivity, therefore the hydrogen and air entering the domain are considered dry gasses. Water in the cathode flow-field is considered at equilibrium with the PA aqueous solution in the membrane, and since the electro-osmotic drag coefficient of water in PA-PBI membranes is nearly zero [6], sorption/desorption of water and electro-osmotic discharge of water [20][21][22][23][24][25] are neglected in the model.

Baseline Performance of MEAs
The model parameters have been initially fitted to match the baseline [9] of a Celtec ® -P1000 MEA in a fuel cell operating at atmospheric pressure and 160 °C with dry hydrogen and air at 1.2 and 2.5 anode and cathode stoichiometric ratios ( Figure 3). Celtec-P1000 MEA was the predecessor to the Advent PBI MEA (Celtec-P1100W).

Baseline Performance of MEAs
The model parameters have been initially fitted to match the baseline [9] of a Celtec ® -P1000 MEA in a fuel cell operating at atmospheric pressure and 160 • C with dry hydrogen and air at 1.2 and 2.5 anode and cathode stoichiometric ratios ( Figure 3). Celtec-P1000 MEA was the predecessor to the Advent PBI MEA (Celtec-P1100W).
Since this study focuses on the cathode catalyst layer operation, when fitting the model parameters to the experimental baseline, particular consideration was given to the parameters affecting the activation and ohmic regions of the polarization curve (parameters in the Butler-Volmer Equations (A23) and (A25)). The best fit to the experimental baseline was obtained for the model parameters shown in Table 1.  Since this study focuses on the cathode catalyst layer operation, when fitting the model parameters to the experimental baseline, particular consideration was given to the parameters affecting the activation and ohmic regions of the polarization curve (parameters in the Butler-Volmer Equation (A23), (A25)). The best fit to the experimental baseline was obtained for the model parameters shown in Table 1. The model values for the Tafel slope and the equilibrium potential in Table 1 represent theoretical values at 160 °C. The model value for the membrane conductivity is significantly lower than the value measured at 160 °C in Reference [9]. The measured value [9] corresponds to a membrane with a doping level of 32 mol PA per PBI repeat unit that results from the sol-gel process before hot-pressing it into the MEA. During hot-pressing some conductivity is lost. Thus, it is reasonable therefore that the fit value for the proton conductivity in the present model is equal to the actual conductivity of a membrane in the MEA. Figure 3 shows a good agreement between the numerical and experimental polarization curves in the activation and ohmic regions, with a small deviation in the highest current density region. This The model values for the Tafel slope and the equilibrium potential in Table 1 represent theoretical values at 160 • C. The model value for the membrane conductivity is significantly lower than the value measured at 160 • C in Reference [9]. The measured value [9] corresponds to a membrane with a doping level of 32 mol PA per PBI repeat unit that results from the sol-gel process before hot-pressing it into the MEA. During hot-pressing some conductivity is lost. Thus, it is reasonable therefore that the fit value for the proton conductivity in the present model is equal to the actual conductivity of a membrane in the MEA. Figure 3 shows a good agreement between the numerical and experimental polarization curves in the activation and ohmic regions, with a small deviation in the highest current density region. This discrepancy is attributed to the mass transport limitations induced by the particular flow-field used in the model, having dimensions different from the undisclosed flow-field dimensions in [9]. The calculation point for the subsequent analysis was selected at 0.5 V (0.9 A/cm 2 ) which corresponds to a relatively high power density and which exhibits good agreement between model and experimental results.

Study of Agglomerate Distribution or Uncoated Regions in the Cathode Catalyst Layer That Would Cause a 10% Drop in Performance
During the catalyst layer fabrication process, the carbon supported platinum (C/Pt)-poly(tetra)fluoroethylene (PTFE) mixture organizes itself into clusters with a bimodal pore-size distribution in which the PTFE binder covers a fraction of catalyst sites. Only the Pt sites covered by PA are accessible by proton-bearing complexes and are electrochemically active. The catalyst layer fabrication process results inherently into a number of defects consisting of agglomerates of the order of 1 µm, in which the Pt sites cannot be accessed by PA. To predict performance variation caused by manufacturing tolerances and defects in the cathode catalyst layer that result from agglomerates or Energies 2020, 13, 1345 6 of 14 uncovered areas, a pristine catalyst layer was compared numerically with a hypothetical one prepared from an identical C/Pt ink (same Pt wt %) with the same catalyst loading (mg Pt/cm 2 ) but exhibiting defects which cause a drop in performance. To achieve this, the effective catalyst area per unit volume (cm 2 /cm 3 ) was changed in the Butler-Volmer Equation (A23) until a 10% drop in performance was obtained at the calculation point. When the defects are much smaller than the characteristic size of the catalyst layer (30 µm) and when they are uniformly distributed, the total volume of defects that cause the performance decay can be calculated as: where a (cm 2 /cm 3 ) is the effective catalyst area per unit volume of the pristine catalyst layer (Table 1), a d is the corresponding value for a catalyst layer with uniformly distributed defects, and V is the catalyst layer volume. Note that while the reference exchange current density, i re f 0 usually cannot be deconvoluted in calculations from the effective catalyst area per unit volume, a, it represents a constant which reduces in Equation (1). Table 2 compares the effective catalyst area per unit volume for the pristine and defected catalyst layers, the total volume of defects, and the number of defects per unit active area that induce a 10% drop in performance at the calculation point. The total volume of defects that cause a 10% performance drop at 0.5 V would represent 39% of the total catalyst layer volume, being much larger than the volume of defects achievable by Advent Technolgies, Inc., with its automated process capabilities, so it represents a massive and easily detectable defect. The number of defects per unit active area in Table 2 was calculated as: where n, A act , t cl , and V 1d represent the number of defects, the catalyst layer active area, the catalyst layer thickness, and the volume of a single defect respectively. The value in Table 2 was calculated based on a 30 µm thick catalyst layer having spherical defects of 1 µm diameter. Figure 4 illustrates the current density distribution along the membrane at 0.5 V and 0.9 A/cm 2 . The current density follows the oxygen reduction reaction (ORR) rate in the cathode catalyst layer, being maximum in the vicinity of the inlet above the channel and minimum in the vicinity of the outlet above the land. Figures 5 and 6 show the current density profiles across the ionomer (PA)-phase of the MEA at locations marked 1 and 2 in Figure 4, for a pristine catalyst at different operating regimes and for a cathode catalyst with defects causing 10% performance decay at 0.5 V. Figures 4 and 5 indicate that the ORR rate depends on the location along the active area and across the cathode catalyst layer thickness. This suggests that the impact of an isolated defect on the performance decay may depend on its location in the catalyst layer.

Reference Exchange Current Density X Effective Catalyst
Area/Unit Volume, (A/cm 2 ) x 10 -4      Figure 4. Open figures correspond to the pristine catalyst layer at different operating regimes; dark triangles correspond to a defected catalyst layer causing 10% performance decay at 0.5 V.
Energies 2020, 13, 1345 8 of 14 Figure 5. Current density profiles along the catalyst layers and membrane in region marked 1 in Figure 4. Open figures correspond to the pristine catalyst layer at different operating regimes; dark triangles correspond to a defected catalyst layer causing 10% performance decay at 0.5 V. Figure 6. Current density profiles along the catalyst layers and membrane in region marked 2 in Figure 4. Open figures correspond to the pristine catalyst layer at different operating regimes; dark triangles correspond to a defected catalyst layer causing 10% performance decay at 0.5 V.
The sensitivity analysis for the location of an isolated defect was performed by selectively canceling the ORR rate in control volumes (250 μm  50 μm  5 μm) along and across the cathode catalyst layer. Since the selected control volumes are larger than a typical defect resulting in the manufacturing process, this sensitivity analysis offers an upper bound for the estimation of the impact of the defect location. Table 3 summarizes the change in the local current density (expressed as a percentage of the local current density in a pristine GDE) when an isolated defect is located at different distances from the membrane-catalyst layer interface. The local current densities were sampled in the membrane along a line normal to the GDE and passing through the defect. The results indicate that the impact on the local current density is higher when the defect is located closer to the interface with the membrane. The local current density decays less than 2% in the presence of an The sensitivity analysis for the location of an isolated defect was performed by selectively canceling the ORR rate in control volumes (250 µm × 50 µm × 5 µm) along and across the cathode catalyst layer. Since the selected control volumes are larger than a typical defect resulting in the manufacturing process, this sensitivity analysis offers an upper bound for the estimation of the impact of the defect location. Table 3 summarizes the change in the local current density (expressed as a percentage of the local current density in a pristine GDE) when an isolated defect is located at different distances from the membrane-catalyst layer interface. The local current densities were sampled in the membrane along a line normal to the GDE and passing through the defect. The results indicate that the impact on the local current density is higher when the defect is located closer to the interface with the membrane. The local current density decays less than 2% in the presence of an isolated defect, regardless of whether it is located in a region of higher or lower reaction rates (region 1 or 2 in Figure 4). It was found that the impact of the isolated defect on the operating (average) current density was negligibly small. Table 3. Change in local current density for the gas diffusion electrode (GDE) with an isolated defect, expressed as a percentage of the local current density in a pristine GDE, at 0.5 V for different locations of the isolated defect.   Figure 7 shows the polarization curve for a PEM fuel cell equipped with the baseline MEA ( Figure 2) and for fuel cells equipped with GDEs having a lower porosity (75% and 50% of the baseline porosity). The power loss for the fuel cell equipped with GDEs having 75% of baseline porosity is Energies 2020, 13, 1345 9 of 14 between 2% at 0.4 A/cm 2 to 20% at 1.3 A/cm 2 , with a power loss of less than 8% at the nominal operating point of 0.9 A/cm 2 (Figure 8). The fuel cell equipped with GDEs having 50% of baseline porosity shows increased mass transport losses with a limiting current density at about 1 A/cm 2 . The power loss induced by the GDE with 50% of baseline porosity is between 5% at 0.4 A/cm 2 to 100% at 1.3 A/cm 2 ( Figure 8).

93
0.3 0.0 Figure 7 shows the polarization curve for a PEM fuel cell equipped with the baseline MEA ( Figure 2) and for fuel cells equipped with GDEs having a lower porosity (75% and 50% of the baseline porosity). The power loss for the fuel cell equipped with GDEs having 75% of baseline porosity is between 2% at 0.4 A/cm 2 to 20% at 1.3 A/cm 2 , with a power loss of less than 8% at the nominal operating point of 0.9 A/cm 2 (Figure 8). The fuel cell equipped with GDEs having 50% of baseline porosity shows increased mass transport losses with a limiting current density at about 1 A/cm 2 . The power loss induced by the GDE with 50% of baseline porosity is between 5% at 0.4 A/cm 2 to 100% 1.3 A/cm 2 (Figure 8).

Conclusions
This work develops a computational fluid dynamics model that predicts performance variation caused by manufacturing tolerances and defects of the cathode GDE of an Advent PBI MEA. The model was used to predict the total volume of defects and the defect distribution that would cause a 10% drop in performance. It was found that a 10% performance drop at the nominal operating regime (0.5 V and 0.9 A/cm 2 ) would be caused by uniformly distributed defects totaling 39% of the catalyst

Conclusions
This work develops a computational fluid dynamics model that predicts performance variation caused by manufacturing tolerances and defects of the cathode GDE of an Advent PBI MEA. The model was used to predict the total volume of defects and the defect distribution that would cause a 10% drop in performance. It was found that a 10% performance drop at the nominal operating regime (0.5 V and 0.9 A/cm 2 ) would be caused by uniformly distributed defects totaling 39% of the catalyst layer volume, or a concentration of more than 46 × 10 6 defects per square centimeter (~0.5 defects/µm 2 ), much larger than the defect concentration achievable by Advent Technologies Inc. with its automated process capabilities.
The study provides an upper bound for the estimation of the impact of the defect location on performance drop. It was found that the impact on the local current density is higher when the defect is located closer to the interface with the membrane. The local current density decays less than 2% in the presence of an isolated defect, regardless of its location along the active area of the catalyst layer.
Sensitivity analysis for GDE porosity shows that the power loss of a fuel cell equipped with GDEs having 75% of baseline porosity to be between 2% at 0.4 A/cm 2 to 20% at 1.3 A/cm 2 , and for a fuel cell equipped with GDEs having 50% of baseline porosity to be between 5% at 0.4 A/cm 2 and 100% at 1.3 A/cm 2 .  In cathode catalyst layer: In anode catalyst layer:

. Momentum Conservation
In anode and cathode channels: In anode and cathode GDLs: In cathode catalyst layer: In anode catalyst layer: