Modeling Electrochemical Performance of the Hierarchical Morphology of Precious Group Metal-Free Cathode for Polymer Electrolyte Fuel Cell

Mass and momentum are conserved in the PM for oxygen, and water vapor using Fick's law and Darcy's law shown in Eq. 29, Eq. 30 and Eq. 33. Liquid water is also modeled in the PM using Darcy's law (Eq. 34). In general, The effective diffusion coefficients, D^eff_{jk, i}, are calculated as shown below:

where τ_i is the tortuosity of the domain i and D_jk is the molecular diffusion coefficient of species j in k. For GDL and MPL, the diffusion coefficient, D_jk, is calculated using Eq. 37, while for CCL domain diffusion coefficient is calculated to represent the hierarchical morphology. In Eq. 37 P_{C, j}, T_{C, j} and M_j are critical pressure, critical temperature and molecular weight of species j.

In our hierarchical model of the cathode, the effective diffusion coefficient used in Eq. 29 is calculated as combination of the diffusivities in the macro- and micro-pore

where D is diffusion coefficient and ε is pore volume fraction. The subscripts micro and macro refer to the micropore and the macropore domains, respectively. The results of our prior pore-scale simulations on high resolution nano-CT images¹⁸ are used for the effective O₂ diffusion coefficient inside the agglomerates, D^eff_micro. In that modeling, we used local pores diffusion coefficients estimated using the Bosanquet formulation for Knudsen transition-regime diffusion. For the larger macropores, Knudsen diffusion is neglected as the pore sizes are generally in the range of 1 μm and larger.

The effective gas permeability, k^eff_i, is a product of the saturated permeability k^sat_i and the relative permeability, k^rel_i.

The relative permeability is dependent on the saturation of the porous media. Here we use the common used cubic dependence relative permeability on saturation, s; for gas is k^rel_i = s_i³ for liquid and k^rel_i = (1 − s_i)³ for gas. The saturated permeability is obtained using Carman-Kozeny relationship³⁰ using representative diameters including those from the nano-CT imaging of the cathodes. The effective permeability for the CCL used in Eqs. 33 and [34] are evaluated using the same form of expression as Eq. 2, but with permeability instead of effective diffusivity parameters. The transport properties in both the micropores and macropores are separately corrected as a function of saturation in each distinct pore size distribution.

The saturation in each domain is determined by Van Genuchten water retention curve model^31,32 (Eq. 38) that relates the saturation to capillary pressure, p_c, defined as

where nw and w are representing non-wetting and wetting phase. The saturations are separately computed in the cathode for the micropore and macropore domains. The pore size distribution obtained from the nano-CT imaging is used to fit a log-normal distribution to obtain a relation between the capillary pressure and saturation.³³ The Van Genuchten fitting parameters p_CB, n and m for GDL and MPL are reported in Supplementary Table S1. Both the micropore and macropore is assumed to be hydrophilic because of Nafion and the graphitized carbon in the electrode. Therefore, with increasing liquid water saturation, the micropores fill first and the liquid water fills the macropores.¹⁸ The water retention curve for the CL is shown in Supplementary Figure S3.

The effective saturation used for estimating relative permeability are corrected for the residual saturations, s_res, that we assume are equal to 0.1:

The water condensation and evaporation are implemented as source terms in Eqs. 29, 33 and 34. The rate of phase change in domain i is represented as

where R^evap_i and R^cond_i is the rate if evaporation and condensation respectively as is shown in Eq. 39 and Eq. 40. In both the rate equations, p_sat is the saturation pressure (Eq. 41) and switch is a mathematical operator defined as

The temperature is obtained by solving energy balance in Eq. 35. Sources of heat include the electrochemical reactions, latent heat of condensation and evaporation, adsorption or desorption of water dissolved in Nafion, and Joule heating due to proton and electron conduction. The effective thermal conductivities, k^eff_{T, i}, in the PM are estimated from volume averaged mixtures of the solid and fluid thermal conductivities as shown in Eq. 36.

Proton and ionomer dissolved water are modeled in the membrane domain and CLs' ionomer binder using Eq. 27 and Eq. 31. The proton conductivity correlation in Eq. 42 includes both water vapor and liquid equilibrated states. The proton conductivity of the CLs is corrected using the equation:

where ε_i and τ_i are the volume fraction and the tortuosity of the Nafion respectively in the domain i. The value of the volume fraction and the tortuosity in the membrane is 1, while these parameters for the cathode are obtained from the nano-CT image analysis (Table S2)¹⁸ or parametrically varied.

Our correlation for the diffusion coefficient of the absorbed water (Eq. 44) is a fit to the data of Zawodzinski et al.³⁴ and is implemented as a continuously differentiable correlation for improved implementation in the computational model. The exchange of water between the vapor and dissolved Nafion phases is implemented as a source term as shown below

where λ_eq is membrane equilibrated water (Eq. 45) that depends on water vapor partial pressure in the gas phase, v is sulphonic acid site density in Nafion, and h_AD is the rate of absorption or desorption (Eq. 47).³⁵

The ACL is considered to be a Pt/C electrode and the HOR is modeled using the Butler-Volmer equation as a rough approximation:

where i^ref_HOR is the reference exchange current density, a_Pt is the electrochemically active Pt surface area in the anode and α is the charge transfer coefficient. The anode overpotential, η_a (Eq. 49), is the potential difference between Pt and Nafion assuming a reversible potential of zero.

To model the PGM-free CCL, the agglomerate model is used to resolve coupling between oxygen diffusion and the ORR within an idealized spherical agglomerate with a closed-form, analytical expression.²⁵ The advantage of the agglomerate model is that it captures transport losses within the microstructure implicitly without have to directly model the agglomerate transport. Following the work of Epting et al.,²⁸ the current density for a size distribution with N bins of diameter is

and the effective current density is a summation of all the diameter weighted by the volume fraction. Figure S1 shows the distribution of the agglomerate diameter obtained from nano-CT for three different Nafion loadings. The effectiveness factor, Eⁱ_r (Eq. 52), is the ratio of actual reaction rate to the ideal reaction rate and depends on the agglomerate's Thiele modulus, ϕ_i (Eq. 53) that represents the ratio of the diffusion rate to reaction rate and is a function of the agglomerate radius, rⁱ_agg. It should be noted that the agglomerate model assumes a spherical agglomerate, however in reality the agglomerates are non-sperical as observed using nano-CT. Instead of the conventional inscribed sphere method to estimate the diameter from CT images, we have utilized the ratio of volume to surface (d = 6 V/A) and the effective diffusion coefficient is estimated from the high-resolution nano-CT images of the agglomerates via computational modeling.¹⁸ In Eq. 11, F,H and represent Faraday's constant, Henry's constant (Eq. 54) and partial pressure of oxygen respectively. f_i and aⁱ_agg are the volume fraction and agglomerate surface area per unit volume, respectively, for an agglomerate with radius rⁱ_agg. Furthermore, δ is the Nafion film thickness, (Eq. 55) is the diffusivity of O₂ in the Nafion, and ε_agg is the electrode's volume fraction of agglomerates not including ionomer film. In Thiele modulus, the effective diffusion in the agglomerate, , is corrected for the saturation in the micropore as follows

where s_microand is the saturation and diffusion coefficient in the micropore respectively.

In Eqs. 11 and 53 k_c is the ORR rate constant defined using a Butler-Volmer equation modified for PGM-free catalyst modeling:

In this approach, we have replaced the typical volumetric exchange current density with R_sJ_s, which is the electrochemically active site density of the cathode multiplied by activity per site. The parameter is the reference oxygen concentration for the R_sJ_s estimation, α_c is the cathodic charge transfer coefficient, and η_c(Eq. 50) is the representative local cathode overpotential.

The flux of protons is zero at the MPL|CL interface and is represented as below

The electric potential on the anode plate is set to the 0 V, while the cathode plate potential is set to the fuel cell voltage, V_FC. To generate the polarization curve V_FC is varied from 1 V to 0.8 V in 0.02 V steps and 0.8 to 0.2 V in steps of 0.05 V.

The relative humidity of both anode and cathode gas feeds are at 100% unless specified. The water vapor concentration at the channel|GDL interface is:

Similarly, the oxygen concentration at the channel|GDL is set as:

where χ = 0.21 if the cathode feed is air and χ = 1 for pure O₂ cathode feed. The flux of gas species across the membrane is zero.

Nafion dissolved water flux at the CL|MPL interface is set to zero.

The gas channel pressure is set to the feed pressure and liquid pressure at the channel|GDL interface is set such that the capillary pressure corresponds to a saturation of 0.0025.

For thermal modeling, the temperatures at the outer surfaces at the top and bottom of the land domains was set to the fuel cell temperature:

A convection boundary condition was applied at the channel interface with the land and the GDL:

where h_{A, C} is the convective heat transfer coefficient and T_gas is temperature of the gas feeds.

Charge, mass and energy conservation are modeled using second order discretization in commercial FEM package (COMSOL Multiphysics v5.1, Comsol Inc., Burlington, MA). The model was run in steady state condition. The potentiostatic parametric sweep was run in COMSOL interface, while the results were exported through MATLAB (Mathworks, Natick, MA). The mesh generated had 34400 triangular elements. The mesh density is higher in the both the CL and membrane domain. The density of mesh is also higher at the CL|MPL and CL|MEM interface.

This section discusses the baseline scenario for validation of the model. The experimental data of the CM-PANI-Fe electrodes with the 35, 50, and 60 wt% Nafion loading was presented previously by Zelenay et al.³⁶ Figure 3 presents both the experimental and model results for CM-PANI-Fe with 35 wt% Nafion loading with different pure O₂ partial pressures. The current density at 0.8 V for the 2 atm partial pressure case was matched between the experimental and model results to fit the R_sJ_s parameter. Comparing with the estimated value of R_sJ_s from the model to the experimental estimation from TOF and active site density by Sahraie et al.¹² has an agreement within a factor of 2, which is reasonably good agreement considering the present challenges in estimating active site properties of PGM-free catalysts. The calculation of this comparison is shown in the supplemental material. In Figure 3, for the 0.3 and 1 atm partial pressure cases the fit parameter, R_sJ_s, was not modified and the polarization curve from the model shows an excellent agreement with the experimental data. This agreement with pure O₂ indicates that the kinetic parameters and the oxygen partial pressure dependence are well approximated by the model.

Zoom In Zoom Out Reset image size

Figure 4 shows the model prediction for CM-PANI-Fe cathode with varying Nafion loading with air as the cathode feed gas. Table V lists the morphological and transport parameters that were changed. In addition, properties computed from the pore and agglomerate size distributions presented in the supplemental information were also updated for each cathode. For the 35 wt% loading, the value of R_sJ_s was not varied from the prior pure O₂ study. Due to the difference in activity observed in higher Nafion loading owing to improved Nafion distribution within the agglomerates,¹⁸ the R_sJ_s, parameter was refit at 0.8 V for the 50 and 60 wt% loading cathode. In other words, the the R_sJ_s parameter was modified to account for the fraction of catalyst that was electrochemically active. In our prior work, we showed evidence of the lower activity of the 35% wt Nafion electrode being due to agglomeration of the Nafion in the macropores and minimal intrusion of the Nafion into the agglomerates. No other parameters were fitted beyond the R_sJ_s parameter and all other parameter changes are those extracted from the morphological and transport analyses of the nano-CT images.

Zoom In Zoom Out Reset image size

As Figure 4 shows, the model captures well the change in the polarization curves with varying Nafion loading. In both the O₂ and air case, the model prediction of the current density in the ohmic region is lower than the experimental values. This is likely due to the model being a steady-state approximation, where the experimental result is a transient process. Thus, the model evaluates the inifinite time water saturation at these low currents, which is likely higher than one would encounted in experiments. Nevertheless, the differences in the approximate limiting current with Nafion loading are well captured.

To understand the differences between the polarization curves with varying Nafion loading, we investigate the spatial distributions of properties predicted in the cathode to resolve the underlying sources of the mass transport resistance. First, we note that the model predicts effectiveness factors close to one in all cases, which means the agglomerate-level transport has only a minor effect when considering the present PGM-free catalyst activity. In addition, the 50 wt% case does not feature a Nafion film in order to be consistent with our prior nano-CT imaging of the Nafion domain.¹⁸ Thus, the origin of the losses is related to the through-plane transport across the cathode thickness rather than local microstructure. Figure 5a shows the liquid water saturation for the 50 wt% Nafion loading electrode at 0.4 V. The hydrophilic CCL is completely flooded at high current densities and reached the upper limit of 0.9, which is reduced from 1 by the residual gas saturation of 0.1. Figure 5b shows the corresponding O₂ distribution that shows a sharp decline in oxygen concentration due to the low effective diffusivity of the flooded pores. Because of the very low oxygen concentration near the membrane, the ORR is inhibited and the region of Peak ORR current density is found approximately one third of the way into the cathode away from the membrane and underneath the channel (see Figure 5c). Since the reaction is primarily occurring 20–40 μm away from the membrane, protons must be conducted over a longer distance, resulting in large ohmic potential drops. Figure 5d shows, the potential drop from the membrane to the location of peak ORR is approximately 150 mV at 0.4 V for that cathode.

Zoom In Zoom Out Reset image size

Comparing Figure 5 with Figure 4, the roughly 150–200 mV ohmic loss through the majority of the cathode in Figure 5 makes up a considerable portion of the mass transport overpotential seen in the polarization curves at 0.4 V in Figure 4. Thus, the oxygen transport resistance manifests itself as an ohmic loss that presents itself in polarization curves as a second ohmic slope, where the mass transport regime appears as steeper linear line without curvature.^18,37 We previously identified this phenomena in experiments with thick Pt/C cathodes by measuring the potential drop through the cathode using microstructured electrode scaffold (MES) diagnostics.³⁷

Considering the electrodes of different Nafion loading, we identify that the differences between the polarization curves in the mass transport regime is due to differences in the flooding of the macropores. First, the 50 and 60 wt% cases have quite similar porosities and 60% wt has much higher effective proton conductivity, but 60% wt presents notably larger losses. Thus, those parameters are not the origin of the difference. Instead, the difference relates to the current at which fully-flooded operation will occur. Nano-CT imaging showed the 60% wt electrode a much lower proportion of larger macropores as can be seen in the pore size distributions (Fig. S1) due to the thick Nafion films surrounding the agglomerates in the 60% wt case. Although the porosity of the macropore domain is similar for 50 and 60 wt%, the 60 wt% porosity resides more with the smaller pores because of the high Nafion loading. Thus, the smaller macropores of the 60% wt case presents early flooding onset and larger losses.

All the further analysis is done for the 50 wt% Nafion loading cathode because of it's good performance in both the activation, ohmic, and mass transport regimes of the polarization curve. Figure 6 shows the voltage gain in the 50 wt% Nafion loading cathode due to the resistance in different domains as a function of current density. We quantify voltage gain by individually setting the various resistance parameters to zero and subtracting the new voltage by the base case voltage at each current density. The effect of liquid water is studied by computationally setting the saturation to zero. The diffusion coefficients in different domains were multiplied by a large factor to make the transport resistance go to zero. Eliminating liquid presents the highest voltage gain, while the oxygen diffusion and proton conductivity in the CL are the next largest voltage gains and are close in magnitude. These voltage gains are similar because setting the O₂ diffusion and proton conductivity offsets the liquid water effect by letting oxygen diffuse through the flooded electrode or by eliminating the potential drop associated with conductivity protons across the cathode. Increasing the oxygen diffusion in the GDL and MPL can also notably reduce the losses at the mass transfer region by up to 250 mV. Diffusion in the Nafion film surrounding the agglomerate and diffusion in the micropores in the agglomerate shows minimal or no effect on the performance. Once the flooding is mitigated, the other resistances would become significant.

Zoom In Zoom Out Reset image size

Different morphological and operating parameters were systematically varied for the 50 wt% Nafion loading electrode under realistic conditions to observe the parameter sensitivity. Figure 7 show the power density's sensitivity to these parameters. Similar to the voltage gain result, the wettability of the CCL has the highest effect on power density. At 0.7 V, increasing the operating pressure and Nafion conductivity increases the power density by 50 mW cm⁻². However, at the voltage corresponding to maximum power, decreasing the thickness increases the power density by approximately 100 mW cm⁻². Other morphological parameters like volume fraction and tortuosity of pore and Nafion domain show an increase of less than 50 mW cm⁻². Oxygen diffusivity in the Nafion has a minimal improvement in power density. The following sections discuss a few parameters in detail.

Zoom In Zoom Out Reset image size

Figure 8 shows the effect of wettability on the polarization curve for the 50 wt% Nafion loading cathode. The base case is the multiphase simulation with hydrophilic wettability. The effect of water saturation was studied by computational setting the saturation to zero in the single-phase case. This case showed significant improvement in the performance for all voltages below 0.8 V. In the mass transfer region, the current density is nearly doubled. Maintaining zero saturation is a practically impossible at these temperature with fully humidified gases and thick electrode. However, fabricating the electrode with increased hydrophobicity is possible.^38,39 The hydrophobic cathode case is shown to have improved mass transport by reduced flooding. In the model, wettability is changed by changing the contact angle in the Young's-Laplace equation to the hydrophobic region (180°). The macropore is changed to hydrophobic while the micropore still is hydrophilic. This is a realistic assumption because the hydrophobic domain has a lower chance of infiltration into the micropore. The improved performance is observed when the through plane transport losses in the macropores becomes dominant, which is below 0.7 V.

Zoom In Zoom Out Reset image size

Effect of the catalyst loading is modeled by varying the cathode catalyst layer thickness. Figure 9 presents the polarization curve for varying thickness/catalyst loadings. While reducing the thickness of the electrode, morphological parameters were constant including the active site density and activity. Reducing the thickness of the CL reduces the current in the activation region, but also reduces the mass transport losses yielding higher limiting currents. If the thickness is reduced by decreasing the porosity, performance in the activation region would remain similar to the base case but a reduced current in the mass transfer region can be expected due to increased oxygen transport resistance with smaller macro-pores.

Zoom In Zoom Out Reset image size

In this section, effect of Nafion loading is studied. The Nafion loading and distribution is informed by the prior nano-CT observation of ionomer distribution.¹⁸ However, in all cases here the Nafion distribution is modeled assuming uniform distribution in the micro- and macro-pore unlike the 35 wt% Nafion loading sample that had feature large Nafion agglomerates, like due to the solid to solvents ratio used in making the electrode ink. The Nafion distribution model assumptions are:

• (1)  
  For loadings below 50 wt%, only the micropores inside the agglomerate are infiltrated with Nafion, and there is no Nafion film formation on the agglomerate.
• (2)  
  Films form around the agglomerates for Nafion loadings greater than 50 wt%.
• (3)  
  Porosities in the micropore and macropore are corrected accordingly to the increase in Nafion volume fraction as well as film thickness.

The Nafion distribution is modeled as follows. The Nafion volume fraction, ε_Nafion, in the agglomerate is estimated according to the expressions:

where ρ_Nafion, L_CCL, and y_Nafion are the density of Nafion and thickness of the cathode catalyst layer respectively. The mass of Nafion, m_Nafion, is calculated as

The film thickness, δ_film, is similarly calculated as

where ε_solid is solid volume fraction and SA is the surface area of the agglomerate obtained from nano-CT imaging.

Figure 10a shows the Nafion volume fraction and the film thickness as a function of Nafion loading. Figure 10b shows the resulting polarization curves for different Nafion loadings. The current density at the activation region increases with increasing Nafion loading due to reduced through-plane ohmic losses since proton transport in the agglomerates and local proton activity is not modeled and are thus not reflected in the activation region current. Observing the mass transfer region, the current density is highest for the Nafion loading of 50 wt%. The 60 wt% Nafion loading shows a significant ohmic-like voltage loss at low potential due to the reduced size and porosity of the macro-pores, which increases the flooding rate.

Zoom In Zoom Out Reset image size

This section discusses the effect of using the distribution of ASD obtained from nano-CT imaging versus a single representative diameter. Figure 11a shows the ASD for the 50 wt% Nafion loading cathode as well as the mean and mode of the distribution. The mean and mode diameter of the distribution are 0.9 μm and 2 μm respectively. In addition, the apparent mean diameter from SEM is roughly 6 μm. The smaller size of the agglomerates from the nano-CT images relative to the SEM images is because the nano-CT uses a volume to surface area ratio to get diameter, which yields smaller diameters for these irregularly shaped particles. The distribution is implemented in the agglomerate model using Eq. 11, whereas the single diameter case was modeled with a single diameter having f_i = 1. Figure 11b shows the resulting polarization curve comparison for the distribution and single diameter. For the current activity, the differences in the polarization curves for the distribution versus the mode case is very low. In the mass transfer region, the performance is notably lower for the SEM's apparent diameter and the mean diameter. Thus, these diameters are limit of where local, agglomerate-level losses are negligible. As ca be seen in Figure 11b, the agglomerate level losses become significant at the 6 μm SEM agglomerate diameter.

Zoom In Zoom Out Reset image size

In this section, we discuss the improvement needed in the density or activity of the active sites to achieve higher power density. The active site's activity or activity site density was varied by changing the R_sJ_s parameter using a multiplatication factor, which we refer to as the normalized site density. The rational for this approach and presentation of the results is based on assuming the active sites remain similar in terms of activity, but the spatial active site density is varied by tuning the catalyst synthesis. The site density can also be increased by increasing the utilization of the metal site. Sahraie et al. estimated that only 20% of the iron metal centers are ORR active in acidic medium.¹² The active iron sites correspond to 1 wt% of iron content present in the catalyst. Figure 12 plots contours of the power density as function of the normalized site density and the cathode thickness. An additional line is ploted showing the thickness yielding the highest power density at a given normalized site density. The optimal thickness of the cathode decreases with increases in the activity or site density. Figure 12a shows that to achieve a proposed interim power density target [Insert ref to strategic analysis AMR presentation] of 475 mW cm⁻² requires a 40 times increase in either the active site activity or the activity per site. The optimal thickness for 40 times increase in volumetric activity at 0.7 V is 47.5 μm while at 0.6 V it is less than 20 μm. For operating at 0.5 V and the voltage corresponding to maximum power, the optimal thickness, shown in Figures 12c and 12d, is below 20 μm.

Zoom In Zoom Out Reset image size

This section discusses the performance and electrode parameters sensitivities of the PGM-free cathodes with increased activity or site density for a 50 wt% Nafion loading cathode. Here, the R_sJ_s parameter is increased by a factor of 10. The morphological parameters were as before except the thickness of the cathode is 35 μm based on the optimal thickness estimation at 0.6 V from Figure 11b. Figure 13a shows the comparison between the electrode with the current activity fitted from experiments and the increased activity. Figure 13b shows the voltage gain due to different resistances for the increased activity catalyst. We quantify voltage gains with the improved catalyst similarly to the method in Figure 6 by individually setting the various resistance parameters to zero and subtracting the new voltage by the base case voltage with the improved activity catalyst at each current density. The voltage gains associated with flooding and through-plane transport still present the highest values. However, the voltage gains associated with oxygen diffusion in the Nafion and agglomerate are now more observable.

Zoom In Zoom Out Reset image size

Figure 14 shows the parametric study of the various morphological parameters for the improved activity catalyst. Even in the improved activity catalyst, the wettability shows the highest performance sensitivity. Unlike the present activity catalyst, other morphological and operating parameters show significant improvement at 0.7 V. Increases in the macropore porosity reduces the performance as the thickness of the CL increases.

Zoom In Zoom Out Reset image size

Figure 15 shows the model prediction for using the distribution and the single representative diameter in the agglomerate model for the improved activity catalyst. The single representative diameter of mode and mean case from the nano-CT distribution have higher mass transport losses than the distribution. The average from the SEM case has lowest performance because of the transport losses due to the larger diameter. Comparing to prior analysis with the present activity, it is apparent that agglomerate size control is important with future increases in active site activity or spatial density. While reducing the size of the diameters may be important, the mode suggests that macropore porosity should be maintained similar to the base case as decreasing the pore volume fraction reduces the through-plane oxygen diffusivity.

Zoom In Zoom Out Reset image size

List of Symbols

Greek

Subscripts and Superscripts