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.
- fuel cells
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/cm2 (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−2 cm. These GDLs are separated by proton-exchange membranes (PEMs) which are thin electrolytes, with widths in the order of 1–5 × 10−3 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 103, involves severe size problems: a regular hexahedral tessellation with 10 elements in the thickness direction implies 109 nodes, namely a very demanding computational dimension . 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 , 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 . 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
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
are segregated by the proton exchange membrane (Figure 1). According to the Nernst equation, the cell’s reversible voltage
The cell open circuit voltage (OCV) is slightly different from
The first subscript of the overpotentials indicates
3.1.1 Exchange current density: activation losses
The activation overpotentials or activation voltage drops, Δ
By using Eq. (6), the Butler-Volmer equation allows to write the current density at the cross-sectional area of each CL as
In this equation, the total equivalent current density
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
3.1.3 Concentration losses and limit conditions
Gas concentration and pressure gradients
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
The linear dependence on
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:
where is the hydrogen concentration at the anodic CL,
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:
Heat transport inside the cell depends on conduction, diffusion, and convection and interacts with thermal capacity in dynamic conditions . An accurate enough estimation of the mean temperature
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.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.
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:
where is the (constant) standard cell potential,
where and are the reference current density and the reference concentration at the anode, respectively. The voltage overpotential 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 , 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 , 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 .
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 . 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 . 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 , 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  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 , 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  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 . In , 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  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−
A typical conductivity value is
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
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,
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:
that is an adynamic circuit element, strongly dependent on the current or voltage.
The dynamic elements
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.
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 .
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  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 . 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−1 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
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
which shows that the distribution of
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 . The distribution of
According to Maxwell’s equations, the current density
Imposing interfaces conditions at the electro-catalyst layers is troublesome since the precise distribution of
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:
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).
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 , 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 . 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 CO2. Numerical analysis shows that gas-phase transport is one of the major issues affecting the fuel cell performance. In , 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 . 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 , 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 , 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 . 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.  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.  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).
In , 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 , 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 , 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.  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.
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) . 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 H2 and CO and can take into account internal reforming by means of a catalytic chemical reaction . 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  to describe reactant flow, transport, and electrochemical reaction in a SOFC. STAR-CD was combined to an electrochemistry module by  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  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 . 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
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 . 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 . A review of deterministic optimization methods used for identification problems in PEMFCs is given in .
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) , particle swarm optimization (PSO) , and differential evolution (DE) . 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 .
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 . 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 ƒ(
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 , uses a detailed multiphysical performance model that employs such parameters and takes into account some physical control quantities . 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 . 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 . 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.
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.
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.
double-layer capacitance (F)
molar concentration (mol m−3)
specific heat capacity (J kg−1 K−1)
diffusion coefficient (m2 s−1)
membrane thickness (m)
standard potential (V)
Faraday constant (A s mol−1)
electric current density (A m−3)
electric current density vector (A m−2)
Boltzmann constant (J K−1)
global thermal exchange coefficient (K W−1)
number of electrons in reaction (−)
molar flow density (mol m−2 s−1)
dissipated power (W)
partial pressure (Pa)
universal gas constant (J mol−1 K−1)
cell voltage (V)
activation energy (eV)
Warburg impedance (Ω)
charge transfer coefficient (−)
electrode thickness (m)
molar entropic variation (J K−1)
mass transfer coefficient (_)
Nafion® hydration (−)
charged particle mobility (m2 V−1 s−1)
electro-osmotic drag coefficient (_)
mass density (kg m−3)
electric conductivity (S m−2)
electric potential (V)
angular frequency (rad s−1)