Modeling approaches for fuel cells

Modeling approaches for fuel cells

CHAPTER Modeling approaches for fuel cells 10 Modeling is of key importance in fuel cell development beyond the current-state-ofthe-art because it ...

2MB Sizes 0 Downloads 47 Views


Modeling approaches for fuel cells


Modeling is of key importance in fuel cell development beyond the current-state-ofthe-art because it is beneficial to understand the mechanisms of various interacting phenomena and effects on the cell performance. It is also hard to measure the local parameters inside fuel cells, particularly inside the small scale functional materials. Modeling and simulation of charge transfer and electrochemical performance are critical to enable optimization of the geometry and performance. For macroscale modeling, the micro- and nano-structure related properties defining a porous media, e.g., the porosity and tortuosity and the specific area available for surface reactions are required. Thus the coupling of models valid at various scales is important for continued progress in the development of fuel cells.

10.1 Introduction In fuel cells various reacting transport phenomena occur at different scales. In the cells and manifold multicomponent and multiphase flow of reactants and products take place and transfer of charges as ions and electrons as well as heat and mass transfer appear at various functional components and sites. Chemical and electrochemical reactions at nano- or microstructured electrodes and electrolytes affect these processes significantly. The electrochemical reactions generate or consume chemical species and produce an electric current and these reactions take place at active sites called triple phase boundaries (TPBs). The reforming reactions of hydrocarbon fuels in solid oxide fuel cells and water phase change in proton exchange membrane fuel cells are coupled with the electrochemical reactions and the transport processes. At the unit cell and component level the various phenomena might be analyzed by methods based on the continuum approach, e.g., CFD (computational fluid dynamics) methods. To analyze and understand the processes in catalyst layers and the active surfaces, microscale models are required. Such models can be based on the density functional theory (DFT) and molecular dynamics (MD) approaches at the microscale and at the mesoscale level by Monte Carlo (MC) methods and Lattice Boltzmann methods (LBM). By such methods the importance of the microscopic structure of the functional materials on the electrochemical reactions at the active sites, the surface chemistry and gas phase chemistry based on elementary reaction kinetics in the porous electrodes can be analyzed. In this chapter, macroscale and micro-/nanoscale modeling approaches are presented for the fuel cell types SOFC Hydrogen, Batteries and Fuel Cells. Copyright © 2019 Elsevier Inc. All rights reserved.



CHAPTER 10 Modeling approaches for fuel cells

(solid oxide fuel cell) and PEMFC (polymer electrolyte membrane fuel cell). In particular the electrodes, electrolyte and unit cells are considered. The macroscopic fuel cell models developed over several decades can be categorized as analytical, semi-empirical or mechanistic. Usually many simplifying assumptions are applied in an analytical approach concerning variable profiles within the cell, which is useful if quick calculations are required for fuel cell systems and does not give an accurate picture of transport processes occurring within the fuel cell components/cells. On the other hand, based on experimental data specific to the applications and operating conditions, semi-empirical fuel cell models combine theoretically derived differential or algebraic equations with empirically determined relationships, when the physical phenomena are difficult to model, or when the theory governing the phenomenon is not well understood. Based on detailed electrochemical, fluid dynamics, species/charge transport and heat transfer relationships, a theoretical fuel cell modeling approach usually employs the basic equations (differential and algebraic ones), to be solved by using some numerical computation method for the strongly coupled electrochemical and transport processes. For proper water and thermal management, this approach includes not only the electrochemical reactions but also thermal- and fluid-dynamic equations. The output of the modeling can provide details of the processes, e.g., fuel cell species distribution and flow pattern, current density and temperature distribution, voltage and pressure drop, etc. The theoretical models can be further divided into multi- or single-domain ones and the equations might be solved separately or simultaneously, see Ref. [1]. On the other hand, the microscopic approaches (e.g., DFT and MD) and the mesoscale ones (e.g., LBM) are more related to deep theoretical knowledge compared to the global models. The detailed chemistry and surface reaction models are able to take into account the effects of the multi-functional materials microscopic structures on the charge-transfer reactions taking place at active sites, the surface chemistry and the gas-phase chemistry based on elementary reaction kinetics (individual chemical reaction steps between intermediates). Such reaction schemes usually consist of more than 10 surface-phase species and around 40 irreversible reactions. Implementation of such large reaction mechanisms in fuel cell design and performing parametric studies using CFD approaches are CPU demanding tasks. More discussion on the microscopic fuel cell modeling development and multi-scale integration technique can be found in, e.g., Refs. [2,3]. Due to limitations of experimental studies at extreme small scales, mathematical modeling is widely used to investigate transport and electrochemical phenomena in CLs. Atomic-scale models (e.g., molecular dynamics, MD) can predict physical properties of materials with ideal, theoretical or proposed microstructures under well defined conditions within the size of several nanometers. However, these are not able to probe the random morphology of the complete CLs due to computational limitations. The approaches beyond the atomic-scale and below the conventional continuum scale use coarse-grained pseudo-particles which can either move on a fixed lattice (lattice-based pseudo-particle models) or continuously in space (off-lattice pseudo-particle models). The lattice Boltzmann method (LBM) has

10.2 Zero-order models of analysis

FIG. 10.1 Approximative ranges for the characteristic time and length scales for various methods of analysis and indication of the need of computational resources. Based on Ref. [16].

been developed to describe the diffusion process in complex pore structures and to study the structure-wettability influence on the underlying two-phase dynamics. The off-lattice pseudo-particle models, e.g., dissipative particle dynamics (DPD), smoothed particle hydrodynamics (SPH) and coarse-grained molecular dynamics (CG-MD), have been applied for exploring the microstructure of CLs. The DPD can be regarded as coarse-grained non-equilibrium molecular dynamics method. It is of the order of 105 times faster than MD, but this substantial increase in computational efficiency comes at the expense of loss of details at the atomic/molecular scale. More systematic and reliable ways of coarse graining molecular dynamics to construct DPD models would be an important advance. The SPH was introduced some 30 years ago to simulate astrophysical fluid dynamics, and now can be used to simulate complex phenomena associated with contact line and contact angle dynamics. However, the full range of scales (from the atomistic to the pore scale) cannot be investigated in a single simulation. The deficiency may be mitigated by simulations based on hybrid multiscale particle-based models. In the CG-MD method, super-atoms are used to represent groups of atoms, and interactions are only defined between these super-atoms. Thus, a reduction of the number of degrees of freedom is achieved. Further details can be found in Refs. [4e15]. Fig. 10.1 illustrates the characteristic time and length scales for the various methods discussed above.

10.2 Zero-order models of analysis In so-called zero-order models of analysis, the transport equations for gas concentration and temperature are not solved but instead voltage-current density



CHAPTER 10 Modeling approaches for fuel cells

polarization curves are considered to provide the performance characteristics. As already presented in Chapters 2 and 9 the voltage is written in terms of the reversible open circuit voltage and various components of the fuel cell overpotentials, i.e., V ¼ E0  hactivation  hohmic  hconcentration


where E0 is the open circuit voltage given by Nernst equation, here re-given for a hydrogen-oxygen fuel cell, i.e., E0 ¼ Ereversible ðTÞ 

pH 2 O RT ln 2F pH2 p0:5 O



The activation loss hactivation is given by the Butler-Volmer equation or the Tafel equation as presented in Chapters 2 and 9. The ohmic voltage loss hohmic occurs due to the resistance of the material to charge transport. The charge transport is composed of the electron transport from the anode to the cathode through the electrodes and interconnect, and the ion transport through the electrolyte. Expressions for the ohmic voltage loss were given in Chapters 2 and 9. The concentration voltage loss hconcentration takes place because the concentrations of the reactants and products at the reaction sites vary. This loss is composed of two components, namely due to the reaction kinetics and concentration effects as in the Nernst equation. Further details can be found in Chapters 2 and 9.

10.3 One-dimensional models of analysis In one-dimensional models the reactant gas variations across the cell, i.e., from the anode and cathode gas channels to the respective electrode-electrolyte interface, are considered. Only variations in one direction are considered. In a hydrogen-oxygen fuel cell, the concentrations of the reactants vary from the bulk values in the gas flow channels due to some resistances to the mass transport. These resistances are (a) convective resistance between the gas stream and the electrode surface, (b) convective and diffusive resistance in the porous electrodes. The heat generated by the reaction kinetics can often be considered as a surface phenomenon while the ohmic heating is a volumetric heat source. Inside the cell, the heat is primarily transported by conduction in the electrolyte, conduction and convection in the porous electrodes, and convection in the gas flow channels. In Chapter 9, such a one-dimensional model for the heat transfer analysis was presented. Water transport through a fuel cell is complex in particular for PEMFCs. For PEMFCs the water is generated at the electrolyte-cathode interface while for SOFCs the water is generated at the anode-electrolyte interface. The water is transported primarily by diffusion through the respective porous electrode layers and by convection in the gas flow channels. For the SOFC electrolytes eventual water transport is by diffusion while for the common Nafion polymer electrolyte of PEMFCs, the water transport is more complex. More details can be found in Ref. [17]. In one-dimensional models for the electrochemical processes, the voltagecurrent density polarization is applied in a similar manner as for zero-order models.

10.4 Multi-dimensional models of analysis

10.4 Multi-dimensional models of analysis Zero-, one- and two-dimensional models for analysis of the fundamental processes occurring in fuel cells are relatively simple and not so demanding on computational resources. Nevertheless, for estimations and qualitative analyses of the performance under various operating conditions such models might be appropriate. On the other hand, as a more accurate, more detailed and quantitative analysis is needed for development and design, complete three-dimensional models are needed. Such models then must include the transport processes of flow, mass, heat and charge in three directions with sufficient resolution of the electrode and electrolyte layers and the flow channels of fuel and oxidant. The purpose of this section is to present some general features of multidimensional models and then present examples to illustrate applications for PEMFCs and SOFCs.

10.4.1 Computational fluid dynamics (CFD) approaches Computational methods started to have a significant impact on analysis of heat transfer and thermal design in the late 1960s. Computational fluid dynamics (CFD) is an interdisciplinary branch of science and engineering with a broad spectrum of applications. Fluid flow, heat transfer, mass transfer, combustion and chemical reactions appear in most aspects of modern life and are of significance in automotive, space and aviation, chemical and process industries as well as in atmospheric science, energy, medicine, micro-and nanotechnology. The development and applications of CFD have been tremendous and it is used as a modeling tool as well as in R & D in many industries nowadays. Besides it continues to be developed for new challenges and being used in basic research at universities. Governing equations The governing partial differential equations of mass conservation, transport of momentum, energy and mass fraction of species etc. can be cast into a general partial differential equation as Ref. [18,19] ! vðrfÞ v v vf þ ðrfuj Þ ¼ G þS (10.3) vt vxj vxj vxj where f is an arbitrary dependent variable, e.g., the velocity components, temperature or species concentration. uj is the velocity vector, G the generalized diffusion coefficient, and S the source term for f. The general differential equation consists of four terms. From the left to the right in Eq. (10.3), they are referred to as the unsteady term, the convection term, the diffusion term and the source term. Numerical solution of the governing equations There are many methods established for numerical solution of the governing equations of fluid flow and heat transfer problems. These are: the finite difference method (FDM) [20], the finite volume method (FVM) [18,19], the finite element



CHAPTER 10 Modeling approaches for fuel cells

method (FEM) [21,22], the control volume finite element method (CVFEM) [23] and the boundary element method (BEM) [24]. In this chapter, only the finite volume method is considered. The finite volume method (FVM) In this method the domain is subdivided into a number of so-called control volumes. The conservation equations are integrated across each control volume. At the center of the control volume a node point is placed and the variables are located at this node. The values of the variables at the control volume faces are obtained by interpolation. The evaluation of the surface and volume integrals is carried out by quadrature formulas. The FVM is very suitable for complex geometries and the method is conservative as long as surface integrals are the same for control volumes sharing boundary. It is frequently used for convective flow and heat transfer. It is also applied in several commercial CFD-codes. Further details can be found in Refs. [18,19]. A brief illustration is presented below. Integration of the general equation across the control volume gives " # ZZZ ZZZ ZZZ vruj f v vf dV ¼ Sf dV (10.4) Gf dV þ vxj V vxj V vxj V Then by applying the Gaussian theorem one has ZZ ZZ ZZZ vf rfuj nj dS ¼ Gf nj dS þ vxj S S

Sf dV



By summing over all the faces of the control volume, this equation is transferred to nf X f¼1

ff C f ¼

nf X

Df þ Sf dV



where the convection flux Cf, the diffusion flux Df and the scalar value of the arbitrary variable f at a face, Ff, have to be determined. Convection and diffusion fluxes To achieve physically realistic results and stable iterative solutions, the way to handle convection-diffusion terms needs to possess properties of conservativeness, boundedness and transportiveness. The upstream, hybrid and power-law discretization schemes all possess these properties and are generally found to be stable but they suffer from numerical or false diffusion in multidimensional flows if the velocity vector is not parallel with one of the co-ordinate directions. The central difference scheme lacks transportiveness and is known to give unrealistic solutions at large Peclet numbers. Higher order schemes like QUICK, van Leer etc. may minimize the numerical diffusion but are less numerically stable. Also implementation of boundary conditions can be somewhat problematic with higher order schemes and the computational demand can be extensive as also additional grid points are needed and the expressions for the coefficients in Eq. (10.6) become more complex.

10.4 Multi-dimensional models of analysis Source term The source term S may in general depend on the variable f and also have a complex mathematical form. In the discretized equation it is desirable to account for such a dependence. Commonly the source term is expressed as a linear function of f. Solution of the discretized equations The discretized equations have the form of Eq. (8) with the f - values at the grid points as unknowns. For boundaries not having fixed f - values, the boundary values can be eliminated by using given or fixed conditions of the fluxes at such boundaries. Gauss elimination is a so-called direct method to solve algebraic equations. For one-dimensional cases the coefficients form a tridiagonal matrix and an efficient algorithm called the Thomas algorithm or the tri-diagonal matrix algorithm (TDMA) is achieved. For two-dimensional and three-dimensional cases, direct methods require large computer memory and computer time. Iterative methods are therefore used to solve the algebraic equations. A popular method is a line-by-line technique combined with a block correction procedure. The equations along the chosen line are solved by the TDMA. Iterative methods are also needed because the equations are non-linear and sometimes interlinked. In many situations, e.g., turbulent forced convection, the change in the value of f from one iteration to another is so high that convergence in the iterative process is not achieved. To circumvent this and to reduce the magnitude of the changes, under-relaxation factors (between 0 and 1) are introduced. Handling pressure in the momentum equations The pressure gradients appear as source terms in the momentum equations and are not known but have to be found as part of the solution. Thus the pressure and velocity fields are coupled and the continuity equation (mass conservation equation) has to be used to develop a strategy. It has been shown that if the velocity components and the pressure are calculated at the same grid points in a straightforward way, some physically unrealistic fields, like checker-board solutions, may arise in the numerical solution. A remedy to this problem is to use staggered grids. The velocity components are then given staggered or displaced locations. These locations are such that they lie on the control volume faces that are perpendicular to them. All other variables are calculated at the ordinary grid points. Another remedy is to use a non-staggered or collocated grid where all variables are stored at the ordinary grid points. A special interpolation scheme is then applied to calculate the velocities at the control volume faces. Most commonly the so-called Rhie-Chow interpolation method is applied, see Ref. [25]. Solution procedures for the momentum equations The oldest algorithm to handle the pressure-velocity coupling is the SIMPLE (semi-implicit-method-pressure-linked-equations) algorithm. A pressure field is guessed and then the momentum equations are solved for this pressure field resulting in a velocity field. Then a pressure correction and velocity corrections are introduced. From the continuity equation an algebraic equation for the pressure correction



CHAPTER 10 Modeling approaches for fuel cells

can be obtained. The velocity corrections are related to the pressure corrections and the coefficients linking the velocity corrections to the pressure correction depend on the chosen algorithm. There are other similar algorithms available today. SIMPLEC (SIMPLEconsistent) and SIMPLEX (SIMPLE-extended) are common. They differ from SIMPLE mainly in the expression for coefficients linking velocity corrections to the pressure correction, see Jang et al. [26]. Another algorithm named PISO (pressure implicit splitting operators), see Issa [27], has become popular more recently. Originally it is a pressure-velocity coupling strategy for unsteady compressible flow. Compared to SIMPLE it involves one predictor step and two corrector steps. Still another algorithm is SIMPLER (SIMPLE e revised). Here the continuity equation is used to derive a discretized equation for the pressure. The pressure correction is then only used to update the velocities through the velocity corrections. Convergence The solution procedure is in general iterative and then some criterion must be used to decide if a solution is converged or not. One method is to calculate residuals R as X R¼ aNB þ b  aP fP (10.7) NB

for all variables. NB means neighboring grid points, e.g., E, W, N, S (east, west, north, south). If the solution is converged, R ¼ 0 everywhere. Practically, it is often stated that the largest value of the residuals R should be less than a certain number. If this is achieved the solution is said to be converged. Number of grid points and control volumes The sizes of control volumes do not need to be uniform nor do the successive grid points have to be equally spaced. It is required that a fine grid is employed where steep gradients appear while coarser grid spacing may suffice where small variations occur. The various turbulence models require certain conditions on the grid structure close to solid walls. The so-called high and low Reynolds number versions of these models demand different conditions. In general it is recommended that the solution procedure is carried out on several grids with different fineness and varying degrees of non-uniformity. Then it might be possible to estimate the accuracy of the numerical solution procedure. Adaptive grid techniques might be beneficial to increase the resolution in some circumstances. Complex geometries To handle complex geometries, methods based on body-fitted or curvilinear orthogonal, non-orthogonal grid systems are employed. Such grid systems may be unstructured, structured or block-structured or composite. Because the grid lines follow the boundaries, boundary conditions can more easily be implemented. There are also some disadvantages with non-orthogonal grids. The transformed equations contain more terms and the grid non-orthogonality may cause

10.4 Multi-dimensional models of analysis

unphysical solutions. Vectors and tensors maybe defined as Cartesian, covariant, contravariant and physical or non-physical coordinate- oriented. The arrangement of the variables on the grid affects the efficiency and accuracy of the solution algorithm. Grid generation is an important issue and today most commercial CFD-packages have their own grid generators but also several grid generation packages, compatible with some CFD e codes, are available. The interaction with various CAD (computer-aided-design) e packages is also an important issue today. Further information on treating complex geometries can be found in Refs. [28,29]. The CFD approach The FVM described briefly above is a popular particularly for convective flow, heat and mass transfer. It is also applied in several commercial CFD-codes. In fuel cells often the fuel and oxidant flow fields in narrow channels are assumed to be laminar. If turbulent flow conditions prevail, so-called turbulence modeling approaches are applied. This sub-section gives a brief introduction to the modeling of turbulent flows. Turbulence models The most commonly applied turbulence models in applications are classified as (1) zero-equation models, (2) one-equation models, (3) two-equation models, (4) Reynolds stress models, (5) algebraic stress models, (6) large eddy simulations (LES). The three first models account for the turbulent stresses and heat fluxes by introducing a turbulent viscosity (eddy viscosity) and a turbulent diffusivity (eddy diffusivity). Linear and non-linear models exist, see Refs. [30,31]. The eddy viscosity is usually obtained from certain parameters representing the fluctuating motion. A popular one-equation model is the Spalart-Allmaras model [32] where a transport equation is solved for the eddy viscosity. In two-equation models, these parameters are determined by solving two additional differential equations. However, these equations are not exact but approximate and involve several adjustable constants. Models using the eddy viscosity and eddy diffusivity approach are isotropic in nature and cannot evaluate non-isotropic effects. Models belonging to this category are the k-ε, and k-u models in high or low Reynolds number versions as well as in linear and non-linear versions. Another popular model is the so-called V2F model introduced by Durbin [33]. It extends the use of the k-ε model by incorporating near-wall turbulence anisotropy and non-local pressure-strain effects, while retaining a linear eddy viscosity assumption. Two additional transport equations are solved, namely one for the velocity fluctuation normal to walls and another one for a global relaxation factor. More recently the shear stress transport k-u model (SST k-u) by Menter [34] has become popular as it uses a blending function of gradual transition from the standard k-u model near solid surfaces to a high Reynolds number version of the k-ε model far away from a solid surface. Accordingly it gives accurate prediction of the onset and the size of separation under adverse pressure gradients. In Reynolds stress equation models (RSM) differential equations for the turbulent stresses (Reynolds stresses) are solved and directional effects are naturally accounted for. Six modeled equations (i.e., not exact equations) for the turbulent



CHAPTER 10 Modeling approaches for fuel cells

stress transport are solved together with a model equation for the turbulent scalar dissipation rate ε. RSM models are very complex and require large computing efforts and for this reason are not widely used for industrial flow and heat transfer applications like heat exchangers or duct flows. Algebraic stress models (ASM) and explicit such (EASM) are economic models to account for the anisotropy of the turbulent stresses without solving the Reynolds stress transport equations. An important idea is that the convective and diffusive terms are modeled or even neglected and then the Reynolds stress equations reduce to a set of algebraic equations. For calculation of the turbulent heat and mass (species) fluxes most commonly a simple eddy diffusivity concept (SED) is applied. The turbulent diffusivities for heat and mass transport are then obtained by dividing the turbulent viscosity by turbulent Prandtl numbers. Such a model cannot account for non-isotropic effects in the thermal field but still this model is frequently used in engineering applications [35]. In the LES model, the time-dependent flow equations are solved for the mean flow and the largest eddies while the effects of the smaller eddies are modeled. The LES model has been expected to emerge as the future model for industrial applications but it still limited to relatively low Reynolds number and simple geometries. Handling wall-bounded flows with focus on the near wall phenomena like heat and mass transfer and shear at high Reynolds number present a problem due to the near-wall resolution requirements. Complex wall topologies also present problems for LES. Nowadays, approaches to combine LES and RANS based methods exist. Such are called hybrid models and the DES, detached-eddy simulation [36,37], is such a one. It should be noted that it is not always evident the flow of the fuel, oxidant and products in fuel cells will be turbulent. Handling wall effects There are two standard procedures to account for wall effects in numerical calculations of turbulent flow and heat transfer. One is to employ low Reynolds number modeling procedures, and the other is to apply the wall function method. The wall functions approach includes empirical formulas and functions linking the dependent variables at the near-wall cells to the corresponding parameters on the wall. The functions are composed of laws of the wall for the mean velocity and temperature, and formulas for the near-wall turbulent quantities. The accuracy of the wall function approach is increasing with increasing Reynolds numbers. In general the wall function approach is efficient and requires less CPU time and memory size but it becomes inaccurate at low Reynolds numbers. When low Reynolds number effects are important in the flow domain the wall function approach ceases to be valid. Then so-called low Reynolds number versions of the turbulence models are introduced and the molecular viscosity appears in the diffusion terms. More recently so-called two-layer models have been suggested to enhance the wall treatment. The transport equation for the turbulent kinetic energy is solved while an algebraic equation is used for, e.g., the turbulent dissipation rate.

10.4 Multi-dimensional models of analysis CFD codes Nowadays the usage of commercial general purpose so-called CFD-codes for simulation of fluid flow, heat and mass transfer in a variety of applications is common in industries, engineering and consulting companies worldwide. Among such codes are: ANSYS-FLUENT, COMSOL and others. Some of these codes have modules for electrochemical equipment. Also many universities and research institutes apply commercial codes besides using their in-house developed codes. More recently open-source codes like OPEN-FOAM have become available with some modules for fuel cell and electrochemical devices.

10.4.2 Porous media approach As mentioned in a previous chapter, porous media appear in fuel cell electrodes and the electrolyte may also be regarded as a porous media. Similarly, the electrodes in batteries are porous media but the transport phenomena are different from those in fuel cell electrodes. In this section a brief summary of the governing equations for multiphase flow, heat and mass transfer in porous media are presented. Mass conservation in phase k ε

 vðrk sk Þ v  _k þ rk uk;j ¼ m vt vxj


where ε is the porosity of the porous medium, sk is the saturation or the volumetric fraction of the void space occupied by phase k, and uk,j is the superficial (or Darcian) velocity vector based on the total cross-section area of the porous medium, including _ k represents an interfacial mass the void space and solid matrix. The source term m transfer rate from all other phases to phase k. If no mass source or mass sink is _ k over all phases is equal to zero. present, the summation of the m For two-phase flow of liquid water-air in PEMFCs, the two-fluid model of the mass conservation becomes:    v rg sg v  _ lg þ rg ug;j ¼ m ε (10.9) vxj vt ε

 vðrl sl Þ v  _ lg þ rl ul;j ¼ m vt vxj

(10.10) Momentum conservation in phase k Here the generalized Darcy’s equation without inertia and macroscopic viscous effects is formulated. uk;j ¼ K

 krk  p  rk g j mk k;j


where K is the permeability of the porous medium, mk the dynamic viscosity of phase k, krk is the relative permeability of phase k and gj is the gravity vector.



CHAPTER 10 Modeling approaches for fuel cells

However, as the flow rate is increased, inertia forces become more significant and the linear relation (10.11) does not hold and instead a non-linear relation is needed. The relation between the pressure gradient pk,j and velocity uk,j is then called the Forchheimer equation and for a one-dimensional case it is written as: 

vp ¼ aV þ bV2 vx


where a and b are constants. Mass conservation of species a in phase k The following equation is valid:   v v  v a  εrk sk Cak þ j þ Jak rk uk;i Cak ¼  vt vxi vxi k;i


where Cak is the mass concentration of species a in phase k, and jak;i is a diffusive flux of species a in phase k due to molecular diffusion, hydrodynamic dispersion or both. The diffusive flux is commonly written in a so-called Fickian form, namely jak;i ¼ εrsk Dak

v  a C vxi k


where Dak is a macroscopic second-order tensor and depends on the molecular diffusion coefficient and the definition of the fluid velocity. The term Jak denotes the species source arising from interphase species transfer, chemical reaction, and phase change. If no species generation occurs due to chemical or biological reactions, the summation of these terms over all phases k becomes zero. For two-phase flows of liquid water-air in PEMFCs, the water transport equations in the gas and liquid read:   v v  v w (10.15) εrg sg Cw rg ug;i Cw j þ Jw g þ g ¼ g vt vxi vxi g;i  v v  ðεrl sl Þ þ rl ul;i ¼ Jw l vt vxi


w where Jw g þ Jl ¼ 0. Energy conservation in phase k The equation below presents the energy conservation in phase k where it is assumed that local equilibrium exists between the phases so that Tk ¼ T.    v v  v vT ðεrk sk hk Þ þ r uk;i hk ¼ (10.17) sk kk þ qk vt vxi k vxi vxi In Eq. (10.17), kk is the effective thermal conductivity of phase k, hk is the phase enthalpy, and qk is the interphase heat transfer rate associated with phase k.

10.5 Example proton exchange membrane fuel cells e PEMFC

The interphase heat transfer rates are related to the external heat source or sink q_ by X qk ¼ q_ (10.18) k Multiphase mixture model Instead of the approach presented so far in Section 10.4.2, a multiphase mixture model for porous media can be established. Then, the mixture’s physical properties are based on the properties of its constituents. The mixture density, velocity, species concentration, enthalpy and pressure are the important properties. Governing equations for mass conservation, momentum conservation, species conservation and energy conservation are set up. The details are not given here but such can be found in Refs. [38,39].

10.4.3 Molecular dynamics based approaches In this section, brief results of an approach to establish the nano/micro-scale structure of the porous anode for, e.g., an SOFC, and then to calculate the thermal-physical properties like thermal conductivity, thermal expansion coefficient and equilibrium lattice constants as well as to predict the thermal behavior and temperature distribution inside the porous anode. To enable this, the following procedures have to be included, namely, (a) AA: all-atom modeling, (b) CG: coarse-graining modeling (c) MD: molecular-dynamics modeling. More details of this can be found in Wang et al. [40] and Sunden [41,42] and accordingly here just a typical result is given in Fig. 10.2. This Fig. shows thermal conductivity versus temperature and it reveals that the theoretical approach captures a similar trend as the experimental values but the precise agreement is not sufficient here. Similar research works have been carried out for polymer electrolyte membrane fuel cells (PEMFC) and for instance some papers by Xiao et al. [43,44] as well as Sunden [41] describe such works.

10.5 Example proton exchange membrane fuel cells e PEMFC In this section, a three-dimensional, non-isothermal, two-phase flow model based on the finite volume method will be briefly presented and an investigation of the effect of the gas diffusion layer (GDL) deformation on transport phenomena of PEMFCs with interdigitated flow channels. The oxygen concentration, temperature, liquid water saturation, volumetric current density, and cell performance of PEMFCs with and without compression are presented and compared. Further details can be found in Li and Sunden [45] and Li [46].


CHAPTER 10 Modeling approaches for fuel cells


Thermal Conductivity (W/(mK))


Simulation Value Literature Value

90 80 70 60 50 300



600 700 800 Temperature/K



FIG. 10.2 CG-MD prediction of thermal conductivity.

10.5.1 Model description Fig. 10.3 provides schematic illustrations of a PEM fuel cell using interdigitated flow channels with and without compression, respectively. Due to the symmetry of the configurations in the y-coordinate direction only, half inlet and outlet channels are included in the computational domain. The geometry of the fuel cell and operating conditions are summarized in Tables 10.1 and 10.2, respectively. The fuel cell operating temperature and pressure are 353 K and 1 atm., respectively. In addition, 100% relative humidity is assumed for both the anode and cathode reactant gases. In the PEMFC mathematical model, the fluid flow is laminar, the ideal gas law is applied for the reactant gases, the CLs are homogeneous and isotropic, the reactant gases cannot diffuse across the membrane and the generated water in the cathode CL is in dissolved phase [50].

10.5.2 Governing equations With the mentioned assumptions, the governing equations of mass, momentum, species, energy, charge, liquid water and dissolved water are presented below [47,50]: Mass conservation equation: v ðruj Þ ¼ Smass vxj


where r and uj are the mixture fluid density and superficial velocity, respectively. Smass represents the source term of the mass conservation equation. The mass changes of hydrogen, oxygen and water vapor due to the processes of electrochemical reactions and phase change are considered in the relevant layers.

10.5 Example proton exchange membrane fuel cells e PEMFC

FIG. 10.3 Schematic illustrations of PEMFCs with interdigitated flow channels, (A) without compression, (B) with compression. The computational domain in a cross sectional plane is indicated.

Table 10.1 Fuel cell geometry dimensions [47e49]. Parameter



Fuel cell length Fuel cell width Flow channel width Flow channel height Rib width Current collector height GDL thickness CL thickness Membrane thickness

50 2 1 1 1 1.5 0.2 0.01 0.05

mm mm mm mm mm mm mm mm mm

Table 10.2 Fuel cell operating conditions. Parameter



Reactant gas Operating pressure, Pa/Pc Stoichiometric ratio, xa/xc Relative humidity, RHa/RHc Operating temperature, Ta/Tc

Hydrogen/water 1 atm 1.5 100% 353 K

Air/water 1 atm 1.5 100% 353 K



CHAPTER 10 Modeling approaches for fuel cells

Momentum conservation equation: v v ðrui uj Þ ¼ vxj vxj

vui m vxj


vp þ Smom vxi


where p and m are the mixture pressure and dynamic viscosity, respectively. Smom denotes the source term of the momentum equation. Species conservation equation: ! v v vYi ðruj Yi Þ ¼ (10.21) rDeff;i þ Si vxj vxj vxj where Yi and Deff,i are the mass fraction and effective diffusivity for the i th species, respectively. Si represents the amount of consumption or generation of the i th species. The effective diffusivity (Deff,i) is obtained by the Bruggeman correlation (see Table 10.4) which accounts for the tortuous path in the GDLs and CLs. Energy conservation equation: ! v v vT ðrcp uj TÞ ¼ (10.22) keff þ ST vxj vxj vxj where cp and keff are the specific heat and effective thermal conductivity, respectively. The irreversible, reversible, ohmic heat generation and phase change terms are all included in the energy equation source term, ST. Charge conservation equation: ! v vfs f þ Ss ¼ 0 (10.23) seff;s vxj vxj s ! v vfm f (10.24) seff;m þ Sm ¼ 0 vxj vxj m where seff,s is the effective electrical conductivity, seff,m the effective protonic conductivity, fs the electrical potential, fm the protonic potential. The Butler-Volmer equation and spherical agglomerate model are adopted to describe the hydrogen oxidation reaction and oxygen reduction reaction in the anode and cathode CLs, respectively. Liquid water transport equation: !   Krl mg v v vS uj ¼ ðr Ds (10.25) rl þ Sl vxj vxj l vxj Krg ml where Krl and Krg are the relative permeability of the liquid and gas phase, respectively. Ds denotes the capillary diffusion coefficient. The source term, Sl, represents the phase change process between water vapor and liquid water, which is determined by the water vapor partial pressure and saturation pressure. s is the liquid water

10.5 Example proton exchange membrane fuel cells e PEMFC

saturation, which is defined as the liquid volume fraction in the pore space and is related to the mixture water concentration. Dissolved water transport equation: ! ! v nd vfm v rm vl  sm Dl (10.26) ¼ þ Sd vxj F vxj Mm vxj vxj where nd and Dl are the electro-osmotic drag coefficient and water diffusivity in the membrane, respectively. The generated water in the cathode CL is in dissolved phase, and the absorption/desorption processes of the ionomer in the anode and cathode CLs take place between water vapor and dissolved water. The parameters and complementary equations used in the mathematical model are summarized in Tables 10.3 and 10.4. In addition, the source terms of the governing equations are given in Table 10.5.

10.5.3 Catalyst layer composition and volume fraction The catalyst layer consists of void space, Pt/C particles, and ionomer phase. The volume fraction of each component is determined by the following expressions [53,64,65]. The volume fraction of the void space (εCL) is given by: εCL ¼ 1  LPt=C  Li


where LPt/C and Li are the volume fractions of Pt/C particles and ionomer phase, respectively. The volume fraction of the Pt/C particles (LPt/C) is determined by: ! mpt 1 1f 1 þ LPt=C ¼ (10.28) f rc tCL rpt where mpt is the platinum loading, tcl the thickness of the CL, rpt the density of platinum, rc the density of carbon. The parameter f is defined as the platinum loading divided by the sum of platinum loading and carbon loading. mpt f¼ (10.29) mpt þ mc where mc is the carbon loading. The volume fraction of the ionomer phase (Li) is expressed by:  h i  LPt=C  r3agg Li;agg þ raggþdi 3  r3agg Li ¼ 3  ragg 1  Li;agg


The ionomer phase has two parts, namely one part inside the agglomerate, ans another part is the ionomer film. When LPt/C, Li, ragg, and Li,agg are given, a unique value of the ionomer thickness di can be determined. It is assumed that the agglomerate is only occupied by Pt/C particles and ionomer phase.



CHAPTER 10 Modeling approaches for fuel cells

Table 10.3 Parameters used in the mathematical model. Parameter



Platinum loading, mpt [51] Platinum density, rpt [49] Carbon loading, mc [49] Carbon density, rc [49] Dry membrane density, rm [49] Membrane equivalent weight, Mm [49] Radius of agglomerate, ragg [49] Volume fraction of ionomer, Li [50] Volume fraction of ionomer in agglomerate, Li,agg [49] Anode reference exchange current density, iref a [49] Cathode reference exchange current density, iref c [52] Anode transfer coefficient, aa [53] Cathode transfer coefficient, ac [53] Reference hydrogen concentration, cref H2 [54] Reference oxygen concentration, cref O2 [54] Hydrogen Henry’s constant, HH2 [50] Oxygen Henry’s constant, HO2 [50] Thermal conductivity of CC/GDL/CL, kCC/GDL/CL [49] Thermal conductivity of membrane, km [49] Electrical conductivity of CC/GDL/CL, ss,CC/GDL/CL [49] Entropy of hydrogen oxidation, DSa [55] Entropy of oxygen reduction, DSc [55] Latent heat of condensation/evaporation, Dhlg [56] Liquid water viscosity, ml [56] Surface tension, s [56] Contact angle of GDL/CL, qGDL/CL [56] Condensation rate, gcon [50] Evaporation rate, gevap [50] Dissolved water phase change rate, g [50] Permeability of CL, KCL [49] Binary diffusivity, DH2 H2 O [57] Binary diffusivity, DO2 H2 O [57] Binary diffusivity, DO2 N2 [57] Binary diffusivity, DH2 ON2 [57]

0.4 2.145  104 0.6 1.8  103 1.98  103 1.1 1 0.4 0.5

mg cm2 kg m3 mg cm2 kg m3 kg m3 kg mol1 mm e e


A m2


A m2

0.5 1 56.4 3.39 4.56  103 0.101325e(666/Tþ14.1) 100/1.7/0.3

e e mol m3 mol m3 Pa m3 mol1 Pa m3 mol1 W m1 K1

0.25 20000/5000/2000

W m1 K1 S m1

0.104 326.36 2.36  106

J mol1 K1 J mol1 K1 J kg1

3.517  104 0.0625 110 /95 100 100 1 1.0e13 9.15  105 2.82  105 2.2  105 2.56  105

Pa s N m1 e s1 s1 s1 m2 m2 s1 m2 s1 m2 s1 m2 s1

Table 10.4 Complementary equations and definitions. Description Deff,i ¼ (1  s)


m2 s1 m2 s1


1Xi Di;m ¼ P n 

Xj Di;j


Binary mass diffusivity [57] Effective thermal conductivity [48] Electrochemical kinetics [58]


Di;j ¼ Di;j ðT0 ; P0 Þ

P0 P

keff ¼ (1  ε)ks þ εkf  PH2 ja ¼ ð1  sÞiref a eff a cref H H2

Over-potential [58] Open circuit voltage [49] Active surface area [53] Proton conductivity [60]



Relative permeability [50] Capillary diffusivity [50]

eaa Fha =RT  eac Fha =RT

104 ðT

 105 TLn

Voc ¼ 1:229  8:456   298:15Þ þ 4:31  mpt  3 2 aeff ¼ tCL 227:79f  158:57f  201:53f þ 159:5  103   1268

W m1 K1 A m3

ha ¼ fs  fm, hc ¼ fs  fm  Voc

sm ¼ ð0:514l  0:326Þe Effective conductivity [49]

m2 s1

T T0

  PH2 P0:5 O2

1 1 303 T

seff;s ¼ ð1  εGDL Þ1:5 ss ; seff;s ¼ ð1  εCL  Li Þ1:5 ss seff,m ¼ L1.5 i sm Krl ¼ Ks3, Krg ¼ K(1  s)3 Ds ¼

3 dPc Ksm ds l


Capillary pressure [50]

Pc ¼ scosðqÞ

Saturation pressure [50]

log10Psat ¼ 2.1794 þ 0.02953(T  273.15)  9.1837  105(T  273.15)2

ε K


1:417s  2:12s2 þ 1:263s3

þ1.4454  10 (T  273.15)

V V m1 S m1 S m1 m2 m2 s1 Pa Pa



10.5 Example proton exchange membrane fuel cells e PEMFC

Effective mass diffusivity [58] Mass diffusivity [59]

Units 1.5 1.5



Description Electro-osmotic drag coefficient [60]

Units nd ¼

l 2:5 22

Dissolved water diffusivity [60]

Dl ¼ 1010 e Equilibrium water content [61] Water activity [50] Oxygen diffusivity in liquid water [62]

( leq ¼


 8 > 2:05l  3:25 1 1 < 303 T 6:65  1:25l > : 2:563  0:33l þ 0:0264l2  0:000671l3

1:41 þ 11:3a  18:8a2 þ 16:2a3 10:1 þ 2:94ða  1Þ

a ¼ XPwvsatP þ 2s DO2 ;w ¼ 7:4  1012

TðjMH2 OÞ mH2 OV0:6 O

DO2 ;i ¼ 2:88  1010 e d2f ε3 16kCK ð1εÞ2

K ¼

Inertial coefficient [63]

b ¼ 2:88  106 ε1:51 K

ð2  l < 3Þ ð3  l < 4Þ ð4 < lÞ

m2 s1



Oxygen diffusivity in ionomer [50] Permeability of the GDL [63]

ða < 1Þ ða  1Þ




m2 s1


1 1 313 T

m2 s1 m2 m1

CHAPTER 10 Modeling approaches for fuel cells

Table 10.4 Complementary equations and definitions.dcont’d

10.5 Example proton exchange membrane fuel cells e PEMFC

Table 10.5 Source terms in the governing equations [50]. Description


Smass ¼ SH2 þ Swv Anode CL Smass ¼ SO2 þ Swv Cathode CL Smass ¼ Swv Anode and cathode GDLs   Smom ¼  Km ! u þ brj! u j! u Anode and cathode GDLs and CLs

kg m3 s1

ja SH2 ¼ 2F MH2

kg m3 s1

Anode CL

jc SO2 ¼ 4F MO2 Cathode CL Swv ¼ Sl  Svd MH2 O Anode and cathode CLs Swv ¼ Sl Anode and cathode GDLs 2 2 a ST ¼ ja ha  TDS 2F ja þ seff;m kVfm k þ seff;s kVfs k þ Sl Dhlg Anode CL 2 2 c ST ¼ jc hc  TDS 4F jc þ seff;m kVfm k þ seff;s kVfs k þ Sl Dhlg Cathode CL 2 ST ¼ seff,mkVfmk membrane ST ¼ seff,skVfsk2 þ SlDhlg Anode and cathode GDLs ST ¼ seff,skVfsk2 Anode and cathode CCs Ss ¼ ja Anode CL Ss ¼ þjc Cathode CL Sm ¼ þja Anode CL Sm ¼ jc Cathode CL Sl ¼ Sphase Anode and cathode GDLs and CLs 8 εð1  sÞ > > MH2 OðPwv  Psat Þ Pwv  Psat < gcond RT Sphase ¼ > > : gevap εs MH OðPwv  Psat Þ Pwv < Psat RT 2 Sd ¼ Svd Anode CL Sd ¼ Svd þ Sl athode CL rm Svd ¼ g M ðleq  lÞ m

jc Sl ¼ 2F

kg m2 s2

kg m3 s1 kg m3 s1 W m3

A m3 A m3 kg m3 s1 kg m3 s1

mol m3 s1

Anode and cathode CLs

Cathode CL

The thickness of the ionomer film (di) is obtained by: # "sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   3 Li 1  Li;agg  Li;agg þ 1  1 di ¼ ragg LPt=C


The total volume of the liquid water is: Vw ¼ sεCL VCL


where s is the liquid water saturation, VCL is the total volume of the CL. The volume of the liquid water covering the individual agglomerate can be obtained by: sεCL (10.33) Vw;i ¼ N where N is the number of agglomerate particles per catalyst layer volume.



CHAPTER 10 Modeling approaches for fuel cells

3LPt=C   4pr3agg 1  Li;agg

The thickness of the liquid water film (dw) is obtained by: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3sεCL 3 dw ¼ ðragg þ di Þ3 þ  ðragg þ di Þ 4pN



10.5.4 Cathode agglomerate model In this analysis, the spherical agglomerate model [65,66] is used to calculate the volumetric current density. Identical agglomerates are uniformly distributed in the cathode CL, and the individual agglomerate is evenly covered by ionomer and liquid water films. !!1 ragg þ di þ dw PO2 1 di dw jc ¼ 4F þ þ HO2 Er kc ð1  εCL Þ ragg aagg;i DO2 ;i aagg;w DO2 ;w (10.36) where PO2 is the partial pressure of oxygen, HO2 the Henry’s constant, Er the effective factor, kc the reaction rate constant, aagg,i the ionomer effective agglomerate surface area, and aagg,w the liquid water effective agglomerate surface area. The effectiveness factor of the spherical agglomerate is defined as:   1 1 1  Er ¼ (10.37) FL tanhð3FL Þ 3FL where the Thiele’s modulus FL for chemical reactions is given as: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ragg kc FL ¼ 1:5 3 Li;agg DO2


The reaction rate constant is determined by: kc ¼

h i iref aa Fhc =RT ac Fhc =RT c aeff þ e  e 4Fð1  εCL Þcref O2


The effective agglomerate surface area is obtained from: aagg;i ¼ aagg;w ¼

3LPt=C εCL  ðragg þ di Þ2 1  Li;agg


3LPt=C εCL  ðragg þ di þ dw Þ2 r3agg 1  Li;agg



10.5 Example proton exchange membrane fuel cells e PEMFC

10.5.5 Determination of porosity of the GDLs after compression A constant and uniform porosity is adopted for simplicity. However, the porosity is varied and non-uniformly distributed due to the deformation resulting from the assembly force. It is assumed that the thickness and porosity are decreased only caused by the change in pore size. Therefore, the GDL porosity is determined by the following expression [47,67]: d0 ε ¼ 1  ð1  ε0 Þ (10.42) d where ε0 is the initial porosity, d0 the initial thickness, d the thickness after compression, ε the porosity after compression. The local porosity is varying with the thickness of the GDL, and the transport parameters associated with porosity are also changed.

10.5.6 On numerical implementation and boundary conditions The commercial software ANSYS FLUENT was used in this particular example. The governing equations of charge, liquid water transport, and dissolved water transport were implemented by using user defined scalar (UDS) equations with a developed under relaxation technique. At the inlet of the flow channels, the mass flow rate, temperature, and species mass fractions are prescribed. The mass flow rates of reactants at both anode and cathode sides are calculated at a reference current density of 1.0 A/cm2 and are given by Eqs. (10.43) and (10.44), see also Chapter 9. The anode and cathode stoichiometric ratios (xa and xc) were given in Table 10.2. In addition, the inlet liquid water saturation is assigned as zero. At the outlet of the flow channels, a pressure-outlet boundary condition is applied. The operating temperature and a constant electric potential, fs ¼ 0, are specified at the anode terminal. At the cathode terminal, the operating temperature and a constant electric potential, fs ¼ Vcell, are applied. A symmetry boundary condition is adopted at the outer surface of the x-z plane. xa MH2 Qa ¼ Iref Am (10.43) 2FYH2 Qc ¼

xc MO2 Iref Am 4FYO2


The SIMPLE algorithm was used to handle the pressure-velocity coupling. A second order upwind discretization scheme was adopted for the momentum, species, energy and UDS equations. A grid independence study was carefully performed to balance the accuracy and computational resources. A mesh system with the total number of control volumes of 208,000 was applied. The performance of the modeling approach was evaluated by comparison with the experimental data reported by Yan et al. [68] for cell voltage versus current density. It was found that the modeling approach was able to reproduce the experimental data well.



CHAPTER 10 Modeling approaches for fuel cells

10.5.7 Some characteristics In this sub-section some results from the above modeling approach are highlighted. The cell performance is commonly presented and compared in terms of the polarization and power density curves. Fig. 10.4 shows the effect of assembly forces on the performance of PEM fuel cells with interdigitated flow fields. It can be seen that the current density is increased with increasing assembly force. For an operating voltage 0.7 V, the current densities for four different cases are 0.2383 A/cm2, 0.2422 A/cm2, 0.2565 A/cm2 and 0.2575 A/cm2, respectively. In addition, the maximum power densities are 0.4227 W/cm2, 0.4263 W/cm2, 0.4464 W/cm2 and 0.4503 W/cm2, respectively. In a similar work carried out by Mahmoudi et al. [69], the cell performance was decreased as the compression ratio was increased. Note that a constant pressure drop between the inlet and outlet channels were adopted in that study. He et al. [70] numerically investigated the influence of pressure drop on the cell performance. It was concluded that the cell performance was substantially improved with increasing pressure drop. The pressure drops between the cathode inlet and outlet flow channels for the four cases at the operating voltage 0.7 V are 146.2 Pa, 428.4 Pa, 1111.3 Pa and 3104.1 Pa, respectively. The pressure drop of the PEM fuel cell with 2.0 MPa assembly force is approximately 21.2 times higher than that of the PEM fuel cell without an assembly force. This is because the mass flow rate is constant for all cases in this study, and the thickness, porosity and permeability of the GDLs under the rib regions and the cross sectional area of the flow channels are decreased due to the deformation of the GDLs caused by the assembly force.

FIG. 10.4 Polarization and power density curves of a PEMFC under different assembly forces.

10.6 Example solid oxide fuel cells e SOFC

FIG. 10.5 Volumetric current density of a PEMFC under different assembly forces.

The volumetric current density distributions at the middle plane of CLs (Line-1) are shown in Fig. 10.5. The minimum and maximum volumetric current densities appear at the outlet channel region and under the rib region for all cases, respectively. The volumetric current density is gradually increased when the assembly force is increased from 0 MPa to 1.5 MPa. When the assembly force is increased from 1.5 MPa to 2.0 MPa, the maximum volumetric current density is increased and the minimum volumetric current density is decreased. This indicates that a higher variation of the volumetric current density is obtained under the assembly force of 2.0 MPa.

10.6 Example solid oxide fuel cells e SOFC In this section, a typical mathematical model for SOFCs including fluid flow, electrochemical processes, heat and mass transfer will be presented. The domains considered are the fuel and oxidant channels, the anode and cathode, the electrolyte and the interconnects. The chemical reactions occur inside the electrodes and water is produced at the anode electrode while oxygen ions are produced at the cathode electrode. SOFCs operate at high temperature so there is no condensation of water at the anode. Thus the important processes taking place are the diffusion of hydrogen and oxygen, the production of water and oxygen ions in the electrodes, the transport of oxygen ions across the electrolyte and the transfer of heat. Further details can be found by Kakac et al. [71] and Sunden and Faghri [72].



CHAPTER 10 Modeling approaches for fuel cells

10.6.1 Transport of fuel and oxidant in channels The governing equations for the fuel and oxidant (commonly air) are identical to those presented in Section 10.5 for the PEM fuel cells, i.e., Eqs. (10.19)e(10.22), expressing conservation of mass, momentum, species balance and energy conservation.

10.6.2 Transport in the porous electrodes The porous electrodes consist of gas and solid phases. The reactions are in most models assumed to occur in the gas phase. The conservation laws have to be written for each phase. The energy equation includes a heat source term due to the chemical reaction rate and in the species conservation equation a source term based on the chemical reaction rate also appears. In the equations below, the porosity ε appears. Conservation of mass v ðrεuj Þ ¼ 0 vxj Conservation of momentum v vp v ðrεui uj Þ ¼ ε þ vxj vxi vxj


vui vuj þ εm vxj vxi

(10.45) #! þ

ε2 mui k


In Eq. (10.46), the Darcy law is introduced to handle the porous media source term. Balance equation for species  v v  ðrεuj ci Þ ¼  (10.47) εJj;i þ Ss;i vxj vxj In this equation, ci is the mass fraction of species i and Jj,i is the diffusive mass flux of species i while Ss,i is an additional source term for species i. Energy conservation " #! v v vT X ðrεci uj Þ ¼  hi Jj;i (10.48) ε keff þ Selectrodes þ Sradiation vxj vxj vxj i Here keff means the effective macroscopic thermal conductivity which is calculated as keff ¼ εkgas þ ð1  εÞksolid


To handle the mass transport in the porous medium commonly three modeling approaches are used. The first one is Fick’s law which has been introduced in earlier chapters. Another one is the Stefan-Maxwell model which is more suitable for multi-component systems and was also given in an earlier chapter. The third is the so-called Dusty Gas model which may be regarded as an extended Stefan-Maxwell equation.

10.6 Example solid oxide fuel cells e SOFC

10.6.3 Transport in the solid electrolyte The transport of O2 ions in the electrolyte is described by taking into account the ion transport due to Ohm’s law and the conservation of charge as jio ¼ seff io Vfio   V$ seff io Vfio ¼ Sio

(10.50) (10.51)

In Eqs. (10.50) and (10.51) fio is the ionic potential, seff io the effective ionic conductivity, jio the ionic current density and Sio the source term from reactions. The heat balance equation in the electrolyte reads   d dTe  A e ke Dx ¼ qeanode þ qecathode þ qconvection; anode dx dx (10.52) þ qconvection; cathode þ qradiation;e þ qgeneration; electrolyte þ qohmic; electrolyte In Eq. (10.52), index e stands for the electrolyte. The left hand term represents conduction in the electrolyte and the right hand side covers the various heat fluxes due to conduction and convection with the electrodes and the reversible and irreversible generated heat fluxes.

10.6.4 Transport in the interconnects The interconnect is used to collect the current from the SOFC and accordingly the electron transport from Ohm’s law and the conservation of charge must be considered. Then one has jel ¼ seff el Vfel   Vf V$ seff el ¼ Sel el

(10.53) (10.54)

In Eqs. (10.53) and (10.54) fel is the electronic potential, seff el the effective electronic conductivity, jel the electronic current density and Sel the source term from reactions. The charge balance on electrons and ions requires V$jio ¼ V$jel


The heat balance for the interconnect is written in similar way as Eq. (10.52), i.e.,   d dTint  Aint kint Dx ¼ qintelectrodes þ qconvection; channels þ qradiation;channels dx dx þ qgeneration; interconnect þ qohmic; interconnect (10.56)



CHAPTER 10 Modeling approaches for fuel cells

10.6.5 Model of the electrochemical processes The SOFC voltage is during operation decreased from its open circuit voltage because of the irreversible losses as already described in previous chapters. Also models how to calculate the various activation or polarization losses have been presented in previous chapters and are not repeated here.

10.6.6 Some results In this subsection, a few results from a fully three-dimensional numerical calculation procedure are presented. This is based on ref. [73] for an intermediate temperature SOFC. The reforming reactions in a reduced temperature SOFC anode are included. The oxidation and shift reactions at the active surface and the porous anode, respectively, are considered. The duct under consideration includes the porous anode, fuel gas flow duct and the solid interconnects. Calculations of the fuel gas species, mass and heat generation, and consumption related to the internal reforming and electrochemical reactions are performed. Varying thermal-physical properties and transport parameters are taken into account. Fig. 10.6 shows a schematic drawing of the considered composite anode duct. In Fig. 10.7 axial velocity contours are presented. Due to the permeation and effects of mass generation or consumption, a uniform distribution and symmetry do not exist. The internal reforming reactions are reflected in the CH4 mass concentration in Fig. 10.8. The temperature distribution along the main flow direction is shown in Fig. 10.9.

10.6.7 Engineering bridges in analysis of multiscale issues The modeling of the nanocomposites in the catalytic layers of fuel cells is a multiscale problem. The continuity of the atomistic and continuum regions of catalyst

FIG. 10.6 Schematics of a composite anode duct.

10.6 Example solid oxide fuel cells e SOFC

FIG. 10.7 Dimensionless axial velocity contours (U/Uin) along the main flow stream.

FIG. 10.8 CH4 mass concentration.

layers necessitates a seamless coupling between these two regions. However, significant conceptual and practical challenges remain how to consistently couple a discrete atomistic description to an averaged continuum mechanical model for mechanical problems. One way to bridge the gap is to introduce a transition region, see Ref. [74]. Due to the large discrepancy between the length and time scales in



CHAPTER 10 Modeling approaches for fuel cells

FIG. 10.9 Temperature distribution along the main flow direction.

the nano- and macro-domains, several researchers have paid attention to the transition region, [75,76]. Ideas for coupling of the MD, DPD and SPH methods with the finite element method have been presented. In engineering modeling and simulations as well as in design and optimization of fuel cells and the ongoing transport, chemical and electrochemical phenomena, it is convenient to apply models or methods based on the continuum formulation approach. Then effective values of the thermal conductivities, gas diffusion coefficients, porosity, tortuosity factors, interfacial and/or volumetric heat transfer coefficients between the fluid and solid surfaces are needed. Small scale particles from micrometer size down to nanometer size are present in the multifunctional porous electrodes in electrochemical cells to enhance the catalytic reactions and the cell performance. The charge transport, the reactions and the heat and mass transfer processes in the cell porous electrodes are strongly affected by the small scale and complex porous structure. The processes are interlinked. Various models have been presented in the literature, based on experimental and numerical approaches of analysis, to calculate effective physical and chemical transport coefficients. Reviews in Refs. [77e79] present evaluations and limitations of such models and recommendations for engineering applications are provided. One of the main goals of multiscale modeling is to drastically reduce the number of degrees of freedom of engineering problems while maintaining accuracy in the regions of interest. So far, the developed methods have not achieved the goal of attaining large realistic system sizes. Even with the use of parallel computing techniques, these methods still require very long time to simulate problems of interest. All single-scaled modeling methods provide results that lead to an understanding of the properties under specific conditions. This information is then passed on to the

10.8 Summary

designers to learn and to be incorporated in future designs. Thus multiscale modeling should be developed further to demonstrate its usefulness for engineering design.

10.7 Softwares For the macroscale modeling approaches using the continuum formulation, many commercial codes (like ANSYS-FLUENT and COMSOL) and their specific fuel cell modules have been considered as well as in-house softwares (at universities and institutes) have been used. Open source codes like Open Foam are also available with some fuel cell modules and used frequently. To handle systems that couple several scales, a theoretical and computational platform including a first principle is needed. This can be built up at least partly by available openesource codes. Here a brief summary of available codes is provided. The Atomic Simulation Environment (ASE) is a common part of a simulation tool developed at DTU, Denmark. It can be used to run molecular dynamic simulations when the atomic numbers and positions are given (all atom modeling), see Ref. [80]. The Visual Molecular Dynamics (VMD) software [81] is a molecular visualization program for display, animation and analysis of large biomolecular systems using 3D-graphics and built-in scripting. The coarse graining builder module in VMD can be employed to, e.g., transform unit structures to CG beads. GROMACS (Groningen Machine for Chemical Simulation) [82] is a molecular dynamics simulation package originally developed at University of Groningen, Netherlands. As the coordinates of the CG beads are at hand, the CG-MD method can be implemented in the GROMACS package. Packmol [83] packs molecules in defined space regions and create a starting point for MD simulations. The packing makes sure that short range repulsive interactions do not disrupt the simulations. LAMMPS (large scale atomic/molecular massively parallel simulator) can be applied to calculate, e.g., the thermal properties, thermal behavior and temperature distribution inside a porous SOFC anode. The application of the Lattice Boltzmann method as a microscale model of an SOFC anode using the programs Palabos with Python and MATLAB together with Paraview for visualization was illustrated in Ref. [16].

10.8 Summary Computational techniques and methods are an evolutionary field from macroscopiclevel to nano-level using continuum and discrete mechanics. It is likely that detailed research methodologies as described in this chapter will continue to contribute to



CHAPTER 10 Modeling approaches for fuel cells

further development of improved effective models which will enhance the prediction ability and promote the introduction of electrochemical cells in engineering sustainable applications. A number of open source computer codes as well as commercial ones are available. The presented information reflects the opportunities to reasonably well model and simulate phenomena and processes at various scales ranging from the catalyst layers, electrodes and membranes to unit cells and stacks.

References [1] D. Cheddie, N. Munroe, Review and comparison of approaches to proton exchange membrane fuel cell modeling, J. Power Sources 147 (2005) 73e84. [2] M. Andersson, J. Yuan, B. Sunden, Review on modeling development for multi-scale chemical-reactions-coupled transport phenomena in SOFCs, Appl. Energy 87 (2010) 1461e1476. [3] M. Khan, B. Sunde´n, J. Yuan, Analysis of multi-phase transport phenomena with catalyst reactions in PEMFC - a review, J. Power Sources 196 (2011) 7899e7916. [4] D. Raabe, Overview of the lattice Boltzmann method for nano- and microscale fluid dynamics in materials science and engineering, Model. Simulat. Mater. Sci. 12 (2004) R13eR46. [5] Y. Dong, D. Bhattacharyya, P.J. Hunter, Experimental characterisation and objectoriented finite element modelling of polypropylene/organoclay nanocomposites, Compos. Sci. Technol. 68 (2008) 2864e2875. [6] M. Fermeglia, S. Pricl, Multiscale molecular modelling in nanostructured material design and process system engineering, Comput. Chem. Eng. 33 (2009) 1701e1710. [7] G. Dorenbos, Y. Suga, Simulation of equivalent weight dependence of nafion morphologies and predicted trends regarding water diffusion, J. Membr. Sci. 330 (2009) 5e20. [8] M. Liu, P. Meakin, H. Huang, Dissipative particle dynamics simulation of fluid motion through an unsaturated fracture and fracture junction, J. Comput. Phys. 222 (2007) 110e130. [9] A.M. Tartakovsky, P. Meakin, T.D. Scheibe, R.M. Eichler West, Simulations of reactive transport and precipitation with smoothed particle hydrodynamics, J. Comput. Phys. 222 (2007) 654e672. [10] F. Jiang, A.C.M. Sousa, Smoothed particle hydrodynamics modeling of transverse flow in randomly aligned fibrous porous media, Transport Porous Media 75 (2008) 17e33. [11] P. Meakin, A.M. Tartakovsky, Modeling and simulation of pore-scale muliphase fluid flow and reactive transport in fractured and porous media, Rev. Geophys. 47 (2009) 47, Paper no. RG3002. [12] L.B. Lucy, A numerical approach to the testing of the fission hypothesis, Astron. J. 82 (1977) 1013e1024. [13] A.M. Tartakovsky, P. Meakin, Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics, Adv. Water Resour. 29 (2006) 1464e1478. [14] F. Muller-Plathe, Coarse-graining in polymer simulation: from the atomistic to the mesoscopic scale and back, ChemPhysChem 3 (2002) 754e769. [15] M.L. Klein, W. Shinoda, Large-scale molecular dynamics simulations of selfassembling systems, Science 321 (2008) 798e800.


[16] H. Paradis, Micro-and Macroscale Modeling of Transport Processes in Solid Oxide Fuel Cells (Ph.D. thesis), Lund University, Lund, 2013. [17] S. Revankar, P. Majumdar, Fuel Cells, Principles, Design and Analysis, CRC Press, 2014. [18] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York, 1980. [19] D.A. Andersson, J.C. Tannehill, R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer, second ed., Taylor and Francis, USA, 1997. [20] G.D. Smith, Numerical Solution of Partial Differential Equations, Oxford University Press, London, 1978. [21] J.N. Reddy, D.K. Gartling, The Finite Element Method in Heat Transfer and Fluid Dynamics, CRC Press, Boca Raton, Fla, 2010. [22] R. W Lewis, K. Morgan, H.R. Thomas, K.N. Seetharamu, The Finite Element Method in Heat Transfer Analysis, J. Wiley and Sons, 1996. [23] M.S. Kandelousi, D.D. Ganji, Hydrothermal Analysis in Engineering Using Control Volume Finite Element Method, Academic Press, Oxford, 2015. [24] L.C. Wrobel, Boundary Element Method e Volume 1 Applications Thermo-Fluids and Acoustics, J. Wiley and Sons, UK, 2002. [25] C.M. Rhie, W.L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA J. 21 (1983) 1525e1532. [26] D.S. Jang, R. Jetli, S. Acharya, Comparison of the PISO, SIMPLER and SIMPLEC algorithms for treatment of the pressure velocity coupling in steady flow problems, Numer. Heat Tran. 10 (1986) 209e228. [27] R.I. Issa, Solution of the implicity discretized fluid flow equations by operator-splitting, J. Comput. Phys. 62 (1986) 40e65. [28] D. McBride, N. Croft, M. Cross, Combined vertex-based-cell-centred finite volume method for flow in complex geometries, in: Third International Conference on CFD in the Minerals and Process Industries, CSIRO, Melbourne, Australia, 2003, pp. 1351e1356. [29] B. Farhanieh, L. Davidson, B. Sunden, Employment of the second-moment closure for calculation of recirculating flows in complex geometries with collocated variable arrangement, Int. J. Numer. Methods Fluids 16 (1993) 525e554. [30] D.C. Wilcox, Turbulence Modeling for CFD, second ed., DCW Industries, Inc., La Canada, California, 2002. [31] P.A. Durbin, T.I.-P. Shih, An overview of turbulence modeling, in: B. Sunden, M. Faghri (Eds.), Modeling and Simulation of Turbulent Heat Transfer, WIT Press, Southampton, UK, 2005, pp. 3e31. [32] P.R. Spalart, S.R. Allmaras, One-equation Turbulence Model for Aerodynamic Flows, 1992. AIAA Paper-92-0439. [33] P.A. Durbin, Separated flow components with k-ε-v2 model, AIAA J. 34 (1995) 659e664. [34] F.R. Menter, Zonal Two-Equation k-u Models for Aerodynamic Flows, 1993. AIAA Paper 93-2906. [35] B.E. Launder, On the computation of convective heat transfer in complex turbulent flows, ASME J. Heat Transfer 110 (1988) 1112e1128. [36] P.R. Spalart, M. Stretlets, S.R. Allmaras, Comments on the feasibility of LES for wings and the hybrid RANS/LES approach, advances in DNS/LES, in: Proceedings of the First AFOSR International Conference on DNS/LES, 1997.



CHAPTER 10 Modeling approaches for fuel cells

[37] P.R. Spalart, V. Venkatakrishnan, On the role of and challenges of CFD in aerospace industry, Aeronaut. J. 120 (2016) 209e232. [38] C.Y. Wang, P. Cheng, Multiphase flow and heat transfer in porous media, Adv. Heat Tran. 30 (1997) 93e196. [39] Y. Wang, Porous media flow fields for polymer electrolyte fuel cells: II. Analysis of channel two-phase flow, J. Electrochem. Soc. 156 (2009) B1134eB1141. [40] B. Sunden, Multiscale modeling approaches of transport phenomena in fuel cells, IHTC15-KN26, in: 15th International Heat Transfer Conference, Kyoto, Japan, 2014. [41] B. Sunden, On Computational procedures of transport phenomena in equipment of energy systems, CHT-15, in: ICHMT Int. Symp. on Advances in Computational Heat Transfer, 2015. CHT-15-188. [42] Y.F. Wang, J. Yuan, B. Sunden, CG-MD investigation of nanostructures and thermal properties of porous anode for SOFC, J. Power Sources 254 (2014) 209e217. [43] Y. Xiao, J. Yuan, B. Sunden, Modeling of micro/meso-scale reactive transport phenomena in catalyst layers of proton exchange membrane fuel cells, Int. J. Low Carbon Technol. 7 (No. 4) (2012) 270e287. [44] Y. Xiao, J. Yuan, B. Sunden, Process based large scale molecular dynamics simulation of a fuel cell catalyst layer, J. Electrochem. Soc. 159 (No. 3) (2012) B251eB258. [45] S. Li, B. Sunden, Effects of agglomerate model parameters on transport characterization and performance of PEM fuel cells, Int. J. Hydrogen Energy 43 (2018) 16279e16292. [46] S. Li, Modeling and Simulations of Multiphysics Phenomena and Performance in Proton Exchange Membrane Fuel Cells (Ph.D. thesis), Lund University, 2018. [47] Y.B. Zhou, K. Jiao, Q. Du, Y. Yin, X.G. Li, Gas diffusion layer deformation and its effect on the transport characteristics and performance of proton exchange membrane fuel cell, Int. J. Hydrogen Energy 39 (2013) 12891e12903. [48] S.A. Li, J. Yuan, M. Andersson, G.N. Xie, B. Sunden, Wavy surface cathode gas flow channel effects on transport processes in a proton exchange membrane fuel cell, ASME J. Electrochem. Energy Conv. Storage 14 (2017). Paper no 031007. [49] S.A. Li, J. Yuan, M. Andersson, G.N. Xie, B. Sunden, Influence of anisotropic gas diffusion layers on transport phenomena in a proton exchange membrane fuel cell, Int. J. Energy Res. 41 (2017) 2034e2050. [50] K. Jiao, X.G. Li, Water transport in polymer electrolyte membrane fuel cells, Prog. Energy Combust. Sci. 37 (2011) 221e291. [51] L. Wang, A. Husar, T.H. Zhou, H.T. Liu, A parametric study of PEM fuel cell performance, Int. J. Hydrogen Energy 28 (2003) 1263e1272. [52] A. Parthasarathy, S. Srinivasan, A.J. Appleby, C.R. Martin, Temperature dependence of the electrode kinetics of oxygen reduction at the Platinum/Nafion interface-a microelectrode investigation, J. Electrochem. Soc. 139 (1992) 2530e2537. [53] N.K.H. Dalasm, M.J. Kermani, D.G. Moghaddam, J.M. Stockie, A parametric study of cathode catalyst layer structural parameters on the performance of a PEM fuel cell, Int. J. Hydrogen Energy 35 (2010) 2417e2427. [54] D.M. Bernardi, M.W. Verbrugge, A mathematical model of the solid polymer electrolyte fuel cell, J. Electrochem. Soc. 139 (1992) 2477e2491. [55] M.J. Lampinen, M. Fomino, Analysis of free energy and entropy changes for half cell reactions, J. Electrochem. Soc. 140 (12) (1993) 3537e3546. [56] X.G. Yang, Q. Ye, P. Cheng, Matching of water and temperature fields in proton exchange membrane fuel cells with non-uniform distributions, Int. J. Hydrogen Energy 36 (2011) 12524e12537.


[57] T. Berning, D.M. Lu, N. Djilali, Three-dimensional computational analysis of transport phenomena in a PEM fuel cell, J. Power Sources 106 (2002) 284e294. [58] H. Meng, Numerical studies of liquid water behaviors in PEM fuel cell cathode considering transport across different porous media, Int. J. Hydrogen Energy 35 (2010) 5569e5579. [59] J. Wang, J. Yuan, J.S. Yu, B. Sunden, Investigation of effects of non-homogeneous deformation of gas diffusion layer in a PEM fuel cell, Int. J. Energy Res. 41 (2017) 2121e2137. [60] T.E. Springer, T.A. Zawodzinski, S. Gottesfeld, Polymer electrolyte fuel cell model, J. Electrochem. Soc. 138 (8) (1991) 2334e2342. [61] T.A. Zawodzinski, T.E. Springer, F. Uribe, S. Gottesfeld, Characterization of polymer electrolytes for fuel cell applications, Solid State Ionics 60 (1993) 199e211. [62] A.J. Wilke, P. Chang, Correlation of diffusion coefficients in dilute solutions, AIChE J. 1 (1995) 264e270. [63] J.T. Gostick, M.W. Fowler, M.D. Pritzker, M.A. Ioannidis, L.M. Behra, In-plane and through-plane gas permeability of carbon fiber electrode backing layers, J. Power Sources 162 (2006) 228e238. [64] M. Secanell, K. Karan, A. Suleman, N. Djilali, Multi-variable optimization of PEMFC cathodes using an agglomerate model, Electrochim. Acta 52 (2007) 6318e6337. [65] A.A. Shah, G.S. Kim, P.C. Sui, D. Harvey, Transient non-isothermal model of a polymer electrolyte fuel cell, J. Power Sources 163 (2007) 793e806. [66] W. Sun, B. Peppley, K. Karan, An improved two-dimensional agglomerate cathode model to study the influence of catalyst layer structural parameters, Electrochim. Acta 50 (2005) 3359e3374. [67] I. Nitta, S. Karvonen, O. Himanen, M. Mikkola, Modelling the effect of inhomogeneous compression of GDL on local transport phenomena in a PEM fuel cell, Fuel Cells 8 (2008) 410e421. [68] Q.G. Yan, H. Toghiani, H. Causey, Steady state and dynamic performance of proton exchange membrane fuel cells under various operating conditions and load changes, J. Power Sources 161 (2006) 492e502. [69] A.H. Mahmoudi, A. Ramiar, Q. Esmaili, Effect of inhomogeneous compression of gas diffusion layer on the performance of PEMFC with interdigitated flow field, Energy Conv. Manag. 110 (2016) 78e89. [70] W.S. He, J.S. Yi, T.V. Nguyen, Two-phase flow model of the cathode of PEM fuel cells using interdigitated flow fields, AIChE J. 46 (2000) 2053e2064. [71] S. Kakac, A. Pramuanjaroenkij, Y.Y. Zhou, A review of numerical modeling of solidoxide fuel cells, Int. J. Hydrogen Energy 32 (2007) 761e786. [72] B. Sunden, M. Faghri (Eds.), Transport Phenomena in Fuel Cells, WIT Press, UK, 2005. [73] J. Yuan, B. Sunden, Analysis of chemically reacting transport phenomena in an anode duct of intermediate temperature SOFCs, ASME J. Fuel Cell Sci. Technol 3 (2006) 89e98. [74] J.M. Wernik, S.A. Meguid, Coupling atomistics and continuum in solids: status, prospects and challenges, Int. J. Mech. Mater. Des. 5 (2009) 79e110. [75] H. Lu, N.P. Daphalaputkar, B. Wang, S. Roy, R. Komandun, Multiscale simulation from atomistic to continuum-coupling molecular dynamics (MD) with the material point methods (MPM), Phil. Mag. 86 (2006) 2971e2994. [76] S.R. Klaus, C. Qiang, Parallel cylindrical water nanochannels in Nafion fuel cell membranes, Nat. Mater. 7 (2008) 75e83.



CHAPTER 10 Modeling approaches for fuel cells

[77] J. Yuan, B. Sunden, On continuum models for heat transfer in micro/nano scale porous structures relevant for fuel cells, Int. J. Heat Mass Transf. 58 (2013) 441e456. [78] J. Yuan, B. Sunden, On mechanisms and models of multi-component gas diffusion in porous structures of fuel cell electrodes, Int. J. Heat Mass Transf. 69 (2014) 358e373. [79] B. Sunden, J. Yuan, Evaluation of models of the effective thermal conductivity of porous materials relevant to fuel cell electrodes, Int. J. Comput. Methods Exp. Meas. 1 (4) (2013) 440e454. [80] [81] [82] [83]