Sensitivity Based Order Reduction of a Chemical Membrane Degradation Model for Low-Temperature Proton Exchange Membrane Fuel Cells

Abstract

The chemical degradation of the perfluorinated sulfonic acid (PFSA) ion-exchange membrane as a result of an attack from a radical species, originating as a by-product of the oxygen reduction reaction, represents a significant limiting factor in a wider adoption of low-temperature proton exchange membrane fuel cells (LT-PEMFCs). The efficient mathematical modeling of these processes is therefore a crucial step in the further development of proton exchange membrane fuel cells. Starting with an extensive kinetic modeling framework, describing the whole range of chemical processes leading to the membrane degradation, we use the mathematical method of sensitivity analysis to systematically reduce the number of both chemical species and reactions needed to efficiently and accurately describe the chemical degradation of the membrane. The analysis suggests the elimination of chemical reactions among the radical species, which is supported by the physicochemical consideration of the modeled reactions, while the degradation of Nafion backbone can be significantly simplified by lumping several individual species concentrations. The resulting reduced model features only 12 species coupled by 8 chemical reactions, compared to 19 species coupled by 23 reactions in the original model. The time complexity of the model, analyzed on the basis of its stiffness, however, is not significantly improved in the process. Nevertheless, the significant reduction in the model system size and number of parameters represents an important step in the development of a computationally efficient coupled model of various fuel cell degradation processes. Additionally, the demonstrated application of sensitivity analysis method shows a great potential for further use in the optimization of models of operation and degradation of fuel cell components.

1. Introduction

Low-temperature fuel cells with proton exchange membranes (LT-PEMFC) are considered as one of the most promising technologies for future energy conversion, especially in the transport sector, where the extended range of fuel cell driven vehicles compared to battery driven vehicles might prove crucial in their market adoption [1]. Despite this clear advantage, the use of fuel cells is hindered by the issues concerning their limited life-time, caused by relatively fast degradation of two main fuel cell components, the catalyst layer and proton exchange membrane (PEM) [2]. Efficient, physically based models of degradation, aimed at a better understanding of the underlying processes and more efficient development of mitigation strategies, are therefore pivotal in the further development of fuel cell technology [3].

In this paper, we focus on the modeling of the degradation in perfluorinated sulfonic acid ionomer membranes (e.g., Nafion^®, most commonly used in LT-PEMFC). A zero-dimensional kinetic framework, describing chemical membrane degradation, consisting of 23 chemical species participating in 23 relevant coupled chemical equations, was recently published by Frühwirt et al. [4]. The basis of the model is the Fenton reaction [5,6,7] between hydrogen peroxide H₂O₂ and iron impurities (Fe${}^{2+}$):

$${\mathrm{Fe}}^{2+}+{\mathrm{H}}_2{\mathrm{O}}_2+{\mathrm{H}}^{+}\to {\mathrm{Fe}}^{3+}+{\mathrm{H}}_2\mathrm{O}+·\mathrm{OH}.$$

(1)

The generated hydroxyl radicals (·OH) subsequently react with other chemical species (H₂O₂, Fe²⁺, Fe³⁺, H₂), leading to the formation of additional radicals, namely hydroperoxyl (·OOH) and hydrogen (·H) [8,9,10,11,12,13,14,15,16,17,18,19,20]. Nevertheless, only the attack of ·OH at the side chains of the ionomer membrane [21,22,23,24,25] (Nafion^®) is regarded in the model, due to thermodynamic and kinetic considerations [26,27,28,29].

A large number of modeled chemical species, however, leads to a difficult calibration of the model, due to the large number of parameters and slow computation times, especially when coupled to spatially resolved degradation models, where the species concentrations need to be tracked in each cell of the computational mesh. These factors limit the applicability of the model when simulating the long-term degradation effects, most relevant in real fuel cell systems. Several simplified models of membrane degradation, with a smaller number of chemical species involved in fewer chemical equations, can be found in the literature addressing this issue [23,24,25,27,30,31,32,33]. The model simplifications, however, are typically done intuitively, without a systematic mathematical approach and analysis, even though the use of mathematical procedures [34] for model reduction is a part of a standard repertoire in chemical engineering, with the combustion of alkanes [35,36] (e.g., n-heptane and iso-octane), the hydrodesulfurization of light gas oil [37], the biogenic production of coalbed methane [38], and biochemical systems such as a yeast glycolysis model [39] being a few examples of its use.

The aim of this paper is therefore a derivation of a computationally-optimized yet adequately accurate chemical membrane degradation model for low-temperature proton exchange membrane fuel cells. Such a model, which is derived from the complete set of chemical reactions, involved in polymer electrolyte membrane (PEM) degradation [4] via systematic reduction, opens an unprecedented capability in high-fidelity, real-time, virtual exploration of membrane degradation. The applicability of the model thus ranges throughout the entire V-development process, starting with the office system layout studies, where very short computational times and a high prediction capability of the model are crucial, due to the unavailability of the experimental data and a large variety of configurations and control strategies that need to be exploited. In addition, the developed model is also fully compatible with the detailed 3D analyses, where it can be integrated as a 0D reactor model in each cell of the computational mesh. Besides its merits in terms of consistency with real phenomena in the fuel cells, and thus a high prediction capability of the model, the use of a single model throughout the entire development process also significantly reduces the effort in preparing the models and increases the consistency of the modeling tool chain.

The model reduction is done by the use of a sensitivity analysis-based model reduction approach [40] to reduce the number of modeled chemical equations in the system to the minimal set that still sufficiently describes the degradation of the FC membrane. The elimination of reactions is analyzed from a physicochemical perspective, explaining and justifying the results of the mathematical procedure. The change in the computational efficiency of the model is tracked throughout the model reduction procedure, using the stiffness analysis of the system of the differential equation to describe the model. The reduced model is tested in different operating conditions, demonstrating the robustness of the simplified model and its applicability in a wide range of fuel cell operating parameters.

2. Mathematical Procedures

2.1. Degradation Model Description

The full model (FM) of the chemical degradation of PEM, recently published by Frühwirt et al. [4], comprising 23 reactions between 23 chemical species, is listed in

Table 1. Overview of the model [4] obtained by mathematical simplification procedures. The rate constants $k=Aexp\left(-E_a/RT\right)$ are given in L mol⁻¹ s⁻¹ and at temperature T = 298.15 K, activation energies E_a in kJ mol⁻¹ and the pre-exponential factors A in L mol⁻¹ s⁻¹. Only the reactions printed in bold are kept in the reduced model. Reactions printed in bold italic are removed for testing in the over-reduced model.

	Reaction		k [L mol${}^{\mathbf{-}\mathbf{1}}$ s${}^{\mathbf{-}\mathbf{1}}$]	A [L mol${}^{\mathbf{-}\mathbf{1}}$ s${}^{\mathbf{-}\mathbf{1}}$]	${\mathit{E}}_{\mathit{a}}$[kJ mol${}^{\mathbf{-}\mathbf{1}}$]
R${}_1$	Fe2+ + H2O2 + H+ → Fe3+ + H2O + ·OH	$r_1$ = $k_1$[Fe2+][H2O2]	66	$1.05\times {10}^8$ [9]	$35.4$ [9]
R${}_2$	Fe3+ + H2O2 → Fe2+ + H+ + ·OOH	$r_2$ = $k_2$[Fe3+][H2O2]	$8.6\times {10}^{-4}$	$8.43\times {10}^{18}$ [11]	126 [11]
R${}_3$	Fe2+ + ·OH + H+ → Fe3+ + H2O	$r_3$ = $k_3$[Fe2+][·OH]	$3.6\times {10}^8$a [20]	$1.37\times {10}^{10}$	9 [20]
R${}_4$	Fe2+ + ·OOH + H+ → Fe3+ + H2O2	$r_4$ = $k_4$[Fe2+][·OOH]	$1.2\times {10}^6$ [12]	$2.74\times {10}^{13}$	42 [12]
R${}_5$	Fe3+ + ·OOH → Fe2+ + H+ + O2	$r_5$ = $k_5$[Fe3+][·OOH]	$8.1\times {10}^4$ [19]	$4.9\times {10}^{10}$	33 [18,19]
R${}_6$	H2O2 → 2 ·OH	$r_6$ = $k_6$[ H2O2]	$6.5\times {10}^{-23}$ [s${}^{-1}$]	${10}^{13}$ [s${}^{-1}$] [10]	201 [10]
R${}_7$	·OH + H2O2 → ·OOH + H2O	$r_7=k_7$[·OH][H2O2]	$3.0\times {10}^7$a [15]	$8.43\times {10}^9$	14 [15]
R${}_8$	·OOH + H2O2 → ·OH + H2O + O2	$r_8=k_8$[·OOH][H2O2]	$3.7$ [8]	$2.71\times {10}^6$	33.5 b [8]
R${}_9$	2 ·OOH → H2O2 + O2	$r_9=k_9$[·OOH]${}^2$	$9.1\times {10}^5$a [13]	$3.70\times {10}^9$	20.6 [17]
R${}_{10}$	2 ·OH → H2O2	$r_{10}=k_{10}$[·OH]${}^2$	$5.3\times {10}^9$a [16]	$1.33\times {10}^{11}$	8 [14]
R${}_{11}$	·OOH + ·OH → H2O +O2	$r_{11}=k_{11}$[·OOH][·OH]	$1.1\times {10}^{10}$a [16]	$3.39\times {10}^{12}$	14.2 [17]
R${}_{12}$	·OH + H2 → H2O + H	$r_{12}=k_{12}$[·OH][H2]	$3.9\times {10}^7$a [16]	$9.14\times {10}^{10}$	19.2 [16]
R${}_{13}$	·H + O2 → OOH	$r_{13}=k_{13}$[·H][O2]	$2.2\times {10}^{10}$a [16]	$1.37\times {10}^{12}$	10.3 [17]
R${}_{14}$	SC−SO3H + ·OH → SC−O· + HOCF2CF2SO3H	$r_{14}=k_{14}$[SC−SO3H][·OH]	$3.7\times {10}^6$ [29]	$6.80\times {10}^{18}$	70 c
R${}_{15}$	SC−O· + 3 ·OH → BB−O· + 6 HF + 3 CO2	$r_{15}=k_{15}$[SC−O·][·OH]	$3\times {10}^7$ [23]	$5.51\times {10}^{19}$	70 c
R${}_{16}$	BB−O· + ·OH → 2 −(CF2)nCOOH + 3 HF	$r_{16}=k_{16}$[BB−O·][·OH]	$8.5\times {10}^7$ [23]	$1.56\times {10}^{20}$	70 c
R${}_{17-23}$	−(CF2)nCOOH + 2 ·OH → −(CF2)n-1COOH + 2HF + CO2	$r_{un}=k_{17}$[−(CF2)nCOOH][·OH]	$1.0\times {10}^6$d [28]	$1.84\times {10}^{18}$	70 c

^a Rate constants were recalculated for $T=298.15$ K. ^b Activation energy is estimated between 6 to 10 kcal mol${}^{-1}$. ^c Average activation energy found in various degradation experiments [27,41]. ^d Upper limit for the attack of ·OH on CF₃COOH [28].

. The model can roughly be separated into three submodels (divided by horizontal lines in Table 1), describing three distinct processes, leading to the chemical degradation of the perfluorinated membrane. The first five reactions (R${}_1$–R${}_5$), which we call the “Fenton” block, describe the redox reactions involving iron ions (Fe²⁺ and Fe³⁺) and hydrogen peroxide H₂O₂, resulting in the formation of ·OH and ·OOH radical species. The second block of eight reactions (R${}_6$–R${}_{13}$) describes the reactions between radical species, hydrogen peroxide, water, hydrogen, and oxygen, thus called the “inter-radical reaction” block. The last block of reactions (R${}_{14}$–R${}_{23}$) describes how the attack of ·OH degrades the perfluorinated membrane and is therefore referred to as the “degradation” block. An overview of the reactions and the involved species (sulfonic acid head group SC−SO₃H, degradation product HOCF₂CF₂SO₃H, the intermediates SC−O· and BB−O·, and backbone fragments −(CF₂)_nCOOH and −(CF₂)_n−1COOH) is schematically presented in Figure 1.

Note, however, that the concentrations of four chemical species which are primarily determined by the fuel cell operation regime, i.e., water H₂O, oxygen O₂, hydrogen H₂, and protons H⁺, are not calculated by the membrane degradation model, thus reducing the number of relevant differential equations of the model to 19.

2.2. Sensitivity Analysis

To determine which of the reactions in the original model are actually relevant and which can be omitted, the model sensitivity was tested on a typical example of a long-term (100 h) degradation test at membrane conditions typically found in a LT-PEMFC during its operation, i.e., T = 90 ${}^{\circ }$C, combined with the iron ion concentration $\begin{array}{c}\hfill {\left[{\mathrm{Fe}}^{x+}\right]}_{\mathrm{tot}}=\left[{\mathrm{Fe}}^{2+}\right]+\left[{\mathrm{Fe}}^{3+}\right]\end{array}$ = 10 ppm and hydrogen peroxide concentration $\begin{array}{c}\hfill \left[{\mathrm{H}}_2{\mathrm{O}}_2\right]\end{array}=1$ mM [4]. Note that the total iron ion concentration ${\left[{\mathrm{Fe}}^{\mathrm{x}+}\right]}_{\mathrm{tot}}$ is conserved during the degradation, due to iron recycling (reactions R${}_1$–R${}_5$), while it is assumed that the same amount of hydrogen peroxide is produced and consumed (steady-state fuel cell conditions) [4].

The time evolution of species concentrations, given by the reactions in Table 1, can be written as a homogeneous, non-linear dynamic system

$$\dot{\mathbf{c}}\left(t\right)=f(\mathbf{c}\left(t\right),\mathit{\theta }),\mathbf{c}(t=0)={\mathbf{c}}_0$$

(2)

where $\mathbf{c}\left(t\right)$ denotes the state vector (consisting of the species concentrations), and $\mathit{\theta }$ the parameter vector (consisting of activation energies $E_{a,i}$ and pre-exponential factors $A_i$, used to calculate the reaction rate constants $k_i=A_{i}exp\left(-E_{a,i}/RT\right)$, where R denotes the general gas constant). In order to evaluate the significance of the individual reactions, the parametric output sensitivity, with respect to their corresponding pre-exponential factors $A_i$ is determined. The sensitivity matrix is given by

$$\mathbf{S}=\left[\begin{array}{c}\mathbf{s}\left(1\right)\\ ⋮\\ \mathbf{s}\left(t_k\right)\\ ⋮\end{array}\right]$$

(3)

with

$$\mathbf{s}\left(t_k\right)=\frac{\partial \mathbf{c}\left(t_k\right)}{\partial {\mathit{\theta }}_{\mathrm{red}}}=\left[\begin{array}{cccc}\frac{\partial c_1\left(t_k\right)}{\partial {\theta }_{1,\mathrm{red}}}& \cdots & \frac{\partial c_1\left(t_k\right)}{\partial {\theta }_{j,\mathrm{red}}}& \cdots \\ ⋮& \ddots & \\ \frac{\partial c_i\left(t_k\right)}{\partial {\theta }_{1,\mathrm{red}}}& \cdots & \frac{\partial c_i\left(t_k\right)}{\partial {\theta }_{j,\mathrm{red}}}& \cdots \\ ⋮& & & \ddots \end{array}\right].$$

(4)

Thereby, ${\mathit{\theta }}_{\mathrm{red}}$ denotes the reduced parameter vector, consisting of the pre-exponential factors $A_i$ for which the sensitivity analysis is carried out, and $t_k$ denotes the (equidistantly sampled) discrete time, indexed by k. Note that the output sensitivity depends on the numerical magnitude of output signals and parameter values. To avoid this, a normalized equivalent is often preferred [42,43]. The elements of the sensitivity matrix (4) are then scaled by

$${\tilde{s}}_{i,j}\left(t_k\right)=\frac{{\theta }_{\mathrm{red},j}}{{\tilde{c}}_i}\frac{\partial c_i\left(t_k\right)}{\partial {\theta }_{j,\mathrm{red}}},$$

(5)

where ${\tilde{c}}_i$ is a suitable scaling factor reflecting at least the order of magnitude of the transient solution of species i. In order to numerically evaluate (5), the partial derivative therein needs to be computed. An obvious choice would be the approximation of the partial derivative via finite differences

$$\frac{\partial c_i\left(t_k\right)}{\partial {\theta }_{j,\mathrm{red}}}\approx \frac{c_i(t,{\theta }_{j,\mathrm{red}}+\Delta {\theta }_{j,\mathrm{red}})-c_i(t,{\theta }_{j,\mathrm{red}})}{\Delta {\theta }_{j,\mathrm{red}}}.$$

(6)

Thereby, $c_i(t,{\theta }_{j,\mathrm{red}}+\Delta {\theta }_{j,\mathrm{red}})$ denotes the solution of species i, obtained by solving system (2) with a parameter vector whose j-th element is altered by a finite step $\Delta {\theta }_{j,\mathrm{red}}$. The choice of the step size $\Delta {\theta }_{j,\mathrm{red}}$ is crucial for a sufficient approximation. Although the approximation error would theoretically converge to zero with a step-size approaching zero, any computer arithmetic has a finite precision and the numerical round-off errors, due to the computation of the difference in (6), become dominant. In practice, the optimal step-size that minimizes the sum of approximation errors and numerical round-off errors is not trivial to obtain. As an alternative, the complex-step derivative has been proposed [44]. An approximation of the partial derivative is then given by

$$\frac{\partial c_i\left(t_k\right)}{\partial {\theta }_{j,\mathrm{red}}}\approx \mathrm{Im}\left(\frac{c_i(t,{\theta }_{j,\mathrm{red}}+i\Delta {\theta }_{j,\mathrm{red}})}{\Delta {\theta }_{j,\mathrm{red}}}\right)$$

(7)

with $i=\sqrt{-1}$ and $Im(·)$ as the operator that evaluates the imaginary part of a complex number. As the computation of differences is avoided in (7), numerical round-off errors are eliminated and the step size can be chosen arbitrarily small.

The significance of parameters with respect to the model output, and therefore the significance of an individual reaction with respect to the overall evolution of concentrations over time, can be concisely analyzed by applying an eigenvalue decomposition to

$${\mathbf{V}}^{-1}\mathsf{\Lambda }\mathbf{V}={\mathbf{S}}^T\mathbf{S}$$

(8)

with $\mathsf{\Lambda }=\mathrm{diag}\left(\left[{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_j,\cdots \right]\right)$ as the eigenvalues in descending order and $\mathbf{V}$ containing the corresponding eigenvectors. For the model reduction, the eigenvector corresponding to the smallest eigenvalue is evaluated. It points directly to the parameter space in the direction of those parameter(s) who, if varied, influence the model output the least. Based on the argument that a reaction whose corresponding pre-exponential factor does not influence the resulting time evolution of species, it can therefore be neglected and removed from the system of differential equations. The errors made in the resulting time traces by doing so were evaluated by calculating the mean-squared error value of the difference between the original result and the result after the removal of the concerning parameter and hence, reaction. This procedure was done iteratively, gradually removing the least significant parameters from the model.

An additional simplification of the model can be made by lumping the reactions between ·OH radicals and the perfluorinated backbone (R${}_{17}$–R${}_{23}$). In the original model [4], the chains of different numbers n of CF₂ groups adjacent to the COOH active group ($0\le n\le 7$) are tracked independently, with the attack of ·OH resulting in the shortening of the chain. Since all these reactions are essentially the same, the number of reactions can be significantly reduced by rather modeling the concentration of the end groups, lumped into a single variable [−COOH] = ${\sum }_{n=1}^7$ [−(CF₂)${}_n$ COOH], and the total concentration of CF₂ groups, [CF₂] = ${\sum }_{n=1}^{7}n$ [−(CF₂)${}_n$ COOH] [23]. Note that the upper limit of $n=7$ in this case does not refer to the total length M of backbone chain between two side chains, but rather to only half of it, since the degradation can proceed along the backbone from two opposite directions. The reactions (R${}_{17}$–R${}_{23}$) are thus reduced to only two differential equations, tracking the concentrations of the moieties [−COOH] and [CF₂]:

$$\begin{array}{cc}\hfill \frac{d}{dt}& \begin{array}{c}\hfill [-\mathrm{COOH}]\end{array}=2r_{16},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \frac{d}{dt}& \begin{array}{c}\hfill \left[{\mathrm{CF}}_2\right]\end{array}=-k_{17}\begin{array}{c}\hfill [-\mathrm{COOH}][·\mathrm{OH}]\end{array},\hfill \end{array}$$

(10)

with the initial values for fresh membrane [−COOH]$(t=0)=0$ and [CF₂]$(t=0)=M\times$[SC−SO₃H], where M denotes the total length of backbone between two side-chains, which is typically $M=13$ for Nafion^® [23].

3. Results and Discussion

3.1. Model Reduction

In order to assess the model error resulting from the removal of reactions at each iteration i, the following mean-squared error (MSE) criterion was used:

$$MS{E}_i=\frac{1}{N}\left|\right|{\mathbf{C}}_{\mathrm{FM}}-{\mathbf{C}}_{\mathrm{RM},i}{\left|\right|}_F,$$

(11)

$${\mathbf{C}}_j=\left[\begin{array}{cccc}{\mathbf{c}}_j\left(1\right)& \cdots & {\mathbf{c}}_j\left(t_k\right)& \cdots \end{array}\right],j\in \left\{\begin{array}{cc}FM,& RM,i\end{array}\right\}$$

(12)

denotes the stacked output matrix for the full model (FM) and the reduced model at the i-th iteration (RM, i) for all time samples of the simulation, and $\left|\right|·{\left|\right|}_F$ denotes the Frobenius norm. The error criterion values are plotted in Figure 2. We see that the removal of the first nine reactions (up to $A_{12}$ in Figure 2) does not change the results more than about ${10}^{-6}$. The removal of the next reaction, (reaction R${}_4$, $A_4$ in Figure 2), changes the results significantly. This point was therefore used as a limit for the model reduction.

The reactions kept in the reduced model (RM) are written in bold font in Table 1. We see that the model reduction eliminates all the reactions in the “inter-radical reaction” block, indicating that the reactions between radicals, hydrogen peroxide, hydrogen, and oxygen are less relevant for the model results, which will be further addressed in Section 3.2. Furthermore, the reaction between Fe³⁺ and ·OOH (R${}_3$) can also be omitted, leaving a reduced set of iron redox reactions. All reactions between radicals and the perfluorinated membrane in the “degradation” block, on the other hand, are also kept in the RM. Note that ·H is only involved in reactions R${}_{12}$ and R${}_{13}$, which are both eliminated in the RM, so the tracking of this species can be discarded. This results (combined with the lumping of backbone groups in a single variable) in a model tracking the concentration of only 12 species in the RM, compared to 19 in the FM.

The validity of the RM was tested by comparing the time traces of most relevant species concentrations, resulting from both the FM and the RM in simulating the degradation of the Nafion^® membrane during 100 h degradation test, with the same parameters as used in the sensitivity analysis in Section 2.2. The results of both simulations are shown in Figure 3, with the results of the FM being plotted as solid lines and the results of the RM as dashed lines.

To verify that further reduction of the model to a smaller number of equations cannot be justified, the over-reduced model (ORM) was constructed by removing further two reactions R${}_4$ and R${}_5$, describing the Fenton recycling of ·OOH radicals. These two additionally removed reactions are marked in bold italic font in Table 1, and the results of this over-simplified model are plotted as a dotted line in Figure 3.

The plotted results indicate that the RM is completely sufficient to produce the same results as the FM in the test simulation, with the only visible difference being in the concentration of ·OOH at the beginning of the simulation, but the results coincide perfectly after the first few seconds of the simulated time. On the other hand, there are significant discrepancies between the results of the ORM and the FM, which is a clear sign that the further reduction of the model is not possible and that the reactions R${}_4$ and R${}_5$ should be kept in the model.

To further verify the plausibility of the RM, the next sections are dedicated to the analysis of its plausibility from a physicochemical perspective, its computational efficiency, and the range of its validity under different boundary conditions.

3.2. Physicochemical Perspective of the Model Reduction

As proven by the mathematically-based procedure of the reaction elimination, the whole block of “inter-radical reactions” (R${}_6$–R${}_{13}$ in Table 1) can be neglected in the reduced model. To further justify this finding and elucidate its significance, the model reduction was considered from the physicochemical perspective of the modeled processes. This is done by considering the concentrations of radical species, calculated from the model, their reaction rate constants, and the comparison with the experimental data found in the published literature.

In operating fuel cell systems, ·OOH, ·OH, and ·H are present at the same time [45], thus it is difficult to investigate the reaction kinetics for the individual radical species (especially ·H). The chemical reactions of the reactive oxygen species (·OOH and ·OH) are well characterized (see ref. [4] for a comprehensive summary), however, the reactivity of ·H towards perfluorinated compounds is not completely understood. Therefore, although the FM yields the transient concentrations of the radicals ([·OOH] $\approx {10}^{-8}$ M; [·OH] $\approx {10}^{-15}$ M; [·H] $\approx {10}^{-17}$ M, which was also found in a previous modeling study [27]), only the ·OH-induced ionomer degradation is included in the FM (R${}_{14}$–R${}_{23}$).

Ghassemzadeh et al. [22] found that ·H abstracts fluorine atoms (to form HF) from the tertiary carbon C-F bonds in both the side and main chain of the ionomer. In the same study, the effects of selectively generated ·H and ·OH radicals (via electron beam irradiation of aqueous solutions of H₂SO₄ and H₂O₂, respectively) on the degradation of Nafion^® 211 were investigated. The authors could show that the fluoride release was almost twice as high, when the membrane is subjected to the attack of ·OH. However, kinetic data could not be extracted, as the absolute relative concentrations of the radicals are not known [22]. As [·H] is about two orders of magnitude lower than [·OH], the rate constant of fluorine abstraction by ·H has to be at least $k_F={10}^2·k_{17}\approx {10}^8$ M${}^{-1}$ s${}^{-1}$ at room temperature (or $k_F\approx {10}^{10}$ M${}^{-1}$ s${}^{-1}$ at T = 90 ${}^{\circ }$C), so that the ·H-induced degradation is in competition with the attack of ·OH on the ionomer. For ·OOH, on the other hand, hydrogen abstractions (which are necessary to induce main chain unzipping, R${}_{17}$–R${}_{23}$) are thermodynamically disfavoured [26], hence it is unlikely, that ·OOH initiates membrane degradation. Nevertheless, ·OOH might react with radical intermediates formed during the degradation.

The second reaction block (R${}_6$–R${}_{13}$) consists only of reactions of the radical species among each other or with other reaction partners like H₂O₂, H₂ and O₂, which are typically present in millimolar concentrations [27]. As a product of usually high rate constants, but small concentrations, the reaction rates of these reactions become very small. This explains why the whole second reaction block is completely removed during the simplification procedure (R${}_6$ is characterized by a negligibly small rate constant, so that this reaction can be completely omitted). During model reduction, reactions involving ·H were completely removed and the role of ·OOH was reduced to its participation in the redox recycling of iron ions (see first reaction block in Table 1).

3.3. Stiffness Analysis of the Reduced Model

The computational speed and efficiency of large systems of coupled differential equations depend on the number of the equations and even more prominently on the stiffness of the system [46]. While the reduction in the number of equations from 19 in the FM to 12 in the RM is obvious, the stiffness of the system needs to be additionally analyzed. This was done by the standard procedure of linearization of the differential Equation (2) into a form

$$\dot{\mathbf{c}}\left(t\right)=\mathbf{A}·\mathbf{c}\left(t\right)$$

(13)

where matrix $\mathbf{A}$ is composed of partial derivatives of function $f\left(\mathbf{c}\right(t),\mathit{\theta })$ over the variables $\mathbf{c}\left(t\right)$, and the calculation of the eigenvalues of $\mathbf{A}$. The inverse values of the eigenvalues can be interpreted as the time scales of different components of the system evolution, and in the case of all eigenvalues having a negative real part, the ratio between the largest and smallest real part is used as an indicator of the stiffness of the system [46].

However, in the case of the analyzed membrane degradation model, several eigenvalues of $\mathbf{A}$ have a vanishing real part or it is very close to zero, which makes the analysis based on the stiffness ratio of little use. Nevertheless, the eigenvalue decomposition of the system still enables a useful insight into its dynamic properties by revealing the reactions and variables governing the step size in the numerical integration. The eigenvalue decomposition of $\mathbf{A}$ using the numerically calculated concentrations $\mathbf{c}\left(t\right)$ from Section 2.2 reveals that the magnitude of the eigenvalue with the most negative real stays almost constant during the model reduction ($\approx {10}^9$ s${}^{-1}$) up to the point when reaction $R_{14}$ is removed. The eigenvector, associated with this eigenvalue, contains a significant component related to the concentrations of species [·OH], [SC−SO₃H], [SC−O·], [BB−O·], [−COOH], indicating that the computational efficiency of the model is to a large extent determined by the dynamics of the side-chain degradation. Since this mechanism is an integral part of any model of chemical degradation and cannot be eliminated during the model reduction, the stiffness of the system and consequently its computational efficiency is not expected to change significantly due to the simplification of the model and the reduction of its size.

3.4. Range of Validity of Reduced Model

To verify that the model reduction is justified for a wider range of parameters, and not only for the ones used for the sensitivity analysis, the FM and the RM were tested with various initial values of iron ion concentration [Fe²⁺] + [Fe³⁺], hydrogen peroxide concentration [H₂O₂], and temperature T. The applicability of the model reduction was evaluated by calculating the $MSE$ (Equation (11)) for the results of the FM and the RM, with the total simulation time $t_{max}$ being either 100 h or the time needed for the modeled membrane side-chain concentration in the FM to be degraded to one percent of its initial concentrations: ([SC−SO₃H]$\left(t_{max}\right)=0.01$[SC−SO₃H]$(t=0)$).

Figure 4 shows the parameter space in which the RM produces results comparable to the FM. As seen in Figure 4a), the concentration of hydrogen peroxide [H₂O₂] can be increased by a factor of 100 and the concentration of iron ions by more than 1000, compared to the reference conditions ([H₂O₂]_ref = 1 mM, [Fe^2+/3+]_ref = 10 ppm), and the discrepancy between the FM and the RM will still remain relatively small ($MSE<{10}^{-2.5}$). It also shows that the increase in the concentrations of the species will increase the discrepancy between the results of the FM and the RM, therefore making the reduction less justified. The impact of the temperature is far less pronounced, which is shown in Figure 4b), where the line of $MSE={10}^{-2}$ is plotted for several different temperatures, indicating only a small change in the quality of the results of the RM. Note that in these limiting cases, the concentrations of hydrogen peroxide and iron ions in the membrane are far above the realistic values expected in a real LT-PEMFC [4].

The robustness of the reduction is further confirmed by plotting the time traces of the relevant species in the perfluorinated membrane during degradation, with model parameters chosen such that significant differences between the results of the FM and the RM in Figure 5 are expected. The simulation in Figure 5 was calculated with a high iron ion [Fe²⁺] + [Fe³⁺] = ${10}^4 \times$[Fe^2+/3+]_ref and high hydrogen peroxide concentration [H₂O₂] = ${10}^3 \times$[H₂O₂]_ref ($MSE\approx {10}^{-2}$), resulting in almost complete degradation of the active side-chains in a few minutes. Even though this scenario is completely unrealistic in real fuel cells, the results of the RM (dashed line in Figure 5) are still very close to the prediction of the FM (solid line in Figure 5). The ORM (dotted line in Figure 5) on the other hand, gives completely wrong results, which indicates the importance of modeling a sufficiently large set of species and chemical reactions. These results indicate that the RM is completely sufficient in simulating and explaining the membrane degradation phenomena in conditions typical for a real fuel cell system, while also producing relatively accurate results in significantly harsher conditions, leading to extremely fast membrane degradation.

4. Conclusions

In this study, we used a sensitivity analysis-based model reduction to systematically reduce the number of chemical reactions and species used in a model of chemical degradation of perfluorinated membranes. The original model, proposed by Frühwirt et al. [4], comprising 23 chemical reactions, coupling 19 species, was reduced to only 8 reactions, coupling 12 species. The systematic analysis of the effects of the model boundary conditions, such as temperature and the iron ion and hydrogen peroxide concentrations, indicates a robustness of the reduced model, which yields results identical to those of the original model for a range of boundary conditions greatly exceeding the ones expected in a real fuel cell system. Furthermore, we show that any further reduction of the model results in significant errors in the model performance, concluding that the proposed reduced model is optimal in describing the chemical degradation of perfluorinated membranes in low-temperature proton exchange membrane fuel cells. The steps taken in the model reduction were justified by physicochemical analysis, explaining that the rates of the chemical reactions, removed from the model, are too small to be relevant, due to the low concentrations of radical species (·OH, ·OOH, and ·H) and their reaction partners. The matching conclusions of the sensitivity analysis and physicochemical analysis show that a good understanding of the chemical background of the analyzed processes is indispensable for model simplification. Since the precise knowledge of chemical kinetics is a crucial ingredient of both methods, further experimental investigation of relevant chemical reactions is strongly encouraged.

The stiffness analysis of the reduced model, however, does not show a significant reduction in its time complexity. Despite this, the simplification of the membrane degradation model is an important step in further development of more complex, interconnected models of fuel cell degradation. On one hand, the significant reduction in the number of differential equations is a crucial step when implementing the model as 0D reactors in individual cells of a larger spatially resolved computational mesh, where the space complexity of the system becomes an important issue. On the other hand, the reduction in the number of modeled chemical reactions eliminates unnecessary parameters from the model and thus greatly simplifies the calibration of the complex interconnected models that couple operation and degradation mechanisms of various fuel cell components. Although several such models, where peroxide and iron ion production and diffusion through the membrane are also considered, already exist [23,24,25,32,33], the use of the reduced model might significantly simplify its parametrization and calibration.

Even more importantly, the paper presents a general and robust procedure for the simplification of chemical membrane degradation models. In the future, this could be used in an improved, extended model, where the degradation mitigating species, such as cerium salts [47,48,49,50] or cerium oxide [32,51,52], are included to study their effect and ensure that the extended model remains as small and efficient as possible.