Inter-laboratory comparison of computational fluid dynamics codes for PEM fuel cell modelling.
Date: 2021-06-14
Institutions: Fuel Cells and Hydrogen 2 Joint Undertaking, Joint Research Centre
Subjects: European Joint Undertaking, benchmarking, energy research, fuel cell, research method, research report, simulation, software, testing

Inter-laboratory  comparison  of 
Computational  Fluid Dynamics  codes   
for PEM fuel cell modelling 

Deliverable: B.2.4 Final report 
on CFD benchmarking  

Bednarek, T, Tsotridis, G 

2021 

EUR 30703 EN 

 
 
 
 
 
 
 
 
 
 
This p ublication is a Technical re port by the Joint Research Ce ntre (JRC), the  European Commission’s science and knowledge service. It 
aims to p rovide evidence-based scientific support to the European p olicymaking process. The scientific output e xpressed does not imp ly a 
p olicy p osition of the European Commission. Neither the European Commission nor any p erson acting on behalf of the Commission is 
re sp onsible  for the use that might be made of this publication. For information on the me thodology and quality underlying the data used 
in this p ublication for which the source  is neither Eurostat nor other Commission service s, users should contact the re ferenced source. The 
de signations employed and the p resentation of material on the maps do not imp ly the exp ression of any opinion whatsoever on the p art 
of the  European Union concerning the le gal status of any country, territory, city or are a or of its authoritie s, or concerning the delimitation 
of its frontie rs or boundaries. 

Contact information  
Name : Tomasz Bednarek 
Addre ss: European Commission, Joint Re search Ce ntre (JRC), We sterduinweg 3, NL-1755 LE Pette n 
Email: Tomasz.Bednarek@ec.europa.eu 
Te l.: +31 224565447 

EU Science Hub 
http s://ec.europa.eu/jrc 

JRC124900 

EUR 30703 EN 

PDF 

ISBN 978-92-76-37770-2 

ISSN 1831-9424 

doi:10.2760/291664 

Luxe mbourg: Publications Office of the European Union, 2021 

© Europ ean Union, 2021 

The  re use p olicy of the European Commission is imp le mented by the Commission Decision 2011/833/EU of 12 December 2011 on the 
re use  of Commission documents (OJ L 330, 14.12.2011, p . 39). Exce pt otherwise noted, the re use of this document is authorised under 
the  Cre ative  Commons Attribution 4.0 International (CC BY 4.0) lice nce  (http s://cre ativecommons.org/licenses/by/4.0/). This me ans that 
re use  is allowed p rovided appropriate cre dit is give n and any changes are indicated. For any use or re production of photos or other 
mate rial that is not owned by the EU, permission must be sought directly from the copyright holders. 

All  conte nt © European Union, 2021 

How to cite  this re p ort: Bednarek T., Tsotridis G., Inter-laboratory comparison of Computational Fluid Dynamics code s for PEM fue l ce ll 
mode lling, EUR 30703 EN, Publications Office of the European Union, Luxembourg, 2021, ISBN 978-92-76-37770-2, 
doi:10.2760/291664 (online), JRC124900. 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  
 
 
 
Contents 

Acknowledgments........................................................................................................................................................................................................................................1 

Abstract.................................................................................................................................................................................................................................................................2 

1  Introduction...............................................................................................................................................................................................................................................3 

2  Objective and methods ..................................................................................................................................................................................................................5 

3  Physical models used for PEM fuel cell modelling............................................................................................................................................... 6 

4  The reference test cases ..............................................................................................................................................................................................................8 

5  Results of the reference test case 2 on catalyst performance................................................................................................................. 9 

6  Conclusions and outlook ........................................................................................................................................................................................................... 12 

References...................................................................................................................................................................................................................................................... 13 

List of figures ..............................................................................................................................................................................................................................................14 

List of tables................................................................................................................................................................................................................................................. 15 

Annexes.............................................................................................................................................................................................................................................................16 

Appendix A 

Details of physical models used in PEM fuel cell simulations ....................................................................16 

Appendix B 

Reference test case 1 for PEMFC single cell............................................................................................................... 27 

Appendix C 

Reference test case 2 for catalyst layer examination ....................................................................................... 37 

iii 

 
 
 
 
 
 
 
 
 
Acknowledgments 

The work described in this report was carried out under the framework contract between the Directorate-
General (DG) Joint Research Centre (JRC) of the European Commission and the Fuel Cells and Hydrogen 2 Joint 
Undertaking (FCH 2 JU). The JRC contractual activities are summarised in the FCH 2 JU  Multi-Annual Work 
Programme1. This report is the Deliverable B.2.4 of the annual rolling plan 20172. 

Authors 

Tomasz Bednarek 

Georgios Tsotridis 

1 available  online at:  
http s://www.fch.europa.eu/sites/default/file s/MAWP%20final%20version_endorsed%20GB%2015062018%20%28ID%203712421%29.p

df, p p  31.  
2 available  online at: 

http s://www.fch.europa.eu/sites/default/file s/FCH2%20JU%202017%20AWP%20and%20Budget_FINAL-20122016-
Cle an%20%28ID%202892681%29_0.pdf, pp 86. 

1 

  
      
 
 
 
 
 
 
                                     
Abstract 

An inter-laboratory comparison of Computational Fluid Dynamics (CFD) codes exercise for Polymer Electrolyte 
Membrane (PEM) fuel cell modelling was performed to assess modelling accuracy. Since PEM fuel cell models 
require a multi-physics approach involving many different phenomena, a simple comparison with experimental 
polarisation curves is not sufficient for the identification of the individual sources of errors the simulation 
software. Therefore, this report presents a methodology based on the comparison of partial simulation results. 

The report introduces first the list physical models available for the simulation of fuel cell phenomena. It 
describes then in details reference numerical test cases. Finally, it provides an example of application showing 
that by this approach, it is possible to verify any simulation software for PEM fuel cells, including commercial 
systems, without access to the source code.  

2 

  
      
 
 
1 

Introduction 

The structure of the report is the following: the objective of the study is presented in Chapter 2, the results of 
the inventory of all the physical models available is summarised in Chapter 3 and further detailed in Annex A. 
The test cases developed by JRC are explained in Chaper 4 and the results are compared and analysed in 
Chapter 5. 

The list of participants invited to contribute to the CFD codes benchmarking exercise is presented in Table 1.  

3 

  
      
 
 
 
 
Table 1 List of invited participants 

Organisation 

Contact  person(s) 

European Commission, Joint Research Centre  
(JRC) 

AVL List GmbH 

Zentrum für Sonnenenergie- und Wasserstoff-
Forschung Baden-Württemberg 
(ZSW) 

Deutsche Zentrum für Luft- und Raumfahrt 
(DLR) 

COMSOL AB 

Robert Bosch GmbH 

Georgios Tsotridis,  
Tomasz Bednarek 

Reinhard Taschl,  
Clemens Fink 

Simon Enz 

Thomas Jahnke 

Henrik Ekström 

Andreas Haeffelin 

Le Commissariat à l’énergie atomique et aux 
énergies alternatives  
(CEA) 

Arnaud Morin,  
Bolahaga Randrianarizafy 

Università degli Studi di Salerno  
(UNISA) 

Aalborg Universitet 

Cesare Pianese,  
Pierpaolo Polverino 

Torsten Berning 

4 

  
      
 
 
 
 
 
2  Objective  and  methods 

The usual way of validation of simulation software consists of its comparison with experimental data. Such a 
comparison shows whether the constitutive laws used in modelling correctly reflect reality. However, for 
complex systems such as fuel cells, the verification of the many internal quantities (i.e. partial pressures of 
reactant species, permeability of diffusion layers, local temperature, resistance of internal components etc.) is 
very challenging. This significantly hinders to analyse, understand and resolve possible divergences between 
the results of simulations and experimental data. More, in some specific cases errors due to programming 
mistakes may cancel each other leading to false conclusion that numerical model correctly reflects real 
experiment. Therefore, the first phase of a simulation software validation should be the verification of the 
applied constitutive laws and interfaces among them, so that the  fuel cell modelling code is free from 
implementation mistakes. 

The overarching objective of the inter-laboratory comparison of CFD codes in PEM fuel cell modelling is to 
establish a set of test cases for validation of the constitutive laws as implemented in the simulation codes. 
Using the proposed reference test cases, a user of simulation software will be able to validate against theory 
any modelling software for PEM fuel cells, including commercial software, even without access to their source 
code. 

In this respect, the JRC started a series of actions to facilitate this objective. The first action performed was an 
inventory of the models for the simulation of fuel cell phenomena used by the participants. A reference 
benchmarking test case 1, covering all aspects of PEM fuel cells, was then proposed by JRC. 

In a following stage, a benchmarking test case 2 was proposed as a modification of first case in order to assess 
the accuracy of the simulation of two specific contributions: the reaction kinetics and the proton conduction in 
ionomer. Therefore, the reference test case 2 validates the modelling of the catalyst layer performance, namely 
the Hydrogen Oxidation Reaction (HOR) in the anode, Oxygen Reduction Reaction (ORR) in the cathode and 
losses due to proton conduction by the ionomer in the catalyst layers (anode and cathode sides). 

The results of PEM fuel cell simulations are compared in terms of: 

—  velocity distribution in the gas channels, 

—  concentration of gas species in both anode and cathode fuel cell compartments, namely: gas channels, gas 

diffusion layers and catalyst layers, 

—  temperature distribution across the active area at catalyst layers and gas diffusion layers, both, anode and 

cathode, 

—  hydration of the membrane, 

—  protonic potential distribution in ionomer (anode and cathode catalyst layers and membrane, 

—  electronic potential at catalyst layers, gas diffusion layers and bi-polar plates, both anode and cathode, 

and  

—  in-plane current density distribution (anode and cathode). 

5 

  
      
 
 
 
 
3  Physical models used  for PEM fuel cell modelling 

The overview of all the physical models used for simulations of fuel cell operation in this benchmarking exercise 
is presented in Table 2. The physical models are grouped per fuel cell component, namely bi-polar plates (BPP), 
gas channels (GCH), gas diffusion layer (GDL), micro-porous layer (MPL), catalyst layer (CL), membrane (MEM). 
Each modelled phenomenon (mass convection, mass diffusion, electric conduction, heat transfer, proton 
conduction, gas-liquid (water vapour-water) phase change and electro-chemical reactions) is described by the 
appropriate  governing  equations  and  corresponding  source  terms  which  account  for  internal  boundary 
conditions between the fuel cell components. 

Details of the physical models used in the simulations are presented in Appendix A. 

6 

  
      
 
 
Table 2 Overview of physical models per fuel cell component used in simulations using CFD codes. Numbers in brackets correspond to equations in Appendix A.  

MEM 

ionomer 

Electron 
transport 

Proton 
transport 

– 

– 

Ele ctric 
e quilibrium of 
H+ (34) 

Ele ctric 
e quilibrium of 
H+ (34) 

– 

– 

–

- 

– 

– 

– 

Transp ort e q. 
(45) 

s
t
n
a
t
c
a
e
r

s
a
G

r
e
t
a
W

O2 

H2 

r
u
o
p
a
v

d
i

u
q
i
l

d
e
v
l
o
s
s
i
d

Anode CL 

carbon 
support 

Lap lacian 
(20) 

– 

– 

– 

– 

– 

pores 

ionomer 

– 

– 

– 

Ele ctric 
e quilibrium of 
H+ (34) 

N-S with Darcy 
(8) 
Sp e cies eq. 
(10) 

N-S with Darcy 
(8) 
Sp e cies eq. 
(10) 

N-S with Darcy 
(8) 
Sp e cies eq. 
(10) 

– 

– 

– 

Cathode CL 

carbon 
support 

Lap lacian 
(20) 

– 

– 

– 

– 

pores 

– 

– 

N-S with Darcy 
(8) 
Sp e cies eq. 
(10) 

N-S with Darcy 
(8) 
Sp e cies eq. 
(10) 

N-S with Darcy 
(8) 
Sp e cies eq. 
(10) 

Darcy (31) 

Darcy (31) 

– 

Darcy (31) 

MPL 

GDL 

pores 

fibres 

pores 

GCH 

BPP 

– 

– 

N-S with 
Darcy (8) 
Sp e cies eq. 
(10) 

N-S with 
Darcy (8) 
Sp e cies eq. 
(10) 

N-S with 
Darcy (8) 
Sp e cies eq. 
(10) 

Darcy (31) 

– 

Lap lacian 
(20) 

– 

– 

– 

– 

– 

– 

– 

– 

– 

– 

N-S with 
Darcy (8) 
Sp e cies eq. 
(10) 

N-S with 
Darcy (8) 
Sp e cies eq. 
(10) 

N-S (1) 
Sp e cies 
e quilib. (2) 

N-S (1) 
Sp e cies 
e quilib. (2) 

N-S with Darc 
(8)y 
Sp e cies eq. 
(10) 

N-S (1) 
Sp e cies 
e quilib. (2) 

Darcy (31) 

N-S (1) 

– 

– 

Lap lacian 
(4) 

– 

– 

– 

– 

– 

– 

carbon 
particles 

Lap lacian 
(20) 

– 

– 

– 

– 

– 

– 

Osmotic and 
back flow 
(37) 

Osmotic (37) 

– 

– 

Osmotic (37) 

– 

– 

Heat 
transfer 

Conduction 
(48) 

Source : JRC, 2017 

Conduction te rm in (18) 

Conve ction 
te rm in (18) 

Conduction te rm in (18) 

Conve ction 
te rm in (18) 

Conduction 
te rm in (18) 

Conve ction 
te rm in (18) 

Conduction 
te rm in (18) 

Conve ction 
te rm in (18) 

Conve ction 
(3) 

Conduction 
(5) 

7 

  
      
 
 
 
 
 
 
 
 
 
4  The reference test cases 

The inter-laboratory comparison of CFD codes effort requires the establishment of a reference computational 
geometry and input parameters. The operating conditions correspond to  those contained within the  EU 
Harmonised Operating Conditions for Automotive Applications [1]. A  typical MEA size of 25 cm2 with 
single serpentine flow field is considered. As agreed with the participants of the exercise, in an effort to simplify 
the  numerical solution and to  minimise the  effects of channel geometry, the “single serpentine cell” is 
approximated by a “single straight channel” of equivalent length. Such a channel configuration decreases the 
effects of flow field shape and simultaneously allows for comparison of results obtained from different 
numerical models with analytical solution results. The numerical values of input data correspond to real 
materials used in fuel cell testing. 

The following modelling assumptions have been made: 

—  The external electrical circuit is modelled by applying appropriate boundary conditions to anode and cathode 
bi-polar plate (BPP) terminals. The applied load corresponds to galvanostatic operation mode, namely fixed 
electric current is applied. 

—  Side walls of the numerical model are assumed as symmetric in order to reduce the effects of boundaries.  

—  Fuel (hydrogen) and oxidant (air) are considered to be ideal gases. Therefore, the diffusion of species 

follows the multicomponent principles applicable to ideal gas mixtures. 

—  Microporous and catalyst layers (MPL and CL) are modelled as single homogeneous sub-domains. The 

catalyst loading is defined by the active surface ratio [m-1]. 

—  Water condensation occurs only in the porous sub-domains (GDL and CL). Therefore, a two-phase simulation 
is performed only in GDL and CL sub-domains. The liquid water that is produced and evacuated to the gas 
channel is homogenised. The liquid water takes the velocity of the surrounding gas, hence the gas flow 
profile is not affected by the existence of water droplets.  

It should be clearly noted that in this study, as in many other similar studies, the open circuit voltage VOC is an 
input to the CFD code. This requires calculation of VOC using the Nernst equation for given operating conditions 
prior to the modelling.  

A description of the input data used in the reference test case 1 is given in Appendix B. 

The reference test case 1 is very demanding for a simulation software as in a single simulation run all crucial 
physical phenomena appear, namely: significant gradients of reactant concentrations, a relatively large volume 
fraction of condensed water and consequential gas density gradients in the anode compartment. Additionally, 
the total length of the fuel cell model (1250 mm) implies a large computational mesh, which is challenging for 
numerical solvers. 

Despite using the same input data and geometry, the results of the reference test case 1 simulation obtained 
by participants of the exercise were significantly different, hence were of very difficult analysis. Due to the 
complexity of the PEM fuel cell numerical model (Section 4) and the large-scale computational domain, it was 
impossible to distinguish differences arising from the use of different implementations of physical models and 
those, which may stem from inevitable numerical errors.  

Therefore, the reference test case 1 was abandoned, and learning from it the reference test case 2 was defined 
allowing for an assessment of the different components contribution to the final loss of modelling accuracy, 
specifically to compare and validate the contribution of the fuel cell reaction kinetics to the results of the 
simulation. In order to minimise the effect caused by the gas inlet humidity, the concentration gradients in the 
reactants and liquid water that are known to affect the reactant diffusion through the GDL, the following 
modifications are introduced: 

—  The length of the channel is reduced to 80 mm in order to reduce computational efforts. 

—  A constant flow of both reactants is set to be equivalent to stoichiometries 10.0 at 1.0 Acm-2 current density 

to obtain a uniform species distribution along with efficient removable of product water. 

—  Fully humidified reactants (100% relative humidity, RH) at the inlets to avoid membrane dry out and ensure 

a sufficiently high and uniform membrane protonic conductivity. 

The details of the reference test case 2 are presented in Appendix C.  

8 

  
      
 
5  Results of the reference test case 2 on catalyst performance 

The usual way of benchmarking CFD-code results is by using reference experimental data. In this study, we 
compare the PEM fuel cell modelling results. This approach provides a verification method of simulation codes, 
which allows for identification of mistakes or bugs in the software. The following example shows how CFD inter-
laboratory comparison exercises helped to improve the modelling software of two participants in the project. 

There are two participants contributing to the exercise with the reference test case 2. Since both participants 
use the same constitutive laws and the same input data, it is expected that the differences in obtained results 
come only from numerical errors (rounding and truncation). Therefore, the obtained results should match each 
other with rather high accuracy. Polarisation I-V curves provided by two participants obtained according to the 
conditions specified in Appendix C, are presented in Figure 1.  

At first glance, it could be concluded from Figure 1A that the I-V curves match each other very well and that 
any differences might be due to numerical errors. However, a close look in Figure 1B shows that there are 
differences reaching 35 mV (2nd important digit) at current density 0.4 A cm-2 between the simulation results 
provided by Participant 1 and Participant 2. It is seen also, that the difference in results increases with increasing 
current density, indicating that reasons of diverging results are not due to typical numerical errors (rounding or 
truncation), but very likely there are some mistakes in implementation of constitutive laws. 

Figure 1 (A) P olarisation curves provided by benchmarking exercise participants according to conditions specified in 
Appendix C. (B) Zoom plot shows low current density region indicated in (A). 

A) 

B)

Source : JRC, 2017 

To  identify the reason for the discrepancy shown in Figure 1B, the polarisation curve was divided into its 
contributions,  namely:  equilibrium  voltage,  activation  overpotential,  ionic  overpotential  and  electronic 
overpotential (Figure 2). The calculated overpotentials were then compared separately.  

A comparison of the activation, ionic and electronic overpotentials is presented in Figure 3. 

9 

  
      
 
 
 
 
 
 
 
 
 
 
Figure 2 Split of polarisation I-V curve provided by P articipant 1 into overpotentials. 

Source : JRC, 2017 

Figure 3 Comparison of simulation results: overpotentials of the polarisation I-V curves presented in Figure 1. 

Source : JRC, 2017 

10 

  
      
 
 
 
 
 
It can be seen in Figure 3 that the activation and electronic overpotentials provided by the participants match 
each other well. However, there is a significant difference in the ionic overpotential. The ionic overpotential 
contains  contributions  arising  from  both  the  anode  and  cathode  catalyst  layers.  Therefore,  the  ionic 
overpotential was split into anode and cathode ionic potentials (potential difference between MEM and GDL) as 
presented in Figure 4.  

Figure 4 Results of simulations provided by participants: (A) anode and (B) cathode ionic overpotentials. 

A) 

B)

Source : JRC, 2017 

In Figure 4 can be seen that differences appear to be only due to the cathode ionic potential. Since the deviation 
is proportional to the current density (Figure 4B) this suggests that the calculation of the ionic conductivity in 
the cathode CL is incorrect, corresponding to Eq. (37), see Appendix A. 

In fact, both participants (Participant 1 and Participant 2) have found a software bug in their codes related to 
Eq. (37) in Appendix A. Participant 2 has identified the bug and corrected the mistake, however Participant 1 
was not able to correct the software bug because it did not have access to the software source code. 

Figure 5 Result of simulations before and after the software update by P articipant 2.  

Source : JRC, 2017 

11 

  
      
 
 
 
 
 
 
 
  
 
6  Conclusions  and  outlook 

The purpose of the present study was to develop a number of test cases in order to compare the results of 
different numerical codes, modelling the performance of a single PEM fuel cell. The results were obtained when 
using the same input physical values, geometry and boundary conditions. To characterise the performance of 
the PEM fuel cell, a polarisation curve test was chosen. The numerical codes were required to evaluate the 
polarisation curve using the reference test case 1, integrating all fuel cell related phenomena. 

The results of the reference test case 1 provided by participants of the exercise were significantly different. Due 
to the complexity of the numerical model it was impossible to identify the differences arising from the employed 
physical models from those associated numerical errors. Therefore, a reference test case 2 was proposed to 
evaluate the contributions of reaction kinetics, minimising the effects of reactants delivery, temperature, 
pressure and material properties of the fuel cell components.  

The evaluation of the modelling results of the reference test case 2 led to the discovery of a number of mistakes 
in  the codes. Hence, the verification of the software was performed in details, accounting internal model 
parameters that could not be measured during real experiment. As a result, the implementation of the physical 
model into simulation software becomes more trustworthy and can be next validated against real experiments.  

The work  described in this report is  just as a  first  step, focussing on the reaction kinetics region of the 
polarisation curve. It should to be completed by similar work dedicated to the other regions of interest, namely 
those limited by ohmic and mass transport. For the ohmic region, additional reference test cases should be 
established to examine the effects of membrane hydration on its protonic resistance, as well as the effects of 
the electronic conductivity of fuel cell components and contact resistances. For examining the third region of 
the polarisation curve where mass transport limitation dominate, another reference test case should investigate 
losses due to reactant diffusion through GDL thickness and the effects of liquid water accumulation.  

12 

  
      
 
 
 
References 

[1] 

[2] 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 

E.L. Cussler, Diffusion. Mass Transfer in Fluid Systems, Cambridge University Press, Cambridge, 2009. 

A.Z.  Weber, R.L. Borup, R.M. Darling, P.K. Das, T.J. Dursch, W. Gu, D. Harvey, A. Kusoglu, S. Litster, M.M. 
Mench, R. Mukundan, J.P. Owejan, J.G. Pharoah, M. Secanell, I. V Zenyuk, J. Electrochem. Soc.  161 (2014) 
F1254. 

R.F. Mann, J.C. Amphlett, B.A. Peppley, C.P. Thurgood, J. Power Sources 161 (2006) 775. 

H. Wu, X. Li, P. Berg, Electrochim. Acta 54 (2009) 6913. 

T.E. Springer, M.S. Wilson, S. Gottesfeld, J. Electrochem. Soc.  140 (1993) 3513. 

A.Z. Weber, J. Newman, Adv. Fuel Cells 1 (2007) 47. 

T.E. Springer, T.A. Zawodzinski, S. Gottesfeld, J. Electrochem. Soc.  138 (1991) 2334. 

T. Berning, M. Odgaard, S.K. Kær, J. Power Sources 196 (2011) 6305. 

13 

  
      
 
 
 
List of figures 

Figure 1 (A) Polarisation curves provided by benchmarking exercise participants according to conditions 
specified in Appendix C. (B) Zoom plot shows low current density region indicated in (A)........................................................ 9 

Figure 2 Split of polarisation I-V curve provided by Participant 1 into overpotentials............................................................10 

Figure 3 Comparison of simulation results: overpotentials of the polarisation I-V curves presented in Figure 1.
 ...................................................................................................................................................................................................................................................................................10 

Figure 4 Results of simulations provided by participants: (A) anode and (B) cathode ionic overpotentials......... 11 

Figure 5 Result of simulations before and after the software update by Participant 2. ....................................................... 11 

14 

  
      
 
 
List of tables 

Table 1 List of invited participants............................................................................................................................................................................................. 4 

Table 2 Overview of physical models per fuel cell component used in simulations using CFD codes. Numbers 
in brackets correspond to equations in Appendix A. ......................................................................................................................................................7 

15 

  
      
 
 
Annexes 

Appendix  A  Details of physical  models used in PEM fuel cell simulations   

Table 3 List of used symbols. 

a 

b 

Ci 

𝑒𝑒𝑒𝑒𝑒𝑒
𝐷𝐷𝑘𝑘
Dgl 

Di 

Dij 

Dλ 

e- 

Ei 

Water activity 

User defined exponent  

Mass concentration of specie i 

Effective diffusivity of specie k in porous media  

Empirical diffusion coefficient for condensation and evaporation of water 

Diffusion coefficient of specie i 

Binary diffusion coefficient 

Diffusion coefficient for water dissolved in ionomer 

Electron 

Activation energy 

EW 

Membrane equivalent weight 

F 

h 

i, i 

i

p
J 

jex 

jref 

K 

k 

Kr 

M 

N 

n 

Faraday constant 

Enthalpy 

Current and its magnitude 

Protonic current  

Volumetric current density 

Exchange current density 

Reference current density 

Permeability  

Effective thermal conductivity 

Relative permeability due to water saturation effect 

Molar mass 

Number of considered species 

Number of electrons taking part in chemical reaction 

16 

  
      
 
 
 
 
 
nd 

p 

p0 

pc 

pi 

pl 

Dimensionless osmotic drag coefficient 

Pressure 

Atmospheric pressure 

Capillary pressure 

Partial pressure of specie i 

Pressure of liquid phase 

psat 

Saturation pressure 

R 

s 

Se 

Universal gas constant 

Liquid water saturation 

Volumetric source term, electric current source 

Sgl, Sgd, Sld 

Mass exchange source tern between different phased of water: gas-liquid, gas-
dissolved in ionomer, liquid-dissolved in ionomer  

Sh 

Sλ 

𝑘𝑘
𝑆𝑆𝑠𝑠
Su 

t 

T 

Volumetric source term, Joule heating 

Volumetric source term, water dissolved in ionomer 

Volumetric source term, source of specie k 

Volumetric source term, body force 

Time 

Temperature  

Half potentials for anode and cathode catalysts 

𝑎𝑎𝑎𝑎,𝑐𝑐𝑎𝑎𝑐𝑐
v 

𝑈𝑈0

Flow velocity 

V 

Diffusion velocity 

x, y, z 

Spatial dimensions 

Xi 

Mole fraction of specie i 

Reference mole fraction of specie i 

𝑟𝑟𝑒𝑒𝑒𝑒

𝑋𝑋𝑖𝑖

Anode and cathode respectively, electron exchange coefficient 

𝛼𝛼𝑎𝑎𝑎𝑎, 𝛼𝛼𝑐𝑐𝑎𝑎𝑐𝑐
β 

Reaction symmetric factor 

ε 

Porosity  

17 

  
      
 
 
 
 
 
 
 
ϵ 

φ 

Γ 

𝛾𝛾𝑎𝑎𝑎𝑎,𝛾𝛾𝑐𝑐𝑎𝑎𝑐𝑐
γ 

γgd, γld 

η 

λ 

Volume fraction of ionomer 

Electric potential 

Factor accounting volume fraction and tortuosity of ionomer 

Anode and cathode, respectively, reactant concentration exponents 

Average surface of a liquid water droplet 

Water  mass  exchange  rate  constants  between  phases:  gas-dissolved  in 
ionomer and liquid-dissolved in ionomer respectively  

Overpotential  

Membrane water content 

Dynamic viscosity, average dynamic viscosity 

𝜇𝜇, 𝜇𝜇̅
θ 

User defined constant 

Density, average density 

Effective electric conductivity 

Electrocatalyst surface area, ECSA 

𝜌𝜌, 𝜌𝜌̅

𝜎𝜎𝑒𝑒𝑒𝑒𝑒𝑒

𝜁𝜁

Gas channels 

A Newtonian fluid is considered for mass transport modelling. The liquid water is accounted as a homogenous 
fine mist in a gas phase. Henceforth, the Navier-Stokes (N-S) equation (momentum diffusion balance) for 
compressible fluids without body forces is used to describe the transport of gas and liquid phases in the gas 
channels, namely: 

, 

(1) 

𝜕𝜕
𝜕𝜕𝑐𝑐 (𝜌𝜌̅𝐯𝐯) + ∇ ⋅ (𝜌𝜌̅𝐯𝐯𝐯𝐯) = ∇ ⋅ (𝜇𝜇̅ 𝐯𝐯) − ∇𝑝𝑝
where v is velocity vector of the respective liquid-gas mixture, p is absolute pressure of the respective liquid-
), 
gas mixture, 
𝜕𝜕

 is the gradient operator of the canonical space variables (

 is the time derivative operator, 

𝜕𝜕𝑐𝑐

 , superscript T denotes transpose and 
 is the divergence operator. Average 
𝑇𝑇
  respectively, of the liquid-gas mixture (humid air and liquid 
density and average dynamic viscosity, 
∇= �
water at the cathode; humid hydrogen and liquid water at the anode) is calculated based on liquid water volume 
fraction s. Namely, average density and average viscosity are given as 

∇
  and 

𝜕𝜕
𝜕𝜕𝜕𝜕 +

𝜕𝜕
𝜕𝜕𝑦𝑦 +

𝜕𝜕
𝜕𝜕𝜕𝜕�

𝜕𝜕
𝜕𝜕𝜕𝜕 ,

𝜕𝜕
𝑦𝑦𝜕𝜕 ,

  and 

∇ ⋅ =

x,y,z

𝜕𝜕𝜕𝜕

𝜌𝜌̅

𝜇𝜇̅

𝜕𝜕

  respectively, where subscripts l and g denote liquid water and gas (air in cathode or hydrogen in 
𝜇𝜇̅ = 𝑠𝑠𝜇𝜇𝑙𝑙 +

anode) respectively.  
(1 − 𝑠𝑠)𝜇𝜇𝑔𝑔
In particular, the proposed homogenized two-phase approach deals with the gas-liquid mixture as an single 
.  No phase change in the gas channels is 
phase fluid, hence its continuity is ensured by 
considered. 

𝜌𝜌̅ = 𝑠𝑠𝜌𝜌𝑙𝑙 + (1 − 𝑠𝑠)𝜌𝜌𝑔𝑔

𝜕𝜕𝜌𝜌�
𝜕𝜕𝑐𝑐 + ∇ ⋅ (ρ� 𝐯𝐯) = 0

Since many different species are involved in diffusion, the Fickian diffusion model, where diffusion flux of 
species is proportional to their concentration gradient neglecting interaction between species, is no longer 
applicable. Therefore the Maxwell-Stefan equations [1] are applied, namely 

18 

  
      
 
 
 
 
 
 
 
(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

𝑁𝑁

𝑋𝑋𝑘𝑘𝑋𝑋𝑙𝑙
𝐷𝐷𝑘𝑘𝑙𝑙 (𝐕𝐕𝑙𝑙 − 𝐕𝐕𝑘𝑘) = ∇𝑋𝑋𝑘𝑘

�
k=1
l≠k 
where 
 denotes number of considered species (oxygen and water vapour in cathode; hydrogen and water 
vapour in anode) Xk,l are molar fractions of species k and l respectively, V is diffusion velocity of considered 
 is presented in [1] depending on molar 
species. The procedure of calculating binary diffusion coefficients 
fraction and molar mass of considered species.  

N

Heat transfer is modelled using the general diffusion equation. Due to relatively low heat conductivity of air 
and dominating convection, heat conduction in the gas channels is neglected. There is no heat source term for 
the gas channels, therefore: 

𝐷𝐷𝑘𝑘𝑙𝑙

𝜕𝜕
�
where h is sensible enthalpy.  
𝜕𝜕𝑐𝑐

+ ∇ ⋅ 𝐯𝐯� 𝜌𝜌̅ℎ = 0 

Bi-polar plates (BPP) 

The electric charge (electrons) in BPP is transferred between the terminal (connector to the external circuit) and 
BPP – GDL contact surface. The electric current model takes the time dependent Laplacian form: 

, 

BPP

𝜕𝜕
BPP
𝜕𝜕𝑐𝑐 + 𝜎𝜎eff
�
where 

2
 ∇

�𝜙𝜙

= 0 
 is Laplacian operator given as 
2

2

2

. 

2

𝜕𝜕
The electric potential and effective electrical conductivity in the bi-polar plate are denoted by 
𝜕𝜕𝜕𝜕
respectively. There are no electric charge sources (S=0) or convective terms (v

0) in BPP. 

2 +

2 +

=

𝜕𝜕𝑦𝑦

∇

∇

𝜕𝜕𝜕𝜕

𝜕𝜕

𝜕𝜕

2

2

The heat transfer in the BPP is modelled by the diffusion equation as follows: 
=

𝜙𝜙

BPP

 and 

, 

BPP
𝜎𝜎eff

BPP

𝐵𝐵𝐵𝐵𝐵𝐵

𝐵𝐵𝐵𝐵𝐵𝐵

𝐵𝐵𝐵𝐵𝐵𝐵

 and k

ℎ) + 𝑘𝑘

2
𝑇𝑇 = 𝑆𝑆ℎ

𝜕𝜕
(𝜌𝜌
where ρ
𝜕𝜕𝑐𝑐
respectively; h is enthalpy of BPP and 
BPP
not dependent on spatial geometry, the k
∇
operator (
electric current and electrical conductivity as follows: 
𝐵𝐵𝐵𝐵𝐵𝐵
𝑆𝑆ℎ

2) appears. The volumetric source term 

∇
 are physical parameters of the BPP material, namely density and thermal conductivity, 
T is the temperature gradient in BPP. Since the thermal conductivity is 
 term is excluded from the divergence operator, hence the Laplacian 
 represents the Joule heating and depends on the local 

BPP

∇

2

𝐵𝐵𝐵𝐵𝐵𝐵

𝑖𝑖
𝐵𝐵𝐵𝐵𝐵𝐵 
𝜎𝜎𝑒𝑒𝑒𝑒𝑒𝑒

𝑆𝑆ℎ
where i is the magnitude of the local electric current i, calculated according to Ohm's law: 

=

. 

𝐵𝐵𝐵𝐵𝐵𝐵
𝐢𝐢 = 𝜎𝜎𝑒𝑒𝑒𝑒𝑒𝑒

𝐵𝐵𝐵𝐵𝐵𝐵
∇𝜙𝜙𝑒𝑒𝑒𝑒𝑒𝑒

Gas diffusion, microporous and catalyst layers (GDL, MPL and CL) 

The CFD approach does not allow for exact modelling of the solid structure of porous layers – carbon fibres for 
GDL, or noble metal catalyst (on carbon support) particles for MPL and CL. According to the Darcy approximation, 
, permeability of the porous media K 
where the flow velocity v is directly proportional to pressure gradient 
and inversely proportional to fluid dynamic viscosity µ and material porosity ε, namely 

∇p

19 

  
      
 
 
 
 
 
 
 
 
 
 
 . 

𝐾𝐾
𝜇𝜇𝜇𝜇 ∇𝑝𝑝

𝐯𝐯 = −
The solid structure of porous media is assumed homogenous and described by a scalar parameter – namely 
the average porosity ε.  The appropriate transport parameters are expressed as functions of the average 
porosity. The superscripts

GDL and CL indicate the parameters corresponding to GDL and CL, respectively. 

Using the Darcy modification, the N-S equation (1) reads  

𝜕𝜕
𝜕𝜕𝑐𝑐 (𝜌𝜌̅𝜀𝜀𝐯𝐯) + ∇ ⋅ (𝜌𝜌̅𝜀𝜀𝐯𝐯𝐯𝐯) = ∇ ⋅ (𝜇𝜇̅𝜀𝜀𝐯𝐯) − 𝜀𝜀∇𝑝𝑝 + 𝑆𝑆𝑢𝑢
where Su is a body force source term considering the viscous resistance of the fluid flowing through the porous 
layer: 

2

𝜇𝜇�
𝐾𝐾 𝜀𝜀

𝑆𝑆𝑢𝑢 = −
where K is the permeability of the porous layer.  

𝐯𝐯

The conservation of gas species in GDL, MPL and CL is modified by homogenous porosity 

𝜕𝜕(𝜀𝜀𝜀𝜀𝑘𝑘)
𝜕𝜕𝑐𝑐
where 

𝑒𝑒𝑒𝑒𝑒𝑒

𝑘𝑘
∇Ck) = 𝑆𝑆𝑠𝑠

+ ∇ ⋅ (𝜀𝜀𝜀𝜀𝑘𝑘𝐯𝐯) + 𝛻𝛻 ⋅ (𝐷𝐷𝑘𝑘
 is a source of species k. 
𝑘𝑘
𝑆𝑆𝑠𝑠

𝑒𝑒𝑒𝑒𝑒𝑒

𝐷𝐷𝑘𝑘

 is the effective diffusion coefficient in porous material: 

𝑒𝑒𝑒𝑒𝑒𝑒

ε
τ

=

𝐷𝐷𝑘𝑘  
𝐷𝐷𝑘𝑘
where τ is tortuosity. However, if the tortuosity factor is not known, it is possible to apply the Bruggemann 
approximation: 

𝑒𝑒𝑒𝑒𝑒𝑒

1.5

𝐷𝐷𝑘𝑘
where Dk is the gas species diffusivity computed according to multicomponent diffusion [1].  

= 𝜀𝜀

𝐷𝐷𝑘𝑘 

The sink terms for reactant species, 
(O2) are consumed. The sink terms for H2 and O2 are given as: 

 , take non-zero values only in the CL, where hydrogen (H2) and oxygen 
𝑘𝑘
𝑆𝑆𝑠𝑠

𝑂𝑂2
𝑆𝑆𝑠𝑠

= −

𝑀𝑀𝑂𝑂2
4𝐹𝐹

𝑖𝑖𝑐𝑐

𝑀𝑀𝐻𝐻2
2𝐹𝐹

𝑖𝑖𝑎𝑎

=  −

𝐻𝐻2
𝑆𝑆𝑠𝑠
where F is the Faraday constant, Mk is molecular weight of species k, and ia and ic are anode and cathode 
currents, respectively. The negative term on the right hand side of equations (14) and (15) denotes consumption 
of reactants. The constants "2" and "4" in the denominators of equations (14) and (15) represent the number of 
electrons taking part in the redox reactions: Hydrogen Oxidation Reaction (HOR) in anode and Oxygen Reduction 
Reaction in cathode, namely 

HOR:

𝐵𝐵𝑐𝑐

   𝐻𝐻2     
𝐵𝐵𝑐𝑐
1
2 𝑂𝑂2   

→  2𝐻𝐻

1
2 (2𝑂𝑂

+1

−2

+ 2𝑒𝑒

−

−

ORR:
where superscripts correspond to electric charge and e- is electron (negative charge). 

− 4𝑒𝑒

→ 

)

(8) 

(9) 

(10) 

(11) 

(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

20 

  
      
 
 
  
  
 
 
 
 
 
 
  
  
Water is produced in dissolved form therefore the water vapour source term corresponds to mass exchange 
between water dissolved in the ionomer and water vapour and evaporation/condensation terms (Sgd and Sgl, 
respectively). 

(18) 

𝐻𝐻2𝑂𝑂

=  −𝑆𝑆𝑔𝑔𝑔𝑔 − 𝑆𝑆𝑔𝑔𝑙𝑙

𝑆𝑆𝑠𝑠
The term Sgd  corresponding to the mass exchange between water vapour and dissolved water is given in 
equation (42). The evaporation/condensation term is given as: 

𝛾𝛾𝑔𝑔𝑙𝑙𝜀𝜀𝑠𝑠𝐷𝐷𝑔𝑔𝑙𝑙 𝑝𝑝 ln�

for 𝑋𝑋𝐻𝐻2𝑂𝑂 𝑝𝑝 ≤ 𝑝𝑝𝑠𝑠𝑎𝑎𝑐𝑐

(evaporation)

(19) 

𝑆𝑆𝑔𝑔𝑙𝑙 =

⎧
⎪

𝑀𝑀𝐻𝐻2𝑂𝑂
𝑅𝑅𝑇𝑇

𝑀𝑀𝐻𝐻2𝑂𝑂
𝑅𝑅𝑇𝑇

⎨
⎪
⎩

� 

𝑝𝑝 − 𝑝𝑝𝑠𝑠𝑎𝑎𝑐𝑐
(1 − 𝑋𝑋𝐻𝐻2𝑂𝑂)𝑝𝑝
𝑝𝑝 − 𝑝𝑝𝑠𝑠𝑎𝑎𝑐𝑐
(1 − 𝑋𝑋𝐻𝐻2𝑂𝑂)𝑝𝑝

𝛾𝛾𝑔𝑔𝑙𝑙𝜀𝜀(1 − 𝑠𝑠)𝐷𝐷𝑔𝑔𝑙𝑙  𝑝𝑝 ln�

� for 𝑋𝑋𝐻𝐻2𝑂𝑂 𝑝𝑝 > 𝑝𝑝𝑠𝑠𝑎𝑎𝑐𝑐
where 
 is mole fraction of water vapour, p is absolute gas pressure and psat is saturation vapour pressure, 
γgl is average surface of liquid water droplet and the Dgl term derives from an empirical formula depending on 
the gas (air or hydrogen) as follows: 

(condensation)

𝑋𝑋𝐻𝐻2𝑂𝑂

⎧0.365 − 4 �
⎪

1.79 − 4 �

𝐷𝐷𝑔𝑔𝑙𝑙 =

⎨
⎪
⎩

2.334

𝑇𝑇
�
343
2.334
𝑇𝑇
�
343

𝑝𝑝0
�
𝑝𝑝
𝑝𝑝0
𝑝𝑝

�

�

�  

for cathode

for anode

where T is gas temperature and p0 is a reference pressure equal to an atmospheric pressure of 0.1MPa. 

The heat transfer in GDLs and CLs is modelled by the diffusion equation modified for a homogenized porous 
medium: 

𝐺𝐺𝐺𝐺𝐺𝐺,𝐶𝐶𝐺𝐺
𝜕𝜕
solid structure + 𝑆𝑆ℎ  
(𝜀𝜀𝜌𝜌̅ℎ)
�������
𝑇𝑇
∇
�����fluid
𝜕𝜕𝑐𝑐
•  where ρGDL,CL is an average (accounting volume of pores) density of the porous layer. Sh is the heat 

𝜕𝜕
�ρGDL,CLℎ�
���������
𝜕𝜕𝑐𝑐
solid structure

+ ∇ ⋅ ((1 − 𝜀𝜀)𝜌𝜌̅𝐯𝐯ℎ)
�����������
fluid

= 𝑘𝑘

+

2

source term, which takes into account the Joule heat produced by transfer of current; 

the heat of water condensation/evaporation; and 

the reaction heat due to the hydrogen oxidation reaction in the cathode catalyst layer. 

2

𝑖𝑖
𝐺𝐺𝐺𝐺𝐺𝐺 + ℎ𝑝𝑝ℎ𝑎𝑎𝑠𝑠𝑒𝑒 
𝜎𝜎𝑒𝑒𝑒𝑒𝑒𝑒

⎧
⎪

2

for GDL

• 

• 

𝑆𝑆ℎ =

𝑖𝑖
𝐶𝐶𝐺𝐺 + ℎ𝑝𝑝ℎ𝑎𝑎𝑠𝑠𝑒𝑒 + ℎ𝑟𝑟𝑒𝑒𝑎𝑎𝑐𝑐𝑐𝑐 − 𝑖𝑖𝑖𝑖
𝜎𝜎𝑒𝑒𝑒𝑒𝑒𝑒

⎨
⎪
⎩

for CL

where hphase is enthalpy change due to phase change (water condensation/evaporation), 
 are 
electrical conductivities of gas diffusion and catalyst layers, respectively, and kGDL,CL is thermal conductivity of 
GDL  or CL. In the cathode catalyst layer, there is another source term hreact corresponding to enthalpy change 
due to OOR while η is the overpotential. 

𝐺𝐺𝐺𝐺𝐺𝐺
𝜎𝜎𝑒𝑒𝑒𝑒𝑒𝑒

𝐶𝐶𝐺𝐺
𝜎𝜎𝑒𝑒𝑒𝑒𝑒𝑒

 and 

The electric charge (electrons) is transferred between the CLs, where the electrochemical reactions occur on the 
triple phase boundary and the electrical terminals of the bi-polar plate. The electric current flow is modelled 
taking into account the homogenized porosity of both gas diffusion and catalyst layers, as follows: 

GDL,CL

+ 𝜎𝜎eff

GDL,CL

2
∇

�𝛷𝛷

= 𝑆𝑆𝑒𝑒 

𝜕𝜕
�
𝜕𝜕𝑐𝑐

21 

(20) 

(21) 

(22) 

(23) 

  
      
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
The source term Se takes non-zero values only in the catalyst layer (where electrochemical reactions occur) and 
is equal to the generated electric current: 

0 
𝑆𝑆𝑒𝑒 = �
𝐽𝐽𝑎𝑎𝑎𝑎,𝑐𝑐𝑎𝑎𝑐𝑐
where J is volumetric current density due to HOR or ORR. Subscripts an and cat correspond to the anode and 
cathode, respectively.  

for GDL and MPL
for CL

The difference between the phase potential of the solid BPP terminal (external boundary condition) and the 
phase potential of the membrane is the main driving force for the electrochemical reactions HOR and ORR.  

The anode and cathode source terms Jan,cat are given in general form by the well-known Butler-Volmer (B-V) 
reaction kinetics equation: 

𝛾𝛾𝑎𝑎𝑎𝑎

𝑒𝑒𝜕𝜕
𝐽𝐽𝑎𝑎𝑎𝑎 = 𝑗𝑗𝑎𝑎𝑎𝑎

�

𝑋𝑋𝐻𝐻2
𝑟𝑟𝑒𝑒𝑒𝑒�
𝑋𝑋𝐻𝐻2

�exp �

𝛼𝛼𝑎𝑎𝑎𝑎𝐹𝐹𝑖𝑖𝑎𝑎𝑎𝑎
𝑅𝑅𝑇𝑇

� − exp �

−(1 − 𝛼𝛼𝑎𝑎𝑎𝑎)𝐹𝐹𝑖𝑖𝑎𝑎𝑎𝑎
𝑅𝑅𝑇𝑇

��

𝛾𝛾𝑐𝑐𝑎𝑎𝑐𝑐

𝑒𝑒𝜕𝜕
𝐽𝐽𝑐𝑐𝑎𝑎𝑐𝑐 = 𝑗𝑗𝑐𝑐𝑎𝑎𝑐𝑐
where: 

�

𝑋𝑋𝑂𝑂2
𝑟𝑟𝑒𝑒𝑒𝑒�
𝑋𝑋𝑂𝑂2

�exp�

−𝛼𝛼𝑐𝑐𝑎𝑎𝑐𝑐𝐹𝐹𝑖𝑖𝑐𝑐𝑎𝑎𝑐𝑐
𝑅𝑅𝑇𝑇

� − exp �

(1 − 𝛼𝛼𝑐𝑐𝑎𝑎𝑐𝑐)𝐹𝐹𝑖𝑖𝑐𝑐𝑎𝑎𝑐𝑐
𝑅𝑅𝑇𝑇

��

—  F and R are Faraday and universal gas constants; 

—  T is the local temperature; 

(24) 

(25) 

(26) 

— 

— 

 and 

 are exchange current densities [A m-3] for anode and cathode catalysts, respectively, given in 
Appendix B, Table 10. Note that Electro-Chemical Surface Area (ECSA, defined as an active catalyst area 
𝑒𝑒𝜕𝜕
𝑗𝑗𝑎𝑎𝑎𝑎
surface per volume) is already accounted for by  

𝑒𝑒𝜕𝜕
𝑗𝑗𝑐𝑐𝑎𝑎𝑐𝑐

(27) 

𝑒𝑒𝜕𝜕
𝑗𝑗𝑎𝑎𝑎𝑎,𝑐𝑐𝑎𝑎𝑐𝑐
where 

𝑟𝑟𝑒𝑒𝑒𝑒
= 𝑗𝑗𝑎𝑎𝑎𝑎,𝑐𝑐𝑎𝑎𝑐𝑐

𝜁𝜁𝑎𝑎𝑎𝑎,𝑐𝑐𝑎𝑎𝑐𝑐

 is reference current density per active catalyst area [A m-2], ζ is ECSA [m-1];  

𝑟𝑟𝑒𝑒𝑒𝑒
𝑗𝑗𝑎𝑎𝑎𝑎,𝑐𝑐𝑎𝑎𝑐𝑐

  are hydrogen and oxygen local concentrations [mole m-3]. The reference values of 
 and 
𝑟𝑟𝑒𝑒𝑒𝑒
 are given in Appendix B, Table 10 and correspond to reactant concentrations at NTP3 conditions;  
𝑋𝑋𝐻𝐻2

𝑋𝑋𝐻𝐻2
𝑟𝑟𝑒𝑒𝑒𝑒
—  Anode and cathode concentration exponents are denoted by γan and γcat, respectively. The values of these 
𝑋𝑋𝑂𝑂2
constants are considered as fitting parameters, however γan=0.5 and γcat=1.0 are mentioned in many 
papers, i.e. [2]; 

 and 

𝑋𝑋𝑂𝑂2

—  αan and αcat are anode and cathode exchange coefficients. The exchange coefficient α is related to reaction 

symmetry factor β, as follows [3]: 

where n is the number of electrons taking part in the reaction: 
𝛼𝛼 = 𝑎𝑎𝑛𝑛

(28) 

(29) 

+

1
1
2 𝐻𝐻2 → 𝐻𝐻
2 𝑂𝑂2 → 𝑂𝑂
According to [3], PEM fuel cell reactions are symmetric (β=0.5), hence the values of anode and cathode 
+ 𝑒𝑒̅ ,                   
exchange coefficients are taken as αan=0.5 and αcat=1.0; and 

− 2𝑒𝑒̅

2−

—  ηan and ηcat are anode and cathode overpotentials, respectively, defined as follows: 

(30) 

0
𝑖𝑖𝑎𝑎𝑎𝑎 = 𝜙𝜙𝑚𝑚𝑒𝑒𝑚𝑚 − 𝜙𝜙𝑠𝑠𝑠𝑠𝑙𝑙 − 𝑈𝑈𝑎𝑎𝑎𝑎

0
,                 𝑖𝑖𝑐𝑐𝑎𝑎𝑐𝑐 = 𝜙𝜙𝑚𝑚𝑒𝑒𝑚𝑚 − 𝜙𝜙𝑠𝑠𝑠𝑠𝑙𝑙 − 𝑈𝑈𝑐𝑐𝑎𝑎𝑐𝑐

3 Normal Te mperature and Pressure: 20oC (293.15 K) and 1 atm (101.325 kPa). 

22 

  
      
 
 
 
 
        
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
                                     
where ϕmem and ϕsol are potentials of the electrolyte and carbon support in the appropriate catalyst layer. 
The half-potentials 

 are defined as follows: 

 and 

0
𝑈𝑈𝑎𝑎𝑎𝑎

0
𝑈𝑈𝑐𝑐𝑎𝑎𝑐𝑐

(31) 

𝑎𝑎𝑎𝑎

𝑅𝑅𝑇𝑇
2𝐹𝐹 ln �
𝑅𝑅𝑇𝑇
4𝐹𝐹 ln �

𝑝𝑝𝐻𝐻2
𝑝𝑝0 �
𝑝𝑝𝑂𝑂2
𝑝𝑝0 � −

= −

𝑈𝑈0
𝑐𝑐𝑎𝑎
𝑈𝑈0
where p denotes pressure, while subscripts H2, O2 and H2O refer to the partial pressure of the gas species, 
p0 is a reference pressure equal to 0.1MPa, psat is the saturation pressure of water vapour. Note, the total 
pressure is the sum of the partial pressure of all species present less the saturation pressure of water. 

𝑅𝑅𝑇𝑇
2𝐹𝐹 ln�

𝑝𝑝𝐻𝐻2𝑂𝑂
𝑝𝑝𝑠𝑠𝑎𝑎𝑐𝑐 �

= 𝐸𝐸0 +

(32) 

E0 corresponds to the reversible potential of fuel cell reactions, defined as: 

Δ𝐺𝐺

𝐸𝐸𝑠𝑠 = −
where ΔG is change of free Gibbs energy, which at 80°C is equal to -2.262 x 105 [J mole-1]. Therefore, at 
an operating temperature of 80°C, the reversible potential E0 is equal to 1.172 [V]. 

2𝐹𝐹

(33) 

The transport of liquid water through GDL, MPL and CL follows the equation: 

𝜌𝜌𝑙𝑙𝐾𝐾𝐾𝐾𝑟𝑟

𝜕𝜕
𝜕𝜕𝑐𝑐 (𝜀𝜀𝜌𝜌𝑙𝑙𝑠𝑠) = ∇ ⋅ �
where, ε is porosity, ρl and μl are liquid water density and viscosity respectively, s is liquid water saturation, 
and 
respectively, Sld is liquid water uptake by the membrane, 

  are absolute permeability of the porous layer and relative permeability due to saturation effect, 
𝐾𝐾
 is an effective diffusion coefficient for liquid water. 

𝜇𝜇𝑙𝑙 ∇𝑝𝑝𝑙𝑙� + 𝛻𝛻 ⋅ �𝐷𝐷𝑚𝑚

∇s�− 𝑆𝑆𝑙𝑙𝑔𝑔 + 𝑆𝑆𝑔𝑔𝑙𝑙

𝑒𝑒𝑒𝑒𝑒𝑒

The left hand side term in equation 34 corresponds to the dynamic inertial term while the right hand side refers 
to the steady state transport of liquid water. The consecutive terms in the right hand side refer to: 

𝑠𝑠
𝐷𝐷𝑚𝑚

𝐾𝐾𝑟𝑟

 which is a liquid water flux due to pressure difference and 

𝜌𝜌𝑙𝑙𝐾𝐾𝐾𝐾𝑟𝑟

 which is a liquid water flux due to saturation gradient. 

— 

— 

𝜇𝜇𝑙𝑙 ∇𝑝𝑝𝑙𝑙
𝑒𝑒𝑒𝑒𝑒𝑒

𝐷𝐷𝑚𝑚

∇s

The liquid water pressure is a sum of gas phase pressure p and capillary pressure pc, namely: 

𝑝𝑝𝑙𝑙 = 𝑝𝑝 + 𝑝𝑝𝑐𝑐
The relative water permeability of porous layer due to liquid water saturation is approximated as:  

𝑏𝑏

𝐾𝐾𝑟𝑟 = 𝑠𝑠
where exponent b is a user-defined constant (assumed b=3.0)[1]. The effective diffusion coefficient in the 
membrane 

 is given in equation (12). 

𝑒𝑒𝑒𝑒𝑒𝑒

𝐷𝐷𝑚𝑚
Membrane 

Transport of H+ ions in the ionomer (membrane and CL) 

Due to the fact that H+ ion transport takes place through the membrane, but between the anode and cathode 
catalyst layers, the ionomer volume fraction in both catalyst layers is taken into consideration. 

The conservation equation for electrical charge in the membrane is expressed as follows: 

eff

∇ ⋅ �𝜎𝜎𝑚𝑚𝑒𝑒𝑚𝑚

∇𝜙𝜙𝑚𝑚𝑒𝑒𝑚𝑚� = 0

23 

(34) 

(35) 

(36) 

(37) 

  
      
 
   
 
 
 
 
 
 
 
   
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  
 
  
   
 
  
(38) 

(39) 

(40) 

  is the effective ionic conductivity of the membrane and 

where, 
  is the membrane potential. The 
effective ionic conductivity is dependent on membrane liquid water content, λ, and temperature, T, and is given 
as follows [4]: 

eff
𝜎𝜎𝑚𝑚𝑒𝑒𝑚𝑚

𝜙𝜙𝑚𝑚𝑒𝑒𝑚𝑚

eff
𝜎𝜎𝑚𝑚𝑒𝑒𝑚𝑚
Ei stands for activation energy for temperature correction. Γ factor is defined as follows: 

(λ, T) = Γ(0.514λ − 0.326) exp �𝐸𝐸𝑖𝑖 �

1
303 −

1
𝑇𝑇��

1
𝜖𝜖𝑎𝑎
𝜏𝜏𝑎𝑎
𝜖𝜖𝑐𝑐

𝛤𝛤 = �

in membrane
in anode catalist

𝜏𝜏𝑐𝑐

in cathode catalist

where, ϵ and τ refer to volume fraction and tortuosity of ionomer in anode (subscript a) and cathode (subscript 
c) catalyst layers, respectively. 

Water transport though the membrane 

In PEM fuel cell water exist in three phases, namely: 

—  water vapour as humidification of anode and cathode gases; 

—  liquid water, considered as fine homogenous mist; and 

—  dissolved water in the membrane ionomer.  

Transport equation of dissolved water in the ionomer 

Dissolved water phase exists in the ionomer in both anode and cathode catalyst layers and in the membrane. 
The transport of dissolved water, both ways due to osmotic drag and back diffusion due to concentration 
gradient, follows the formula [4]: 

∂
∂t �𝜖𝜖𝑖𝑖𝑀𝑀𝐻𝐻2𝑂𝑂
where: 

𝜌𝜌𝑙𝑙
𝐸𝐸𝐸𝐸  𝜆𝜆� + ∇ ⋅ �𝐢𝐢

𝑎𝑎𝑑𝑑
𝐹𝐹 𝑀𝑀𝐻𝐻2𝑂𝑂� = ∇ ⋅ �𝑀𝑀𝐻𝐻2𝑂𝑂𝐷𝐷𝜆𝜆∇𝜆𝜆� + 𝑆𝑆𝜆𝜆 + 𝑆𝑆𝑔𝑔𝑔𝑔 + 𝑆𝑆𝑙𝑙𝑔𝑔

—  λ is a dimensionless water content defined as the number of water molecules per number of charge sides 

(-SO3

-H+)[2,5,6];  

—  ϵi is a dimensionless volume fraction of the ionomer (equal 1 for the membrane) [4];  

—  ρi is density of the dry ionomer [kg m-3]; 

—  EW is membrane equivalent weight (weight per one acid site -SO3

-) [kg mol-1]; 

— 

 is molecular weight of water [kg mol-1];  

—  nd is dimensionless osmotic drag coefficient;  

𝑀𝑀𝐻𝐻2𝑂𝑂

—  Dλ is diffusion coefficient of water content [mol s-1];  

—  Sλ is dissolved water production rate in cathode catalyst due to ORR reaction [kg s-1], equals 0 in anode;  

—  Sgd is mass change rate between gas and dissolved phases [kg s-1];  

—  Sld is mass change rate between liquid and dissolved phases [kg s-1]; and  

—  i is ionic current [A] calculated as 

. 

Source terms on the right hand side of equation 40 are defined as follows: 

𝑒𝑒𝑒𝑒𝑒𝑒
𝐢𝐢 = −𝜎𝜎𝑚𝑚𝑒𝑒𝑚𝑚

∇𝜙𝜙𝑚𝑚𝑒𝑒𝑚𝑚

24 

  
      
 
  
  
 
 
  
𝐻𝐻2𝑂𝑂
𝑀𝑀𝑤𝑤
2𝐹𝐹 𝑗𝑗𝑐𝑐𝑎𝑎𝑐𝑐

𝑆𝑆𝜆𝜆 =

𝜃𝜃

𝑆𝑆𝑔𝑔𝑔𝑔 = �1 − 𝑠𝑠

𝐻𝐻2𝑂𝑂 𝜌𝜌𝑖𝑖

�𝛾𝛾𝑔𝑔𝑔𝑔𝑀𝑀𝑤𝑤

𝐸𝐸𝐸𝐸 (𝜆𝜆𝑒𝑒𝑒𝑒 − 𝜆𝜆) 

𝜃𝜃

𝐻𝐻2𝑂𝑂 𝜌𝜌𝑖𝑖

𝛾𝛾𝑙𝑙𝑔𝑔𝑀𝑀𝑤𝑤

𝐸𝐸𝐸𝐸(𝜆𝜆𝑒𝑒𝑒𝑒 − 𝜆𝜆) 

𝑆𝑆𝑙𝑙𝑔𝑔 = 𝑠𝑠
where jcat is current density in the cathode catalyst layer, s is liquid water saturation, θ is user defined constant, 
γgd and γld are gas to dissolved and liquid to dissolved mass exchange rate constants for water. The equilibrium 
water content is given as follows [4]: 

𝜆𝜆𝑒𝑒𝑒𝑒 = 0.3 + 6𝑎𝑎(1 − tanh(𝑎𝑎 − 0.5)) + 0.69(𝜆𝜆𝑎𝑎=1 − 3.52)𝑎𝑎
where a stands for water activity. 

0.5

�1 + tanh �

𝑎𝑎−0.89
0.23 �� + 𝑠𝑠(𝜆𝜆𝑠𝑠=1 − 𝜆𝜆𝑎𝑎=1)

𝜕𝜕𝐻𝐻2𝑂𝑂 𝑝𝑝
𝑝𝑝𝑠𝑠𝑎𝑎𝑐𝑐   

𝑎𝑎 =
where 
Parameters 
𝑥𝑥𝐻𝐻2𝑂𝑂

 is mole fraction of water vapour, p is absolute gas pressure and psat is saturation vapour pressure. 

 and 

 in equation 44 are user-defined constants equal to 16.8 and 9.2, respectively [5]. 

Diffusivity of the water content in the membrane is given as follows [4]: 

𝜆𝜆𝑠𝑠=1

𝜆𝜆𝑎𝑎=1

𝜌𝜌𝑖𝑖
𝐸𝐸𝐸𝐸 4.1 ⋅ 10

−10

𝜆𝜆
25�
�

0.15

𝐷𝐷𝜆𝜆 = 
The value of the osmotic-drag coefficient, nd, is given as [7]: 

�1 + tanh �

𝜆𝜆−2.5
1.4  ��

2.5
22 𝜆𝜆

𝑎𝑎𝑔𝑔 =
However, other researchers provide a bi-linear osmotic drag coefficient relation to water content, λ, or simply 
a constant value equal to nd=1.0 [8]. 

Transport of liquid water due to concentration, pressure and temperature difference 

The transport of liquid water through the membrane follows the equation: 

𝜌𝜌𝑙𝑙𝐾𝐾𝐾𝐾𝑟𝑟

𝜌𝜌𝑙𝑙𝑠𝑠� = ∇ ⋅ �

𝜕𝜕
𝑙𝑙𝑖𝑖𝑒𝑒𝑢𝑢𝑖𝑖𝑔𝑔
𝜕𝜕𝑐𝑐 �𝜖𝜖𝑚𝑚𝑒𝑒𝑚𝑚
where, 
water density and viscosity, respectively, 
membrane, respectively, and Sld is the uptake of liquid water by the membrane defined in equation (43). 
𝐾𝐾𝑟𝑟

 is the volume fraction of membrane that could be occupied by liquid water, ρl and μl are liquid 
  are absolute and relative water permeability of the 

𝜇𝜇𝑙𝑙 ∇𝑝𝑝𝑙𝑙� + 𝑆𝑆𝑙𝑙𝑔𝑔

The liquid water pressure is a sum of gas phase pressure, p, and capillary pressure, pc, namely: 

𝑙𝑙𝑖𝑖𝑒𝑒𝑢𝑢𝑖𝑖𝑔𝑔
𝜖𝜖𝑚𝑚𝑒𝑒𝑚𝑚

  and 

𝐾𝐾

𝑝𝑝𝑙𝑙 = 𝑝𝑝 + 𝑝𝑝𝑐𝑐
The relative water permeability of the membrane is a function of membrane water content:  

𝑀𝑀𝐻𝐻2𝑂𝑂
𝜌𝜌𝑖𝑖
𝑀𝑀𝐻𝐻2𝑂𝑂
𝜌𝜌𝑖𝑖

𝜆𝜆𝑠𝑠=1+

𝐸𝐸𝐸𝐸
𝜌𝜌𝑖𝑖
𝐸𝐸𝐸𝐸
𝜌𝜌𝑖𝑖

𝜆𝜆+

𝐾𝐾𝑟𝑟 = �

2

𝜆𝜆
𝜆𝜆𝑠𝑠=1�

25 

(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

  
      
 
  
  
  
  
  
  
  
 
  
  
 
  
 
Heat transfer 

The heat transfer in the membrane is modelled as follows: 

(51) 

2

ℎ

𝑇𝑇 = 𝑆𝑆𝑚𝑚𝑒𝑒𝑚𝑚

ℎ� + 𝑘𝑘𝑚𝑚𝑒𝑒𝑚𝑚∇
 and k

𝜌𝜌𝑖𝑖
𝜕𝜕
�
where ρ
𝜕𝜕𝑐𝑐
𝐸𝐸𝐸𝐸
respectively, h is enthalpy of the membrane and 
is  not dependent on spatial geometry, the k
Laplacian operator (
depends on the local protonic current, i, and protoic conductivity as follows: 

 are physical parameters of the membrane material, namely density and thermal conductivity, 
T is the temperature gradient. Since the thermal conductivity 
 term is excluded from the divergence operator hence the 
∇
, represents the Joule ohmic heating and 
BPP

2) appears. The volumetric source term, 

mem

i

𝐵𝐵𝐵𝐵𝐵𝐵
𝑆𝑆ℎ

∇

2
𝑖𝑖𝑝𝑝
𝑒𝑒𝑒𝑒𝑒𝑒  
𝜎𝜎𝑚𝑚𝑒𝑒𝑚𝑚

=

ℎ
𝑆𝑆𝑚𝑚𝑒𝑒𝑚𝑚
The local protonic current, ip, and its magnitude, ip, is proportional to the protonic potential gradient 
is calculated according to Ohm's law  

 and 

(52) 

(53) 

∇𝛷𝛷𝑚𝑚𝑒𝑒𝑚𝑚

𝑒𝑒𝑒𝑒𝑒𝑒
𝐢𝐢𝐩𝐩 = 𝜎𝜎𝑚𝑚𝑒𝑒𝑚𝑚

∇𝜙𝜙𝑚𝑚𝑒𝑒𝑚𝑚 

References 

[1] 

E.L. Cussler, Diffusion. Mass Transfer in Fluid Systems, Cambridge University Press, Cambridge, 2009. 

[2] 
A.Z. Weber, R.L. Borup, R.M. Darling, P.K. Das, T.J. Dursch, W. Gu, D. Harvey, A. Kusoglu, S. Litster, M.M. 
Mench, R. Mukundan, J.P. Owejan, J.G. Pharoah, M. Secanell, I. V Zenyuk, J. Electrochem. Soc.  161 (2014) F1254. 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 

R.F. Mann, J.C. Amphlett, B.A. Peppley, C.P. Thurgood, J. Power Sources 161 (2006) 775. 

H. Wu, X. Li, P. Berg, Electrochim. Acta 54 (2009) 6913. 

T.E. Springer, M.S. Wilson, S. Gottesfeld, J. Electrochem. Soc.  140 (1993) 3513. 

A.Z. Weber, J. Newman, Adv. Fuel Cells 1 (2007) 47. 

T.E. Springer, T.A. Zawodzinski, S. Gottesfeld, J. Electrochem. Soc.  138 (1991) 2334. 

T. Berning, M. Odgaard, S.K. Kær, J. Power Sources 196 (2011) 6305. 

26 

  
      
 
 
 
 
 
 
 
 
Appendix  B  Reference  test case 1 for  PEMFC single cell 

Computational geometry 

To simplify the numerical solution, the “single serpentine cell” is approximated by a “single straight channel” of 
equivalent channel length. The cross-sectional area of the anode and cathode flow channels is 1.5 mm2 – the 
channel width is 1 mm and the height of the channel is 1.5 mm. The width of the bi-polar plate ribs is also 1 
mm. The length of the gas channels is 1250 mm, which corresponds to the total length of the gas channels in 
a serpentine configuration representing a 25 cm2 active area. 

Due to vertical symmetry planes of the gas channels and bi-polar plate, only half of the gas channel and half 
of the rib are modelled. The geometry of the numerical model is presented in Figure 6. The computational model 
covers half of the fuel cell active area, namely 12.5 cm2. 

The gas diffusion layers on the cathode and anode sides are the same and of 200 µm thickness according to 
SGL®GDL24BC data sheets. The catalyst layers and the microporous layers are 20 µm thick on both sides. The 
membrane is assumed to be 20 µm in thickness. 

Figure 6 Computational domain. 

y 

z 

x 

Source : JRC, 2017 

27 

  
      
 
 
 
                           
 
 
 
 
Model assumptions 

The following working assumptions were made: 

—  The flow of reactant gases (air and hydrogen) within the fuel cell is counter flow; 

—  The open circuit voltage VOC is estimated on the basis of the Nernst equation. If this is not possible, VOC 

should be assumed as a constant input value; and 

—  The electrical boundary conditions at the anode and cathode bi-polar plate terminals are expressed as 

current density related to the catalyst area. 

External boundary conditions 

The fuel cell model boundary conditions are assumed as follows (Figure 7). 

—  The top and bottom surfaces of the bi-polar plates are assumed to be the fuel cell electric terminals. The 
potential on the anode bi-polar plate terminal is set to 0 V while the current or potential is applied to the 
cathode  bi-polar  plate  terminal,  depending  on  the  set  operation  mode,  namely  galvanostatic  or 
potentiostatic. 

—  A nominal operating temperature is applied to both terminals. 

—  Symmetry boundary conditions are applied on the fuel cell model symmetry plane.  

—  The inlet flow rates of oxidant (air) and fuel (hydrogen) are calculated according to the stoichiometry 

regime.  

—  The inlet pressure is equal to the operating pressure. However, if it is not possible to apply inlet pressure, 

the outlet pressure could be specified instead. 

—  There are no flux conditions (species, temperature and electric charge) on the solid surfaces at other 

external walls of the fuel cell model. 

External boundary conditions are presented in Figure 7 and Table 4 while operating conditions are given in Table 
4. Subsequent tables show the physical properties of fuel cell component materials, air and hydrogen. 

28 

  
      
 
 
 
 
Figure 7 External boundary conditions. 

Source : JRC, 2017 

29 

  
      
 
 
                        
 
 
Table 4 External boundary conditions. 

Location 

Boundary condition 

Inlet of the anode gas channel 

Mole fractions of fuel species: 

(hydrogen) 

𝑄𝑄𝐻𝐻2 = 𝑄𝑄𝐻𝐻2�����

Mass flow rate of fuel (hydrogen + water vapour): 

Temperature of inlet gas: 
��������,
𝜀𝜀[𝐻𝐻2] = 𝜀𝜀[𝐻𝐻2]

����������������
𝜀𝜀[𝐻𝐻2𝑂𝑂] = 𝜀𝜀𝑎𝑎𝑎𝑎𝑠𝑠𝑔𝑔𝑒𝑒[𝐻𝐻2𝑂𝑂]

Outlet of the anode gas channel 

Outlet pressure: 

𝑇𝑇 = 𝑇𝑇�

Inlet of the cathode gas channel 

(air) 

𝑝𝑝 = 𝑝𝑝̅ 𝑎𝑎𝑎𝑎𝑠𝑠𝑔𝑔𝑒𝑒
Mass flow rate of humid air (air+water vapour): 
Mole fraction of air species: 

�����
𝑄𝑄𝑎𝑎𝑖𝑖𝑟𝑟 = 𝑄𝑄𝑎𝑎𝑎𝑎𝑟𝑟

Temperature of inlet gas: 
�������,
𝜀𝜀[𝑂𝑂2] = 𝜀𝜀[𝑂𝑂2]

������������������
𝜀𝜀[𝐻𝐻2𝑂𝑂] = 𝜀𝜀𝑐𝑐𝑎𝑎𝑐𝑐ℎ𝑠𝑠𝑔𝑔𝑒𝑒[𝐻𝐻2𝑂𝑂]

Outlet of the cathode gas channel  Outlet pressure: 

𝑇𝑇 = 𝑇𝑇�

𝑝𝑝 = 𝑝𝑝𝑐𝑐𝑎𝑎𝑐𝑐ℎ𝑠𝑠𝑔𝑔𝑒𝑒
����������
Zero current flux normal to the symmetry plane: 

Zero velocity normal to the symmetry plane: 

∂j
∂x = 0

Zero species concentration gradients normal to the symmetry 
plane: 

Vx = 0

Symmetry planes, orange surfaces 
in Figure 7 

∂C[H2O]
Zero temperature gradient, normal to the symmetry plane: 
∂x

∂C[O2]
∂x

∂C[H2]
∂x

= 0,

= 0,

= 0

∂T
∂x =

Zero current flux normal to the boundary: 
0

End of model, green surfaces in 
Figure 7 

No slip velocity constraint: 
Zero species concentration gradient normal to the boundary: 

∂j
∂z = 0
T

V = [0 0 0]

,    

,    

∂C[H2]
Zero temperature gradient normal to the boundary: 

∂C[H2O]

∂C[O2]

∂z = 0

∂z = 0

∂z = 0

Anode bi-polar plate terminal 

Cathode bi-polar plate terminal 

Source : JRC, 2017 

Potential equal zero: 

Temperature: 

𝛷𝛷 = 0

𝑇𝑇 = 𝑇𝑇�
Prescribed electric current density: 

Temperature: 

𝑇𝑇 = 𝑇𝑇�

𝐣𝐣 = 𝐣𝐣̅

∂T
∂z = 0

30 

  
      
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Table 5 Operating Conditions 

Operating conditions 

Symbol 

Unit 

Setting 

Tcell 

Tfuel,in 

K 

K 

353 

353 

RH fuel,in 

% 

50 

Pfuel,in 

kPa 

250 

Parameter 

Nominal cell 
operating 
temperature 

Fuel gas inlet 
temperature 

Fuel gas inlet 
humidity RH 

Fuel gas inlet 
pressure 
(absolute) 

Fuel inlet 
stoichiometry 

λH2 

– 

Oxidant gas inlet 
temperature 

Tair,in 

Oxidant gas inlet 
humidity 

RH Ox,in 

% 

30 

Oxidant inlet 
pressure 
(absolute) 

Oxidant inlet 
stoichiometry 

POx,in 

kPa 

230 

λOx 

– 

1.3 

353 

1.5 

0.2 

Minimum current density 
for stoichiometry operation 

Jmin 

A cm-2

Source : JRC, 2017 

E
D
O
N
A

E
D
O
H
T
A
C

The simulation of the fuel cell operating at a current density of 0.2 Acm-2 and higher, follows the stoichiometry 
regime, i.e. the amount of supplied oxidant (air) and fuel (hydrogen) is directly proportional to the applied current 
density. For low current densities of up to 0.2 Acm-2, the constant flow rates of the supplied gases are assumed 
to be equal to the flow rate corresponding to the current density 0.2 Acm-2.  

The inlet mas flow rate of the supplied gases is calculated according to the formula: 

𝑀𝑀𝜕𝜕𝐽𝐽𝐽𝐽𝜆𝜆𝜕𝜕
𝑧𝑧𝜕𝜕𝐹𝐹𝜙𝜙𝜕𝜕

𝑄𝑄𝜕𝜕 =
where subscript x corresponds to hydrogen or air. Mx is the  molar mass, J is the applied current density, A is 
the active area of the catalyst layer, zx is the coefficient corresponding to the number of atoms taking part in 
the reactions, F is Faraday constant, ϕx is the volumetric content of hydrogen or oxygen in the supplied gases. 
The values of the quantities are given in Table 6. 

31 

  
      
 
 
 
 
 
 
 
 
 
 
Table 6 P roperties of reactant gases for calculation of the inlet mass flow rates 

air 

hydrogen 

Mx  
[kg mol-1] 

28.186e-3 

3.509e-3 

z 

4 

2 

ϕx  [%] 

19.66 

90.57 

Source : JRC, 2017 

The composition of inlet gases (mole fractions of gas species) is defined for the nominal operating fuel cell 
temperature. The 30% RH air, the oxidant gas, is composed as follows [mole fractions]: 
[H2O]=0.0615 and 
[mole fractions].  

 [N2]=0.7419, while the fuel (RH 50%) is composed of 

[O2]=0.1966, 
 [H2O]=0.0943 
𝜀𝜀̅

 [H2]=0.9057 and 
𝜀𝜀̅
𝜀𝜀̅

𝜀𝜀̅

𝜀𝜀̅

Physical properties 

The physical constants and parameters of the fuel cell used in the current computations are listed in Tables 7 
to 13. 

Table 7 P hysical properties of inlet gases at operating conditions 

Parameter 

Unit 

H2 

O2 

N2 

H2O 
(vapour) 

Specific heat 

J kg-1 K-1 

14.43e3 

0.928e3 

1.041e3 

1.88e3 

Thermal conductivity 

W m-1 K-1 

0.253 

0.031 

0.030 

0.025 

Viscosity 

Kg m-1 s-1 

8.41e-6 

2.34e-5 

2.01e-5 

1.30e-5 

Molecular weight 

kg mol-1 

2.02e-3 

31.99e-3 

28.01e-3 

18.01e-3 

Source : JRC, 2017 

Table 8 P roperties of bi-polar plates 

Parameter 

Unit 

Value 

Specific heat 

J kg-1K-1 

708 

Thermal conductivity 

W m-1K-1 

10.0 

Electrical conductivity 

Ohm-1m-1 

1.0e5 

Density 

Source : JRC, 2017 

kg m-3 

2100 

32 

  
      
 
 
 
 
 
 
Table 9 P roperties of anode and cathode Gas Diffusion Layers 

Parameter 

Unit 

Value 

Under          the 
channel 

Under 

the rib 

Specific heat 

J kg-1K-1 

708 

Thermal conductivity* 

W m-1K-1 

0.83 

Electrical conductivity* 

Ohm-1m-1 

184 

422 

200 

70 

Density* 

Thickness 

Porosity 

Gas permeability 

Contact angle 

kg m-3 

µm 

% 

m2 

deg 

* e ffe ctive value, porosity is taken into account 

Source : JRC, 2017 

708 

0.95 

195 

850 

200 

50 

6.74e-14 

2.50e-14 

112 

112 

Table 10 P roperties of anode and cathode Microporous Layers 

Parameter 

Unit 

Value 

Specific heat 

J kg-1K-1 

708 

Thermal conductivity* 

W m-1K-1 

1.5 

Electrical conductivity* 

Ohm-1m-1 

2.0e3 

Density* 

Thickness 

Porosity 

Gas permeability 

Contact angle 

kg m-3 

1260 

µm 

% 

m2 

deg 

20 

40 

3.33e-15 

112 

* e ffe ctive value due to porosity is taken into account 

Source : JRC, 2017 

33 

  
      
 
 
 
 
 
 
 
 
Table 11 P roperties of anode and cathode Catalyst Layers 

Parameter 

Unit 

Value 

Specific heat 

J kg-1K-1 

680 

Thermal conductivity* 

W m-1K-1 

1.7 

Electrical conductivity* 

Ohm-1m-1 

5.0e3 

kg m-3 

1350 

µm 

% 

% 

m2 

deg 

Density* 

Thickness 

Porosity 

Volume fraction of ionomer 

Ionomer tortuosity 

Gas permeability 

Contact angle 

Catalyst loading 

ECSA** 

Anode 

Cathode 

Anode 

Cathode 

mg cm-2 

m2 m-3 

* e ffe ctive value, porosity is taken into account 
** EC SA – Ele ctro-Chemical Surface Are a 

Source : JRC, 2017 

20 

38 

25 

1.5 

2.0e-15 

112 

0.2 

0.4 

1.7e7 

3.0e7 

34 

  
      
 
 
 
 
 
Table 12 P roperties of membrane. 

Parameter 

Unit 

Value 

Specific heat 

J kg-1K-1 

2000 

Thermal conductivity* 

W m-1K-1 

0.2 

Equivalent weight 

Thickness 

- 

µm 

1100 

20 

Protonic conductivity 

S m-1 

Osmotic drag coefficient 

Water 
membrane 

diffusivity 

in 

m2 s-1 

Source : JRC, 2017 

𝜅𝜅 = (0.5139𝜆𝜆 − 0.326)𝑒𝑒𝑥𝑥𝑝𝑝 �1268�

1
303

−

1
𝑇𝑇

��

𝑎𝑎𝑔𝑔𝑟𝑟𝑎𝑎𝑔𝑔 =

2.5
22

𝜆𝜆

𝐷𝐷𝜆𝜆 = 

𝜌𝜌𝑖𝑖
𝐸𝐸𝐸𝐸

−10

4.1 ⋅ 10

𝜆𝜆
�
�
25

0.15

�1 + tanh �

𝜆𝜆 − 2.5
1.4

 ��

Table 13 Electrochemical parameters corresponding to Butler-Volmer reaction kinetics 

Parameter 

Unit 

Value 

Anode (H2)                        Cathode (O2) 

1.0e4 

8.65 

1 

1 

Eq. (3) 

Exchange current 
density* 

Reference specie 
concentration 

Concentration 
exponent 

Exchange coefficient 

Open circuit voltage 

Half-cell potentials 

V 

V 

Am-3 

1.0e9 

mole m-3 

41.1 

0.5 

0.5 

Eq. (1) 

Eq. (2) 

*   The  e xchange current density takes into account the EC SA as we ll. 

Source : JRC, 2017 

35 

  
      
 
 
 
 
 
 
 
 
 
 
 
 
 
𝐸𝐸𝑠𝑠 = −
𝑎𝑎𝑎𝑎

Δ𝐺𝐺

2𝐹𝐹

(1) 

(2) 

𝑅𝑅𝑇𝑇
2𝐹𝐹 ln�
𝑅𝑅𝑇𝑇
4𝐹𝐹 ln �

𝑝𝑝𝐻𝐻2
𝑝𝑝0 �
𝑝𝑝𝑂𝑂2
𝑝𝑝0 � −

=  −

𝑈𝑈0
𝑐𝑐𝑎𝑎
𝑈𝑈0
) is natural logarithm, T is temperature, ΔG is change of free 
where R and F are gas and Faraday constants, ln(
Gibbs energy, p denotes pressure, while subscripts H2, O2 and H2O refer to partial pressures of gas species, p0 
is a reference pressure equal to 0.1MPa, psat is the saturation pressure of water vapour. 

𝑅𝑅𝑇𝑇
2𝐹𝐹 ln�

𝑝𝑝𝐻𝐻2𝑂𝑂
𝑝𝑝𝑠𝑠𝑎𝑎𝑐𝑐�

= 𝐸𝐸0 +

(3) 

⋅

36 

  
      
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Appendix  C  Reference  test case 2 for  catalyst layer examination 

INTRODUCTION 
The polarisation IV curve can be split into three regions of interest as shown in Figure 8, namely: 

—  Region I; Kinetics Region where electrochemical kinetics play a major role.  

—  Region II; Ohmic Region where the ohmic resistance of PEMFC components is the main source of voltage 

losses. In this region in particular, the membrane models will be evaluated. 

—  Region III; Transport Limitations Region where liquid water and properties of diffusion media are dominant.  

Figure 8 Regions of interest 

Source : JRC, 2017 

In region I, electrochemical kinetics and catalyst layer properties play the dominant role in determining the 
polarisation curve. Therefore, in benchmarking test case 2, the geometry, boundary and operating conditions 
are adjusted accordingly to place emphasis on the modelling of the electrochemical kinetics. 

37 

  
      
 
 
 
 
 
 
Geometry 

In an effort to simplify the numerical solution and to minimize the effects of gas flows on the fuel cell kinetics 
the “single straight channel” is considered. The cross-sectional area of the anode and cathode flow channels is 
1.5 mm2 with a channel width 1.0 mm and channel height of 1.5 mm. The width of the bi-polar plate ribs is 
also 1.0 mm. The length of the gas channels is equal to 80.0 mm, representing a 1.6 cm2 active area. The length 
of the fuel cell model and operating conditions (Table 14) were assumed to minimize gradients of pressure, 
concentration of reactants and temperature within the fuel cell active area.  

Due to vertical symmetry planes of the gas channels and bi-polar plate, only half of the gas channel and half 
of  the  rib  are  modelled. The geometry of  the numerical model is presented in  Figure 9.  Therefore, the 
computational model covers half of the fuel cell active area, namely 0.80 cm2. 

The gas diffusion layers on the cathode and anode sides are the same having 200 µm thickness – according to 
SGL®GDL24BC data sheets. The catalyst layers and the microporous layers are 20 µm thick on both sides. The 
membrane is assumed to be 20 µm in thickness. 

Figure 9 Computational domain for Region I – electrochemistry kinetics. 

Source : JRC, 2017 

38 

  
      
 
 
 
 
 
 
Common Operating Conditions 

The Common Operating conditions for electrochemical kinetics evaluation of the PEM fuel cell are listed in Table 
14. 

Table 14 Operating Conditions for benchmarking test case 2. 

Parameter 

Nominal cell 
operating 
temperature 

Fuel gas inlet 
temperature 

Fuel gas inlet 
humidity RH 

Fuel gas inlet 
pressure 
(absolute) 

Fuel inlet 
stoichiometry 

Oxidant gas inlet 
temperature 

Oxidant gas inlet 
humidity 

Oxidant inlet 
pressure 
(absolute) 

Operating conditions 

Unit 

Setting 

Justification 

OC 

OC 

% 

80 

80 

100 

bar 

2.5 

Typical temperature for PEM fuel 
cell operation 

The same as cell operating 
temperature 

To avoid drying out of the 
membrane 

Ambient pressure to avoid 
pressure corrections in 
electrochemical parameters 

Constant flow 
equivalent to 
stoichiometry 10.0 
@1 Acm-2 

To satisfy uniform distribution of 
species 

80 

100 

The same as operating 
temperature 

To avoid drying out of the 
membrane 

– 

OC 

% 

bar 

2.3 

To avoid pressure corrections 

E
D
O
N
A

E
D
O
H
T
A
C

Oxidant inlet 
stoichiometry 

– 

Constant flow 
equivalent to 
stoichiometry 10.0 
@1 Acm-2 

To satisfy uniform distribution of 
species 

Current density 

Acm-2 

OCV* – 1.0 

Only kinetics region of polarization 
curve 

* if O C V conditions cannot be simulated, a curre nt density of 5 mAcm-2 should be applied. 

Source : JRC, 2017 

39 

  
      
 
 
 
 
 
 
 
 
 
Gas composition and mass flow rates of supplied air and hydrogen for the 0.8cm2 active area model are 
presented in Figure 9 and are listed in Table 15. 

Table 15 Mass flow rates and gas composition of supplied gases 

Parameter 

Setting 

E
D
O
N
A

E
D
O
H
T
A
C

Total mass flow rate 

kg s-1 

1.539E-7 * 

H2 gas species 

Mole fraction 

Water vapour 

Mole fraction 

0.811 

0.189 

Total mass flow rate 

kg s-1 

3.314E-6 * 

O2 gas species 

Mole fraction 

H2O gas species 

Mole fraction 

0.167 

0.205 

* mass flow rates are equivalent to 0.8cm-2 active area. 

40 

  
      
 
 
 
 
 
 
 
 
 
 
GETTING  IN TOUCH  WITH  THE EU 

In person 

All  ove r the  European Union there are hundreds of Europe Dire ct information ce ntres. You can find the address of the ce ntre 
ne arest you at: http s://europa.eu/european-union/contact_en 

On the phone or by email 

Europ e Dire ct is a service that answers your questions about the  European Union. You can contact this service : 

- by fre e p hone: 00 800 6 7 8 9 10 11 (ce rtain operators may charge for the se calls), 

- at the  following standard number: +32 22999696, or 

- by e le ctronic mail via: http s://europa.eu/european-union/contact_en 

FINDING  INFORMATION  ABOUT  THE  EU 

Online 

Information about the  European Union in all the official languages of the EU is available  on the Europa website at: 
http s://europa.eu/european-union/index_en 

EU publications 
You can download or order free and p rice d EU publications from EU Bookshop at: http s://publications.europa.eu/en/publications. 
Multip le  cop ies of fre e publications may be obtained by contacting Europe Direct or your local information ce ntre (see 
http s://europa.eu/european-union/contact_en). 

  
      
 
 
 
 
K
J
-
N
A
-
3
0
7
0
3
-
E
N
-
N

doi:10 .2760/291664 

IS BN 9 78-92-76-37770-2