Excerpt from the Proceedings of the COMSOL Conference 2007, Boston

COMSOL Modelling and Simulation of PEM Fuel Cell’s

Flow Channels
Eric Robalinho1,Edgar F Cunha1, Alexandre B Andrade1,
Martha L M Bejarano1, Marcelo Linardi1, Efraim Cekinski2
1
Instituto de Pesquisas Energéticas e Nucleares (IPEN / CNEN - SP)
Av. Professor Lineu Prestes 2242
05508-000 São Paulo, SP
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
2
Instituto de Pesquisas Tecnológicas (IPT - SP)
Av. Prof. Almeida Prado, 532
05508-901 São Paulo, SP
[email protected]

Abstract: Fuel cells are one of the most accurately solve the equations over the solution
important devices to obtain electrical energy domain. Process simulations of diffusion in
from hydrogen. The Computational Fluid porous media and conduction/convection were
Dynamic (CFD) is a very useful tool to explore studied. Once all of the global iterations had
the connection between the transport of reactants been completed, the experimental data were
and products within the fuel cell and overall cell compared with modelling results. The graphite
performance. CFD simulates hydrogen and plate B exhibited better performance than the A
oxygen gases flow channels to reduce the costs one, that was validated by numeric methods.
of bipolar plates’ production through Such behavior was attributed to the higher
optimization of the mass transport. In this work, channel area in contact with gas diffusion
two configurations of gas flow channels were electrode which provides higher active area
studied by simulation with COMSOL software available for the gas reaction. Besides the area,
and a preliminary experimental test to verify the higher gas velocities in flow fields decrease the
model applied in these configurations was gas stagnation points along the channels.
performed. Two single Proton Exchange
Membrane Fuel Cells – PEMFC – with 144 cm² Keywords: Numerical model, CFD simulation,
of active area were manufactured. The first one PEMFC, hydrogen technology.
(A) was a commercial graphite plate supplied by
ELECTROCELL. The other one (B) was 1. Introduction
projected and constructed at IPEN. The
serpentine flow field was chosen. The Fuel cells and hydrogen technology represent
dimensions of flow field of prototype A were: 2 the most promising alternative pathway for
mm width and 2 mm depth channels and 2 mm automotive and stationary applications. The
ribs. In this case each part of the serpentine flow PEMFC offers low to zero emission for subWatt
field was composed by 6 channels and 5 ribs. to MWatt power generation, meaning
The dimensions of flow field of prototype B applications in transportation, industries and
were: 1.5 mm width and 1.5 depth channels and portable supplies units [1, 2]. In this paper we
0.5 mm ribs. For this prototype, each part of focus on the role of computational tools in order
serpentine flow field was composed by 12 to verify some experimental results and
channels and 11 ribs. A mesh study was demonstrate the better performance of PEMFC
performed on the 3D models to determine the constructed from optimizations techniques.
minimum number of elements required to
 Excerpt from the Proceedings of the COMSOL Conference 2007, Boston

The PEMFC is an electrochemical device 2. Methodology

and consists of two principal parts: (i) an
electrolyte between two electrodes, an anode and Two single cells of 144 cm² of active area
a cathode with two diffusion layers, called were manufactured, called prototype 1 and
Membrane Electrode Assembly (MEA) and (ii) prototype 2. The serpentine flow field was
gas flow field plates, called bipolar plates, in chosen for channel geometry. The configuration
case of the stack. The plates have many functions of flow field in the prototype 1 was: 2 mm width
in a fuel cell [3]: distribute reactant gases and 2 mm depth channels and 2 mm ribs. In this
(hydrogen and air or oxygen) uniformly, collect case each part of the serpentine flow field was
and conduct electrical current and remove heat composed by 6 channels and 5 ribs. The flow
and water from the electrodes. The bipolar plates field of prototype 2 was: 1.5 mm width and 1.5
make the cell robust and give a rigid structure to depth channels and 0.5 mm ribs. For this
support the impacts of portable and automotive prototype, each part of serpentine flow field was
applications. The cost of the bipolar plates composed by 12 channels and 11 ribs. In Figure
corresponds up to 45% of the total cost and 1 (a) and (b) the bipolar plates of prototype 1 and
about 80% of the total weight in a stack fuel cell prototype 2 were showed, respectively. For these
[4, 5]. flow field configurations, the prototype 2 has
The search for reliable computational models 50 % more channel area than prototype 1 (0.110
is a challenge because it involves several m2 and 0.073 m2 respectively).
transport phenomena: multi-component, multi- The best-known membrane material for PEM
phase and multi-dimensional flow processes, fuel cell is Nafion™ membrane by Dupont,
electrochemical reactions, convective heat and which has high proton conductivity. Besides the
mass transport in flow channels, diffusion of membrane, the MEA is composed by two
reactants through porous electrodes, transport of electrodes as catalytic layers, using platinum
water through the membrane and transport of supported in carbon black as electro catalysts;
electrons through solid matrix. and two gas diffusion layers (GDL) like carbon
The Computational Fluid Dynamic (CFD) is paper or carbon cloth. Fuel cell polarization
a very useful tool to simulate hydrogen and measurements were carried out galvanostatically
oxygen gases flow channels configurations, with the single cell at 70 °C, using oxygen and
reducing the costs of bipolar plates’ production hydrogen with saturated water vapor (Millipore
and optimizing mass transport [6-8]. quality – Elix 3 model) at 85 °C, atmospheric
pressure and hydrogen volumetric flux
8.33 10-6 m3 s-1.

(a) (b)

Figure 1: Flow field channels from prototype 1 (a) and prototype 2 (b).
 Excerpt from the Proceedings of the COMSOL Conference 2007, Boston

3. Numerical Models

In this paper we concentrate efforts in two

goals. The first work was the analysis of four
partial geometries consisting of the inlet and two
channels with the porous media. Then, in the
second place, the work was to simulate the two
graphite plates in order to investigate hydrogen
flow, pressure drops and temperature
distribution. Implementations of the modeling
allow us optimizing the design of the plates and
better understanding the fuel cell dynamics. Figure 4: Partial 2 mm square geometry.
The Nonlinear solver for stationary analysis
was set automatically by COMSOL, and the
Spooles direct solver was selected for linear
conditions [9]. The Spooles uses the multifrontal
method and direct LU factorization of the sparse
matrix A when solves systems of the form Ax=b.
There are several solver methods available at
COMSOL Multiphysics, including the iterative
ones that are the natural choice for problems
with many degrees of freedom. However the
iterative solvers are less stable than direct solvers
and they do not always converge.
Figure 5: Partial 1.5 mm square geometry.

The geometries used were showed in

Figures 2-5. The red region represents the gas
inlet, the blue one is the distribution place and
the green is the initial part of two channels. The
area above the channels (0.25 mm thickness
transparent rectangle) is the GDL. The triangular
and the trapezoidal geometries have 1.5 mm
depth. In the Figures 2-5 the x, y and z axis are
showed.
At the Table 1 the different meshes and
corresponding solutions time were compared.
Figure 2: Partial triangular geometry. The convergence in these cases was obtained in
3 to 7 iterations and the relative tolerance at
nonlinear solver was set on 1.0E-6. The
refinement method was longest, with the best
quality option. The mesh in each simulation is
defined by the Number of Degrees of Freedom
(NOF’S), the Number of Elements (NE) and the
Minimum Element Quality (MEQ). As one can
see at Table 1, the refinement of the mesh means
much more solution time and availability of
memory.
The flow distribution is described by two
equations: Navier-Stokes (Eq.1.) and continuity
Figure 3: Partial trapezoidal geometry. (Eq.2.),
 Excerpt from the Proceedings of the COMSOL Conference 2007, Boston

ρ
∂u
∂t
( T
)
−∇.η ∇u +(∇u) + ρ (u.∇)u +∇p = 0 (1)

∇ .u = 0 (2)

where ρ denotes the density [kg m-3], u the

velocity vector [m s-1], η the viscosity [N s m-2],
and p the pressure [Pa]. The density and
viscosity of the modeled gases were calculated in
function of the pressure and temperature using
models offered for the CFD software. The model
uses normal/pressure boundary conditions for
inlet and convective flux for outlet. The
condition for porosity (porous media) was found
in Navier-Stokes module, which was set 1.0E-12
for permeability of GDL.

Table 1: Different meshes and corresponding NOF’S, NE, MEQ and SOLUTION TIME.

MESH NOF’S NE MEQ SOLUTION

TIME (s)
Coarser 7351 957 0.2829 9.1
Coarse 10808 1464 0.2990 14.6
TRIANGULAR Normal 16111 2298 0.2500 27.4
Fine 30187 4498 0.3756 98.7
Finer 83234 13081 0.3558 1158.3
Coarser 8305 996 0.3339 12.7
Coarse 13399 1689 0.3566 29.1
TRAPEZOIDAL Normal 21845 3064 0.3656 113.6
Fine 45315 6560 0.2509 551.8
Finer 119717 18594 0.3401 5239.1
Coarser 6729 769 0.3080 9.2
Coarse 10425 1273 0.2944 17.7
SQUARE (2mm) Normal 19148 2563 0.3498 56.2
Fine 29260 4144 0.3696 149.8
Finer 83845 12718 0.3422 470.1
Coarser 7496 937 0.3230 12.7
Coarse 14962 1867 0.2823 42.3
SQUARE (1.5mm) Normal 31808 4509 0.3195 241.5
Fine 56461 8369 0.3868 908.5
Finer 158145 24959 0.2715 17592.2
 Excerpt from the Proceedings of the COMSOL Conference 2007, Boston

4. Results & Discussions

Pressure lines graphics were plotted for each

profile studied. For one predetermined channel
length, y =3.50mm, six pressure lines were
plotted at the GDL interior and one at the
channel outlet (0.05mm from the upper surface
of the channel), in order to evaluate the pressure
gradient there. The pressure gradient behavior of
the studied models provides homogeneity of the
reactants gas flux data through diffusion layer. A
homogeneous pressure gradient means higher
gas availability for the fuel cell reactions.
According to the graphics shown, a huge Figure 7: Pressure lines of trapezoidal profile.
stagnation region (dead zone) is noticed at the
top of rib between the 2mm square profile
channels. At Figure 8 this effect is noticed from
the inexistence of a high pressure gradient at
2.50mm < x < 3.50mm range (every pressure
lines were close to zero).
The 1.5 square profile presents stagnation,
but it is local surrounding x=1.75mm (Figure 9).
The trapezoidal and triangular profile
channels presented similar responses about
pressure gradient throughout the geometry, with
no stagnation zones.

Figure 8: Pressure lines of 2mm square profile.

Figure 6: Pressure lines of triangular profile.

Figure 9: Pressure lines of 1.5mm square profile.

In the Figures 10-13 the z-velocity slice

profiles for all channels configurations are
showed.
At velocity slice graphic (Figure 12) of 2mm
square profile, there is a remarkable region of
low velocity in z, confirming the obtained results
by pressure lines.
 Excerpt from the Proceedings of the COMSOL Conference 2007, Boston

Although both triangular and trapezoidal

profiles presented pressure gradient (Δp ≈ 0.55
kPa at GDL) below of the 1.5mm square profile
(Δp ≈ 0.65 kPa at GDL), the z-direction velocity
profile stood around 0.2 m.s-1 throughout x-axe;
the square profiles, instead, presented areas or
spots of large velocity decreasing (hitting 0.1
m.s-1 at diffusion layer outlet, z=2.25mm at 2mm
square profile and z=1.75mm at 1.5 square
profile).
In order to check the studied models
correspondence within fuel cell experiments,
both square profiles (graphite plates) were used.
The polarization curve measures the fuel cell
performance in terms of potential, under Figure 12: z-velocity slice of 2mm square
profile.
different values of current. The value of 0.6V
(work potential), is set to compare two different
fuel cells performance.

Figure 13: z-velocity slice of 1.5mm square

profile.

Figure 10: z-velocity slice of triangular profile.

At Figure 14 the polarization curves of
graphite plates (prototype 1 and 2), were
presented. The values for work operation were:
166.7mA cm-2 and 470.6mA cm-2, for the cases 1
and 2, respectively.
The graphite plate B had better performance
than the A one, that was validated by numeric
methods. Such behavior was attributed to the
higher channel area in contact with gas diffusion
electrode which provides higher active area
available for the gas reaction. Besides the area,
higher gas velocities in flow fields decrease the
gas stagnation points along the channels.

Figure 11: z-velocity slice of trapezoidal profile.

 Excerpt from the Proceedings of the COMSOL Conference 2007, Boston

1000
900
800
700
Potential [mV]

600
500
400
300
Potential Prot 1
200
Potential Prot 2
100
0
0 100 200 300 400 500 600 700 800
-2
Current Density [mA cm ]

Figure 14: Polarization curves of prototype 1 and prototype 2.

5. Conclusions Challenges and opportunities, Energy (2006),

doi:10.1016/j.energy.2006.08.007.
The CFD techniques help the fuel cell’s 2. E.F. Cunha, A.B. Andrade, E. Robalinho,
reserches to explore the connection between the M.L.M. Bejarano, E. Cekinski, M. Linardi,
transport of reactants and products and overall Modeling and simulation of PEM fuel cell’s flow
cell performance. CFD simulates hydrogen gas channels using CFD techniques, International
flow channels to reduce the costs of bipolar Nuclear Atlantic Conference – INAC, 2007.
plates’ production through optimization of the 3. E. Middelman, W. Kout, B. Vogleaar, J.
mass transport. New geometries are proposed Lenssen, E. de Waal, Bipolar plates for PEM fuel
and tested at laboratory. cells, Journal of Power Sources, Volume 118,
Several models’ simulations of diffusion in pp.44-46 (2003).
porous media and conduction/convection 4. A. Kummar, R.G. Reddy, Materials and design
processes were studied. Then, the experimental development for bipolar/end plates in fuel cell,
data were used to visualize the differences in Journal of Power Sources, Volume 129, pp.62-
performance between PEMFC A and B. The 67 (2004).
graphite plate B exhibited better performance 5. W. Quian, D. P. Wilkson, J. Shen, H. Wang, J.
than the A one, that was validated by numeric Zhang, Architecture for portable direct liquid
methods. Such behavior was attributed to the fuel cells, Journal of Power Sources, Volume
higher channel area in contact with gas diffusion 154, pp.202-213 (2006).
electrode which provides higher active area 6. P.T. Nguyen, T. Berning, N. Djilali,
available for the gas reaction. Besides the area, Computational model of PEM fuel cell with
higher gas velocities in flow fields decrease the serpentine gas flow channels, Journal of Power
gas stagnation points along the channels. Sources, Volume 130, pp.149-157 (2004).
7. S. Karvonen, T. Rottinen, J. Saarinen, O.
6. Acknowledgements Himanen, Modeling of flow field in polymer
electrolyte membrane fuel cell, Journal of Power
This work was supported by IPEN/CNEN-SP Sources, Volume 161, pp.876-884 (2006).
at PROCEL laboratories. 8. G. Squadrito, O. Barbera, I. Gatto, G.
Giacoppo, F. Urbani, E. Passalacqua, CFD
7. References analysis of the flow field scale-up influence on
the electrodes performance in a PEFC, Journal
1. N. Djilali, Computational modeling of of Power Sources, Volume 152, pp.67-74 (2005).
polymer electrolyte membrane (PEM) fuel cells: 9. Comsol Multiphysics Advanced Features
Course, Version: December 2005, COMSOL 3.2.