How to Build Simple Models of PEM Fuel Cells for Fast Computation

Abstract

Hydrogen is one of the leading candidates in the search for an alternative to fossil hydrocarbon fuels. The spread of these technologies requires a real-time control of generator performances. Artificial intelligence (AI) and mathematic tools can make smarter the smart grid. The electrochemical modeling can be coupled successfully with artificial intelligent approach, if these models can be quickly computed with a large numerical stability. This chapter shows a methodology to build this kind of modeling work. Thanks to a simplified but physically reasonable model of PEM fuel cell, we will show that the reactant access (oxygen) or water management (a product of the reaction) and the reaction rate can be easily described with low computing time consuming. In addition, the artificial neural network could be trained with a reduced amount of data generated by these cell models.

1. Introduction

In the current context of the spread of renewable energies, these are by nature variable, therefore subject to both daily and seasonal intermittencies. Electrochemical devices have been successful in proving their applicability in terms of energy storage (power to gas) [1]. Controlling in real time, predicting the performance is the advantage that electrochemical generators can offer.

In addition, electricity consumption and production must, at every moment, be in perfect adequacy with the demand of the users. However, this demand is variable and cannot be completely regulated. The production must be able to adapt instantly to the demand, to preserve the stability of the network. Thus, exchanging information between the production and storage sites becomes a major issue. This will result in a strategy of predictive actions in a power grid strongly constrained by intermittent sources of energy. The interpretation and exchange between blocks of storage and energy supplies will be the key to a decentralized energy production and smart grid [2, 3]. According to Ramchurn et al. [4] the grand challenge for artificial intelligence is to put the smarts into the Smart Grid.

Electrochemical modeling can provide “smart” tools for smart grid. The establishment of a mathematical model of an electrochemical generator is handled by the scientific culture of the researcher who establishes it. Therefore, a great deal of subjectivity appears in any modeling work, the approach of a mechanic/energy specialist, an electrochemist, or a physicochemist will differ mainly in the basic assumptions of modeling (model simplifications). However, whatever the cultural origin of the modeler, the numerical resolution of a multiphysical problem makes it possible to assure three major functions in the phase of development [5]:

• Assistance with the understanding of experimental results

• Study and optimization of design

• The prediction of performances

The modeling of an electrochemical system involves the mathematical expressions of the physical phenomena that take place there (a priori). Obviously, all model representations only offer a fragmentary assessment of the real systems. These various representations are distinguishable by their scales of time and space. However, the notion of adapted or appropriate modeling remains subjective. Indeed, as described in the literature about fuel cell models [6, 7], each of the approaches has limitations of description or prediction, and their main interest is to highlight one specific process. Despite this subjectivity, the model must prove its validity.

The validity of the model could be named external, i.e., related to theories, concepts, assumptions, and experimental data. Thus, the model is theoretically valid if it accepts theories or models already validated. In addition, if the model well matches to its potential of scientific explanation (the state of the art), one will qualify its heuristic validity. However, building a model cannot be done without solving it in all its intended range. Consequently, it is also necessary to define criteria of internal validity, which are criteria of evaluation of the model independent of the theories, results, and hypotheses. The algorithm (solver) must be appropriate, and the evaluation errors must remain within “valid” limits.

External validations that may be acceptable include empirical validity (the model corresponds to the available data) or pragmatic validity (the model satisfies the intended use).

A fuel cell is a nonlinear and strongly coupled dynamic system. It is a multi-input multi-output system based on multiphase flow, electrochemical reactions, and heat transfer. For example, the control strategies of PEMFC can be built on a prediction of the future output of the system to compute the current control action [8]. In practice, the current control action is obtained by solving online an optimization problem. The aim of the optimization problem is to find the optimum of a cost function that minimizes the mean squared difference between predicted outputs and target values.

To spare computation time due to computing of multiphysic fuel cell models, artificial intelligence (AI) techniques are useful as alternate approaches to conventional multiphysic modeling: e.g., artificial neural network (ANN) simulator could be employed to predict the fuel cell behavior [9, 10, 11, 12]. The ANN could be trained with a reduced amount of data generated by a validated cell model [13]. Once this network is trained, it can predict different operational parameters of the fuel cell reducing the computation time [14]. This strategy has many possibilities [15]: spectroscopic analysis, prediction of reactions, chemical process control, and the analysis of electrostatic potentials. The ANN is trained to learn the internal relationships from data. These data may be taken from the real process even if there are noisy. The database should have a significant size and contain the maximum combination of inputs-outputs covering the studied domain. Therefore it is possible to generate database from model [14].

As shown below, electrochemical modeling can be coupled successfully with artificial intelligent approach; thus in the following subsections, we provide basic mathematical models of fuel cells. These models can be quickly computed with a large numerical stability.

2. Main electrochemical phenomenon

An electrochemical cell is characterized by the I-V (current–voltage) behavior: the current that passes across the cell to the applied cell voltage. The I-V relation depends on various physical phenomena and is fundamental to achieve efficiently electrochemical conversion. When current is drawn using computational tools, the current density may not be uniformly distributed on the electrode surfaces. The performance and lifetime of electrochemical cells, such, is often improved by a uniform current density distribution. Therefore, it is necessary to optimize the current distribution. The electric current is a flow of electric charges: through the electrolyte between the anode and the cathode in the ions form and within the wires and current collectors/electrode materials in the electrons form. When the overall current of the cell is equal to zero, for example, on the disconnection of an electrode from the power supply or the load, the cell voltage is equal to UOCV the open circuit voltage:

$$I=0;U_{\mathit{\text{cell}}}=U_{\mathit{OCV}}=E_{a,0}-E_{c,0}$$E1

where $E_{a,0}\mathit{\text{and}}{E}_{a,0}$ are the potential of each electrodes at OCV. It is shown that the $U_{\mathit{OCV}}$ is related to the difference of free enthalpy differences of each reaction, involving the number n of electrons exchanged and the Faraday constant:

$$U_{\mathit{OCV}}=\frac{∆G_i}{\mathit{nF}}$$E2

where $∆G_i$ is the Gibbs energy of the species i with specified cell temperature (T) and cell pressure (P). At a non-zero current, the cell voltage of an electrochemical reactor (electrolyzer) is greater than $U_{\mathit{OCV}}$, and electrochemical generator (fuel cell or battery) is smaller than $U_{\mathit{OCV}}$ due to the various irreversibilities standing in the electrochemical conversion. The electrode potentials differ from the equilibrium values $E_{a,0}\text{and}{E}_{a,0}$, and this difference is called overvoltage:

$${\eta }_a=E_a-E_{a,0}$$E3

$${\eta }_c=E_c-E_{c,0}$$E4

According to this description, an anode overvoltage is positive, while the cathode overvoltage is negative in all cases. The overvoltages depend on the current density at the electrode, depending on the involved electrochemical reactions, the electrode materials, and several operating conditions: concentration species, flow rate, etc. The current density is an extensive quantity that can be defined in any point of the electrochemical device, i.e., at the surface of the electrodes and through the electrolyte. In addition Ohm’s law expresses the current density according to the local potential gradient, using the conductivity as follows:

$$\overrightarrow{I}=-\sigma \overrightarrow{\nabla E}$$E5

where $\sigma$ is the ionic conductivity of an electrolyte. In the case of a one-dimensional system where the electrolyte is confined in a finite space between two surface electrodes S and separated from each other by the distance e, this potential difference (called ohmic drop ${\eta }_{\mathit{ohm}}$) can be expressed as follows:

$${\eta }_{\mathit{ohm}}=-\frac{e}{\mathit{\sigma S}}I$$E6

where J is the total current through electrochemical device: J is the average current density I multiplied by the electrode surface S. For other geometries, it will be necessary to calculate ${\eta }_{\mathit{ohm}}$ from the relation (Eq. 5) by integration. The ohmic drop in the electrolytic solution of an electrochemical cell is one of the parameters to be evaluated to optimize the cell efficiency. This limitation is also called primary current distribution. Thus, the cell potential is calculated from following expression:

$${\mathrm{U}}_{\mathrm{\text{cell}}}={\mathrm{U}}_{\mathrm{OCV}}-{\mathrm{\eta }}_{\mathrm{a}}+{\mathrm{\eta }}_{\mathrm{c}}-{\mathrm{\eta }}_{\mathrm{ohm}}$$E7

Inside the nonaqueous electrochemical device (Figure 1), the primary current distribution is well controlled, because a solid electrolyte is confined in a finite space between two surface electrodes and only the misalignment of both electrodes could affect the current distribution through the electrolyte. In nonaqueous case, the optimization endeavor is devoted to secondary (electrochemical activation) and tertiary (mass transport limitation) current distributions. In this context, fuel cells are the best example to scrutinize energetic balance of electrochemical devices.

Figure 1.

Schematic description of a MEA.

The main advantage associated with fuel cells is that they are not limited by Carnot efficiency. Besides, no moving parts are required to convert thermal energy into mechanical energy. The energy release from interatomic bonds of the reactants is converted efficiently into electrical energy. In this document, we consider proton-exchange membrane fuel cells (PEMFCs): a PEMFC consists of a polymer electrolyte sandwiched between two electrodes to form a membrane electrode assembly (MEA), placed between two graphite bipolar plates, which feed the device with gases and cool it down. At the anode, fuel H₂ is oxidized, liberating electrons and producing protons. The electrons flow to the cathode via an external circuit, where they combine with the proton and the dissolved oxidant O₂ to produce water and heat:

$${\mathrm{H}}_2\leftrightarrows {\text{2H}}^{+}+{\text{2e}}^{-}$$E8

$${\text{2H}}^{+}+{\text{2e}}^{-}+\mathrm{½}{\mathrm{O}}_2\leftrightarrows {\mathrm{H}}_2\mathrm{O}$$E9

The efficiency of fuel cell is easily computed according to:

$$\mathrm{\psi }\frac{{\mathrm{U}}_{\mathrm{OCV}}}{{\mathrm{U}}_{\mathrm{\text{cell}}}}$$E10

Proton transfer from the anode to the cathode via the membrane closes the electrical circuit. A careful water management is required to ensure that the membrane remains fully hydrated in order to improve the ionic conductivity and to avoid the electrode flooding (which occurs when an excess of liquid water restrains the active species access to the active layer). The gas diffusion electrodes (GDEs) are made up of two distinct areas: an inactive area (backing layer) and an active area (active layer) which is a place of the electrochemical and chemical reactions (Figure 1). Therefore, in the first approach it is possible to observe separately secondary current distribution (only in the active layer) and tertiary current distribution (only in the backing layer).

PEMFC and all solid electrolyte cells are the model devices in order to sketch the primary (in the polymeric membrane), the secondary (in active layer), and tertiary (in backing layer) current distributions. In addition, operating cell potential and each overvoltage could be compared to thermodynamic potential (Gibbs energy) to access enthalpic balance.

3. One-dimensional modeling of electrochemical membrane cell

A single cell can be described schematically as an assembly of several layers constituting four distinct areas. A fuel cell is a composite structure of anode, cathode, and electrolyte. Good electrochemical performance of the cell requires effective electrocatalysts. On both sides of the cell, the interconnect plates cumulate three functions: current collector, gas feeding via the gas channels, and thermal control thanks to cooling water channels. Flow field is used to supply and distribute the fuel and the oxidant to the anode and cathode electrocatalysts, respectively. The distribution of flow over the electrodes should ideally be uniform to try to ensure a uniform performance of each electrode across its surface. Thus, it is possible to develop a single-cell 1D model. Although most of the models used are one-dimensional, they correctly predict the electrochemical behavior of membrane electrode assembly.

In order to build a 1D model, a particular attention is required on water flux. The water managements is the key issue in single-cell modeling. In order to predict cell performance, the single-cell model must take into account gas diffusion in the porous electrodes, water diffusion, and electroosmotic transport through the polymeric membrane. Ramousse et al. [16] have developed one-dimensional coupled charge and mass transfer model in the electrodes. The three main types of model can be employed to evaluate the overvoltages at both electrodes:

• The Tafel law results from experimental observations: this simple and widely used description is based on the phenomenological Eq. (12).

• The classical catalyst layer model exhibited in the previous section (Appendix 1).

• The agglomerate model takes into account the existence of gas pore and the ohmic drop in the agglomerate. This model results to considering a cylindrical model of pores inside the electrode [17, 18].

In these three models of active layer, the authors have computed the same values of kinetic coefficients $b_i$ and $i_{o\left(a,c\right)}$ (Table 1).

Parameter	Values
γ: roughness factor (−)	100
δ: gas backing layer thickness (m)	230 × 10^−6
LCA: catalyst Layer thickness (m)	10 × 10^−6
δ2: size of the agglomerate (m)	3 × 10^−6
Lm: membrane thickness (m)	125 × 10^−6
ε^GDL: gas backing layer porosity	0.8
ε^CL: active layer porosity	0.5
j0c: exchange current density (A/cm^2)	4 × 10^−7
j0a: exchange current density (A/cm^2)	1 × 10^−2
bc: Tafel slope (mV/dec)	120
ba: Tafel slope (mV/dec)	30
νe: total number of exchanged electrons	4
$D_{H_2/H_{2}O}^{\mathit{eff}}$: effective diffusion coefficients (m^2/s)	1.63 × 10^−4
$D_{O_2/H_{2}O}^{\mathit{eff}}$: effective diffusion coefficient (m^2/s)	3.20 × 10^−5
$D_{O_2/N_2}^{\mathit{eff}}$: effective diffusion coefficients (m^2/s)	2.41 × 10^−5
$D_{H_{2}O/N_2}^{\mathit{eff}}$: effective diffusion coefficients (m^2/s)	3.35 × 10^−5
D^CL: diffusion coefficient in the CA (m^2/s)	10^−9 × (ε^CA)^1.5
EW: equivalent weight (kg/mole)	1.1
ρdry: dry Nafion density, (kg/m^3)	2020
τ0: osmotic coefficient	2.5/22
Dm: effective diffusion coefficients of water in the membrane (m^2/s)	3 × 10^−10

Table 1.

Common parameter values of PEMFC cells.

The simulations presented below describe, in steady state, one-dimensional mass transfer in the whole cell and charge and mass transfers in the electrodes. Then, only a numerical solution can be accessed. Thanks to numerical computation the cell overpotential η_a/c and the ohmic drop through the membrane η_ohm can be simulated as a function of current density.

The total current is then calculated by integrating the local current density all over the MEA surface:

$$I={\iint }_{S}i\mathit{dS}$$E11

Figure 2 exhibits the evolutions of the cathodic and anodic overpotentials, calculated in the same operating conditions, as functions of the current density. There are great differences between the overpotentials predicted by the three models, which emphasize the important influence of the active layer on fuel cell performances. The porous catalyst layer model (agglomerate model) and nonporous catalyst layer models require only intrinsic parameters: active layer thickness, platinum load (roughness factor), and kinetic coefficients (exchange current density). The polarization curves obtained with the porous and nonporous models exhibit Tafel behavior. Due to the slow diffusion of active species in the homogeneous solid phase, the nonporous model tends to overestimate the overpotentials. Particular attention should be paid to the use of the Tafel law. The effective platinum area being higher than the MEA section and the exchange current density i0 should be multiplied by a geometric parameter such as the roughness factor γ (Eq. 12) to provide corrected Tafel mode l as follows:

Figure 2.

Influence of electrode description on overvoltages à T = 60°C. (a) Anodic overvoltage; c, cathodic overvoltage; 1, Tafel law corrected with the roughness factor; 2, porous catalyst layer model (agglomerate model); 3, nonporous catalyst layer model.

$${\mathrm{\eta }}_{\mathrm{a}/\mathrm{c}}={\mathrm{b}}_{\mathrm{a}/\mathrm{c}}log\frac{\mathrm{i}}{\mathrm{\gamma }{\mathrm{j}}_{0.\mathrm{a}/\mathrm{c}}}$$E12

In Figure 2, the anodic prediction is incoherent because the anodic overvoltage is negative. Therefore, the Tafel law is only adapted to high current density, far from thermodynamic equilibrium.

However, Singh et al. [19] suggested that a multidimensional model could improve the description, since gas composition and temperature vary along the feeding channels. In addition, these phenomena become critical for large cell areas used to obtain great current intensities. The flow field allows gas to flow along the length of the electrode while permitting mass transport to the electrocatalyst normal to its surface. One of the simplest flow field designs consists of a series of narrow parallel rectangular channels where fuel or oxidant is fed at one end and removed from the opposite side. Figure 3 shows a schematic illustration of the PEMFC that consists of a membrane sandwiched between two gas diffusion electrodes, this assembly being pressed between two current collectors. Within a PEMFC, reactant depletion and cathodic water production occur along the length of the fuel cell due to the electrochemical processes. Gas supply to the catalytic active site takes place in the porous gas diffusion electrode that contains a GDL consisting of hydrophobic gaseous pores and a reaction region (CL). Namely, slow diffusion in the GDL and CL can induce oxygen depletion in the case of air cathode. The membrane, commonly Nafion®, acts as both a separator and an electrolyte. Water transport across the membrane results from the electroosmotic water dragging with proton migration from anode to cathode and water “back-diffusion” that reduces the concentration gradient. During the operation, these effects can be responsible of either membrane dehydration or flooding on the porous electrode.

Figure 3.

Scheme of a proton exchange membrane fuel cell and the equivalent electrical circuit for the whole PEMFC for DC and AC solutions.

A finite volume method using a computational grid [20] can be used to solve mass, charge, energy, momentum balances including transport through porous media, and chemical and electrochemical reactions within the porous electrodes in a gas diffusion electrode model. Freshly, Zhang et al. [21] have developed a two-phase multidimensional model to properly handle mass transport. The results of these computations reveal that liquid water transport inside and across the polymeric membrane plays an important role in PEM fuel cell water distributions. In addition, it was shown that increasing the contact angle at GDL/channel interface is found to be able to improve the water removal process in channels. However, the resolution of fuel cell models requires important calculation times to obtain the detailed variations of flow field, species concentrations, temperature, liquid saturation, and electronic and ionic phase potentials.

4. Simple multidimensional modeling

4.1 Pseudo 2D modeling

To cut down these difficulties, the pseudo 2D as 1D + 1D approach implies that (i) meander-like channel is replaced with the straight one and (ii) water and oxygen concentrations in the channel provide “boundary conditions” for local transport of these species across the cell. PEMFC cell model, which couples 1D transport across the cell with 1D description of the flow in the feed channel, has been proposed [22, 23]. An approach is to model the flow as an ideal gas in a straight channel of cross-section and ignore the influence of viscosity. Dohle et al. [24] have proposed a pseudo two-dimensional model of the cathode compartment of a PEMFC. The model is based on the continuity equations for the gases and assumes constant pressure and plug flow conditions. The influence of hydrodynamics was not considered. Regardless of the latter point, the model shows some interesting effects for a simple, parallel flow field channel configuration. The mole fraction change along the gas channel length results from the normal molar flux N_i^GDL in the y-direction through the gas diffusion layer (GDL) allowing to the electrochemical processes in the active layer of the GDL and the corresponding current density changes. Note that this normal flux also depends on z. The mass balance for a single-phase species i (= O₂, H₂, N₂) may be written as:

$$\frac{\mathrm{d}{\mathrm{F}}_{\mathrm{i}}}{\mathrm{dz}}={\mathrm{N}}_{\mathrm{i}}^{\mathrm{GDL}}$$E13

where F_i is respectively the molar flow rate for the component Xi. Both oxygen and hydrogen fluxes are proportional to current density, while the nitrogen flux is zero because nitrogen is neither consumed nor produced in the fuel cell:

$${\mathrm{N}}_{\mathrm{i}}^{\mathrm{GDL}}=\frac{{\mathrm{\nu }}_{\mathrm{i}}\mathrm{i}}{{\mathrm{\nu }}_{\mathrm{e}}\mathrm{F}}$$E14

For the specific case of PEMFC, water produced at the cathode can flow through the gas diffusion layer to reach the gas channel and transport through the membrane:

$${\mathrm{N}}_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{GDL}}+{\mathrm{N}}_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{m}}=\frac{\mathrm{i}}{2\mathrm{F}}$$E15

Even if water transport processes into the membrane are not well understood, the phenomenological model of water transport in the membrane takes into account water diffusion and water electroosmotic fluxes [25].

4.2 Pseudo 2D modeling: closed-form expression

Dohle et al. [24] have proposed an analytical solution of the oxygen profile along the channel. This model assumed that oxygen reduction is the determining rate step and neglects water transport. The oxygen concentration profile is given by:

$$\left[X_{O_2}\right]={\left[X_{O_2}\right]}_0\left(1-\left(1-\zeta \right)\frac{z}{L_{\zeta }}\right)$$E16

where L_ζ is the characteristic length of oxygen consumption.

$$L_{\zeta }=4F\frac{h{v}_0{\left[X_{O_2}\right]}_{\ast }}{\gamma j_{0,c}}{\left(\frac{{\left[X_{O_2}\right]}_0}{{\left[X_{O_2}\right]}_{\ast }}\right)}^{1-\zeta }\mathit{exp}\left(-\frac{2.3{\eta }_c}{b_c}\right)$$E17

where h is the channel height and v₀ is the inlet flow velocity, assuming a kinetic Tafel law for oxygen reduction (Eq. (12)). Note that for γ = 1, the oxygen profile distribution becomes:

$$\left[X_{O_2}\right]={\left[X_{O_2}\right]}_0\mathit{exp}\left(-\frac{z}{L_1}\right)$$E18

The mean current density in the cell is defined as an average over the channel length L:

$$i=i_{\mathit{lim}}\left[1-{\left(1-\left(1-\zeta \right)\frac{z}{L_{\zeta }}\right)}^{\frac{1}{1-\zeta }}\right]$$E19

The limiting current density (i_lim) represents the limiting current density, which is attained when all the available oxygen is consumed:

$$i_{\mathit{lim}}=4F\frac{h{v}_0{\left[X_{O_2}\right]}_0}{L}$$E20

The cathode overpotential is then given by:

$${\eta }_c=\frac{b_C}{2.3}\mathit{ln}\left[\frac{K_{\zeta }}{1-\zeta }{\left(1-\left(1-\frac{\overline{i}}{i_{\mathit{lim}}}\right)\right)}^{1-\zeta }\right]$$E21

where:

$$K_{\zeta }=4F\frac{h{v}_0{\left[X_{O_2}\right]}_{\ast }}{\gamma j_{0,c}L}{\left(\frac{{\left[X_{O_2}\right]}_0}{{\left[X_{O_2}\right]}_{\ast }}\right)}^{1-\zeta }$$E22

Based on the previous physical models, the experimental polarization curves can be approximated. The fitting parameters are generally the exchange current density, the Tafel slope, and the limiting current density. Figure 4 displays the effect of the channel length. The decrease in channel length L increases the limiting current density (I_lim). The decrease in L not only increases I_lim but improves cell performance in the whole range of current densities. Physically, for lower L, oxygen concentration in the channel is more uniform, and this improves the overall performance per unit cross-section area.

Figure 4.

Voltage–current curves for the three indicated values of channel length L (cm) [24].

This single-cell model exhibits that the reactive gas diffusion becomes a limiting step. However, the water management and electrochemical kinetics play a critical role. For a given cell voltage U_cell assumed to be constant along the gas channel, the current density distribution Iz must take in to account the membrane conductivity σ_m,z, which depends on membrane water content. Therefore, accurate numerical approaches should be proposed in order to calculate the electrical characteristics (current density, ohmic resistance, overpotential) as well as gas concentration distributions along the gas channel.

4.3 Continuous stirred tank reactor approach to fuel cell

A first rigorous 3D modeling of PEM single cell using CFD programs based on commercial finite volume technic solver has been developed by He et al. [26], to solve the fully coupled governing equations. The model assumes that liquid film is formed on the electrode surface during liquid water condensation and computes the diffusion flux, electroosmotic drag force, and water back-diffusion in order to assess to water management. According to Zhang and Jiao [27], despite the multiplicity of 3D CFD models available in the literature, a satisfactory 3D multiphase CFD model which is able to simulate the detailed gas and liquid two-phase flow in channels and reflect its effect on PEM fuel cell performance precisely still was not real, because it is difficult to solve coupling physics and computation amount is real barrier. Nandjou et al. [28] have proposed a pseudo 3D model to reduce the computational cost (i.e., the transport equations are formulated using a pseudo 3D approximation and coupled to an analytical electrochemical model at the catalyst layers/membrane interface). Nevertheless, the experimental measurement cannot easily emphasize or “validate” the computing results. Krewer et al. [29] have showed that 3D CFD results could be compared to experimental results via assuming that the experimental RTD of the inlet and outlet pipes can be used as input signal for the CFD simulations. Diep et al. [30], thanks to RTD experiments, have shown that a reduction in the unit cell model dimensionality to 1 + 1 or 1 + 0 D based on scaling arguments and contrasts with higher dimensional 3D models is accurate. Figure 5 exhibits a comparison between the model and the outlet RTD. Good agreement between laminar flow 1+ 0 D-based model and experimental data was found.

Figure 5.

Comparison between tracer step outlet RTD model fuel cell measurements and one-dimensional numerical model computations (cathode) [30].

On the other hand, to simplify the hydrodynamic description (modeling) or experimental observations, it is possible to employ the ideal reactor models. The continuous stirred tank reactor (CSTR) model does not require a long computing time. Benziger et al. [31] have employed this method to explore the auto-humidification of the PEMFC. In Figure 6, the anode and cathode chambers in the PEM fuel cell are modeled as two stirred tank reactors (STR). The previous models can only be considered as approximate, and more refined approach shall be preferred, in particularly for larger fuel cells.

Figure 6.

STR PEM fuel cell model [31].

Boillot et al. [32] have demonstrated that the bipolar plate is similar to a series of continuous stirred tank reactors. These authors have studied various models with or without exchange zones. They have observed that best fitting was obtained by the simple model of CSTR in series and the optimal number of CSTR (J) varied from four to six in the range investigated (Figure 7).

Figure 7.

RTD measurements and modeling of the column flow pattern (signals in arbitrary units) [32].

Deseure [33] proposed series of CSTR to solve the mass balance equation inside the gas channel (Figure 8). The electrode is divided into k elementary units. With such a hypothesis, gas transport in the gas channel is assumed to be infinitely fast, leading to homogeneous gas composition for each CSTR. In particular, the partial pressure of a component in the elementary unit (k) is equal to that of the outlet and to that of the elementary unit (k + 1) inlet. The resolutions of mass balances are numerically obtained thanks to the continuity of the flux between the GDL, the CL, and the membrane.

Figure 8.

Schematic description of the PEMFC half cell (cathode) [33].

The electrode is divided into k elementary units, where the gas transport in the gas channel is assumed to be infinitely fast leading to homogeneous gas composition for each CSTR. Partial pressure of a component in the elementary unit (k) is equal to the outlet one and equal to the elementary unit (k + 1) inlet one. Like the 1D + 1D approach, all elementary units have the same voltage U due to the high electrical conductivity of bipolar plate and backing layer. Therefore, partial pressures of each component in each elementary unit vary and make necessary to solve mass balances and current density in each elementary unit at the same time. The current–voltage expression used is the one proposed previously in Section 2.1. The molar balance equation of each CSTR with electrode geometric area S^k for each component in gas phase in steady state is written as:

$$F_{H_{2}O}^{k,\mathit{\text{in}}}-F_{H_{2}O}^{k,\mathit{\text{out}}}={\xi }^k{S}^k{N}_{H_{2}O}^{\mathit{GDL},k}$$E23

$$F_{O_2}^{k,\mathit{\text{in}}}-F_{O_2}^{k,\mathit{\text{out}}}=S^k{N}_{O_2}^{\mathit{GDL},k}$$E24

$$F_{N_2}^{k,\mathit{\text{in}}}-F_{N_2}^{k,\mathit{\text{out}}}=0$$E25

$$S^T=\sum _k{S}^k$$E26

The term ξ^k represents the ratio of gas molar flux divided by liquid molar flux that is estimated using a numerical optimization (1 ≤ ξ^k ≥ 0). The global mass balance for each CSTR is given by:

$$F_T^{k,\mathit{\text{in}}}-F_T^{k,\mathit{\text{out}}}=S^k\left(N_{O_2}^{\mathit{GDL},k}+N_{H_{2}O}^{\mathit{GDL},k}\right)$$E27

and the gas flux of consumption and production are expressed as:

$$N_{H_{2}O}^{\mathit{GDL},k}=\frac{i^k}{2F}-N_{H_{2}O}^{m,k}$$E28

$$N_{O_2}^{\mathit{GDL},k}=-\frac{i^k}{4F}$$E29

To satisfy saturated vapor pressure of water, the following expression is required:

$$P_{\mathit{sat}}=e^{\left(13.66-\frac{5096.23}{\left(T+273.15\right)}\right)}$$E30

The resolution of mass balance is numerically obtained thanks to the continuity of flux between the backing layer, the active layer, and the membrane. The molar flux density of oxygen and nitrogen are nil at membrane/electrode layer interface, and the molar flux density of water are given by the solution of the equation (Eq. (15)). Thus, the sum of diffusion and electroosmotic fluxes yields in steady state a first-order differential equation:

$$N_{H_{2}O}^{m,k}=\frac{{\tau }_0{i}^k}{F}\left[{\lambda }_a^k+\frac{{\lambda }_c^k-{\lambda }_a^k}{1-e^{k_m^k{L}_m}}\right]$$E31

$$k_m^k=\frac{\mathit{EW}{\tau }_0{i}^k}{{\rho }_{\mathit{dry}}{D}_{m}F}$$E32

${\lambda }_c^k$ and ${\lambda }_a^k$ stand for polymer water content at the membrane/active layer interface respectively at cathode and anode. The thermodynamic equilibrium between water vapor in the backing layers and liquid water in the polymer has been assumed using Eq. (30).

The ionic conductivity of polymeric membrane could be accessed thanks to the Neubrand [34] equation as follows:

$${\sigma }_{H^{+}}=e^{\left(-E_A\left(\frac{1}{T}-\frac{1}{353}\right)\right)}\left(0.0013{\lambda }^3+0.0298{\lambda }^2+0.2658\lambda \right)$$E33

$$E_A=2640e^{\left(-0.6\lambda \right)}+1183$$E34

${\lambda }_c^k$ and ${\lambda }_a^k$ are determined using sorption equilibrium of [35]:

$$\lambda =0.3+10.8\left(\frac{{\mathit{Py}}_{H_{2}0}}{\mathit{\text{Psat}}}\right)-16{\left(\frac{{\mathit{Py}}_{H_{2}0}}{\mathit{\text{Psat}}}\right)}^2+14.1{\left(\frac{{\mathit{Py}}_{H_{2}0}}{\mathit{\text{Psat}}}\right)}^3$$E35

The set of these equations could be solved by iterative methods and have needed only these inlet conditions ${\left[X_{O_2}\right]}_{\ast },{\left[X_{N_2}\right]}_{\ast },{\left[X_{N_2}\right]}_{\ast }$ and$F_{O_2}^{\mathit{\text{in}},0},F_{N_2}^{\mathit{\text{in}},0},F_{H_{2}O}^{\mathit{\text{in}},0}$, where $F_{O_2}^{\mathit{\text{in}},0}+F_{N_2}^{\mathit{\text{in}},0}+F_{H_{2}O}^{\mathit{\text{in}},0}=F_T^{\mathit{\text{in}},0}$. The conventional way to investigate the influences of the mass flux on electrochemical behavior uses the inlet molar flux, which are proportional to the current density and the stoichiometry at the cathode. This flux is currently expressed using the stoichiometry oxygen coefficient $r_{O_2}$, and this expression is given by:

$$F_{O_2}^{\mathit{\text{in}},0}=r_{O_2}{S}^T\frac{I}{4F}$$E36

with I as the average current density of i^k given by Eq. (17). In order to initialize the computing, the initial inlet molar flux density is obtained considering only one STR according to Benziger et al. [31] investigations. This model does not take into account the liquid water apparition; according to this assumption the values of $r_{O_2}$ are calculated in order to obtain a gaseous condition. A simple mass balance could estimate the minimal air flux or $r_{O_2}$ value for non-apparition of liquid water on the cathode:

$$r_{O_2}\left(T\right)=0.21\left[\frac{2+\frac{P}{P-P_{\mathit{sat}}\left(T\right)}}{\frac{P}{P-P_{\mathit{sat}}\left(T\right)}-\frac{{\mathit{HR}}_{C}P}{P-{\mathit{HR}}_C{P}_{\mathit{sat}}\left(T\right)}}\right]$$E37

This mass balance equation considers that the produced water by oxygen reduction is fully exhausted by the gas flux. However, this relation does not take into account the water flux provided from anode side: for this reason, the stoichiometry oxygen coefficient is proposed at high value. Indeed, Figure 4 shows the limit of liquid water apparition as a function of water activity (HR_C) of inlet gas. In any simulated cases, the selected values of $r_{O_2}$ allow to avoid the flooding. To insure no saturation by water in all studied cases, the feeding flux in simulations could be proposed at high level such as presented in Figure 9.

Figure 9.

The minimal stoichiometry oxygen coefficient $\left({\mathit{r}}_{{\mathit{O}}_{\mathit{2}}}\right)$ to obtain a single-phase condition as a function of INLET HR_C with t = 60°C without membrane water exchange.

4.4 Results of continuous stirred tank reactor approach to fuel cell

The aim of these simulations is to emphasize the dependence of the locale electrochemical behavior on overall electrical performance. Three kinds of humidification mode are investigated: the auto-humidification, the humidification of inlet gas, and the humidification via the membrane. The steady-state equations of each elementary unit are solved thanks to the “simplex” numerical optimization method. First, the effect of water and ion transport within the membrane has been investigated; Figure 9 shows the minimal stoichiometry oxygen coefficient r_O2 to obtain a single-phase condition without membrane water exchange. Therefore it is possible to scrutinize the average ohmic drop within the membrane versus humidification process. The ohmic behavior appears as linear contribution of overall polarization curves. Figure 10 shows thee ratio of ideal conductivity of fully hydrated Nafion® membrane (Eq. (35)) divided per average ionic conductivity trough the membrane according to the set parameter of gathered in Table. Considering a series of two CSTR with r_O2 = 2.2, it was observed that the auto-humidification or the humidification of inlet gas involves similar ion conductivity performance. Of course, a humidification of inlet gas at 45°C improves the effective membrane conductivity, but the impact of water production inside catalyst layer is more significant.

Figure 10.

Membrane ratio conductivity as a function of hydration process and current density considering a series of two CSTR, r_O2 = 2.2 and T = 60°C.

In the cases of dry inlet condition which corresponds to a PEM fuel cell operating with pure hydrogen and oxygen, the stoichiometric effects have been simulated by Han and Chung [36]. The results of these computations have shown decreasing membrane ratio conductivity with increasing stoichiometry oxygen coefficient $r_{O_2}$ and the main effect of water production. Therefore as highlighted by Nguyen and White [37], the membrane ratio conductivity is expected to depend on gas distribution. Figure 11 shows a decreasing effective conductivity with increasing number of CSTR series: it is worth mentioning that piston flow distribution (20 CSTR) increases the ohmic drop and the static mixing (1 or 2 CSTR) improves membrane conductivity.

Figure 11.

Membrane ratio conductivity as function of CSTR series approaches and current density, considering r_O2 = 2.2, T = 60°C with dry inlet conditions.

In Figure 12 the opposite effect is observed for r_O2 = 2.2; gas access is limited when static mixing is developed (two CSTR). Then, cathodic overvoltage is higher considering only 2 CSTR than using a series of 20 CSTR. Nevertheless, ionic access to the catalyst areas is modeled here, and for a large amount of gas flow (r_O2 = .11), piston flow distribution increases the cathodic overvoltage.

Figure 12.

Simulation of cathodic overpotential with varying CSTR approach and stoichiometry oxygen coefficient, considering r_O2 = 2.2, T = 60°C with dry inlet conditions.

At this stage, it is important to study the distribution of current density. In Figure 13, the simulations show the effects of oxygen gas access. The gas distribution is the main control on current densities along the gas channel. At the cathode inlet, dry gas feeds the fuel cell; the water content in the cathode stream increases with decreasing number of CSTR because the amount of water transported back to the anode by back-diffusion is lower than convective water flux. For higher water content in the cathode side, diffusion limitations are increased. Consequently the local current density decreases, and the net water flux across the membrane also decreases, resulting in a lower depletion rate of water from the anode gas stream, a lower production rate of water in the cathode, and a lower depletion rate of hydrogen and oxygen.

Figure 13.

Simulation of local current density distribution along the channel (r_O2 = 2.2 for dry conditions at T = 60°C with |η_c| = 0.55 V).

However, near the inlet of the cell, as the current density is high and cathodic water concentration is low, the electroosmotic drag coefficient is higher, and ionic access limits the reaction rate. This observation is possible in case of piston flow (20 CSTR). In Figure 14, at low stoichiometric value, a depletion of reactants can be observed close to the outlet, thus generating a cathodic overpotential increase and a dramatically current density decrease in this part. It is interesting to note that water molar fraction increases at the cathode side in connection to the electroosmotic drag of water, from anode to cathode, and the insufficient “back-diffusion” from the cathode to anode. Water transport and production at the cathode result in increasing water content at the cathode side along the gas channel while it decreases at the anode side.

Figure 14.

Simulation of local oxygen (■) and water (×) molar fraction distribution along the channel (r_O2 = 2.2 for dry conditions at T = 60°C with |η_c| = 0.55 V).

This study has drawn attention to the relative advantages of using a CSTR description. Moreover liquid water phase only exists at the cathode side, thus leading to the polymer drying out at the anode and consequently increasing ohmic resistance. In the future, it will be important to include the flooding or partial flooding effects that were not taken into account in this study.

5. Conclusion: single-cell modeling results and limits

We have shown, through a simplified but physically reasonable model of PEM fuel cell steady-state multiplicity, caused by reactant access (oxygen) or water management (a product of the reaction) and the reaction rate. Of course, catalyst loading is critical but 0-D models do not have enough accuracy to perform realistic predictions. Closed-form model (pseudo 2D) can provide relevant simulations but does not take into account water management. The succinct analogy with CSTR can provide fruitful analysis of water management through the fuel cell. This approach was settled with a view to determining the impact of gas distribution in the gas channel using a fluid dynamics observation (i.e., RTD) coupled with the usual electrochemical model. Water production and removal are analogous to heat production and removal. Therefore this analogy (energy balance) can be added to the present electrochemical pseudo 2D approach of mass balance (CSTR-PEM).

Operating conditions that can threaten the life of the PEM cell are not easy to detect; using simple models for fast computation with A.I. shall avoid the problematic lack of large experimental databases. Predictive models and A.I. can take up the challenge of “Smart Grid”.

A.1 The classical catalyst layer model

Under steady-state conditions, the dimensionless equations relating the concentrations of species (Eq. (17)) are obtained by setting $\partial /\partial t=0$:

$$\frac{{\partial }^2{X}_i}{\partial Y^2}+a_e\frac{{\nu }_i{j}_{0,a/c}{L}^2}{\left|{\nu }_e\right|F{D}_i^{\mathit{eff}}{\left[X_i\right]}_{\ast }}\mathit{exp}\left(\pm 2.3\frac{{\eta }_{a/c}}{b}\right)X_i=0$$EA1

By introducing the dimensionless parameter $X_i=\left[X_i\right]/{\left[X_i\right]}_{\ast }$ and $Y=y/L\mathit{CL}$, the Thiele’ modulus (mL) appears as follows:

$${\left(\mathit{mL}\right)}^2=a_e\frac{\left|{\nu }_i\right|j_{0,a/c}{L}^2}{\left|{\nu }_e\right|F{D}_i^{\mathit{eff}}{\left[X_i\right]}_{\ast }}\mathit{exp}\left(\pm 2.3\frac{{\eta }_{a/c}}{b}\right)=0$$EA2

Note that this dimensionless parameter mL enables us to compare diffusion resistances with regard to kinetic resistance on the active layer performances. The stationary problem has an analytical solution in the regime in which the effect of charge transport can be neglected. The steady-state formulation of the mass balance equation within the CL (Eq. (33)) leads to an expression of the geometrical current density i in the CL versus the overpotential:

$$i=\gamma j_{0,a/c}\mathit{exp}\left(\pm 2.3\frac{{\eta }_{a/c}}{b}\right)\frac{{\left[X_i\right]}_0}{{\left[X_i\right]}_{\ast }}\frac{\mathit{tanh}\left(\mathit{mL}\right)}{\mathit{mL}}=0$$EA3

where γ (=a_e × L_CL) is the roughness factor. This model predicts two distinct regions arising from control by kinetics and the other due to the control of both kinetic and diffusion. Under kinetics control, for low current density corresponding to mL <<1, concentration depletion in the active layer does not occur, and the current density remains equal to the kinetic one:

$$i=\gamma j_{0,a/c}\mathit{exp}\left(\pm 2.3\frac{{\eta }_{a/c}}{b}\right)\frac{{\left[X_i\right]}_0}{{\left[X_i\right]}_{\ast }}=0$$EA4

Conversely, under kinetics and diffusion control corresponding to mL>> 1, since tanh(mL) converges to unity, the current density can be approximated by [38]:

$$i=\sqrt{\frac{\left|{\nu }_e\right|F{D}_i^{\mathit{eff}}{c}_{i,0}\gamma j_{0,a/c}}{\left|{\nu }_i\right|L}}\mathit{exp}\left(\pm 2.3\frac{{\eta }_{a/c}}{2b}\right)\frac{{\left[X_i\right]}_0}{{\left[X_i\right]}_{\ast }}=0$$EA5

A catalyst layer controlled by both kinetics and diffusion is analogous to the classical chemical engineering problem of porous catalyst pellet with a first-order reaction controlled by mass transport limitation of reactant characterized by Thiele’s modulus. This explanation also holds for fuel cell model where the active layer may be considered to be a single pellet or agglomerate. When both kinetic and diffusion in the active layer are controlling the electrochemical process, it is well-known that limiting current behavior does not occur, but rather a doubling of the effective Tafel slope on a steady-state polarization curves and a first-order dependence on the concentration [38]. The shapes of the concentration profiles depend strongly on parameter mL [39]. The model allows us to understand easily how the diffusion limitation within the catalyst layer acts on the DC responses. Under kinetic control (mL < <1), the diffusion is fast enough to supply the reacting species up to the catalyst surface in the GDE, and the concentration profile remains close to a constant value in the catalyst layer. Conversely, with increasing diffusion limitation (i.e., increasing parameter mL), an increasing concentration gradient is predicted. Figure A1 shows the influence of the diffusion limitation within the active layer on the polarization curve. Under kinetic and diffusion control in the active layer, a doubling of the apparent Tafel slope is observed in the DC response. Simulation was performed by considering parameters summarized in Table 1.

Figure A1.

Simulated i - |η| curve for a catalyst layer.

$$\mathrm{b}=2\times 120\mathrm{mV}/\mathrm{dec}.$$

$$\mathrm{b}=120\mathrm{mV}/\mathrm{dec}.$$

List of symbols

a_i

the specific surface area (m⁻¹)

b

Tafel slope (V dec⁻¹)

C_dl

double-layer capacitance (F m⁻²)

c_T

total gas concentration (mol m⁻³)

d

pore diameter (m)

D_i

diffusion coefficient (m² s⁻¹)

E₀

equilibrium potential (V)

EW

dry membrane weight per mole of sulfonate group (kg mol⁻¹)

F

Faraday’s constant (C mol⁻¹)

F_i

molar flow rate for the component X_i (mol s⁻¹)

h

enthalpy (J mol⁻¹)

I

current (A)

i

current density (A m⁻²)

j

faradaic current density (A m⁻²)

L_CL

CL thickness (m)

M_i

molecular weight of i (kg mol⁻¹)

N_i

flux density of each dissolved species (mol m⁻² s⁻¹)

R_i

homogeneous reaction rate (mol m⁻³ s⁻¹)

r_i

heterogeneous reaction rate (mol m⁻² s⁻¹)

T

temperature (K)

t

time (s)

V

potential (V)

v

bulk velocity (m s⁻¹)

[X_i]

total gas concentration of the species (mol m⁻³)

x_i

mole fraction of the species X_i

Z

impedance (Ω)

z_i

charge (−)

ΔS

entropy (J mol⁻¹ K⁻¹)

δ

GDL thickness (m)

ε

porosity (−)

Φ

electrostatic potential

Γ_i

surface concentration of a solute species i (mol m⁻²)

γ

roughness factor

η

overpotential (V)

λ

water content of the membrane (−)

ν_e

exchange number of electron (−)

ν_i

stoichiometric coefficient of species i (−)

ρ

density (kg m⁻³)

σ

conductivity of the solution (S m⁻¹)

τ

tortuosity (−)

τ₀

electroosmotic drag coefficient (−)