MEV


 

 

J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

Journal of Mechatronics, Electrical Power, 
and Vehicular Technology 

 
e-ISSN:2088-6985 
p-ISSN: 2087-3379 

 

 

 

www.mevjournal.com 
 

© 2016 RCEPM - LIPI All rights reserved. Open access under CC BY-NC-SA license. doi: 10.14203/j.mev.2016.v7.7-20. 
Accreditation Number: (LIPI) 633/AU/P2MI-LIPI/03/2015 and (Ministry of RTHE) 1/E/KPT/2015. 

A CFD MODEL FOR ANALYSIS OF PERFORMANCE, WATER AND 
THERMAL DISTRIBUTION, AND MECHANICAL RELATED FAILURE IN 

PEM FUEL CELLS 
 

Maher A.R. Sadiq Al-Baghdadi* 
Department of Mechanical Engineering, Faculty of Engineering, University of Kufa, 

Najaf, Kufa, Iraq 
 

Received 25 January 2016; received in revised form 03 May 2016; accepted 08 May 2016 
Published online 29 July 2016 

 

Abstract 
This paper presents a comprehensive three-dimensional, multi-phase, non-isothermal model of a Proton 

Exchange Membrane (PEM) fuel cell that incorporates significant physical processes and key parameters affecting 
the fuel cell performance. The model construction involves equations derivation, boundary conditions setting, and 
solution algorithm flow chart. Equations in gas flow channels, gas diffusion layers (GDLs), catalyst layers (CLs), 
and membrane as well as equations governing cell potential and hygro-thermal stresses are described. The algorithm 
flow chart starts from input of the desired cell current density, initialization, iteration of the equations solution, and 
finalizations by calculating the cell potential. In order to analyze performance, water and thermal distribution, and 
mechanical related failure in the cell, the equations are solved using a computational fluid dynamic (CFD) code. 
Performance analysis includes a performance curve which plots the cell potential (Volt) against nominal current 
density (A/cm2) as well as losses. Velocity vectors of gas and liquid water, liquid water saturation, and water content 
profile are calculated. Thermal distribution is then calculated together with hygro-thermal stresses and deformation. 
The CFD model was executed under boundary conditions of 20°C room temperature, 35% relative humidity, and 1 
MPA pressure on the lower surface. Parameters values of membrane electrode assembly (MEA) and other base 
conditions are selected. A cell with dimension of 1 mm x 1 mm x 50 mm is used as the object of analysis. The 
nominal current density of 1.4 A/cm2 is given as the input of the CFD calculation. The results show that the model 
represents well the performance curve obtained through experiment. Moreover, it can be concluded that the model 
can help in understanding complex process in the cell which is hard to be studied experimentally, and also provides 
computer aided tool for design and optimization of PEM fuel cells to realize higher power density and lower cost. 

 
Key words: CFD; PEM; fuel cell; multi-phase; hygro thermal stress. 

 
I. INTRODUCTION 

Water management is one of the critical 
operation issues in proton exchange membrane  
fuel cells (PMFCs). Spatially varying 
concentrations of water in both vapour and liquid 
form are expected throughout the cell because of 
varying rates of production and transport. Water 
emanates from two sources i.e. the product water 
from the oxygen-reduction reaction in the 
cathode catalyst layer and the humidification 
water carried by the inlet streams or injected into 
the fuel cell [1]. One of the main difficulties in 
managing water in a PEMFC is the conflicting 
requirements of the membrane and the catalyst 

gas diffusion layer (GDL). On the cathode side, 
excessive liquid water may block or flood the 
pores of the catalyst layer, the GDL or even the 
gas channel, thereby inhibiting or even 
completely blocking oxygen mass transfer. On 
the anode side, as water is dragged toward the 
cathode via electro-osmotic transport, 
dehumidification of the membrane may occur, 
resulting in deterioration of protonic conductivity. 
In the extreme case of complete drying, local 
burnout of the membrane can result. Devising 
better water management is a key issue. 

Thermal management is also required to 
remove the heat produced by the electrochemical 
reaction in order to prevent drying out of the 
membrane, which in turn can result not only in * Corresponding Author. Tel: +96-47719898955 

E-mail: mahirar.albaghdadi@uokufa.edu.iq 

http://dx.doi.org/10.14203/j.mev.2016.v7.1-6
tel:%2B9647719898955


M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

8 

reduced performance but also in eventual rupture 
of the membrane [2]. Thermal management is 
also essential for controlling the water 
evaporation or condensation rates. The difficult 
experimental environment of fuel cell systems 
has stimulated efforts to develop models that 
could simulate and predict multi-dimensional 
coupled transport of reactants, heat, and charged 
species using CFD methods. A comprehensive 
model should include equations and numerical 
relations needed to fully define fuel cell behavior 
over the range of interest. 

Early multidimensional models described gas 
transport in the flow channels, gas diffusion 
layers, and the membrane [3-5]. Iranzo et al. [6] 
developed a multi-phase, two-dimensional model 
to describe the liquid water saturation and flood 
effect, and have been studied transport limitations 
due to water build up in the cathode catalyst 
region. Djilali [7] developed a CFD multiphase 
model of a PEMFC. This model provides 
information on liquid water saturation and flood 
under various conditions, however it does not 
account for water dissolved in the ion-conducting 
polymer to calculate water content through the 
membrane. Hu et al. [8] and Fouquet [9] 
developed an isothermal, three-dimensional, two-
phase model for a PEMFC. Their model 
describes the transport of liquid water within the 
porous electrodes and water dissolved in the ion-
conducting polymer. The model is restricted to 
constant cell temperature. 

The need for improving lifetime of PEMFCs 
necessitates that the failure mechanisms be 
clearly understood, so that new designs can be 
introduced to improve long-term performance. 
Weber and Newman [10] developed one-
dimensional model to study the stresses 
development in the fuel cell. Bograchev et al. 
[11] studied the hygro and thermal stresses in the 
fuel cell caused by step-changes of temperature 
and relative humidity. However, their model is 
two-dimensional. In addition, constant 
temperature was assumed at each surface of the 
membrane. 

An operating fuel cell has various local 
conditions of temperature, humidity, and power 
generation across the active area of the fuel cell. 
Nevertheless, no models have yet been published 
to incorporate hygro-thermal stresses in three-
dimension. Therefore, in order to acquire a 
complete understanding of the damage 
mechanisms in the membranes, mechanical 
response under steady-state hygro-thermal 
stresses should be modelled and studied under 
real cell operation conditions. 

 

II. MODEL DESCRIPTION 
This paper presents a comprehensive three-

dimensional, multi-phase, non-isothermal model 
of a PEMFC. It accounts for both gas and liquid 
phase. It includes the transport of gaseous species, 
liquid water, protons, energy, and water dissolved 
in the ion-conducting polymer. The model 
features an algorithm to improve prediction of the 
local current density distribution. It takes into 
account convection and diffusion of different 
species, heat transfer, and electrochemical 
reactions. It also incorporates the effect of hygro-
thermal stresses.  

More specifically, this paper describes the 
development of the model, the determination of 
properties for use in the model, the validation of 
the model using experimental data, and the 
application of the model to explain observed 
experimental phenomena. Formulas governing 
the process inside a PEMFC have already been 
disclosed in the previous papers [20-24]. In order 
to make it self-content to some extent, some 
important equations are re-presented. 

 
A. Computational Domain 

In order to save computing resources and 
shorten simulation times, computational domain 
is limited to one straight flow channel. It consists 
of gas flow channels, and membrane electrode 
assembly (MEA) as shown in Figure 1. 

 
B. Model Equations 
1) Gas Flow Channels 

In the fuel cell channels, the gas-flow field is 
obtained by solving the steady-state Navier-
Stokes equations. The continuity equations for 
the gas phase and the liquid phase inside the 
channel is given by: 

 
Figure 1. Computational domain 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

9 

( ) 0g g gr ρ∇⋅ =u  (1) 
( ) 0=⋅∇ lllr uρ  (2) 

The following momentum equations are 
solved in the channels, and they share the same 
pressure field. 

( )
( )[ ]Tggggg

ggggg

Pr uu

uuu

∇⋅∇+







⋅∇+∇−

=∇−⊗⋅∇

µµ

µρ

3
2  (3) 

( )

( )[ ]Tlllll
lllll

Pr uu

uuu

∇⋅∇+







⋅∇+∇−

=∇−⊗⋅∇

µµ

µρ

3
2  (4) 

The mass balance is described by the 
divergence of the mass flux through diffusion and 
convection. Multiple species are considered in 
the gas phase only, and the species conservation 
equation in multi-component, multi-phase flow 
can be written in the following expression for 
species i: 

( )
0

                                                     

  
1 =



















∇
+⋅

+∑ 






 ∇
−+






 ∇

+∇−
⋅∇ =

T
T

Dyr

P
P

yx
M
M

yy
M
M

Dyr

T
igigg

N

j
jjjj

j
ijigg

uρ

ρ  (5) 

where the subscript i denotes oxygen at the 
cathode side and hydrogen at the anode side, and 
j is water vapour in both cases. Nitrogen is the 
third species at the cathode side. 

The Maxwell-Stefan diffusion coefficients of 
any two species are dependent on temperature 
and pressure. They can be calculated according to 
the empirical relation based on kinetic gas theory 
[12]: 

21

23131

375.1 1110











+






















∑+








∑

×
=

−

ji

k
kj

k
ki

ij MM
VVP

T
D

 (6) 

In this equation, pressure is in [atm] and the 
binary diffusion coefficient is in [cm2/s]. The 
values for ( )∑ kiV  are given by Fuller et al. [12]. 
The temperature field is obtained by solving the 
convective energy equation: 

( )( ) 0=∇−⋅∇ TkTCpr ggggg uρ  (7) 
The gas phase and the liquid phase are 

assumed to be in thermodynamic equilibrium; 
hence the temperature of the liquid water is the 
same as the gas phase temperature. 
2) Gas Diffusion Layers (GDL) 

The physics of multiple phases through a 
porous medium is further complicated here with 
phase change and the sources and sinks 

associated with the electrochemical reaction. The 
equations used to describe transport in the GDL 
are given below. Mass transfer in the form of 
evaporation ( )0>phasem  and condensation 
( )0<phasem  is assumed, so that the mass balance 
equations for both phases are: 

( )( ) phasegg msat =−⋅∇ uερ1  (8) 
( ) phasell msat =⋅∇ uερ.  (9) 

The momentum equation for the gas phase 
reduces to Darcy’s law, which is based on the 
relative permeability for the gas phase. The 
relative permeability accounts for the reduction 
in pore space available for one phase due to the 
existence of the second phase [7]. 

The momentum equation for the gas phase 
inside the GDL becomes: 

( ) PKpsat
g

g ∇−−= µ
1u  (10) 

Two liquid water transport mechanisms are 
considered; shear, which drags the liquid phase 
along with the gas phase in the direction of the 
pressure gradient, and capillary forces, which 
drive liquid water from high to low saturation 
regions [7]. Therefore, the momentum equation 
for the liquid phase inside the GDL becomes: 

sat
sat
PKP

P
KP c

l

l

l

l
l ∇∂

∂
+∇−=

µµ
u  (11) 

The functional variation of capillary pressure 
with saturation is prescribed following Leverett 
[7] who has shown that: 

𝑃𝑐 = 𝜏 �
𝜀
𝐾𝑃
�
1 2⁄

(1.417(1 − 𝑠𝑎𝑡) −
2.121−𝑠𝑎𝑡2+1.2631−𝑠𝑎𝑡3 (12) 

The liquid phase consists of pure water, while 
the gas phase has multi components. The 
transport of each species in the gas phase is 
governed by a general convection-diffusion 
equation in conjunction which the Stefan-
Maxwell equations to account for multi species 
diffusion: 

( ) ( )

( )
phase

T
igig

N

j
jjjj

j
ijig

m

T
T

Dysat

P
P

yx
M
M

yy
M
M

Dysat
=



















∇
+⋅−

+∑ 






 ∇
−+






 ∇

+∇−−
⋅∇ =

εερ

ερ

u1                                                      

1
1

 (13) 

In order to account for geometric constraints 
of the porous media, the diffusivities are 
corrected using the Bruggemann correction 
formula [3]: 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

10 

5.1ε×= ij
eff
ij DD  (14) 

The heat transfer in the gas diffusion layers is 
governed by the energy equation as follows: 

( )( )( )
( ) evapphasesolid

geffggg

HmTT

TkTCpsat

∆−−

=∇−−⋅∇

εεβ

εερ ,1 u  (15) 

where the term ( ( )TTsolid −εβ ), on the right 
hand side, accounts for the heat exchange to and 
from the solid matrix of the GDL. β is a modified 
heat transfer coefficient that accounts for the 
convective heat transfer in [W/m2] and the 
specific surface area [m2/m3] of the porous 
medium [3]. Hence, the unit of β is [W/m3]. 

The gas phase and the liquid phase are 
assumed to be in thermodynamic equilibrium, i.e., 
the liquid water and the gas phase are at the same 
temperature. The potential distribution in the 
GDLs is governed by: 

( ) 0=∇⋅∇ φλe  (16) 
In order to determine the magnitude of phase 

change inside the GDL, expressions are required 
to relate the level of over- and undersaturation as 
well as the amount of liquid water present to the 
amount of water undergoing phase change. In the 
present work, the procedure of Djilali [7] was 
used to determine the magnitude of phase change 
inside the GDL. 
3) Catalyst Layers (CLs) 

The catalyst layer is treated as a thin interface, 
where sink and source terms for the reactants are 
implemented. Due to the infinitesimal thickness, 
the source terms are implemented in the last grid 
cell of the porous medium. The sink terms for 
oxygen and hydrogen are given by: 

c
O

O iF
M

S
4

2
2

−=
; 

a
H

H iF
M

S
2

2
2

−=
 

(17) 

The production of water is modelled as a 
source terms, and hence can be written as: 

c
OH

OH iF
M

S
2

2
2

=  (18) 

The generation of heat in the cell is due to 
entropy changes as well as irreversibilities due to 
the activation overpotential [13]: 

( )
ccact

e
i

Fn
sT

q  , 







+

∆−
= η  (19) 

The local current density in the CLs is 
modelled by the Butler-Volmer equation [14, 15]; 
















−+























= cact

c
cact

a
ref
O

Oref
coc RT

F
RT

F
C

C
ii ,,, expexp

2

2 η
α

η
α

 (20) 
















−+























= aact

c
aact

a
ref
H

Href
aoa RT

F
RT

F
C

C
ii ,,

21

, expexp
2

2 η
α

η
α

 (21) 

4) Membrane 
The balance between the electro-osmotic drag 

of water from anode to cathode and back 
diffusion from cathode to anode yields the net 
water flux through the membrane: 

( )WWOHdW cDF
i

MnN ∇⋅∇−= ρ
2

 (22) 

The water diffusivity in the polymer can be 
calculated as follow [6]: 

















−×= −
T

DW
1

303
1

2416exp103.1 10  (23) 

The variable Wc  represents the number of 
water molecules per sulfonic acid group (i.e. mol 
H2O/equivalent SO3

-1). The water content in the 
electrolyte phase is related to water vapour 
activity via [8, 9]: 

( )
( ) ( )

( )3                                                         8.16
31                                      14.10.14
10        0.3685.3981.17043.0 32

≥=
≤<−+=
≤<+−+=

ac
aac
aaaac

W

W

W  (24) 

The water vapour activity given by: 

sat

W
P

Px
a =  (25) 

For heat transfer purposes, the membrane is 
considered a conducting solid, which means that 
the transfer of energy associated with the net 
water flux through the membrane is neglected. 
Heat transfer in the membrane is determined by: 

( ) 0=∇⋅⋅∇ Tkmem  (26) 
The potential loss in the membrane is due to 

resistance to proton transport across membrane, 
and determined by: 

( ) 0=∇⋅∇ φλm  (27) 
5) Cell Potential 

Electrical energy is obtained only when a 
current is drawn, but the actual cell potential 
(Ecell) is decreased from its equilibrium 
thermodynamic potential (E) because of 
irreversible losses. The cell potential is obtained 
by subtracting all over potentials (losses) from 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

11 

the equilibrium thermodynamic potential as the 
following expression [16, 17]: 

Diffmemohmactcell EE ηηηη −−−−=  (28) 

The equilibrium potential is determined using 
the Nernst equation [18]: 

𝐸 = 1.229 − 0.83𝑥10−3(𝑇 − 298.15) +
4.3085𝑥10−5𝑇�𝑙𝑛�𝑃𝐻2� +

1
2
𝑙𝑛�𝑃𝑂2�� (29) 

The anode and cathode activation 
overpotentials are calculated from Butler-Volmer 
equation. The ohmic overpotentials in GDLs and 
protonic overpotential in membrane are 
calculated from the potential equations, 
respectively. 

The anode and cathode diffusion 
overpotentials are calculated as follows: 











−=

cL

c
cDiff i

i
F

RT

,
, 1ln2

η

; 










−=

aL

a
aDiff i

i
F

RT

,
, 1ln2

η  (30) 

GDL

OO
cL

CFD
i

δ
22

2
, =

; GDL
HH

aL
CFD

i
δ

22
2

, =  (31) 

 
C. Hygro-Thermal Stresses in Fuel Cell 

As a result of the changes in temperature and 
moisture, all the PEM, GDL, and bipolar plates 
will experience expansion and contraction. 
Because of the different thermal expansion and 
swelling coefficients between these materials, 
hygro-thermal stresses are expected to be 
introduced. In addition, the non-uniform current 
and reactant flow distributions in the cell can 
result in non-uniform temperature and moisture 
content of the cell, which could in turn, 
potentially causing localized increases in the 
stress magnitudes. 

Using hygrothermo elasticity theory, the 
effects of temperature and moisture as well as the 
mechanical forces on the behavior of elastic 
bodies have been addressed. An uncoupled 
theory is assumed, for which the additional 
temperature changes brought by the strain are 
neglected. The total strain tensor is determined 
using the following expression [11]: 

STM ππππ ++=  (32) 

where Mπ  is the contribution from the 
mechanical forces and ,T Sπ π  are the thermal 
and swelling induced strains, respectively. 

The thermal strains resulting from a change in 
temperature of an unconstrained isotropic volume 
are given by: 

( )fT TT Re−℘=π  (33) 
The swelling strains caused by moisture 

change in membrane are given by: 

( )fmemS Reℜ−ℜ= π  (34) 
Assuming linear response within the elastic 

region, the isotropic Hooke's law is used to 
determine the stress tensor: 

πσ G=  (35) 

where G is the constitutive matrix. The effective 
stresses according to von Mises, 'Mises stresses', 
are given by: 

( ) ( ) ( )
2

2
13

2
32

2
21 σσσσσσσ

−+−+−
=v

 (36) 

where 1σ , 2σ , 3σ  are the principal stresses.  
 

D. Boundary Conditions 
In order to reduce computational cost, the 

advantage of the geometric periodicity of the cell 
is taken. Hence symmetry is assumed in the y-
direction, i.e. all gradients in the y-direction are 
set to zero at the x-z plane boundaries of the 
domain. With the exception the channel inlets 
and outlets, a zero flux condition is applied at all 
x-boundaries (y-z planes). 

The inlet values at the anode and cathode are 
prescribed for the velocity, temperature, and 
species concentrations (Dirichlet boundary 
conditions). The gas streams entering the cell are 
fully humidified, but no liquid water is contained 
in the gas stream. The inlet velocities of air and 
fuel are calculated based on the desired current 
density according to: 

chinc

cin

inO
MEAccin AP

RT
x

A
F
I

u
11

4 ,
,

,
,

2

ζ=  (37) 

china

ain

inH
MEAaain AP

RT
x

A
F
I

u
11

2 ,
,

,
,

2

ζ=  (38) 

At the outlets of the gas-flow channels, only 
the pressure is being prescribed as the desired 
electrode pressure, for all other variables, the 
gradient in the flow direction (x) is assumed to be 
zero (Neumann boundary conditions). At the 
external surfaces in the z-direction, (top and 
bottom surfaces of the cell), temperature is 
specified and zero heat flux is applied at the x-y 
plane of the conducting boundary surfaces. 
Combinations of Dirichlet and Neumann 
boundary conditions are used to solve the 
electronic and protonic potential equations. 
Dirichlet boundary conditions are applied at the 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

12 

land area. Neumann boundary conditions are 
applied at the interface between the gas channels 
and the GDLs to give zero potential flux into the 
gas channels. Similarly, the protonic potential 
field requires a set of potential boundary 
condition and zero flux boundary condition at the 
anode CL interface and cathode CL interface 
respectively. The mechanical boundary 
conditions are noted in Figure 1. The initial 
conditions corresponding to zero stress-state are 
defined, all components of the cell stack are set 
to reference temperature 20°C, and relative 
humidity 35% [11, 19]. In addition, a constant 
pressure of 1 MPa is applied on the surface of 
lower plate, corresponding to a case where the 
stack is clamped with springs to control the force. 

 
III. SOLUTION ALGORITHM 

The governing equations were discretized 
using a finite volume method and solved using a 
general-purpose computational fluid dynamic 
(CFD) code. Stringent numerical tests were 
performed to ensure that the solutions were 
independent of the grid size. A computational 

mesh of 150,000 computational cells was found 
to provide sufficient spatial resolution.  

Figure 2 shows flow diagram of the solution 
algorithm used in this study. It begins by 
specifying a desired current density of the cell for 
calculating inlet flow rates at the anode and 
cathode sides. An initial guess of the activation 
over potential is obtained from the desired 
current density using the Butler-Volmer equation, 
then follows by computing the flow fields for 
each phase for velocities u, v, w, and pressure P. 
Once the flow field is obtained, the mass fraction 
equations are solved for oxygen, hydrogen, 
nitrogen, and water. Scalar equations are solved 
last in the sequence of the transport equations for 
the temperature field in the cell and potential 
fields in the GDLs and the membrane. 

Hooke's law with total strain tensor is solved 
to determine the stress tensor. The local current 
densities are solved based on the Butler-Volmer 
equation, then local activation over potentials can 
be calculated from the Butler-Volmer equation. 
They are updated after each global iterative loop. 
Convergence criteria are performed on each 
variable and the procedure is repeated until 
convergence. The properties are updated after 
each global iterative loop based on the new local 
gas composition and temperature. Source terms 
reflect changes in the overall gas phase mass 
flow due to consumption or production of gas 
species via reaction and due to mass transfer 
between water in the vapor phase and water that 
is in the liquid phase or dissolved in the polymer.  

The set of equations was solved iteratively, 
and the solution was considered to be convergent 
when the relative error in each field between two 
consecutive iterations was less than 10-6. The 
number of iterations required to obtain converged 
solutions dependent on the nominal current 
density of the cell, the higher the load the slower 
the convergence. 

 
Figure 2. Flow diagram of the solution procedure used 

 
Figure 3. Effect of liquid water build up on cell performance 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

13 

IV. RESULTS AND DISCUSSION 
The values of the electrochemical transport 

parameters for the base case operating conditions 
are listed in Table 1. The geometric and the base 
case operating conditions are listed in Table 2. 
Since this model accounts for all major transport 
processes and the domain comprises all the 
elements of a complete cell, no parameters 
needed to be adjusted. Performance curves with 
and without phase change as well as experimental 
data are shown in Figure 3 for the base case 

conditions. The experimental data is provided by 
Wang et al. [14]. Comparison of the two curves 
demonstrates that the model represents better cell 
potential (V) versus current density (A/cm2) 
curve than single phase model. The effects of 
liquid water accumulation become apparent even 
at relatively low values of current density. This 
result validated the developed model. Results 
with phase change for the cell operates at 
nominal current density of 1.4 A/cm2 under the 
base operating conditions are analyzed and 

Table 1. 
Electrode and membrane parameters for base case operating conditions 

Parameter [Ref.] Symbol Value Unit 
Electrode porosity [16] ε  0.4 - 
Electrode electronic conductivity [5] eλ  100 mS /  

Membrane ionic conductivity (humidified Nafion117) [16] mλ  17.1223 mS /  

Transfer coefficient, anode side [17] aα  0.5 - 

Transfer coefficient, cathode side [15] cα  1 - 

Cathode reference exchange current density [3] refcoi ,  1.81e-3 2/ mA  

Anode reference exchange current density [3] refaoi ,  2465.6 2/ mA  

Electrode thermal conductivity [4] effk  1.3 KmW ./  

Membrane thermal conductivity [4] memk  0.455 KmW ./  

Electrode hydraulic permeability [14] kp  1.76e-11 2m  
Entropy change of cathode side reaction [13] S∆  -326.36 KmoleJ ./  

Heat transfer coefficient between solid and gas phase [3] β  4e6 3/ mW  

Protonic diffusion coefficient [16] +HD  4.5e-9 sm /2  
Fixed-charge concentration [16] fc  1200 3/ mmole  

Fixed-site charge [16] fz  -1 - 

Electro-osmotic drag coefficient [18] dn  2.5 - 

Droplet diameter [7] dropD  1e-8 m  

Condensation constant [7] C  1e-5 - 
Scaling parameter for evaporation [7] ϖ  0.01 - 
Electrode Poisson's ratio  GDLℑ  0.25 - 

Membrane Poisson's ratio [11] memℑ  0.25 - 

Electrode thermal expansion [11] GDL℘  -0.8e-6 K1  

Membrane thermal expansion [19] mem℘  123e-6 K1  

Electrode Young's modulus [11] GDLΨ  1e10 Pa  

Membrane Young's modulus [19] memΨ  249e6 Pa  

Electrode density [11] GDLρ  400 
3mkg  

Membrane density [19] memρ  2000 
3mkg  

Membrane humidity swelling-expansion tensor [11] mem  23e-4 %1  

 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

14 

discussed in this section. The selection of 
relatively high current density is due to illustrate 
the phase change effects, where it becomes 
apparent between single and multi-phase model 
in the mass transport limited region. 

The velocity fields inside the GDLs are shown 
in Figures 4 and 5. The liquid water saturation 
inside the GDL are shown in Figure 6, and the 
profile of water content in the membrane is 
shown in Figure 7. These results are newly 
disclosed here and they have not been presented 
in the previous papers [20-24]. From Figure 4 
and 5 it is clear that the pressure gradient induces 
bulk gas flow from the channels into the GDL 
while the capillary pressure gradient drives the 
liquid water out of the GDL into the flow 
channels. The velocity of the liquid phase is 
lower than the gas phase, and the highest liquid 
water velocity occurs at the corners of the 
channel/GDL interface.  

The liquid water oozes out of the GDL, 
mainly at the corners of the GDL/channel 
interface. From Figure 6 it can be observed that 
condensation occurs mainly in two areas inside 
the cathodic GDL i.e. at the catalyst layer the 
molar water vapour fraction increases due to the 
oxygen depletion, and at the channel/GDL 
interface, where the oversaturated bulk flow 
condenses out. Condensation occurs throughout 

Table 2. 
Parameters for base case conditions 

Parameter Symbol Value Unit 
Channel length L 0.05 m  
Channel width W 1e-3 m  
Channel height H 1e-3 m  

Land area width landW  1e-3 m  

GDL thickness GDLδ  0.26e-3 m  

Membrane thickness (Nafion 
117) mem

δ  0.23e-3 m  

Catalyst thickness CLδ  
0.028e-
3 

m  

H2 reference mole fraction 
ref
Hx 2  0.84639 - 

O2 reference mole fraction 
ref
Ox 2  0.17774 - 

Anode pressure aP  3 atm  

Cathode pressure cP  3 atm  

Inlet fuel and air temperature cellT  353.15 K 
Relative humidity of inlet fuel 
and air  

ψ  100 % 

Air stoichiometric flow ratio cξ  2 - 

Fuel stoichiometric flow ratio aξ  2 - 
 

 
(a) 

 
(b) 

Figure 4. Velocity vectors inside the cathode GDL: (a) Gas 
phase velocity vectors; (b) Liquid water velocity vectors 

 
(a) 

 
(b) 

Figure 5. Velocity vectors inside the anode GDL: (a) Gas 
phase velocity vectors; (b) Liquid water velocity vectors 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

15 

the anodic GDL due to hydrogen depletion. 
Similar to the cathode side, the liquid water can 
only leave the GDL through the build-up of a 
capillary pressure gradient to overcome the 
viscous drag, because at steady state operation, 
all the condensed water has to leave the cell.  

The liquid water saturation is distributed 
through the entire GDL with maximum 
saturations found under the land areas. The high 
spatial variation of the saturation demonstrates 
again the three-dimensional nature of transport 
processes in PEMFC. Clearly, liquid water 
saturation depends strongly on the specified 

 
Figure 6. Liquid water saturation inside the GDLs 

 
Figure 7. Water content profile through the MEA 

 
Figure 8. Water vapour molar fraction distribution in the cell at a nominal current density of 1.4 A/cm2 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

16 

capillary pressure, and the permeability of the 
GDL becomes the central parameter. In the 
membrane, primary transport is through electro-
osmotic drag associated with the protonic current 
in the electrolyte, which results in water transport 
from anode to cathode, and diffusion associated 
with water-content gradients in the membrane. 
From Figure 7 it can be noticed that the influence 
of electro-osmotic drag and back diffusion are 
apparent.  

Water vapour molar fraction distribution is 
shown in Figure 8. It remains almost constant 
throughout the GDL. In the absence of phase 
change, this would not be the case, as the 
nitrogen and water vapour fraction would 
increase as the oxygen fraction decreases. 
Reactant gas distribution in the cell cannels and 
the GDLs demonstrate similar distribution with 
those under operation at current density of 0.8 
A/cm2 [24]. They are not plotted in this paper. 
The molar fraction of oxygen is highest under the 
channel, and exhibits a three-dimensional 
behavior with a fairly significant drop under the 

land areas, and a more gradual depletion towards 
the outlet. The molar hydrogen fraction is almost 
constant inside the GDL due to the higher 
diffusivity of the hydrogen. 

The local current density distribution is shown 
in Figure 9. This result resembles the result 
reported earlier when the fuel cell was operated 
at nominal density current of 1.2 A/cm2 under  
temperature of 60°C and 90°C [20]. The 
differences exist at its magnitude. The maximum 
value in this study is 1.5 while that in the 
previous study was 1.4. The local current 
densities have been normalized by divided 
through the nominal current density (i.e. Iic ). 
In the model it has a much higher fraction of the 
total current, generated under the channel area. 
This is due to the effects of liquid water inside 
the GDL which decreases its permeability to 
reactant gas flow and lead to the onset of pore 
plugging by liquid water. This can lead to local 
hot-spots inside the membrane electrode 
assembly. This leads to a further drying out of the 
membrane and increasing the electric resistance, 

 

Figure 9. Dimensionless local current density distribution ( Iic ) at the cathode side catalyst layer 

 
Figure 10. Temperature distribution inside the fuel cell 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

17 

which in turn leads to more heat generation and 
can lead to a failure of the membrane. It is 
important to keep the current density relatively 
even throughout the cell. Therefore, the model is 
capable of identifying important parameters for 
the wetting behavior of the GDLs and can be 
used to identify conditions that might lead to the 
onset of pore plugging. The temperature 
distribution inside the cell has important effects 
on nearly all transport phenomena, and 
knowledge of the magnitude of temperature 
increases might help preventing failure.  

Figure 10 shows temperature distribution 
throughout the x-axis inside the cell at a nominal 
current density of 1.4 A/cm2. Temperature 
distribution only on the y-z plane at x=10 mm 
has been published earlier [23]. From Figure 10 it 
can be observed that the increase in temperature 
can exceed several degrees Kelvin near the inlet 
area, where the electrochemical activity is highest. 
The temperature peak appears in the cathode CL. 
In general, the temperature at the cathode side is 
higher than at the anode side, this is due to the 
reversible and irreversible entropy production. 

Figure 11 shows von Mises stress distribution 
inside the membrane during the cell operating at 
base case conditions and a nominal current 
density of 1.4 A/cm2. The similar distribution has 
been published when the cell was operated at 
nominal density current of 1.2 A/cm2 [21]. 
Compared to the previous results Figure 11 
shows larger stresses. The maximum stress 
occurs, where the temperature is highest, which is 
near the cathode side inlet area. 

The Membrane-Electrode-Assembly (MEA) 
is the core component of PEMFC and consists of 
membrane with the GDLs including the catalyst 
attached to each side. Von Misses stress 
distribution (contours plots) in the cell on the y-z 
plane at x=10 mm at nominal current density of 
1.4 A/cm2 has been published [23]. In this paper 

the close up of the same distribution is plotted in 
Figure 12 (scale enlarged 300 times). Because of 
the different thermal expansion and swelling 
coefficients between GDLs and membrane 
materials, hygro-thermal stresses and 
deformation are introduced. It induces bending 
stresses which can contribute to delaminating 
between the membrane and the GDLs. 

The distribution of total displacement and 
deformed shape for MEA at base conditions on 
the y-z plane at x=10 mm are similar to those at 
nominal current density of 1.2 A/cm2 [21]. 
However, the magnitude is different. At nominal 
current density of 1.4 A/cm2 larger displacement 
can be seen. It is related to the temperature where 
the temperature is highest in the centre of the 
channel and coincide with the highest reactant 
concentrations. The deformation under the land 
areas is much smaller than under the channel 
areas due to the clamping force effect. This result 
may explain the occurrence of cracks and pin 
holes in the membrane under steady–state 
loading during regular cell operation.  

Fuel cell losses originate from sources such as 
activation overpotential, ohmic overpotential in 
GDL, membrane overpotential, and diffusion 
overpotential. When the cell operates at nominal 
current density of 0.8 A/cm2 all of these losses 
have been published [24]. Here, only activation 
losses are plotted in Figure 13. Table 3 shows 
comparison of the maximum values of these 
losses (in Volt) between the two cases of 

 
Figure 11. Mises stresses distribution inside the membrane 

Table 3. 
Comparison of maximum losses 

Source of losses Current density:  
0.8 (A/cm2) 

Current density: 
1.4 (A/cm2) 

Activation 0.398 (volt) 0.387 (volt) 
Ohmic 0.060 (volt) 0.086 (volt) 
Membrane 0.180 (volt) 0.250 (volt) 
Diffusion 0.065 (volt) 0.040 (volt) 
 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

18 

operating nominal current density. The activation 
overpotential loss profile at the cathode side CL 
correlates with the local current density. The 
current densities are highest in the centre of the 
channel and coincide with the highest reactant 
concentrations. It can be seen that the loss is 
directly related to the nature of the 
electrochemical reactions. The reaction 
propagates at the rate demanded by the current. 

 
V. CONCLUSION 

The developed model has been validated by 
comparing its results to experimental data. It has 
been proved that the model represents better cell 
potential (V) versus current density (A/cm2) 
curve than single phase model. In the cathodic 
and anodic GDLs, the velocity of the liquid phase 
is lower than the gas phase, and the highest liquid 
water velocity occur at the corners of the 
channel/GDL interface. The liquid water 
saturation is distributed through the entire GDL 
with maximum saturations found under the land 
areas. The molar water vapour fraction remains 
almost constant throughout the GDL. The molar 
fraction of oxygen is highest under the channel, 
and exhibits a three-dimensional behavior with a 
fairly significant drop under the land areas, and a 

more gradual depletion towards the outlet. The 
molar hydrogen fraction is almost constant inside 
the GDL due to the higher diffusivity of the 
hydrogen. The model has been used to analysis 
important variables at nominal current density of 
1.4 A/cm2. When compared to previous results at 
nominal current density of 1.2 A/cm2, the results 
in this paper shows larger local current density, 
larger von Mises stresses, and larger 
displacement. When compared to previous results 
at nominal current density of 0.8 A/cm2, the 
results in this paper show smaller activation 
losses, larger ohmic losses, larger membrane 
losses, and smaller diffusion losses. The model is 
shown to be able to: (1) understand complex 
phenomena which cannot be studied 
experimentally, (2) identify limiting steps and 
components, and (3) provide a computer-aided 
tool for design and optimization of PEMFC with 
much higher power density and lower cost.  

 
REFERENCES 
[1] E. E. Kahveci and I. Taymaz, 

"Experimental investigation on water and 
heat management in a PEM fuel cell using 
response surface methodolog," 

 
Figure 12. Von Misses stress distribution in MEA 

 
Figure 13. Activation overpotential loss distribution 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

19 

International Journal of Hydrogen Energy, 
39(20), pp. 10655-10663, 2014. 

[2] A. P. Sasmito et al., "Numerical 
evaluation of various thermal 
management strategies for polymer 
electrolyte fuel cell stacks," International 
Journal of Hydrogen Energy, 36, pp. 
12991-13007, 2011. 

[3] H. Meng, "Numerical studies of liquid 
water behaviors in PEM fuel cell cathode 
considering transport across different 
porous layer," International Journal of 
Hydrogen Energy, 35(11), pp. 5569–5579, 
2010. 

[4] N. K. H. Dalasm et al., "A parametric 
study of cathode catalyst layer structural 
parameters on the performance of a PEM 
fuel cell," International Journal of 
Hydrogen Energy, 35(6), pp. 2417–2427, 
2010. 

[5] S. W. Perng, and H. W. Wu, "Non-
isothermal transport phenomenon and cell 
performance of a cathodic PEM fuel cell 
with a baffle plate in a tapered channel," 
Applied Energy, 88(1), pp. 52–67, 2011. 

[6] A. Iranzo et al., "Numerical model for the 
performance prediction of a PEM fuel 
cell. Model results and experimental 
validation," International Journal of 
Hydrogen Energy, 35(20), pp. 11533–
11550, 2010. 

[7] N. Djilali, "Computational Modelling of 
PEM Fuel Cells: Challenges and 
Possibilities," Energy, 32(4), pp. 269-280, 
2007. 

[8] M. Hu et al.,  “Three dimensional, two 
phase flow mathematical model for PEM 
fuel cell. Part I. Model development,” 
Energy Conversion Manag, 45(11-12): 
1861–1882, 2004. 

[9] C. Fink, and N. Fouquet, "Three-
dimensional simulation of polymer 
electrolyte membrane fuel cells with 
experimental validation," Electrochimica 
Acta, 56(28), pp. 10820–10831, 2011. 

[10] A. Webber, and J. Newman, “Theoretical 
Study of Membrane Constraint in 
Polymer-Electrolyte Fuel Cell,” AIChE J., 
50(12), 3215–3226, 2004. 

[11] D. Bograchev et al., "Stress and plastic 
deformation of MEA in running fuel cell," 
International Journal of Hydrogen 
Energy, 33, pp. 5703–5717, 2008. 

[12] M. A. R. S. Al-Baghdadi, "Performance 
comparison between planar and tubular-
shaped ambient air-breathing PEM fuel 
cells using three-dimensional 

computational fluid dynamics models," 
Journal of Renewable and Sustainable 
Energy; 1(2), pp. 023105-1-023105-15, 
2009. 

[13] M. A. R. S. Al-Baghdadi, "A CFD 
analysis of transport phenomena and 
electrochemical reactions in a tubular-
shaped ambient air-breathing PEM micro 
fuel cell," The Hong Kong Institution of 
Engineers Transactions (HKIE 
Transactions), 17(2), pp. 1-8, 2010. 

[14] M. A. R. S. Al-Baghdadi, “PEM Fuel 
Cells - from Single Cell to Stack,” 
International Energy and Environment 
Foundation (IEEF), 2015, ISBN-13: 
9781505885644. 

[15] M. A. R. S. Al-Baghdadi. "Novel design 
of a compacted micro-structured air-
breathing PEM fuel cell as a power source 
for mobile phones," International Journal 
of Energy and Environment; 1(4), pp. 
555-572, 2010. 

[16] M. A. R. S. Al-Baghdadi, "Analysis of 
transport phenomena and electrochemical 
reactions in a micro PEM fuel cell," 
International Journal of Energy and 
Environment, 5(1), pp. 1-22, 2014. 

[17] M. A. R. S. Al-Baghdadi, "Analysis of 
transport phenomena and electrochemical 
reactions in a micro PEM fuel cell with 
serpentine gas flow channels," 
International Journal of Energy and 
Environment, 5(2), pp. 139-154, 2014. 

[18] M. A. R. S. Al-Baghdadi, "Analysis of 
transport phenomena and electrochemical 
reactions in a micro PEM fuel cell with 
nature inspired flow field design," 
International Journal of Energy and 
Environment; 6(1), pp. 1-16, 2015. 

[19] Product Information, DuPont™ Nafion® 
PFSA Membranes N115, N117, N1110 
Perfluorosulfonic Acid Polymer. 
NAE101, 2016. 

[20]  M. A. R. S. Al-Baghdadi, and Haroun 
A.K.S. Al-Janabi, "Parametric and 
optimization styudy of a PEM fuel cell 
performance using three-dimensional 
computational fluid dynamics model," 
Renewable Energy (2007), pp. 1077-1101, 
32. 

[21]  M. A. R. S. Al-Baghdadi, and Haroun 
A.K.S. Al-Janabi, "Effect of operating 
parameters on the hydro-thermal stresses 
in proton exchange membrane of fuel 
cells," International Journal of Hydrogen 
Energy, 32, pp. 4510-4522, 2007. 



M.A.R.S. Al-Baghdadi / J. Mechatron. Electr. Power Veh. Technol 07 (2016) 7-20 
 

 

20 

[22] M. A. R. S. Al-Baghdadi, and H. A. K. S. 
Al-Janabi, "Modeling optimizes PEM 
fuel cell performance using three-
dimensional multi-phase computational 
fluid dynamics model," Energy 
Conversion and Management, 48, pp. 
3102-3119, 2007. 

[23] M. A. R. S. Al-Baghdadi, "A CFD study 
of hygro-thermal stresses distribution in 
PEM fuel cell during regular cell 
operation," Renewable Energy, 34, pp. 
674-682, 2009. 

[24] M. A. R. S. Al-Baghdadi, "Performance 
comparison between air flow-channel and 
ambient air-breathing PEM fuel cells 
using three-dimensional computation fluid 
dynamics models," Renewable Energy, 
34, pp. 1812-1824, 2009. 

 
NOMENCLATURE 
a  water activity 

MEAA  geometrical area of MEA (m
2) 

chA  cross sectional area of flow channel (m
2) 

C  local concentration (mole/m3) 
refC  reference concentration (mole/m3) 

Cp  specific heat capacity (J/kg.K) 

D  diffusion coefficient (m2/s) 
E  equilibrium thermodynamic potential (V) 

cellE  cell operation potential (V) 
F  Faraday’s constant = 96487 (C/mole) 
I  cell operating current density (A/m2) 
i  local current density (A/m2) 

ref
oi ,  reference exchange current density 

effk  effective electrode thermal conductivity 
(W/m.K) 

memk  membrane thermal conductivity (W/m.K) 
kp  hydraulic permeability (m2) 

phasem  mass transfer kg/s) 

M  molecular weight (kg/mole) 

WN  net water flux across the membrane (kg/m
2.s) 

nd  electro-osmotic drag coefficient 

en  number of electrons transfer 
P  pressure (Pa) 

cP  capillary pressure (Pa) 
q  heat generation (W/m2) 

R  universal gas constant = 8.314 (J/mole.K) 
s  specific entropy (J/mole.K) 
sat  saturation 
T  temperature (K) 
u  velocity vector (m/s) 

ix  mole fraction 

iy  mass fraction 

aα  charge transfer coefficient, anode side 

cα  charge transfer coefficient, cathode side 
ε  porosity 
η  over potential (V) 

eλ  electrode electronic conductivity (S/m) 

mλ  membrane ionic conductivity (S/m) 
µ  viscosity (kg/m.s) 
ρ  density (kg/m3) 
τ  surface tension (N/m) 
ξ  stoichiometric flow ratio  
σ  stress (N/m) 
π  strain 
ℜ  relative humidity (%) 

 


	I. Introduction
	II. Model Description
	A. Computational Domain
	B. Model Equations
	1) Gas Flow Channels
	2) Gas Diffusion Layers (GDL)
	3) Catalyst Layers (CLs)
	4) Membrane
	5) Cell Potential

	C. Hygro-Thermal Stresses in Fuel Cell
	D. Boundary Conditions

	III. Solution Algorithm
	IV. Results and Discussion
	V. Conclusion
	References
	Nomenclature