Distributed and Lumped Parameter Models for Fuel Cells

Abstract

The chapter presents a review of modeling techniques for three types of fuel cells that are gaining industrial importance, namely, polymer electrolyte membrane (PEMFC), direct methanol (DMFC), and solid oxide (SOFC) fuel cells (FCs). The models presented are both multidimensional, suitable for investigating distributions, gradients, and inhomogeneities inside the cells, and zero-dimensional, which allows for fast analyses of overall performance and can be easily interfaced with or embedded in other numerical tools, for example, for studying the interaction with static converters needed to control the electric power flow. Thermal dependence is considered in all models. Some special numerical approaches are presented, which allow facing specific problems. An example is the Proper Generalized Decomposition (PDG) that allows overcoming the challenges arising from the extreme aspect ratio of the thin electrolyte separating anode and cathode. The use of numerical modeling as part of identification techniques, particularly by means of stochastic optimization approaches, for extracting the material parameters from multiple in situ measurements is also discussed and examples are given. Merits and demerits of the different models are discussed.

1. Introduction

Several types of fuel cells (FCs) are under development for small to mid-size applications, both mobile and stationary. The modeling of FCs is the subject of a vast body of literature, with contributions coming from the fields of chemistry, material science, and engineering. This paper focuses on the authors’ experience and provides references for further reading and for the derivation of more models. In particular, the distributed and lumped modeling of direct methanol fuel cells (DMFCs), proton exchange membrane fuel cells (PEMFCs), and solid oxide fuel cells (SOFCs) is addressed. A section is devoted to the numerical optimization of such devices and to the identification of the parameters appearing in their equations by means of stochastic optimization algorithms. The two approaches, that is, distributed and lumped modeling, derive from opposite necessities, namely on one hand, the necessity of studying local details of the physical phenomena and on the other hand the necessity of having computationally efficient tools for large-scale simulations and integration of systems.

2. Modeling challenges

The performance of a fuel cell depends on the materials used, on the sizes of the cell components, on their geometry and arrangement, and on the combination and interactions of all these factors, involving interdisciplinary effects of, notably, interface, mass charge and heat transport, electrochemistry, catalysis, and materials science. Consequently, identifying an optimized stack configuration is a very demanding task that can require long and expensive experimental programs. In such a framework, models are very helpful in exploring possible system behavior and addressing the search for optimal solutions. However, due to the diversity and complexity of such phenomena, which occur at multidimensional level and on a wide range of length and time scales, model analyses have to necessarily be carried out by computer simulations. The numerical modeling of fuel cells, often dubbed CFCD (computational fuel cell dynamics), deals with multidimensional mass transport phenomena, electrochemical kinetics, and transport of charge (electrons and ions), in complex temperature-dependent relationships that are strongly coupled to each other. Namely, they are strongly coupled nonlinear problems and their solutions require advanced iterative algorithms capable of efficiently ensuring converged and accurate solutions. The analyses can be conveniently formulated in terms of electric potentials of the electronic and electrolyte phases, whose equations are strongly nonlinear and coupled to each other at electrochemical kinetics level. Techniques have been developed for efficiently solving these potential equations. Fuel cells present a stratified structure, made of thin layers of different materials where interactions occur. These domains with very different dimensions along coordinate axes give rise to additional computational challenges. Electrochemical activity sites generate current densities that can reach values in the order of 1 A/cm² (referred to the cross-sectional area) so that a current of 500 A requires membrane cross sections of 20 × 20 cm or more, whereas typical thicknesses of gas diffusion layers (GDLs) are in the order of 2–3 × 10⁻² cm. These GDLs are separated by proton-exchange membranes (PEMs) which are thin electrolytes, with widths in the order of 1–5 × 10⁻³ cm.

PEMs allow ion transport while preventing the passage of electrons, a key feature of electrochemical devices in order to force electrons to flow in the external electric circuits. Ion-conductive membranes allow the flow of protons when sufficiently hydrated. In optimal conditions, the proton conductivity can reach values as high as 0.2 S/cm at 100°C, which is a fairly good value for ions, but a poor one if compared with the electronic conductivity of metals, so that in order to reduce the inherent voltage and power losses such electrolytes must be very thin. At the same time, the relative fluctuations of the thickness due to manufacturing must be small, to ensure the membrane performance to be very uniform. The numerical simulation of such a domain, with aspect ratio exceeding 10³, involves severe size problems: a regular hexahedral tessellation with 10 elements in the thickness direction implies 10⁹ nodes, namely a very demanding computational dimension [1]. Commercial CFD codes commonly use un-structured meshes, but this fact does not alter the dimensional complexity of the problem. As a further concern, depending on the transient time scale, a large number of time steps may be needed in order to accurately compute time dynamics. These features raise very challenging problems that can only be faced with parallel computing and multiprocessor computers. Analyzing a fuel cell behavior requires the full characterization of the materials used, that is, the determination of their chemical, physical, thermal, and electrical parameters [2], which are involved in: mass, heat, and charge transport in the electrolyte; thermodynamics and electrokinetics in the catalyst layers; mass, charge, and heat transport in the diffusion layers; and mass, heat, and charge transport in the bipolar plates and their flow channels. In fact, these parameters are needed in the equations (Nernst equation, Butler-Volmer equation, Darcy’s equation, Fourier’s law, Ohm’s law, etc.) constituting FC models [3]. Since these equations are strongly nonlinear and coupled, the models present a strong vulnerability to parameter uncertainties. A number of diagnostic methods, such as cyclic voltammetry, thin-film rotating ring-disc electrode (CV-TF-RRDE), electrochemical impedance spectroscopy (EIS), and broadband electrical spectroscopy (BES), can provide accurate ex situ measurements. Nevertheless, the use of the measurement data thus obtained to model fuel cells in working conditions presents several setbacks. On the other hand, parameter values fully suitable for operative conditions can be achieved by means of in situ measurements, but few are available and they are often difficult to carry out, allowing to determine a limited number of parameters. Examples are EIS, neutron radiography, and voltammetric and chronoamperometric approaches in the “driven-cell” mode. To cope with material characterization and parameters extractions, a number of sophisticated numerical techniques that resort to optimization methods for extracting more parameters at once from a large number of experimental data have been developed in recent years. Stochastic methods have proven to be particularly successful to meet such target. Efficient modeling suitable to provide a comprehensive understanding of fuel cell dynamics thus involves: physicochemical model identification, advanced numerical algorithms, materials parameters extraction, and model validation against detailed distribution data. This last aspect arises since data of global nature, such as I-V curves, can result inadequate to capture the material parameters with their correlations.

3. Analytical models

3.1 Proton exchange membrane fuel cells

PEMFCs have been subject to a vast body of studies aimed at modeling and some of the most important issues are described hereafter. In PEMFCs, the hydrogen oxidation reaction (HOR) occurring at the anode catalyst layer (CL) and the oxygen reduction reaction (ORR) at the cathode catalyst layer

$$\begin{array}{llll}\text{anode}:& H_2& \to & 2H^{+}+2e^{-}\\ \text{cathode}:& \frac{1}{2}{O}_2+2H^{+}+2e^{-}& \to & H_{2}O\end{array}$$E1

are segregated by the proton exchange membrane (Figure 1). According to the Nernst equation, the cell’s reversible voltage E varies with temperature T and gas pressures $p_{H_2},p_{O_2}$ (or equivalently, with concentrations $c_{H_2},c_{O_2}$) [4]:

Figure 1.

Sketch of a PEMFC section with anode and cathode flow channels of bipolar plates (BPs), diffusion layer (DLs), catalyst layers (CLs), and proton exchange membrane (PEM). Convective (in BPs) and diffusive (in DLs) fluid flows are sketched (courtesy of Journal of Power Sources).

$$E=E^0+\Delta E_s\left(T\right)+\Delta E_c\left(T,p_{H_2},p_{O_2}\right)=E^0+\frac{1}{\mathit{nF}}\underset{T^0}{\overset{T}{\int }}\mathrm{\Delta }\widehat{s}\left(T\right)\mathit{dT}+\frac{T}{f_e}ln\left[\left(\frac{c_{H_2}}{c_{H_2}^0}\right){\left(\frac{c_{O_2}}{c_{O_2}^0}\right)}^{0.5}\right]$$E2

where E⁰ = 1.229 V is the value in standard temperature and pressure conditions, ΔE_s is the entropic variation due to $\mathrm{\Delta }\widehat{s}$, and ΔE_c is the term related to the variation of gas pressures and hence of gas concentrations. By introducing the “bulk” (undisturbed) concentrations ${\overline{c}}_{H_2},{\overline{c}}_{O_2}$, ΔE_c can be split into two terms:

$$\mathrm{\Delta }{E}_c=\frac{T}{f_e}ln\left[\left(\frac{{\overline{c}}_{H_2}}{c_{H_2}^0}\right){\left(\frac{{\overline{c}}_{O_2}}{c_{O_2}^0}\right)}^{0.5}\right]+\frac{T}{f_e}ln\left[\left(\frac{c_{H_2}}{{\overline{c}}_{H_2}}\right){\left(\frac{c_{O_2}}{{\overline{c}}_{O_2}}\right)}^{0.5}\right]=\mathrm{\Delta }{E}_{\mathit{co}}+\mathrm{\Delta }{E}_{\mathit{cl}}$$E3

Here f_e = nF/R, where R is the gas constant, F is the Faraday constant, and n is the number of electrons transferred in the reaction. ΔE_co arises from the differences between the actual bulk concentrations and the standard-condition values, particularly in a no-load state. ΔE_cl springs from the differences between the concentrations at the triple-phase boundaries (TPBs, which are the sites of electrochemical reaction in the catalyst layers—CLs) and the bulk concentrations in the CLs, namely from the concentration gradients which arise to sustain the species molar flows needed in load conditions. In order to allow an accurate modeling over a wide temperature range, ΔE_s can be calculated by integration. Combining the above and excluding the gradient-dependent term ΔE_cl, we have:

$$E_{\mathit{OC}}=E^0+\frac{1}{\mathit{nF}}\underset{T^0}{\overset{T}{\int }}\mathrm{\Delta }\widehat{s}\left(T\right)\mathit{dT}+\frac{T}{f_e}ln\left[\left(\frac{{\overline{c}}_{H_2}}{c_{H_2}^0}\right){\left(\frac{{\overline{c}}_{O_2}}{c_{O_2}^0}\right)}^{0.5}\right]$$E4

The cell open circuit voltage (OCV) is slightly different from E_OC because some gas crossover through the membrane occurs in every condition, including no-load, producing small concentration gradients and hence a minimal ΔE_cl. On the other hand, when the cell is operated at a steady-state load with an electric current density j, its voltage V differs markedly from E_OC, due to much larger concentration gradients, and hence much larger ΔE_cl, and to overpotentials η_kh (voltage drops ΔV_kh, in electrical engineering terms, [3]):

$$V\left(j\right)=E_{\mathit{OC}}-{\eta }_{\mathit{aa}}-{\eta }_{\mathit{ac}}-{\eta }_{\mathit{ca}}-{\eta }_{\mathit{cc}}-{\eta }_m$$E5

The first subscript of the overpotentials indicates a, activation losses, or c, concentration losses, whereas the second subscript stands for a, anode, or c, cathode. The single subscript m indicates the membrane (PEM), where ohmic losses occur. All these overpotentials are strong functions of j.

3.1.1 Exchange current density: activation losses

The activation overpotentials or activation voltage drops, ΔV_aa and ΔV_ac, are an effect of the electrochemical kinetics that appear when the species react at the anodic and cathodic CLs. ΔV_aa and ΔV_ac increase with the rate of charge density separation ∂_tρ_e (namely the proton and electron creation at the anode and recombination at the cathode—∂_t represents the partial time derivative), which, in steady-state operation equates the current density at the TPB, j_TPB. ΔV_aa and ΔV_ac are typically modeled with the Butler-Volmer equation [4]. Due to the particular porous structure of the CLs, the area A_TPB of the TPBs where j_TPB is produced is much larger than the active cell cross-sectional area A (A_TPB/A can be larger than 10³) and, when modeling is devoted to analyzing the cell electric performance, the current density is preferably reported to the cross-sectional area, as:

$$j_0=\frac{A_{\mathit{TPB}}}{A}{j}_{0\mathit{TPB}}$$E6

By using Eq. (6), the Butler-Volmer equation allows to write the current density at the cross-sectional area of each CL as

$$j_t=j_0\left(\frac{c_R}{{\overline{c}}_R}{e}^{\alpha f_e\mathrm{\Delta }{V}_a/T}-\frac{c_P}{{\overline{c}}_P}{e}^{-\left(1-\alpha \right)f_e\mathrm{\Delta }{V}_a/T}\right)$$E7

In this equation, the total equivalent current density j_t accounts for the effect of hydrogen crossover on the overpotentials, that can be modeled as an equivalent internal loss current to be added to the cell internal current density and is not available at current collectors, but contributes to the activation overpotentials. j₀ is the exchange current density. The accurate evaluation of its values at the anode and cathode half-reactions is important, because they strongly affect the cell performance and round-trip efficiency. To take into account the effects of the temperature on j₀, an Arrhenius-like dependence with T can be considered [5]. As a consequence of the low temperatures at which PEMFCs operate and of the exponential dependence of j_t on ΔV_a, the activation losses are the major responsible factors for voltage drops at low current densities. In addition, ΔV_ac is typically one order of magnitude larger than ΔV_aa, so that the cathodic activation losses are the dominant effect at low current densities [3].

3.1.2 Concentration gradients

When the cell operates in load condition, the inflow of reagents and outflow of products are necessary to sustain the electrochemical reactions at the CLs. In turn, these species flows are ensured by convective mass flow in the transport channels of the bipolar plates (BPs) and by the diffusive mass transport in the diffusion layers (DLs, Figure 1). These flows are sustained by concentrations and pressure gradients of reagents and products, namely by values c and p at the CLs different from the bulk (and inlet) values $\overline{c}$ and $\overline{p}$ [6]. In order to model the concentration gradients ∇c, Fick’s first law N = –D∇c is often used, which invokes the medium diffusivity D and the gas molar flow vector N, related in turn to the current density vector j through the Faraday constant F and the charge carriers n as j = nFN, for both hydrogen at anode and oxygen at cathode. Since the diffusivity depends on temperature, the effect of the latter on concentration gradients can also have important role in the FC performance.

3.1.3 Concentration losses and limit conditions

Gas concentration and pressure gradients ∇c and ∇p produce the term ΔE_cl of the Nernst equation, reducing the electromotive force (emf) with respect to the no-load value E_OC and constitute the concentration losses which dominate the cell’s performance at high current densities. ΔE_cl can be split in the anodic and cathodic concentration voltage drops as:

$$-\mathrm{\Delta }{E}_{\mathit{cl}}={\eta }_{\mathit{ca}}+{\eta }_{\mathit{cc}}={\kappa }_{\mathit{ca}}\frac{T}{f_e}ln\left(\frac{{\overline{c}}_{H_2}}{c_{H_2}}\right)+{\kappa }_{\mathit{cc}}\frac{T}{2f_e}ln\left(\frac{{\overline{c}}_{O_2}}{c_{O_2}}\right)$$E8

where κ are the mass transfer coefficients. These two terms provide the anodic and cathodic current density limits $j_{\mathit{La}}$ and $j_{\mathit{Lc}}$, namely the theoretical values of the current densities which cause the CL concentrations and pressures to vanish at the TPBs and the half-reactions stop, when the fuel cell starvation occurs. Since the smaller limit current density occurs at the cathode, due to the lower diffusivity of oxygen compared to hydrogen, $j_{\mathit{Lc}}$ sets the device’s current density limit.

3.1.4 Membrane ohmic losses

The voltage drop in the membrane is proportional to the current and to its thickness and inversely proportional to its conductivity σ which depends on temperature and hydration, which is the ratio λ = c_w/c_as (with c_as = 1970 mol m⁻³) between water and sulfonic acid concentrations that varies in the range 0–22 for typical Nafion^® membranes. The dependence of conductivity on hydration can be expressed as:

$$\sigma \left(\lambda \right)={\sigma }_0\left(T\right)\mathit{B\lambda }$$E9

The linear dependence on λ via the dimensionless coefficient B is the adaptation of an empirical model [7] aimed at avoiding a negative value of σ at lower λ. The temperature dependence can be expressed with the model of Vogel-Tamman-Fulcher [8]. Although λ varies along the PEM’s thickness according to back-diffusion and electro-osmotic drag [9], the average between the PEM boundary values λ_a and λ_c can be used, consistently with a linear profile between λ_a and λ_c. These values depend on the water activities a_wa and a_wc of the reacting gases at the CLs, and are computed with an empirical polynomial [7], which depends on the water vapor saturation pressure p_ws which is also a function of the temperature [3].

3.1.5 Crossover

Fuel crossover consists of hydrogen that, instead of reacting at the anode, migrates through the PEM and reacts with oxygen at the cathode, without producing electric power. This is a major side effect that affects the FC performance and efficiency and depends on two mechanisms, diffusion, and electro-osmotic drag. These two contributions of hydrogen mass flow can be modeled as equivalent current densities so that the resulting equivalent crossover current density is:

$$j_{\mathit{co}}=\frac{{\mathit{nFD}}_{{\mathit{mH}}_2}}{d_m}{c}_{H_2}+\mathit{\text{nξλj}}$$E10

where $c_{H_2}$ is the hydrogen concentration at the anodic CL, D_mH2 is the hydrogen diffusivity (with an Arrhenius-like temperature dependence), d_m is the membrane thickness and ξ is a dimensionless electro-osmotic drag coefficient (giving a maximum value ν_w = ξλ = 2.5 at full hydration λ = 22, as reported). Crossover hydrogen is the main cause of the difference between the open circuit emf E_OC and the observed OCV V(0) [10]. It also causes a loss of stored energy that reduces round-trip efficiency. Also, oxygen crosses the PEM, but, since its diffusivity is much lower than hydrogen’s [11], this effect is usually neglected.

Dissipations occurring inside the cell produce thermal gradients which affect the temperature-dependent parameters. The main loss phenomena are Peltier heating (thermodynamic heat generation), losses due to the electrochemical kinetic activity at the anode and cathode CLs, and Joule losses in the PEM, so that, inside the cell, the dissipated power per unit area is:

$$P_{\mathit{\text{diss}}}=j\left(-\frac{\mathit{T\Delta }\widehat{s}}{\mathit{nF}}+{\eta }_a\left(j\right)\right)+\frac{d_m}{\sigma }{j}^2$$E11

Heat transport inside the cell depends on conduction, diffusion, and convection and interacts with thermal capacity in dynamic conditions [3]. An accurate enough estimation of the mean temperature T inside the cell with respect to the room and gas inlet temperature T_r can be obtained by:

$$T=T_r+k_T{P}_{\mathit{\text{diss}}}$$E12

with k_T a global thermal exchange coefficient. This expression can be reformulated as a function of the current density as:

$$T=T_r+k_{t1}j+k_{t2}{j}^2$$E13

where k_t1 and k_t2 are properly fitted parameters. In the numerical implementation of such models, consistent analytical expressions can be used without introducing approximation if the electric current density j is chosen as the independent variable to compute all voltage terms. In order to deal with the non-invertible Butler-Volmer equation, a look-up table can be conveniently used.

3.2 Methanol fuel cells

DMFCs suffer from two fundamental problems: (i) the sluggish kinetics of the methanol electro-oxidation reaction and (ii) the high degree of permeation of the methanol through the membrane (crossover). Analytical and numerical models are necessary for better understanding the interactions between mass transfer and electrochemical phenomena, and for optimizing the power output and runtime to interface electronics [3].

3.2.1 DMFC analytical models

A schematic of a typical DMFC inside a cell stack is sketched in Figure 2. It consists basically of an anode flow channel (AFC), an anode diffusion layer (ADL), an anode catalyst layer (ACL), a proton exchange membrane (PEM), a cathode catalyst layer (CCL), a cathode diffusion layer (CDL), and a cathode flow channel (CFC). Analytical models of DMFCs account for the following phenomena: electrochemical reactions occurring at catalyst layers, protonic conduction, and methanol crossover across the PEM, diffusion of reactants in porous media layers, fluid flow inside channels, and coupling between heat and mass transfer. Small DMFCs for portable electronics usually include a fuel reservoir with no pumps for feeding the cell and make use of atmospheric oxygen. Simplifying assumptions are typically made on both model geometry and physics to obtain the current-voltage characteristic in closed form. In 2-D models, methanol/oxygen concentrations and electric potentials depend also on the direction (y-axis) orthogonal to the species flow direction (x-axis). In 1-D models, this dependency is neglected by assuming constant concentrations along the y-axis.

Figure 2.

DMFC schematic with reactant and species flows (a, anode; c, cathode; pem, proton exchange membrane; fc, flow channel; dl, diffusion layer; cl, catalyst layer).

The following assumptions are usually made in analytical DMFC models: (i) constant physical parameters; (ii) one-phase model (vapor phase is neglected); (iii) ideal and diluted solutions; (iv) homogeneous electrochemical reactions in the electrodes; and (v) negligible overpotential variation across catalyst layers and along the y-axis.

The basic equation for analyzing the cell voltage as a function of current density is:

$$V\left(J\right)=E_0+\frac{\mathrm{\partial E}}{\mathrm{\partial T}}∆T-{\eta }_{\mathrm{a}}\left(J\right)-{\eta }_c\left(J\right)-\frac{{\mathrm{\delta }}_{\mathrm{m}}}{{\mathrm{\sigma }}_{\mathrm{m}}}J-R_{s}J$$E14

where $E_o$ is the (constant) standard cell potential, E is the fuel cell standard potential, T is the temperature, ${\sigma }_m$ and ${\delta }_m$ are the electric conductivity and the thickness of the membrane, η is the activation overpotential, J is the current density, and $R_s$ is the overall contact area specific resistance (ASR). This resistance per unit cross section, typically assumed constant, accounts for all resistances between gas diffusion layers and current collectors. The current density at the anode in Eq. (14) is related to the activation voltage overpotential at the anode ${\eta }_a$. In the simplifying assumption of the Tafel equation, it can be expressed as:

$$J_a=J_{a,\mathit{ref}}\frac{C_{\mathit{ac}}}{C_{\mathit{ac},\mathit{ref}}}\mathit{exp}\left(\frac{{\alpha }_{a}F}{\mathit{RT}}{\eta }_a\right)$$E15

where $J_{a,\mathit{ref}}$ and $C_{\mathit{ac},\mathit{ref}}$ are the reference current density and the reference concentration at the anode, respectively. The voltage overpotential ${\eta }_a$ is required to overcome the activation energy of the electrochemical reaction, and causes an energy loss. A relationship similar to Eq. (15) holds for the cathode side, where the atmospheric oxygen is provided at the catalyst layer through the gas diffusion layer and the flow channel.

In the one-dimensional model proposed in [12], the cell performance is assessed by describing mass transport in the porous electrode structures, including methanol crossover, and the potential and concentration distributions in electrodes. Multiphase flow is accounted for in [13], where liquid-vapor mixtures at anode and cathode compartments are considered. A two-dimensional analytical model of a DMFC, which describes electrochemical reactions on the anode and cathode and main transport phenomena in the fuel cell including methanol crossover, diffusion of reactants in porous media layers, and fluid flow in the reactants distributor, is presented in [14].

Previous models, mathematically and physically based, cannot be used directly in CAD software for electronic circuits, which is used to design the DMFC-circuitry interface. By using the one-dimensional assumption, a lumped circuit model for simulating the DMFC runtime is derived in [15]. The methanol consumption at the reservoir is simulated by mass conservation equation and the charge/discharge electrode dynamics on the short time scale is simulated by a fictitious equivalent capacitance. A dynamic nonlinear circuit model for passive DMFC, including water mass flow and membrane hydration, is presented in [16]. By synthesizing physics equations in circuit form, this model is able to take into account mass transport, current generation, electronic and protonic conduction, methanol adsorption, and electrochemical kinetics. Adsorption and oxidation rates, which affect the cell dynamics, are modeled by a detailed two-step reaction mechanism. The fuel cell discharge and methanol consumption are computed by combining mass transport and conservation equations. As a result, the runtime of a DMFC can be predicted from current and initial methanol concentration.

3.2.2 DMFC semi-empirical models

Most DMFC performance models combine differential and algebraic equations with empirically determined parameters. Simpler empirical expressions can be derived from such models, allowing designers and engineers to predict the fuel cell performance as a function of different operating conditions (such as pressure, temperature, or fuel concentration). In [17], a semi-empirical approach to account for the limiting current behavior is proposed. The general expression for the polarization curve is based on analytical formulations of the catalyst layer reaction and includes the overpotential due to transport limitation in diffusion backing layer and the effect of methanol crossover. A simple semi-empirical model for evaluating the cell voltage of a DMFC is provided in [18] by using the semi-empirical Meyers-Newman rate equation for assessing the anode overpotential and a Tafel-type kinetic model for the cathode. Ohmic losses are accounted for by considering constant and uniform membrane conductivity. Such semi-empirical models are useful to describe the steady-state fuel cell behavior. Dynamic models can resort as well to semi-empirical relationships in order to limit the model complexity in time-domain simulations, notably when implementing models in real-time controllers. In [19], a completely different approach for extracting equivalent parameters is used. Instead of fitting current-voltage equations, a nonlinear equivalent circuit is derived from impedance fuel cell models. The equivalent circuit is composed of nonlinear electrical circuit elements that include resistors, capacitors, and inductors in order to simulate the DMFC performance in a wide operative range, including steady-state and transient conditions. A model of a micro-DMFC battery is developed in [20] to capture all pertinent dynamic and steady-state electrical performance parameters, including capacity and its dependence on current and temperature, open circuit voltage, methanol-crossover current, polarization curve and its dependence on concentration, internal resistance, and time-dependent response under various loading conditions. The main advantage of this model is that it is able to predict the DMFC runtime for a given usable energy capacity of the methanol reservoir. A 1-D analytical model for simulating both DMFC static and dynamic operations is interfaced to a PSO (particle swarm optimization) stochastic algorithm in order to maximize the battery duration while minimizing methanol crossover [21]. In [22], a semi-analytical model of a passive-feed DMFC with non-isothermal effects and charge conservation phenomenon is proposed. The 1-D model is suitable for predicting the current-voltage curve by resolving iteratively both (bi-phase) mass transfer, electrochemical, and heat transfer equations, starting from imposed cell current, ambient temperature, and methanol feed concentration. The phenomenon of cathodic mixed potential, which is related to methanol crossover is accounted for as well by the DMFC model. Adaptive operations for voltage stabilization are proposed in [23] by using a simplified semi-empirical model combined to an on-demand control system. A rapid response of DMFC (less than 5 s) to variations of operating parameters is experimentally observed, ensuring the applicability of the proposed adaptive control strategy.

3.3 Solid oxide fuel cells

A solid oxide fuel cell (SOFC) consists of an anode and a cathode with a ceramic electrolyte between them that transfers oxygen ions. Oxidant reduction occurs in the cathode catalyst layer, oxygen ions are transported through the electrolyte and oxidation of the fuel occurs at the anode catalyst layer. A SOFC typically operates between 700°C and 1000°C, that is, temperatures at which the ceramic electrolyte becomes sufficiently conductive of oxygen ions. In fact, the electrolyte ionic conductivity is a strongly increasing function of the operating temperature. Moreover, since the operating temperature accelerates electrochemical reactions, precious metal catalysts are not required to promote the reactions and cheaper materials such as nickel can be used as catalysts. In addition, the SOFC can be fed with conventional hydrocarbon fuels, reforming being performed inside the cell. SOFCs can be planar or tubular, with the former having gained increasing success because of easier manufacturing and higher performance.

Among the ceramic electrolytes used in SOFCs, the most important is yttria-stabilized zirconia (YSZ), which is the most conductive. The anode is generally made of nickel/yttria-stabilized zirconia cermet and the cathode is an LSM layer chemically expressed as La1−xSrxMnO3. The open circuit potential of a SOFC is given by the Nernst equation, whereas the activation overpotentials in both electrodes are high, so that the electrochemical kinetics of the both electrodes can be approximated by the Tafel equation, with the concentration dependence of exchange current density given by [24]. Zirconia (zirconium oxide, Zr O₂) has a crystalline structure and is stabilized in the allotropic cubic form YSZ with the addition (doping) of Yttria molecules (yttrium oxide, Y₂O₃). Moreover, Yttria significantly increases oxygen holes, namely an Y³⁺ ion replaces a Zr⁴⁺ inducing an oxygen hole to maintain electrical neutrality and such oxygen holes facilitate the transport of O⁺ ions in the lattice, by means of jumps. Doping with yttrium therefore dramatically increases the concentration of charge carriers (holes) and conductivity is proportionally improved. However, if the concentration of Yttria increases too much, the holes interact with each other thus reducing their mobility. Maximum conductivity is obtained with doping concentrations of about 8%: in such conditions, hole concentration is about 10² times greater than the zirconia intrinsic concentration. The probabilistic dependence of the charge carrier mobility on the energy yields an exponential dependence of the conductivity on temperature:

$$\sigma \left(T\right)={\sigma }_o{e}^{-W_a/\mathit{kT}}$$E16

A typical conductivity value is σ ≅ 0.02 S/cm at T = 1000 K, with activation energy W_a = 0.5–1.2 eV. This is a rather high conductivity value but not the highest among other types of FC electrolytes and, consequently, in SOFCs, ionic conduction losses are relatively higher.

4. Lumped parameter models

FC models must strike a balance between opposite requirements. On the one hand, they should be extremely rich in order to be able to represent the complete behavior of the cell by capturing the tridimensional distribution and time evolution of the physical quantities inside the cell. Such performance is achievable with a multiphysics, three-dimensional model, described by partial differential equations (PDEs) and characterized by a large number of physical parameters. On the other hand, they should be sufficiently simple to be included in other numerical procedure, such as a stochastic optimization loop, that is, the model should be numerically computable in a very short time, and they should be characterized by a relatively small number of parameters in order to avoid the curse of dimensionality issue. In such cases, a zero-dimensional stationary model is preferable since it avoids PDE numerical discretization and their inherent computational burden, allowing to run the algorithm on standard PCs.

4.1 Circuit models

The fuel-cell balance equations can be arranged so as to correspond to lumped equivalent circuits. Indeed, several of them can be identified, around the same basic concept, depending on the required level of accuracy and the behaviors to be highlighted.

In Figure 3, E is the open circuit voltage (the reversible voltage, i.e., the ideal electromotive force provided by the Nernst equation that depends on temperature and concentrations), R_ion is the electrolyte ionic resistance, R_a and R_c are the equivalent resistances of the electrode kinetics (at anode and cathode), equal to the activation overpotentials divided by the electrode-generated current. Depending on the analysis purpose, both the incremental and differential values can be used:

Figure 3.

Simple lumped equivalent circuit of a fuel cell.

$$R=\frac{{\eta }_c}{i\text{'}}\text{or}R=\frac{d{\eta }_c}{\mathit{di}\text{'}}$$E17

The latter resistance characterizes transient responses to small variation of the cell’s working point. Expanding it in the neighborhood of a working point for which the Tafel approximation holds provides:

$$R=\frac{d}{\mathit{di}}\left(-\frac{\mathit{RT}}{\alpha nF}ln{j}_0+\frac{\mathit{RT}}{\alpha nF}lni\right)=\frac{\mathit{RT}}{\alpha nFi}=\frac{\mathit{RT}}{\alpha nF{i}_0{e}^{\alpha nF\eta /\mathit{RT}}}$$E18

that is an adynamic circuit element, strongly dependent on the current or voltage. C_dl is the double layer capacity, related to TPB interfacial charge behavior, J_l are leakage current densities, due to various phenomena and particularly to the reactant crossover through the electrolyte, which make the output electric current smaller than the electrode-generated current (they are controlled source elements since such currents depend on the cell useful currents and on concentrations), Z_w is an element accounting for the losses due to concentration gradients. This is a dynamic element because it synthesizes the diffusion of reactant concentration at the electrodes, which are governed by the diffusion dynamics equations. In case of pulsating variations at an angular frequency ω around a steady-state working point, this is a Warburg nonlinear impedance that can be expressed as:

$$Z_w=\frac{\mathit{RT}}{A{\left(\mathit{nF}\right)}^2\overline{c}\sqrt{D}\sqrt{\omega }}tanh\left(\delta \sqrt{\frac{\mathit{j\omega }}{D}}\right)e^{j\frac{\pi }{4}}$$E19

where δ is the electrode thickness. By taking the limit ω → 0 it reduces to:

$$Z_w\left(0\right)=\frac{\mathit{RT}}{A{\left(\mathit{nF}\right)}^2\overline{c}}\frac{\delta }{D}$$E20

The dynamic elements C_dl and Z_w limit the response dynamics of the FC. In other words, an FC is not able to face arbitrary power variations and when these are too rapid it is necessary to support the FC with a storage device capable of faster response. This solution is adopted in automotive FCs, which are coupled with a lithium-ion battery or a supercapacitor, sized to meet the fastest transients and capable of reversible operations, enabling regenerative braking, but unable to accumulate as much energy as the H₂ tank of the FC.

A more sophisticated equivalent circuit of a DMFC is presented in Figure 4 in the case of a DMFC. Here, the dependences of the leakage currents on the reagent gradients and on the electrolyte currents have been separated and the concentration losses have been represented as voltage sources controlled by the methanol and oxygen gradients. Three separated equivalent circuits have been added, to represent the water, methanol, and oxygen behavior. Each of these includes capacitances, which account for the volumes of species accumulating at electrolyte, catalyst, and diffusion layers, whereas the controlled current sources describe drain and source at the catalyst layers and current-driven crossover.

Figure 4.

Multiphysics equivalent circuits capable of simulating electrical, chemical, and mass transport interactions (courtesy of IEEE Trans. on Industrial Electronics).

In any case, equivalent circuits provide an approximated description of the complex reactions and transport events occurring inside the FC. They have the merit to highlight separately single events (such as the ohmic losses in the electrolyte), thus allowing to easily study the effects of the variation of single quantities/parameters on the FC behavior. The main drawback of circuit models is that complex interactions and nonlinearity are not simulated in detail. More importantly, they can be implemented into circuit simulation software, to study the electrical interface of the FC with the power management electronics and system supervisor that is devoted to provide FC control, in order to study the overall system dynamics [16].

5. Distributed parameter models

The knowledge of how physical fields (electric field, current density, flow, velocity, temperature, species concentrations) are distributed within internal components constitutes a pivotal aspect in FC analysis and design, since gradients and irregularities hamper the achievement of optimal performance, but they can not be gripped by zero-dimensional, lumped models. On the other hand, distributed models have to cope with additional challenges deriving from the huge number of grid points needed for a complete tessellation of the multilayer 3-D domain, resulting in issues of “curse of dimensionality” which can hardly be faced without resorting to supercomputers. Parallel computing with domain decomposition can overcome this challenge if less powerful computers are used, by assigning one subdomain to each processor and implementing the few interactions between subdomains. The two electric potentials, for the electronic and electrolyte phases, are coupled by the surface overpotentials at the catalyst layers, where reaction kinetics is modeled by the Butler-Volmer equation. In the case of a one-dimensional formulation, Newton’s method is an efficient algorithm to integrate the equations, with LU factorization used at each interaction. Nevertheless, in the case of 2-D and 3-D formulations, the sparse Jacobian matrix produced by Newton’s method is too large to be efficiently handled. In this case, Gauss elimination with a generalized minimal residual subroutine (GMRES) preconditioned with a Gauss-Seidel block and a multigrid algorithm has proven to be more suitable to face the non-symmetric Jacobian matrix.

5.1 Proton exchange membrane fuel cells

Typical multiphysics coupled models include, among others, proton conduction, water and fuel transport, joule dissipation, and thermal diffusion. The models, typically discretized with the finite element method (FEM), pose significant numerical challenges. Some commercial simulation tools like COMSOL^® Multiphysics allow the solution of general time-dependent systems of partial differential equations [25] and are therefore very useful tools for this class of problems. For the computation of the fluid-dynamic field, particularly in the case of turbulent motion at high Reynolds numbers, the finite volume method is also used. Ansys Fluent^® is a commercial package based on this method particularly efficient in modeling fuel cells. PEM fuel cells, as the name implies, are based on proton-conducting polymeric membranes. The most commonly used material for their realization is persulfonated polytetrafluoroethylene, commercialized as Nafion^® by Chemours. This material has a structure similar to the one of PTFE, but is functionalized with sulfonic acid groups providing charge sites for proton conduction [7]. If the membrane is properly hydrated, protons can form hydronium complexes which once freed from sulfonic acid groups can move through the membrane. In these conditions, that is, proper hydration, proton conduction strongly depends on the water content and the temperature of the membrane, and can reach values as high as 20 S m⁻¹ at 100°C.

5.1.1 Electrical conductivity model

As briefly mentioned above, protonic conduction in Nafion^® strongly depends on the temperature, since the mechanism is based on charged particles jumping from site to site, with a rate described by the diffusivity D. This statistical parameter depends on the activation barrier energy which exhibits an exponential dependence on the temperature T according to the law:

$$D=D_o{e}^{-W_{\mathit{ai}}/kT}$$E21

where D_o is a diffusivity reference value, W_ai is activation barrier energy, and k is Boltzmann’s constant. The charged particle mobility μ is proportional to D according to the Einstein relation:

$$\mu =\frac{\left|z\right|FD}{\mathit{RT}}$$E22

with |z| the ion charge number, so that the proton conductivity σ = ρ_c μ can be written as:

$$\sigma =\frac{{\left(zF\right)}^{2}c{D}_o}{\mathit{RT}}{e}^{-W_{\mathit{ai}}/\mathit{kT}}={\sigma }_o{e}^{-W_{\mathit{ai}}/\mathit{kT}}$$E23

being ρ_c = |z|, Fc the charge density, and c the molar concentration.

Apart from the temperature, the proton conductivity also depends on the water content in the membrane. A common modeling approach to represent such dependence is based on the hydration λ, that is, the ratio between the number of water molecules and the number of charge sites available for proton conduction. In the specific case of Nafion^®, such ratio can be rewritten in a modified form using the water concentration c_w and the sulfonic acid concentration c_as, that is, λ = c_w/c_as. In absence of more sophisticated models, a linear dependence of conductivity on hydration can be assumed:

$$\sigma =\alpha \lambda =0.5139\lambda S/m$$E24

where λ is derived from a correlation empirically derived for Nafion^® [7]. Combining the above, the factorized expression of σ(λ, T) is obtained:

$$\sigma \left(\lambda ,T\right)=\alpha \lambda e^{\frac{W_{\mathit{ai}}}{k}\left(\frac{1}{303}-\frac{1}{T}\right)}$$E25

with W_ai/k = 1268 K for ions hopping. This is a more sophisticated model than the one expressed by Eq. (9). The conductivity σ influences the scalar potential φ according to the charge conservation equation in quasi-static conditions:

$$\nabla \cdot \sigma \left(\lambda ,T\right)\nabla \phi =0$$E26

which shows that the distribution of φ depends on λ and T.

5.1.2 Hydration model

A critical issue in modeling PEMFCs consists in providing an accurate description of the hydration effects, which rules proton conductivity [26]. The distribution of λ in the membrane can be computed resorting to specific equations at the surfaces and in the bulk. For the membrane bulk, two mechanisms are taken into account, namely electro-osmotic drag and back-diffusion, giving rise to the following expression of the water molar flow:

$${\mathit{N}}_w={\mathit{N}}_{\mathit{\text{we}}}+{\mathit{N}}_{\mathit{wd}}=\mathit{\xi \lambda }\frac{\mathit{J}}{\left|z\right|F}-D_w{c}_{\mathit{as}}\nabla \lambda$$E27

where N_i is the ionic molar flow vector and J the current density vector, and D_w is the water diffusivity in the membrane. This equation is nonlinear because D_w itself depends on λ and also on T. Such dependences can be expressed by factorizing the statistical mechanics exponential dependence on T with a polynomial regression obtained from experimental data:

$$D_w\left(\lambda ,T\right)={\sum }_{i=0}^3{d}_i{\lambda }^i{e}^{\frac{W_{\mathit{aw}}}{k}\left(\frac{1}{303}-\frac{1}{T}\right)}$$E28

where W_aw/k = 2416 K for water in Nafion^® and d₀ = 2.563671 × 10⁻⁶, d₁ = −0.33671 × 10⁻⁶, d₂ = 0.0264 × 10⁻⁶, and d₃ = 0.000671 × 10⁻⁶ for D_w (in cm² s⁻¹). The dynamics of N_w is ruled by Fick’s second law, which can be written in terms of λ:

$$\nabla \cdot {\mathit{N}}_w+c_{\mathit{as}}{\partial }_t\lambda =0$$E29

Letting Eq. (27) in Eq. (29) and assuming |z| = 1 for protons provide the following diffusion equation:

$$\nabla \cdot D_w\nabla \lambda -{\partial }_t\lambda -\nabla \cdot \frac{\xi \lambda \mathit{J}}{c_{\mathit{as}}F}=0$$E30

According to Maxwell’s equations, the current density J can be expressed in terms of $\phi$, so that Eq. (30) becomes:

$$\nabla \cdot D_w\left(\lambda ,T\right)\nabla \lambda -{\partial }_t\lambda -\nabla \cdot \frac{\xi \lambda }{c_{\mathit{as}}F}\sigma \left(\lambda ,T\right)\nabla \phi =0$$E31

Imposing interfaces conditions at the electro-catalyst layers is troublesome since the precise distribution of λ at such surfaces cannot be determined with sufficient accuracy. However, hydration can be expressed in terms of a more easily representable quantity, that is, the water vapor activity (that is, the relative humidity). The relationship between these two physical quantities can be modeled by resorting to an empirical relationship [27], which in turn can be expressed as a function of the water vapor pressure by another empirical model.

5.1.3 Thermal model

Since most of quantities appearing in the conductivity and hydration models depend on temperature, it is also necessary to take into account the transient heat equation:

$$\rho c_p{\partial }_{t}T-\nabla \cdot k\left(\lambda \right)\nabla T-\sigma \left(\lambda ,T\right){\left|\nabla \phi \right|}^2=0$$E32

where ρ is the hydrated Nafion^® density, c_p its specific heat, and k(λ) the thermal conductivity. According to data given in [28], k = 0.12 + 0.81λ (W m⁻¹ K⁻¹) can be assumed. The last term at the left-hand side represents Joule’s losses.

5.1.4 Coupled multiphysics model

The complete model to be solved is assembled from the above equations together with proper boundary (time-dependent Dirichlet and homogeneous Neumann) and initial conditions. The specific characteristics of the set of partial differential equations together with the high aspect ratio of the geometry (i.e., the very small thickness of the membrane compared to its extension in the plane) lead, after discretization with FEM, to a badly conditioned system of linear equations. The system tends to be quite large so that direct solvers may become inapplicable and it is also difficult to precondition so that iterative solvers tend to converge slowly or to fail altogether. Due to the strong nonlinearity of the complete coupled problem, a further numerical challenge concerns the nonlinear solver, typically the Newton-Raphson (NR) algorithm, which usually converges only, if at all, with strong under-relaxation.

The numerical solution of the final system may require substantial effort in correctly setting linear and nonlinear solver parameters to achieve convergence. Extreme care is needed in the choice of the drop tolerance if GMRES coupled to an ILU pre-conditioner with threshold is applied to solve the linear system arising in Newton’s method. Numerical experiments have shown that the selection of the drop tolerance needed to obtain convergence is highly problem-dependent so that the only way to reliably obtain a solution was to apply an efficient direct solver (PARDISO^®). This however limited the maximum admissible problem size on the available hardware. Furthermore, the fully coupled nonlinear system of equations can typically not be solved in its complete assembled form by the standard NR technique, even with strong under-relaxation. The problem was therefore solved with a further iterative loop, that is, a so-called “segregated solution,” in which one or more blocks of equations are fed into the next one in a simple iteration until convergence. In our specific case, the electrical model formed the first block, and its solution was inserted into the thermal model; then, the solution of both was used in the diffusion model. An important technological problem in the construction of PEM fuel cells is the reproducibility of the production process in the case of large membranes. In particular, variations in the thickness can result in hot spots which can cause aging effects and may lead to membrane failure. In order to investigate such effects, a lens-shaped compression of 2 mm radius and 50 μm depth (1/4 of the total PEM thickness of 200 μm) was studied. Results show an increase in current density in the order of 100%, leading to strong localized overheating (Figure 5).

Figure 5.

Local increase of current density due to a manufacturing defect (courtesy of IEEE Trans. on Magnetics).

5.2 Direct methanol fuel cells

Spatially resolved analyses of DMFCs allow studying their limitations, like sluggish kinetics and methanol crossover. These models resort to a detailed description of the real device geometry and materials, and of local nonlinear physics phenomena such as heat, momentum, multicomponent mass transport, and electrochemical processes. Early distributed parameter models were multi-domain ones, in the sense that the problem variables were defined in separated domains by introducing appropriate internal boundary conditions. As an example, in [12], the equations describing the concentration and potential distribution within the electrode were solved numerically using the finite difference method (FDM) and Newman’s BAND algorithm for the resulting simultaneous nonlinear equations. After the introduction of computational fluid dynamics (CFD) in the simulation of fuel cells, mostly single-domain models have been developed. The main advantages over the multi-domain approach is that internal BCs and continuity conditions at each domain interface are not required, thus simplifying the model geometry construction and speeding up problem set-up into a commercial CFD code. Single-domain CFD models can be classified into two-dimensional (2-D) and three-dimensional (3-D) models, depending on the simplifying assumptions on geometry and on boundary conditions. 2-D models generally provide a strong reduction in terms of computational cost, but their solution is less accurate compared to 3-D models.

5.2.1 DMFC two-dimensional models

A two-phase, multicomponent-flow 2-D model of a DMFC that accounted for capillary effects in both anode and cathode backings was developed [29]. In addition to electrochemical reactions, this model takes into account diffusion and convection of both gas and liquid phases in backing layers and flow channels. The effect of mixed potential related to methanol oxidation at the cathode, as a result of methanol crossover caused by diffusion, convection, and electro-osmosis, is simulated as well. Multiphysics equations are solved after discretization by the finite volume method (FVM). Numerical results concerning polarization curves are validated by experimental measurements. The main contribution of this work is the two-phase flow modeling of the anode, that is, the gas phase at the anode saturated with water and methanol and the liquid phase saturated with CO₂. Numerical analysis shows that gas-phase transport is one of the major issues affecting the fuel cell performance. In [30], a similar two-phase, two-dimensional model is presented. A capillary pressure function is used in order to simulate the methanol adsorption of backing materials. In addition, detailed multistep reaction models for both ORR and methanol oxidation as well as the Stefan-Maxwell formulation for gas diffusion are proposed. The effect of methanol and water crossover trough the membrane is accounted for in [31]. The two-phase mass transport in the anode and cathode porous regions is formulated based on the classical multiphase flow in porous media. A micro-agglomerate model, that is able to reflect the effect of the microstructure of the catalyst layer on cell performance, is proposed. The resulting polarization curves and methanol crossover rates at different concentrations are in very good agreement with experimental data. In [32], a realistic passive liquid-feed DMFC in transient charge/discharge conditions is simulated. The main contribution of this work is that effects of feed methanol concentration in the reservoir and current density on both mass transport and performance are investigated. Analyses show that when the initial feed concentration in the reservoir decreases, methanol crossover is minimized, but the fuel cell runtime is shortened. Recent works, for example [33], provide a detailed description of the whole DMFC system, including reservoirs. By developing a transient multiphase model of a passive cell, the effects of operating current density, voltage, micro-porous layer, and methanol feeding condition are comprehensively investigated for the whole operating processes, that is, with the fuel tank evolving from full to empty. Results highlight that for all operating conditions, it is necessary to operate at moderate current density or voltage to limit the methanol crossover and ensure the energy conversion efficiency. A 2-D multiphase non-isothermal mass transfer model for the DMFC is presented in [34]. The model includes the reaction of methanol and oxygen at the anode and cathode and the diffusion of every component involved in the DMFC, such as water, oxygen, and methanol at the diffusion layer and methanol crossover. It is shown that a maximum output power can be achieved for optimal temperature and concentration values.

5.2.2 DMFC three-dimensional models

A few papers are reported on 3-D two-phase DMFC models, which can capture the species distributions and the transport limitations along any direction inside the DMFC. Ref. [35] proposes a 3-D, two-phase, multicomponent model. Catalyst layers are incorporated in the computational domain instead of being modeled as zero-thickness interfaces. This model includes the effects of the second phase on the reduction of active catalyst surface areas and the mixed potential effects due to methanol crossover. The amount of carbon dioxide obtained from 3-D models indicates that the porosity of the anode diffusion layer plays a very important role in the DMFC performance. With a low porosity, the produced carbon dioxide cannot be removed effectively from the catalyst layer, thus reducing the active catalyst surface area as well as blocking methanol from reaching the reaction area. Ref. [31] extends the 2-D two-phase mass transport model for liquid-feed DMFCs to a fully 3-D model. The two-phase mass transport in the anode and cathode porous regions is formulated based on the multiphase flow theory in porous media without defining the mixture pressure of gas and liquid and assuming a constant gas pressure in the porous regions. The interaction between the phases in this 3-D model is captured by taking into account the effect of non-equilibrium evaporation/condensation at the phase interface, as opposed to the assumption of other models of thermodynamic equilibrium condition between the phases. From 3-D analysis, it can be observed that methanol concentration in the diffusion layer is higher in the channel than under the ribs, demonstrating that the flow-field channels cause methanol to be distributed unevenly over the reaction area (Figure 6).

Figure 6.

Methanol concentration in the anode and membrane (unit: m³/s) [30] (courtesy of Electrochimica Acta).

In [36], a commercial flow solver (i.e., Ansys Fluent^®) is used to solve at the same time flow, species, and charge transport equations. 3-D simulations are carried out in order to explore mass transport phenomena occurring in DMFCs for portable applications as well as to reveal an interplay between the local current density and methanol crossover rate. In [37], 3-D modeling is then extended to transient conditions. The authors note that cathode processes, for example oxygen and water transport coupled to electrochemical reaction, are inherently transient so that an unsteady-state model gives more accurate prediction than a steady-state model. Numerical simulations indicate that the cathode catalyst layer porosity has major effects on oxygen transfer and water removal. A three-dimensional multiphase model of DMFC is developed in [38], in which the Eulerian-Eulerian model is adopted to treat the gas and liquid two-phase flow in channel. By 3-D simulation, cell performance is found to be severely affected by accumulation of carbon dioxide mainly at the anode channel and by high-temperature operations. Ref. [39] shows that three-dimensional models are suitable for analyzing DMFC stacks with flowing electrolyte. A multiscale approach is therefore proposed in order the reduce the computational cost arising from 3-D modeling of the entire stack geometry. By this solution strategy, fully 3-D flow fields, backing layers, and membranes are numerically solved, whereas electrochemical reactions are analytically simulated.

Figure 7.

Identification by means of a hybrid PSO-DE algorithm run simultaneously on all polarization curves. Dashed lines: experimental, continuous lines: optimized models (courtesy of Journal of Power Sources).

It should be finally noted that multiphysics models coupling electrochemical reactions, methanol, water, and heat transport are still under investigation due to their high complexity.

5.3 Solid oxide fuel cells

Three different designs are used in planar SOFCs depending on their operating temperatures. High-temperature SOFCs (around 1000°C) usually present an electrolyte-supported structure, with thin electrodes (e.g., 50 μm) supported by a thick electrolyte (above 100 μm) [1]. The high temperature ensures so high a conductivity that the electrolyte resistance remains within acceptable values. In low-temperature SOFCs (though not less than 600°C), thinner electrolytes are used (e.g., 20 μm) and the cell is supported by either anode or cathode (300–1500 μm) with the other electrode being thinner (e.g., 50 μm). Understanding the multiphysics behavior is indispensable in identifying optimal design and operation of such SOFCs and a multiphysics numerical model is required at this purpose. Both Fluent and STAR-CD have been applied with success to this task, providing precious information on the internal distribution of reactant and product, current density, temperature, and stress and, more generally, on the detailed operation of a SOFC. Models can include different fuels such as H₂ and CO and can take into account internal reforming by means of a catalytic chemical reaction [40]. Moreover, transport phenomena are not as complex as in a PEFC and DMFC.

Due to the high working temperatures, their gradients contribute to stress formation, which is a major technical issue of SOFC. Consequently, early modeling studies were aimed at predicting the current and temperature distributions, whereas flow and multicomponent transport were typically simplified. The subsequent use of CFD models has allowed more detailed three-dimensional multiphysics analyses. A Fluent-based CFD model has been developed by [41] to describe reactant flow, transport, and electrochemical reaction in a SOFC. STAR-CD was combined to an electrochemistry module by [42] to simulate a SOFC. Results showed that the co-flow reactant distribution at anode and cathode had the most uniform temperature distribution and the smallest thermal gradients. As a drawback, this approach treated the electrodes-electrolyte as a solid component, neglecting mass diffusion that is important in the case of thick electrodes. The effect of mass transport in a thick electrode has been analyzed by [43] in a two-dimensional study that couples mass and heat transport and included methane/steam reforming by means of a catalytic chemical reaction. Basically, in three-dimensional analysis, the governing equations for conservation of mass, momentum, species, thermal energy, electric charge, and electrochemical kinetics in anode and cathode of a SOFC are the same as those of a PEMFC, except that water transport through the electrolyte is not considered [44]. Orthogonal meshes are typically used in modeling planar SOFCs. A five-channel geometry can be modeled with 80 × 80 × 35 mesh, resulting in 224,000 cells, on which the model can converge in some 300 iterations. Simulations of this kind reveal how current density increases in correspondence of the channels in the current collectors, where the electrodes receive more reactants. The effect of the flow channel design on the cell performance can be analyzed in detail with such numerical tools. Due to high operating temperature, experimental validation of the numerical model is particularly challenging and few works are reported in the literature on this topic.

6. Optimization and identification

6.1 Optimization

A number of nonlinear deterministic optimization methods (DOMs) have been applied to PEMFCs in the last decade, proving successful in dealing with specific tasks. Least squares methods have been applied to the estimation of single material parameters as well as parameters evolution under degradation events [2]. The gradient method has been exploited in the search for optimal designs and parameters evolution, such as cathode configuration optimization, geometric optimization, and flow field serpentine optimization [45]. A review of deterministic optimization methods used for identification problems in PEMFCs is given in [46].

Deterministic methods are known for their efficiency, that is, speed of convergence, but their applicability may be hindered, depending on the specific algorithm, by lack of flexibility in handling arbitrary constraints, sensitivity to noise in the objective function, possible need of function derivatives, and premature convergence to local minima. On the other hand, stochastic optimization methods (SOMs), in spite of their comparatively low efficiency, typically overcome the above-mentioned shortcomings of deterministic methods. Although the convergence to the global optimum for SOMs is only asymptotically guaranteed, there is abundant numerical evidence that very good solutions can be obtained for many problems, including FC design problems. It is worth noticing that a crucial feature of FC circuit models is that they avoid partial differential equations, thus resulting in numerical formulations with relatively low computational costs, which make them ideally suited for SOMs. It should also be noted that since optimization problems related to PEMFCs are characterized by highly nonlinear device models, the resulting objective functions subject to minimization have many local minima, which need to be all identified in the search for the global one. Therefore, for this type of problems, stochastic optimizers may end up being almost as efficient as deterministic ones. A further advantage, which however is also shared by some deterministic methods, is that stochastic optimizers can also deal with non-differentiable or fully discrete, optimization problems. In recent years, the application of stochastic methods for the solution of FC optimization problems has been constantly increasing, and interesting results have been reported, for example, with genetic algorithms (GA) [47], particle swarm optimization (PSO) [48], and differential evolution (DE) [49]. The above-mentioned methods are all population-based, that is, they explore several candidate solutions concurrently, which makes them ideally suited for parallelization. SOMs can also be combined among them or hybridized with DOMs in order to tailor their behavior to the specific optimization problem. Multiobjective stochastic approaches have also been recently investigated [50].

6.2 Identification

Since fuel cells present a stratified structure of thin layers made of different materials, analyzing their behavior requires the full characterization of these materials, that is, the determination of their chemical, physical, thermal, and electrical parameters. The identification of these parameters is crucial for guiding the research for advanced functionalized materials. These parameters are also needed in FC models, used in the fast exploration of different operating scenarios and in the research of optimized structural design and operating conditions [51]. The systems of equations involved (Nernst equation, Butler-Volmer equation, Darcy’s equation, Fourier’s law, Ohm’s law, etc.) are strongly nonlinear, making the models extremely sensitive to parameter variations and uncertainties.

Unfortunately, they are hard to measure in real operating conditions and their identification still remains challenging when dealing with direct measurements. Careful ex situ measurements can be performed by means of a number of diagnostic techniques; however, the transferability of results to operative fuel cells raises a number of issues. Conversely, in situ measurements can provide meaningful operational values, but very few, often complicated and cumbersome, techniques are available to determine a limited number of parameters. A different approach consists in tackling the identification of multiple parameters by using a very large body of experimental data collected at different of physical conditions (e.g., temperature, pressure, and humidity). However, this approach suffers from the well-known “curse of dimensionality,” that is, the problem becomes exponentially harder to solve as the number of parameters increases. This challenging problem can be approached with an optimization procedure (i.e., the search of the minimum of a function ƒ(x)). When using optimization algorithms for model parameter identification, x is the n-dimensional vector of the unknown parameters to be identified and ƒ(x) consists in an error function that assesses the difference between the output of the parameter-based model and the measurements. The parameter identification problem is a constrained one, that is, the domain A where x values are defined is supplemented with a number of constraints and the problem is also typically burdened by model nonlinearity, as is the case of an FC model, which results in the non-convexity of ƒ and consequent local minima. Moreover, large problems lead to high computational cost. Given the problems of “curse of dimensionality,” presence of local minima, and computational costs, smart strategies are needed to find good solutions, if not the absolute best one, which actually may be impossible to identify.

In the last decade, stochastic methods have been applied to the study of FC parameter identification problems and their use has been strongly increasing in the last 3 years. Studies reported in literature typically aim at using stochastic methods in order to obtain a good fit of PEMFC polarization curves and usually resort to simplified empirical PEMFC models. Such models tend to be interpolatory in nature and contain enough parameters (e.g., 5–7) to allow for a good fit. Given their nature, fitting a set of empirical parameters to match a given polarization curve is a relatively easy task for most optimization procedures, but the usefulness of the obtained results is rather limited. The relatively small number of parameters also helps in avoiding so-called duplicity problems, that is, multiple distinct solutions achieving the same minimal values of ƒ. However, duality is not crucial, because empirical parameters have no direct physical meaning. A more ambitious challenge consists in identifying several physical parameters of the materials of a PEMFC by means of an optimization approach. An algorithm of this kind, built over an early investigation on the capability of stochastic methods to deal with FCs [52], uses a detailed multiphysical performance model that employs such parameters and takes into account some physical control quantities [53]. A straightforward use of a stochastic optimizer with a large number of unknowns and a weakly constrained nonlinear problem can result in duplicity. In order to overcome it, a possible strategy consists in splitting the overall identification problem into a sequence of distinct identification sub-problems, each having a lower number of unknowns, and thus suffering less from the “curse of dimensionality.” This approach basically relies on isolating a group of equations, dominated by some of the unknowns only, and has already been applied successfully to fuel-cell problems [54]. This behavior emerges, for example, since the parameters related to the activation losses and to the concentration losses, which are strongly nonlinear, prevail over the ohmic losses, which are linear at given hydration and temperature, and their identified values tend to vanish [55]. This behavior can be exploited by separately considering the typical three parts of the polarization curves. The experimental data obtained at low current density can be used to identify the parameters related with activation losses, while those obtained at high current density can been used to identify the parameters related with the concentration losses. Finally, experimental data at intermediate current density values can be used to identify the parameters related with the ohmic losses, which dominate in the central part of the polarization curve (Figure 7).

Some final considerations emerge: first of all, the accuracy and reproducibility the experimental data and the experimental conditions must increase with growing number of unknown parameters. This may be hard to obtain in the case of the polarization curves of PEMFCs, which depend on many factors related to both the samples under test and the experimental conditions, some of which are hard to control. Such difficulties can be mitigated by collecting more curves in the same nominal operating conditions and performing a statistical selection of the data. Moreover, enriched experimental plans may allow to identify also some parameters that have a minor effect on the polarization curve, for example, a set of experimental data obtained at different back pressures may allow the identification of the anodic exchange current density, which is otherwise masked by the larger effect of the cathodic one. On the other hand, a failure of the identification procedure, that is, a poor fit of the polarization curves with varying operating conditions, typically hints at weaknesses in the model.

7. Conclusion

The decarbonization of power grids and the development of electric vehicles and renewable resources are promoting researches on advanced fuel cells for both stationary and mobile applications. Several programs are running that aim at reducing costs and improve performance and duration. These researches can be strongly supported by numerical models which are the core of a computer-aided engineering approach capable of reducing tentative experimentation by selecting those solutions which result more competitive on the basis of analytical/numerical computations, both in terms of device design and its operation. Effective models must take into account all relevant chemical-physical-electrical quantities, their interdependence, and their evolutions. These models play an interactive game with diagnostic and measurement issues, because their reliability depends on their fitting to the cases under investigation and this occurrence involves the determination of the many parameters used in nonlinear equations. Design optimizations based on numerical procedures are strong tools not only in identifying device materials and geometries capable of competitive performance, but also in determining the whole FC system, including static converters and system supervisors.

Fuel cell computational modeling confirms to be an important topic of applied research, involving multiphysics, multiscale problems which are still challenging for the researchers working on this subject, both at the scientific and industrial level.

Acknowledgments

This work was supported by MAESTRA 2011—“From Materials for Membrane-Electrode Assemblies to Electric Energy Conversion and Storage Devices,” the 2011 University of Padua Strategic Project cod. STPD11XNRY_002.

Acronyms

A

area (m²)

a/c

anode/cathode (−)

C_dl

double-layer capacitance (F)

c_i

molar concentration (mol m⁻³)

c_p

specific heat capacity (J kg⁻¹ K⁻¹)

D

diffusion coefficient (m² s⁻¹)

d_m

membrane thickness (m)

E₀

standard potential (V)

F

Faraday constant (A s mol⁻¹)

j

electric current density (A m⁻³)

J

electric current density vector (A m⁻²)

k

Boltzmann constant (J K⁻¹)

k_T

global thermal exchange coefficient (K W⁻¹)

n

number of electrons in reaction (−)

N

molar flow density (mol m⁻² s⁻¹)

P_diss

dissipated power (W)

p_i

partial pressure (Pa)

R

universal gas constant (J mol⁻¹ K⁻¹)

T

temperature (K)

V

cell voltage (V)

W_a

activation energy (eV)

Z_w

Warburg impedance (Ω)

α

charge transfer coefficient (−)

δ

electrode thickness (m)

$\mathrm{\Delta }\widehat{s}$

molar entropic variation (J K⁻¹)

η

overpotential (V)

κ

mass transfer coefficient (_)

λ

Nafion^® hydration (−)

μ

charged particle mobility (m² V⁻¹ s⁻¹)

ξ

electro-osmotic drag coefficient (_)

ρ

mass density (kg m⁻³)

σ

electric conductivity (S m⁻²)

φ

electric potential (V)

ω

angular frequency (rad s⁻¹)