Research of Proton Exchange Membrane Fuel Cell Modeling on Concentration Polarization under Variable-Temperature Operating Conditions

Abstract

In this study, a concentration overvoltage model that focuses on describing variable-temperature operating condition properties for PEMFCs is established. Sensitivity analysis and a quantification study of oxygen transport resistance are carried out based on the oxygen transport resistance model and measurement data. By analyzing the influence of temperature on cathode oxygen transport resistance, the key structural parameters of the cathode oxygen transport resistance models are estimated, and the parameter modification method of fuel cell limiting current density under variable temperatures is proposed. Based on the polarization curve test experiments under variable-temperature conditions, it is demonstrated that the newly developed concentration overvoltage model reduces the relative error of simulation for a low Pt loading fuel cell in the high current region by 2.97% and 10.06% at 60 °C and 80 °C, respectively. The newly established concentration overvoltage model of a PEMFC solves the problem that the parameter of limiting current density is set without considering the influence of fuel cell temperature fluctuation, which leads to the poor simulation accuracy of the concentration overvoltage model in the high current region.

1. Introduction

Under the common vision of human efforts to achieve carbon neutrality, proton exchange membrane fuel cells (PEMFCs) possess the potential for large-scale application in future power generation and energy storage systems due to their zero-emission and high energy efficiency features [1]. As one of the alternatives to traditional internal combustion engine vehicles, with long endurance, fast refueling capability and adaptability to ambient temperature, PEMFC vehicles have ushered in a rapid development stage [2]. However, a series of technical bottlenecks, such as high cost of catalyst, high complexity of power system integration and short PEMFC lifetime, still put PEMFC vehicles at a disadvantage in competition with traditional internal combustion engine vehicles [3]. At the present stage, applying PEMFCs as an auxiliary power source and developing a compact PEMFC system to reduce the vehicle cost has become a favored technical solution [4] used in, for example, the BAIC Motor Corporation Limited CF53FCV extended range PEMFC vehicle, Mercedes-Benz GLC F-Cell and Audi A7 Sportback H-tron Quattro. Coping with high load conditions, a compact PEMFC system needs to run in a high-power (high-current) operating state. The voltage output of a PEMFC in the high-current range is dominated by concentration polarization [5]. Although the PEMFC thermal management system restricts the stack within a certain range of operating temperatures, the temperature fluctuation of the PEMFC stack due to cold starting or instantaneous power load change may still influence the voltage output [6]. At present, there is no unified analytical scheme of concentration overvoltage model suitable for characterizing PEMFCs under variable operating conditions (such as stack temperature, inlet pressure, humidity, etc.) [7].

In 1986, a concentration overvoltage model for describing mass transport in the high-current region of PEMFCs based on temperature, partial pressure of reaction gas and morphological characteristic assumptions of proton exchange membranes was proposed in Physical Chemistry by an Oxford University scholar named Atkins. Owing to its theoretical significance, this model has been widely used to describe the concentration overvoltage of PEMFCs [8]. In 1989, a simplified semi-empirical model for describing the PEMFC volt-ampere characteristics was first proposed by Srinivasan from Los Alamos National Laboratory in the United States [9]. The model of open-circuit voltage was obtained by Gibbs free energy theory, the activated overvoltage was expressed by the Tafel equation summarized from experimental rules and the ohmic overvoltage was represented by the slope of the linear region of a polarization curve. The model could not give the analytical mechanism interpretation of the concentration overvoltage or ohmic overvoltage. In 1995, Kim introduced a concentration overvoltage term described by an exponential function on the basis of the Srinivasan model to describe the overvoltage caused by mass transport in high-current regions [10]. The same year, Amphlett proposed a semi-empirical voltage model considering multiple PEMFC operating parameters, which accurately expressed the PEMFC-activated overvoltage through the parameterization formula [11,12]. In 2000, Mann and Amphlett gave an analytical expression of ohm overvoltage in PEMFCs by analyzing the topology of a proton exchange membrane [13]. This model is recognized as the semi-empirical model that can best reflect the polarization characteristics of PEMFCs under different operating conditions at present. The model was widely used in the study of stack parameter identification, since there is no evidence of a relationship between Kim’s exponential function parameters and experimental mass transport variables for describing concentration overvoltage. In 1999, Squadrito proposed an improved empirical formula of concentration overvoltage based on the Atkins model [14]. The model requires fewer identification parameters; however, different stack operating conditions have a remarkable impact on the model precision in fitting concentration overvoltage. In 2004, Kulikovsky started with the Srinivasan and Kim models and took the influence of oxygen supply into consideration by applying stoichiometry of activation overvoltage and concentration overvoltage of PEMFCs [15]. In 2009, Jiang Hui Yin established a PEMFC model that considered all areas of PEMFCs and could characterize the influence of membrane water content variation on PEMFC performance [16]. In 2018, Chakraborty analyzed the nonequilibrium state error of the Nernst voltage for PEMFCs [17].

Low platinum (Pt) loading and high performance requirements have impelled the researchers to explore the influence mechanism of the membrane electrode assembly (MEA) structure and material composition on overpotential of PEMFCs. [18]. The essence of concentration polarization is the phenomenon that cathode reaction gas encounters mass transport resistance in the process of passing through a gas diffusion layer, microporous layer and catalyst layer, resulting in an insufficient catalytic reaction of reaction gas [19]. Model simulation and experimental studies on mass transport resistance for the gas diffusion layer, microporous layer and catalyst layer of a PEMFC become a promising research direction to explore the possibility of improving the catalytic efficiency, reducing mass transport resistance of reactant gas and improving the performance of PEMFCs in the high-current region [20].

Nonoyama quantified the oxygen transport resistance in the PEMFC structural layer under different operating conditions (temperature, pressure, humidity) by carrying out an oxygen transport resistance measurement experiment. The simplified models of each oxygen transport subresistances were presented [21]. Hydrophobic particle materials were added to the cathode catalyst layer by Wan Zhaohui to analyze the influence on the oxygen transport resistance of the catalyst layer. The transport resistances were measured by limiting current testing [22]. In view of the difficulty in determining the key structural parameters of the catalyst layer, Epting used typical PEMFC electrodes with nano-CT imaging to measure the diameter distribution and analyzed the effect of agglomeration size distribution on concentration overvoltage of PEMFCs. The data source of concentration overvoltage was not mentioned in the literature [23]. Li Shan considered the structure of the catalyst layer and the effect of liquid water on a spherical agglomerate model to determine the electrochemical kinetics of the cathode catalyst layer. The effects of the agglomerate radius $r_{agg}$ and the volume fraction of the ionomer $L_{i,agg}$ on PEMFC local mass transport characteristics were studied by using the established agglomerate model [24].

The limiting current is defined as the current at which oxygen diffuses through the GDL and CL layers to the active site of the catalyst and is consumed instantaneously, and the voltage drops sharply to zero [21]. This physical quantity can indirectly represent the oxygen transport resistance in the gas diffusion layer and catalyst. In order to obtain a relatively stable measurement value of the fuel cell limiting current density, it is necessary to control the oxygen supply of the cathode to a certain flow, that is, to control the oxygen supply with a constant flow rate rather than a fixed stoichiometric ratio.

Many scholars have carried out relevant studies by taking advantage of the fact that there is zero oxygen concentration on the surface of catalysts under limiting current conditions. For example, Reshetenko proposed an experimental method to determine the mass transport coefficient of oxygen in the gas phase and ionomer. The method was based on the mathematical model of limiting current density distribution and diluted gases of different molecular weights, and gave the conclusion that the total transport coefficient had a linear relationship with the molecular weight of the diluent gas [25]. Spingle and Schuler [26] developed a hydrogen pump limiting current device to analyze the mass transport loss and gas diffusion for the catalyst layer and microporous layer. Attuluri’s experimental data on the PEMFC polarization curve at variable temperatures show that the limiting current density is positively correlated with the PEMFC temperature in general; however, this paper does not give the mechanism explanation [27].

High-performance PEMFCs for automobiles are important platforms for the application of hydrogen energy, and the breakthroughs in PEMFC technology will strongly support the large-scale application of hydrogen energy as a sustainable energy source. Based on analyzing mass transport rules of PEMFCs under the condition of variable temperature, it will be a potential modeling direction to combine the concentration overvoltage model of a PEMFC with a mass transport resistance mechanism. This study aims to explore the temperature dependence of the cathode mass transport resistance of PEMFCs in the high-current region and to optimize the performance of the model fitting for concentration overvoltage under variable temperature conditions. The novelty of this study is that we focused on modeling analysis of the ionomer/water membrane permeation resistance, which dominates the oxygen transport resistance under low Pt loading conditions of the catalyst, and reconstructed the ionomer/water membrane permeation resistance model based on the key structural parameters of the catalyst layer at the microscopic level. Meanwhile, the experimental methods for the separation of oxygen transport resistance in PEMFCs are systematically summarized for the first time. The importance of this study is to reveal the law of the PEMFC oxygen transport resistance in the high current density interval as affected by the difference in operating temperature, which provides a theoretical analysis basis for the development of high-performance PEMFCs and simulation research of cold-start operating conditions for a PEMFC stack.

2. Model Description

2.1. Concentration Overvoltage Model

At present, the most widely used concentration overvoltage electrochemical model of a PEMFC is a semi-empirical model proposed by Atkins [8], and the subsequent models proposed by Amphlett, Squadrito and others are developed based on this model. The concentration overvoltage of PEMFCs can be calculated by using the Nernst equation during concentration variation of the reaction gas at the electrode site [7].

The reversible voltage of PEMFCs can be expressed according to the Nernst equation.

$$E_{Nernst}=E^0+\frac{RT}{2F}\ln (\frac{P_{H_2}·P_{O_2}^{\frac{1}{2}}}{P_{H_{2}O}})$$

(1)

The concentration overvoltage derivation process is detailed in Section S2 of the Supplementary Materials.

$$V_{con}=\frac{RT}{F}\left(\frac{1}{4}+\frac{\mathsf{\Upsilon }}{\mathsf{\alpha }}\right)\ln (1-\frac{i}{i_{\lim }})$$

(2)

2.2. Limiting Current Density

As a highly sensitive parameter, $i_{\lim }$ can significantly affect the simulation result of the concentration overvoltage model in the high-current region. Correa [28] conducted a sensitivity analysis on the $i_{\lim }$ parameter based on BCS Technology stack experiment data, which also confirmed the conclusion above. Attuluri and Ozen, through experiments, proved the conclusion that appropriately increasing the operating temperature can significantly increase the $i_{\lim }$ of PEMFCs and reduce concentration overvoltage [29]. In the past, according to the actual operation condition of PEMFCs, researchers generally set a fixed value or range estimate value for the $i_{\lim }$, ignoring the influence of PEMFC temperature fluctuation on $i_{\lim }$ [30,31,32,33]. Therefore, in dealing with the situation of temperature fluctuation, simulation accuracy will decline in the high-current region of the classic concentration overvoltage model. The problem can be proved in Riascos’s simulation polarization curve image, which indicates convergence tendency on PEMFC concentration overvoltage (high-current region) under different temperature conditions. This trend is obviously inconsistent with the actual PEMFC polarization curve [34]. Therefore, it is necessary to introduce temperature as an influencing factor to modify the $i_{\lim }$ (Equation (S5), Supplementary Materials) parameter in the classical concentration overvoltage model.

Atkins characterized the relationship between the molar conductivity of the electrolyte and the ionic diffusion coefficient by combining the Nernst–Einstein Equation (3), where the molar conductivity $\mu =2uF$, and $u$ is the ionic mobility.

$$D=\frac{RT\mu }{z^2{F}^2}$$

(3)

Substituting Equation (3) into (Equation (S5), Supplementary Materials), the limiting current formula described by molar conductivity is obtained. ${\delta }_{Nernst}$ is the thickness of Nernst diffusion layer.

$$i_{\lim }=\frac{cRT\mu }{zF{\delta }_{Nernst}}$$

(4)

Due to the complexity of the multilayer and multimaterial structure of the PEMFC MEA, it is difficult to use some fuzzy, immeasurable or undefinable boundary parameters such as ionic conductivity or Nernst diffusion layer thickness to characterize it.

PEMFC concentration overvoltage is dominated by the oxygen transport resistance of the cathode. The mass transport resistances with different cathode materials, structures and components are significantly different [35,36,37,38,39]. Oxygen transport resistance can be expressed by Formula (5) [40].

$$R_{total}=\frac{\Delta C}{J}$$

(5)

According to equation (Equation (S5), Supplementary Materials), $i_{\lim }$ is negatively correlated with oxygen transport resistance. Nonoyama et al. found in their study that the sensitivity of each oxygen transport subresistance was different under the influence of PEMFC temperature [21]. As a result, coupling subresistances of cathode mass transport into the analytic expression $i_{\lim }$, with oxygen an transport resistance measurement experiment to confirm the temperature effect of mass transport, putting analytical expression $i_{\lim }$ into the classic PEMFC concentration overvoltage model, will be a meaningful research direction.

2.3. Mass Transport Resistance

At present, it is known that PEMFC cathode oxygen transport resistance consists of three kinds of mass transport subresistances, which are molecular diffusion resistance $R_{MD}$, Knudsen diffusion resistance $R_{CL,Knud}$ and ionomer/water membrane permeation resistance $R_{CL,ion}$. Due to the structural difference of each layer of the MEA, the location of the mass transport resistance is determined by the structural characteristics [41]. Molecular diffusion resistance mainly exists in the gas diffusion layer (support layer, microporous layer) of the MEA, in which the microporous layer is dominant, and the catalyst layer takes up a very small proportion. Knudsen diffusion resistance mainly exists in the catalyst layer, and a small part exists in the microporous layer. The permeation resistance of the ionomer mainly exists in the catalyst layer. Some studies have suggested that the outer boundary of the agglomerates that produce ionomer permeability resistance is surrounded by ionomer and water film (Figure 1). In this paper, the ionomer/water film on the surface of the agglomerates is considered as the same medium. Then the permeation characteristics of oxygen through the ionomer/water film can be effectively described.

In terms of $R_{MD}$ calculation, Wilke and Fuller–Schetler–Giddings equations are applied to calculate the molecular diffusion coefficient, where Equation (6) is used to calculate the binary diffusion coefficient of oxygen in a single gas component.

$$D_{A:B}=\frac{0.001T^{1.75}}{P(v_A^{1/3}+v_B^{1/3})}\sqrt{\frac{1}{M_A}+\frac{1}{M_B}}$$

(6)

As cathode gas usually contains nitrogen, vapor and other components, the multicomponent gas diffusion coefficient for oxygen must be calculated. According to Equation (6), binary diffusion coefficients of oxygen in nitrogen and water vapor are calculated respectively. Then, taking the coefficients into the multicomponent gas diffusion coefficient calculation formula, the calculation formula for the oxygen diffusion coefficient in the air $D_{O_2:Air}$, are (7).

$$D_{O_2:Air}={\left({\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne A\end{array}}^{Nc}\frac{x_j}{D_{A:j}}}\right)}^{-1}={\left(\frac{W_{N_2}·M_{Ca,in}}{W_{Ca,in}·M_{N_2}·D_{O_2:N_2}}+\frac{W_{Va}·M_{Ca,in}}{W_{Ca,in}·M_{Va}·D_{O_2:V_a}}\right)}^{-1}$$

(7)

In Equation (7), nitrogen and vapor mass flow derivation process details can be seen in Supplementary Materials (S4. Humidification and mass flow rate conversion modules). Since the gas diffusion layer structure will affect molecular diffusion resistance $R_{MD}$, in the actual analytic expression fitting process on $R_{MD}$, the coefficient should be modified with the actual $R_{MD}$ measurement result of the PEMFC. Then the effective molecular diffusion coefficient $D_{DM,O_2}^{eff}$ is introduced. The $R_{MD}$ is the equation dividing the thickness of the diffusion layer $h_{DM}$ by the effective molecular diffusion coefficient $D_{DM,O_2}^{eff}$, which is correlated with the molecular diffusion coefficient $D_{O_2:Air}$.

$$R_{MD}=\frac{h_{DM}}{D_{DM,O_2}^{eff}}\propto \frac{1}{D_{O_2:Air}}$$

(8)

The diffusion coefficient of Knudsen diffusion $D_{Knud,O_2}$ can be expressed as (15).

$$D_{Knud,O_2}=\frac{2r_{Knud}}{3}\sqrt{\frac{8RT}{\pi M_A}}$$

(9)

Knudsen diffusion mainly exists in the catalyst layer. Considering the heterogeneity and complexity of the catalyst, the catalyst effective thickness $h_{CL}^{eff}$ and effective diffusion coefficient $D_{CL,O_2}^{eff}$ are introduced to modify Knudsen diffusion resistance. The expression of Knudsen diffusion resistance $R_{CL,Knud}$ for the catalyst layer is (10).

$$R_{CL,Knud}=\frac{h_{CL}^{eff}}{D_{CL,O_2}^{eff}}\propto \frac{h_{CL}^{eff}}{D_{Knud,O_2}}$$

(10)

At present, there are still controversies about the catalyst layer ionomer binder and the Pt/C particle morphology. The interface model, homogeneous model and aggregate model are the most commonly used models to describe the catalyst layer structure [42]. This paper intends to introduce the agglomerate model as the architectural basis for the ionomer permeation resistance modeling. The permeation resistance of the ionomer $R_{CL,ion}$ is expressed as Equation (11).

$$R_{CL,ion}=\frac{{\delta }_{ion}^{eff}}{A_{ion}^{eff}}\frac{H_{ion,O_2}}{D_{ion,O_2}RT}\propto \frac{{\delta }_{ion}^{eff}}{A_{ion}^{eff}}\frac{1}{{\psi }_{ion,O_2}T}$$

(11)

${\psi }_{ion,O_2}$ is the permeation coefficient; this parameter is the ratio of the diffusion coefficient $D_{ion,O_2}$ to Henry’s constant $H_{ion,O_2}$. Henry’s constant was introduced into the analytical Formula (11) to describe the oxygen permeation phenomenon at the interface of the ionomer film. The variation in oxygen diffusion concentration between the gas phase and the ionomer phase in the catalyst layer was effectively characterized by $H_{ion,O_2}$.

2.4. Parameter Sensitivity of Oxygen Transport Resistances under Operating Conditions

The mathematical model of oxygen transport resistances in the section above indicates significant diversity in terms of temperature sensitivity (

Table 1. Sensitivity of transport resistances toward operating temperature.

Transport Resistance	Temperature Sensitivity
$R_{MD}$	$T^{-1.75}$
$R_{CL,Knud}$	$T^{-0.5}$
$R_{CL,ion}$	$T^{-1}{\left(e^{17200/RT}\right)}^{-1}$

); $R_{CL,ion}$ is significantly affected by temperature. For the factor of pressure, $R_{MD}$ is positively correlated with the change of cathode pressure, and $R_{CL,Knud}$ and $R_{CL,ion}$ are not affected by pressure variation. $R_{CL,ion}$ was significantly affected by the condition of variable humidity. The type of balance gas also affects the oxygen transport resistance, for example, the diffusion coefficient of oxygen in different balance gases (such as $N_2$ and $He$), resulting in different $R_{MD}$. In addition, some key structural parameters of the MEA, such as the thickness of the gas diffusion layer, the thickness of the catalyst layer $h_{CL}$, the thickness of ionomer film ${\delta }_{ion}^{eff}$ and the radius of the agglomerate $r_{agg}$, also have significant effects on oxygen transport resistance. According to the analysis above, under the operating condition of variable temperature (60–90 °C), the $R_{CL,ion}$ changes drastically compared with $R_{CL,Knud}$ and $R_{MD}$.

2.5. Method for Measuring Oxygen Transport Resistance

Due to the complexity and difficulty of observing the microstructure characteristics of the MEA, the structural parameters are difficult to obtain in general. In order to accurately fit each oxygen transport resistance with a mathematical model, it is necessary to appropriately modify the key structural parameters of oxygen transport resistance. Before modification, oxygen transport subresistance values should be measured. Some relevant measurement methods to separate and measure the oxygen transport subresistances under specific operating conditions have been proposed. The core idea to quantify oxygen transport subresistances is taking advantage of the sensitivity difference for cathode pressure, temperature, humidity and the type of balance gas. The experimental measurements of oxygen transport resistance methods are summarized and evaluated in (Supplementary Materials Table S1).

This paper summarizes the typical experimental methods of oxygen transport resistance measurement from the aspects of experimental variables, theory and implementation difficulties, for the first time. The four experimental methods can be used cooperatively according to existing experimental equipment and experimental conditions in order that oxygen transport resistance can be separated efficiently and accurately. Finally, the key structural parameters of each subresistance model are modified by using the measurement results of oxygen transport resistances under operating conditions. The key parameters to be corrected include $r_{agg}$, $r_{Knud}$, $h_{DM}$, etc.

2.6. Structural Parameter Analysis and Optimization for Ionomer Permeation Resistance

According to the temperature sensitivity analysis above, the penetration resistance of the ionomer can significantly affect the PEMFC output voltage in the high-current region in terms of different temperatures. Due to the complexity and irregularity of the aggregate structure, effective thickness ${\delta }_{ion}^{eff}$ and effective area $A_{ion}^{eff}$ of the ionomer film can only be estimated and measured in a certain accuracy range by imaging tools such as a scanning electron microscope [43,44,45]. Nonoyama et al. estimated the two parameters above at same time based on the experimental measurement data $R_{CL,ion}$ (setting the value ranges of ${\delta }_{ion}^{eff}$ and $A_{ion}^{eff}$ respectively, selecting the appropriate ${\delta }_{ion}^{eff}$ and $A_{ion}^{eff}$ values to be used in conjunction with the formula of $R_{CL,ion}$, so as to fit the measurement results of $R_{CL,ion}$). However, the direct simultaneous estimation of ${\delta }_{ion}^{eff}$ and $A_{ion}^{eff}$ lacks the estimation basis, and the estimated value is likely to be far beyond the actual size range of the topological structure for the ionomer film [46].

In this paper, catalyst component parameters are used to effectively characterize the effective thickness ${\delta }_{ion}^{eff}$ and effective area $A_{ion}^{eff}$ of ionomer film. The ionomer film thickness is expressed by the following relationship (18). For the detailed derivation, please refer to Li Shan’s literature [24]. In that paper, catalyst component parameters are used to effectively characterize ${\delta }_{ion}^{eff}$ and $A_{ion}^{eff}$. The ionomer film thickness is expressed by following relationship.

$${\delta }_{ion}^{eff}=r_{agg}\left[\sqrt[3]{\frac{L_i\left(1-L_{i,agg}\right)}{L_{Pt/C}}-L_{i,agg}+1}-1\right]$$

(12)

The idea of deducing the effective area of ionomer film per unit of MEA area $A_{ion}^{eff}$ is to share the effective area of the ionomer film to the surface of each agglomerate, which is the sum of the agglomerates’ effective area $S_{agg,i}$ per unit catalyst volume. $A_{ion}^{eff}$ can be represented by Formula (13).

$$A_{ion}^{eff}=\frac{S_{agg,i}·N·{\epsilon }_{CL}·A_{CL}·h_{CL}}{A_{CL}}=\frac{3L_{Pt/C}·{\epsilon }_{CL}·h_{CL}}{{r_{agg}}^3\left(1-L_{i,agg}\right)}{\left(r_{agg}+{\delta }_{ion}^{eff}\right)}^2$$

(13)

where $S_{agg,i}$ is the area of a single aggregate, $N$ is the number of aggregates, $A_{CL}·h_{CL}$ is the volume of unit catalyst layer and ${\epsilon }_{CL}$ is the pore volume fraction of catalyst layer.

According to Formula (13), by introducing catalyst component parameters, ${\delta }_{ion}^{eff}/A_{ion}^{eff}$ can be expressed as (14).

$${\delta }_{ion}^{eff}/A_{ion}^{eff}=\frac{{r_{agg}}^2\left[\sqrt[3]{L_i\left(1-L_{i,agg}\right)/L_{Pt/C}-L_{i,agg}+1}-1\right]}{3L_{Pt/C}·{\epsilon }_{CL}·h_{CL}}$$

(14)

Substituting (14) into (11), then the newly constructed ionomer film resistance formula $R_{CL,ion}$ is obtained (15).

$$R_{CL,ion}=\frac{{r_{agg}}^2\left[\sqrt[3]{L_i\left(1-L_{i,agg}\right)/L_{Pt/C}-L_{i,agg}+1}-1\right]}{3L_{Pt/C}·{\epsilon }_{CL}·h_{CL}}·\frac{H_{ion,O_2}}{D_{ion,O_2}RT}$$

(15)

Compared with Nonoyama’s model, Equation (15) uses the catalyst layer component parameters (Pt/C particle volume fraction $L_{Pt/C}$, volume fraction of ionomer $L_i$ and other measurable parameters) to effectively connect the two isolated catalyst structural parameters ${\delta }_{ion}^{eff}$ and $A_{ion}^{eff}$, which is worthy of putting forward to analyze the $R_{CL,ion}$ detailed mechanism with more structural parameters.

3. Electrochemical Model of PEMFC

This paper is based on the electrochemical model of a PEMFC proposed by Amphlett, in which the Nernst reversible voltage is:

$$E_{Nernst}=1.229-0.85\times {10}^{-3}(T-298.15)+4.3085\times {10}^{-5}T⌈\ln \left(P_{H_2}\right)+0.5\ln \left(P_{O_2}\right)⌉$$

(16)

Based on the Tafel equation obtained from the electrode reaction kinetics, the activation overvoltage is:

$$V_{act}=\frac{RT}{{\mathsf{\alpha }}_{0}F}\ln (\frac{i}{i_0})$$

(17)

The ohmic overvoltage is:

$$V_{ohmic}=i\left(R_m+R_c\right)$$

(18)

$R_m$ is membrane equivalent resistance, which dominates in ohmic overvoltage, and $R_c$ is electronic impedance.

$$R_m=l·A_{CL}^{-1}·\left(\frac{181.6\left[1+0.03\left(\frac{i}{A}\right)+0.062{\left(\frac{T}{303}\right)}^2{\left(\frac{i}{A}\right)}^{2.5}\right]}{\left[\lambda -0.634-3\left(\frac{i}{A}\right)\right]\exp \left[4.18\left(\frac{T-303}{T}\right)\right]}\right)$$

(19)

From Equations (2), (8), (10) and (11), the newly established analytical formula of concentration overvoltage can be derived (20), in which the mathematical description related to oxygen transport resistance has not been expanded in detail. For detailed mathematical descriptions, please refer to (7), (9) and (15) above.

$$V_{con}=\frac{RT}{F}\left(\frac{1}{4}+\frac{\mathsf{\Upsilon }}{\mathsf{\alpha }}\right)\ln \left(1-i\frac{RT}{4FP{x}_{O_2}}\left(\frac{h_{DM}}{D_{DM,O_2}^{eff}}+\frac{h_{CL}^{eff}}{D_{DM,O_2}^{eff}}+\frac{{\delta }_{ion}^{eff}}{A_{ion}^{eff}}\frac{H_{ion,O_2}}{D_{ion,O_2}RT}\right)\right)$$

(20)

4. Result and Discussion

In this paper, the electrochemical mathematical one-dimensional model of a PEMFC is developed in the Matlab2021b environment; the model realizes oxygen transport resistance numerical simulation and PEMFC V-I characteristic curve simulation analysis. Nonoyama’s and Zhaohui Wan’s experimental operation conditions and MEA parameters were used in this paper (details in Supplementary Materials Tables S2 and S3). The experimental data of PEMFC cathode oxygen transport subresistances under the condition of different temperatures in Nonoyama’s and Wan’s literature are extracted and analyzed. The ultimate purpose is to apply the modified subresistances’ key structural parameters into the newly established PEMFC concentration overvoltage model, so as to analyze the simulation result in the high-current region.

4.1. Oxygen Transport Resistance Numerical Simulation

In Figure 2 and Figure 3, $R_{MD}$ is separated by switching the balance gas method while $R_{CL,Knud}$, $R_{CL,ion}$ are separated by variable temperature operating conditions in MEA1. In Figure 4, $R_{MD}$ is separated by intercept method, while $R_{CL,Knud}$, $R_{CL,ion}$ are separated by variable temperature operating conditions in MEA2. It can be seen from Figure 2 and Figure 3 that $R_{MD}$ is significantly different in balance gases ($N_2$, $He$), which reflects distinct diffusion difference in different balance gases. According to Tables S2 and S3 (Supplementary Materials) and comparatively analyzing Figure 3, Figure 4 shows that the Pt loading level can significantly affect the magnitude of $R_{CL,ion}$. At lower Pt loading conditions, $R_{CL,ion}$ could even dominate $R_{total}$, which indicates that Pt loading level is an important factor to describe the difficulty of oxygen molecules being catalyzed. Meanwhile, $R_{CL,ion}$ can be seen to be more significantly affected by temperature than $R_{MD}$ and $R_{CL,Knud}$ by observing its slope with temperature variation.

4.2. Experiment

The fuel cell MEAs were prepared according to the parameters in (Supplementary Materials Tables S2 and S3), and the fuel cell polarization curves were collected on an HTS-125 Hephas fuel cell test bench at variable working temperature. The effective area of the single fuel cell used for testing with carbon plate structure and serpentine flow channels is 2 cm².

The polarization curve test conditions are as follows: H₂/air flow rate, 500 cm³/min and 650 cm³/min for MEA1; 1700 cm³/min and 4600 cm³/min for MEA2; H₂/air humidity 100/70% RH; H₂/air back pressure 150/150 kPa (relative pressure). Fuel cell test temperatures are 60, 70 and 80 °C.

Although the problem of flooding in the high-current region should be considered, the literature experiment and our experiment adopted the measures of low oxygen concentration, high flow and small-reaction-area membrane electrode, which effectively avoided the problem of flooding for the single fuel cell [47].

4.3. The Effect of Variable-Temperature Operating Conditions on PEMFC Concentration Polarization

Both the simulation results of concentration overvoltage models based on $i_{\lim }$ parameter correction and the classical concentration overvoltage model in the high-current region have been compared with experiment data in Figure 5 and Figure 6. It can be seen from Figure 5 and Figure 6 that when $i_{\lim }$ is set at a constant value without considering the influence of temperature, the polarization curve of the classical PEMFC voltage model in the temperature range of ±10 °C tends to converge to the point of fixed current density in the high-current region, which is inconsistent with the experimental data of the PEMFC polarization curve under the actual temperature change condition. The $i_{\lim }$ should be positively correlated with the working temperature of the PEMFC between 60 and 80 °C. Through relative error analysis, it is found that the relative errors of the concentration overvoltage model based on parameter modification are reduced by 1.94% and 4.58% when fitting MEA2 at 60 °C and 80 °C. The relative errors of the fitting effect for MEA1 decrease by 2.97% and 10.06% for the modified model at 60 °C and 80 °C. There is a more reasonable result for the parameter modification concentration overvoltage model on describing the variation of $i_{\lim }$ affected by temperature fluctuation in the high-current region.

5. Conclusions

Reducing catalyst Pt loading is the PEMFC development direction. The permeation resistance of ionomer film will increase sharply with the decrease in Pt loading, which will significantly improve the total oxygen transport resistance of the PEMFC and even make the permeation resistance of the ionomer film become the dominant mass transport resistance in the high-current region. The limiting current density is an important physical quantity to characterize the total oxygen transport resistance of PEMFCs. The high sensitivity of ionomer film permeation resistance to temperature makes the factor of temperature fluctuation cannot be neglected when setting the PEMFC voltage model parameter of limiting current density. In this paper, the influence laws of temperature on the oxygen transport subresistances of cathodes are analyzed in detail. Experimental methods of oxygen transport resistance measurement are summarized comprehensively. The analytic expression of ionomer permeation resistance was optimized by reconstructing the representation of effective ionomer film area. Combined with the cathode structure parameters determined by oxygen transport resistance measurement data, a modification analytic expression for limiting current density parameter of a PEMFC under variable temperature has been proposed. Finally, a concentration overvoltage model of a PEMFC has been established under the condition of temperature variation. The newly established concentration overvoltage model of a PEMFC solves the problem of the limiting current density parameter setting with the traditional method of fixed or range estimation value, which does not consider the influence of PEMFC temperature fluctuation, leading to the poor simulation accuracy of concentration overvoltage in the high-current region. This paper provides theoretical support for the design and modeling of high-performance and high-power density PEMFCs.