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 Nafion117) [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