Excerpt from the Proceedings of the COMSOL Conference 2008 Hannover

Modeling Polybenzimidazole/Phosphoric Acid Membrane Behaviour in

a HTPEM Fuel Cell
C. Siegel*1,2, G. Bandlamudi1,2 and A. Heinzel1,2
1
Zentrum für BrennstoffzellenTechnik (ZBT) gGmbH, Duisburg, Germany
2
University of Duisburg-Essen, Institut für Energie- und Umweltverfahrenstechnik, Germany
*Corresponding author: Zentrum für BrennstoffzellenTechnik (ZBT) gGmbH, Carl-Benz-Straße 201, D-47057
Duisburg, Germany, [email protected]

Abstract: Phosphoric acid (H3PO4) doped in order to achieve the requested high proton
polybenzimidazole (PBI) membranes are conductivity [1].
commonly used in today’s high-temperature
polymer-electrolyte-membrane (HTPEM) fuel Much modeling and simulation work has been
cell technology. COMSOL Multiphysics is used published for low-temperature PEM fuel cells
to model and simulate the three-dimensional, but only a few computational fluid dynamics
single-phase, non-isothermal overall cell studies are currently available in literature
behaviour (e.g., species mass fraction covering HTPEM fuel cells. Cheddie and
distribution, solid-phase and fluid-(gas)-phase Munroe [1-3] published several studies
temperature distribution) at different operating accounting for various operating conditions and
points. The sol-gel PBI/H3PO4 membrane layouts. These models were solved using
behaviour is modeled using an Arrhenius COMSOL Multiphysics in one-, two, and three-
approach, accounting for the activation energy dimensions. Peng et al. [4,5] presented three-
and the membrane doping level. Simulation dimensional, steady-state and transient studies
results are compared to experimental using control volumes. Our group recently
investigations performed at the Zentrum für published a two-dimensional study with an
BrennstoffzellenTechnik (ZBT) gGmbH in emphasis on overall quantities behaviour [6].
Duisburg, Germany and a good agreement is The sol-gel PBI/H3PO4 membrane was modeled
found. In terms of numerical analysis, using an empirical relationship. Moreover, a
computational times, solution procedures and the two-equation energy conservation was
model’s overall convergence behaviour are introduced to specifically locate hot-spots.
briefly discussed. Experimental data was compared with simulation
results and a good agreement was found.
Keywords: High-temperature polymer-
electrolyte-membrane (HTPEM) fuel cell, Fuel In this study, the former developed two-
cell modeling, Polybenzimidazole (PBI), dimensional model is extended to the third
Phosphoric acid (H3PO4), Membrane behaviour. dimension to get a more realistic view of the
internal quantities and membrane behaviour at
1. Introduction different operating points. Cell optimization
potential is shortly discussed. Further, solution
Polymer-electrolyte-membrane (PEM) fuel procedures and convergence behaviour for
cells are expected to play a key role in different implemented solvers are highlighted.
forthcoming power-delivering devices. Beside
the well known advantages of the low 2. HTPEM model set-up
temperature PEM there are some major
drawbacks that are directly related to the 2.1 Computational domain
operating conditions (e.g., 80°C, high gas
humidification). Consequently, much attention Fig.1 shows the computational domain as
has been paid to high-temperature polymer- well as the used computational mesh. It consists
electrolyte-membrane (HTPEM) fuel cells. of two single gas channels to continuously feed
Phosphoric acid (H3PO4) doped the air and hydrogen in counterflow
polybenzimidazole (PBI) membranes have been configuration. The channels are machined into
investigated over the last decade and are our in-house, high-temperature stable, highly
commonly used in today’s HTPEM technology conductive (electrical and thermal) bipolar plates
(polyphenylene sulfide (PPS) based). The total gases through the membrane, viii) all material
thickness of a pristine, uncompressed high parameters are isotropic and homogeneous, ix)
temperature membrane electrode assembly the PBI/H3PO4 membrane is considered as a
(MEA) is expected to be 900-1000·10-6 [m]. In system of PBI, H3PO4 and H2O only (other forms
our model, the 5-layer MEA consists of two gas of acid are neglected), x) the initial amorphous
diffusion layers (GDL – light grey – ±320·10-6 phase H3PO4 concentration is 80-85 [wt.%], xi)
[m]) to equally distribute the incoming gases no interaction between the water and the
over the entire reaction layer surface. Within amorphous phase H3PO4, and xii) no water
both reaction layers (RL – black – ±50·10-6 [m]), transfer through the membrane.
the well known electrochemical half-cell
reactions take place. 2.3 Governing equations - Subdomains

2 ⋅ H 2 → 4 ⋅ H + + 4 ⋅ e− (1) COMSOL Multiphysics is used to solve this

complex HTPEM fuel cell model. The
O2 + 4 ⋅ H + + 4 ⋅ e − → 2 ⋅ H 2Ov , l (2) conservation equations are solved sequentially
(and grouped) for the variables velocity vector u
The sol-gel PBI/H3PO4 membrane (orange – [m·s-1], pressure P [Pa], species mass/mole
±100·10-6 [m]) separates the anode from the fraction ωi [-], electrical φ s [V] and protonic
cathode side and transports protons via the phase potential φm [V], and solid Ts [K] and
Grotthuss proton switching mechanism [9].
fluid-(gas)-phase temperature Tf [K]. All
conservation equations are strictly used in their
O2, N2 inlet conservative form and are very similar to the
(Gas channel)
ones given in [6].
5-layer MEA
The incompressible Navier-Stokes application
mode is used to describe the gas flow through
the channels (3)

∇⋅u = 0
ρ ⋅ u ⋅ ∇u =
( (
∇ ⋅ − P ⋅ I + η ⋅ ∇u + (∇u )T )) (3)

+F

and through the porous media (Brinkman

equations) (4).
PBI/H3PO4 H2 inlet z
membrane y x
(Gas channel)
∇ ⋅u = 0
Figure 1. Left: Computational domain including gas η
channels, bipolar-plates, and the 5-layer MEA; Right: ⋅u =
k p ,i
Structured mesh (4)
⎛ 1
(
∇ ⋅ ⎜⎜ − P ⋅ I + ⋅η ⋅ ∇u + (∇u )T
ε
)⎞⎟⎟
2.2 Assumptions ⎝ i ⎠
+F
The following assumptions are used in this
model: i) steady-state operating conditions, ii) In (3) and (4), ρ is the density [kg·m-3], η the
single-phase water flow since the fuel cell dynamic viscosity [Pa·s-1], kp,i [m2] the
operates at higher temperatures, iii) laminar gas permeability, and εi [-] the porosity of the porous
flow within all gas channels, iv) fully developed media.
gas inlet conditions, v) all gases are treated as
ideal, vi) all electrochemical reactions are The governing equations for species (Maxwell-
gaseous phase reactions, vii) no crossover of Stefan diffusion and convection application
mode), charge (conductive media DC application conductivity can accurately be described using
mode), and energy conservation (conduction and an Arrhenius approach (6) [8,9].
convection respectively conduction application
mode) are the same as in [6], extended to the ⎛ ΔE (k i , X ) ⎞
third dimension and are not repeated herein. σ 0 (ki , X ) ⎜−
⎜ R ⋅Ts ⎟⎠
⎟
(6)
σ= ⋅ e⎝
Moreover, all source term couplings are Ts
identically to the ones in [6].
The pre-exponential term σ0 [S·K·m-1] is
2.4 PBI/H3PO4 membrane considerations assumed to be independent from cell operating
temperature and decreases with higher doping
A highly doped sol-gel PBI/H3PO4 levels X [7]. In (7), C [mol·m-3] is the
membrane of approximately 30-35 mol of H3PO4 concentration of the mobile species. Among
per PBI repeat unit finally consists of more than
others, the concentration changes with X and σ0
±85 [wt.%] H3PO4 [7]. Since two H3PO4
is expected to reach the value of concentrated
molecules are bonded to PBI (crystalline region), H3PO4 [7].
X-2 molecules remain free and tend to form an
amorphous phase within the sol-gel PBI/H3PO4 ΔS + ΔS f
matrix. ⎛ z2 ⋅ F 2 ⎞
σ 0 = ⎜⎜ ⎟ ⋅ α ⋅ υ0 ⋅ d 2 ⋅ C ⋅ e
⎟
R (7)
⎝ R ⎠
In this model, the amorphous phase volume
fraction is calculated based on densities [3].
The activation energy ΔE [J·mol-1] is also
−1 expressed in terms of X and is the sum of the
⎛ 4.81 ⎞ enthalpies ΔH+ΔHf. Only minor valuable data
VH 3 PO4 = ⎜ + 1⎟ (5)
⎝ X −2 ⎠ are available in open literature for sol-gel
membranes [9-12]. Consequently, both terms (σ0
Fig.2 represents VH 3 PO4 [Vol.%] as a function of and ΔE) are empirically correlated to the values
the membrane doping level X [-], showing a of concentrated H3PO4 (85 [wt.%]) (see Fig.3 –
saturation behaviour at very high doping levels. VH 3 PO4 = 1) and described using the following
equation.
to the amorphous phase volume fraction [-]
Amorphous phase volume fraction [Vol.%]

1 0.25
σ 0 , ΔE = k1 ⋅ ln( X ) + k 2
Additional contribution: Doping level

0.9 (8)
0.8 0.2
0.7 40
0.6 0.15 1E+05
Pre-exponential factor σ0 [S·K·cm-1]

Sol-gel PBI/H3PO4
Activation energy ΔE [kJ·mol-1]

0.5
membrane
0.4 0.1 30
Sol-gel PBI/H3PO4 9E+04
0.3
membrane
0.2 0.05
20
0.1 6E+04
0 0
0 5 10 15 20 25 30 35 40 45 50 55 60 65 70
Doping level X [-] 10 3E+04

Figure 2. Amorphous phase volume fraction for

0 0
different membrane doping levels; Additional 0.6 0.7 0.8 0.9 1
contribution: Doping level to the amorphous phase Amorphous phase volume fraction [Vol.%]
volume fraction
Figure 3. Empirically correlated values for activation
It is widely accepted that the free (amorphous energy; Pre-exponential factor for different membrane
phase) H3PO4 mainly contributes to high doping levels
membrane conductivity via a Grotthuss proton
switching mechanism [8,9]. Moreover, the strong
temperature dependency of the membrane
2.5 Porous media considerations and physical- scripting while using the computed results as
chemical properties initial guesses within the next solution step. The
charge conservation was solved first. In a second
This model considers a woven-type GDL step, the mass/momentum conservation was
(e.g., E-tek - ELAT®) with an assumed porosity solved. Finally, the energy conservation was
of 0.7 [-] (uncompressed). The RL consists of i) solved while iteratively updating previous
a GDL weight fraction ii) a carbon results.
black/platinum weight fraction of 30 [wt.%]
(e.g., Vulcan XC-72), and iii) a PBI/H3PO4 Depending on the cell voltage, calculations were
weight fraction of 30 [wt.%]. The anode RL performed using the PARDISO direct solver
platinum loading is supposed to be 0.01 [kg·m-2] (±8,400 seconds for base case operating
and the cathode RL loading 0.0075 [kg·m-2]. The conditions) or the segregated group solver (group
calculated weight fractions are used to correlate 1: P, u; group 2: ωi; group 3: φ s , φm , Ts, Tf)
e.g., the electrical/protonic conductivities and the (±4,850 seconds for base case operating
diffusion coefficients of porous media using a conditions) on an 8GB RAM, Quad-Core
Bruggemann correlation. machine (64-Bit). The total number of degrees of
freedom was roughly 400,000. The convergence
Similar to [6], the reaction layer permeability is criterion was set to 1·10-3 – 1·10-6.
calculated taking a packed bed of spherical
particles into account. Further, all material and 1E+04
Charge balance Mass/Momentum/Species
gas properties are carefully declared with respect Balance (update)
to the higher operating temperatures.
(model subsystem / group) ErrEst [-]

1E+02 Energy balance

(update)
Relative error estimate

2.6 Governing equations – Boundary

conditions 1E+00
0 20 40 60 80 100 120

All boundary and initial conditions are carefully 1E-02

defined according to the experimental set-up.
Gas composition, temperature (fluid-(gas)-
phase), and velocity are declared at the gas 1E-04

channel inlets. At the outlets, pressure and PARDISO direct solver

Segregated group solver
convective flux boundary conditions are used. 1E-06
# of iteration steps [-]
The cell operating temperature (solid-phase) and
cell voltage are defined at the bipolar plate
Figure 4. Model subsystem / group convergence
boundaries. At internal boundaries, symmetry,
behaviour over the number of iteration for base case
continuity or thermal insulation boundary operating conditions at Ucell = 0.6 [V]: PARDISO
conditions are used. solver – black, Segregated group solver – grey)

3. Computational techniques 4. Simulation results

The dimensions of the geometry were 0.002 Fig.5 shows the oxygen mass fraction at
[m] x 0.004915 [m] x 0.01 [m]. A structured different internal subdomains boundaries within
mesh was used for all computations and the cell. The cell voltage is 0.6 [V]. The mass
consisted of 8,820 hexahedral elements. The fraction is higher under the gas channel than
minimum element quality was 0.0178 (element under the both current collectors and is
volume ratio 0.06). The gas channel was continuously consumed in z-direction due to
discretization using 6 x 6 x 15 elements. The electrochemical reactions.
GDL, RL, and PBI/H3PO4 membrane
subdomains were discretized using 5 elements in
y-direction (see Fig.1).

This highly coupled system of partial differential

equations was solved sequentially using solver
 0.23
and z-direction. It is proposed that hydrogen
easily reaches electrochemical active sites.

GDL to bipolar plate contact area

0.22
The average cathode side RL current density is
0.21 shown in Fig.7. The distribution is mainly
dependent on the species mass fractions and the
0.2 local overpotentials. Consequently highest
values are located under the bipolar plate,
0.19
possibly due to ohmic losses.
0.18
Fig.8 shows the distribution of the solid- and
fluid-(gas)-phase temperatures. In [6] it has
already been noticed that a severe influence
Figure 5. Oxygen mass fraction [-] at the internal between both variables exists, especially when
boundaries PBI/H3PO4 – RL (left); RL – GDL entering the cell at low gas temperatures. The
(middle); GDL – bipolar plate (right) solid-phase heat flux vectors indicate the heating
of the cell, keeping it at 160°C. The fluid-(gas)-
phase heat flux vectors show that heat is slightly
absorbed form the solid-phase (e.g., GDL solid
0.95
structure) by the gas streams. Consequently, the
GDL to bipolar plate contact area

0.9
temperature of both gases slightly rises along the
channel length. In this model set-up, Ts equals Tf
0.85 within the RL.

0.8 Cathode side

420
Ts = Tf within anode and cathode RL

0.75
400
0.7
380

Figure 6. Hydrogen mass fraction [-] at the internal 360

boundaries PBI/H3PO4 – RL (left); RL – GDL 340

(middle); GDL – bipolar plate (right)
320
3,800
300
3,700 Anode side
3,600

3,500
Figure 8. Solid-phase temperature and heat flux
vectors (left); Fluid-(gas)-phase temperature and heat
3,400 flux vectors (right) (Ucell = 0.6 [V]) (Slice plot at
3,300 z/zmax = 0.25)
3,200
The sol-gel PBI/H3PO4 membrane conductivity
3,100
distribution is given in Fig.9. The cold gases
3,000 entering the cell may cool the cell at local spots,
2,900 leading to lower proton conductivity values
implying higher ohmic (protonic) losses over the
membrane, again influencing the local current
Figure 7. Average cathode side RL current density density distribution. The average values are
[A·m-2] (Ucell = 0.6 [V]) slightly under 22 [S·m-1] (Ts = 160°C) and 25
[S·m-1] (Ts = 170°C) which is comparable to the
The anode side H2 mass fraction distribution is values given in [7]. Comparing these values with
much smoother, showing small gradients in x-, the conductivity values of H3PO4 [10] for similar
conditions (160°C – 80 [wt.%] → ±65 [S·m-1] 1.35 [-] (H2) and 2.5 [-] (air) (ambient pressure at
and 170°C – 80 [wt.%] → ±69 [S·m-1])) returns a both outlets).
ratio of ±0.3 [-], indicating that the transport
process through the sol-gel PBI/H3PO4
membrane is possibly not as continuous as in
H3PO4, even at high doping levels.

25

24

23

22

21

20

19
Hydrogen flow

18
Air flow

17
Figure 10. Left: Specially designed HTPEM fuel cell;
Right: Cell during operation at ZBT
16

Fig.11 compares the simulated at measured IV-

Figure 9. Sol-gel PBI/H3PO4 membrane conductivity
and power density curves. It should be noted that
distribution for 160°C (left) and 170°C (right) (Ucell =
0.6 [V])
the exchange current density i0 [A·m-2] is used as
a tuning parameter within a reasonable range
The voltage losses over the membrane are given in [12].
calculated using the membrane thickness lm [m] 1 16
and the following equation
14

im 0.75 12
U Ω, mem = ⋅ lm (15)
Cell voltage Ucell [V]

σm

Cell power P [W]

10

0.5 8
The IV-curve in Fig.12 is mainly dominated by
6
activation and ohmic losses. For 160°C, the
calculated ohmic voltage losses returned ±0.04 0.25 4

[V] at Ucell = 0.6 [V]. 2

0 0
5. Experimental investigations 0 2000 4000 6000
Current density J [A/m2]

The Zentrum für BrennstoffzellenTechnik Figure 11. Measured (black) and simulated (grey) IV-
(ZBT) gGmbH in Duisburg, Germany, is and power curves for a HTPEM fuel cell (Ts = 160°C,
working for several years in the field of HTPEM Tf = 21°C)
fuel cell research and development.
Consequently, all tests were run at their facilities 6. Conclusions
under precisely defined operating conditions (see
Fig.10). This work presented a complete three-
dimensional HTPEM fuel cell models using a
The base case cell operating temperature was sol-gel PBI/H3PO4 membrane. Most material
controlled at 160°C using a heated aluminum parameters were based on the latest publications
end-plate. Gases (H2 – ±0.25% relative humidity, in the field and from data sheets. It was shown,
air – ±4% relative humidity) were feed at room that the cathode side flow-field layout is one
temperature using stoichiometric flow rates of major optimization parameter in order to
increase and equally distribute the oxygen near
the RL. Further, the gas-(fluid)-phase 10. Chin, D.-T., Chang, H.H., On the
temperature should be considered as an conductivity of phosphoric acid electrolyte, J.
important operating parameter for cell Appl. Electrochem., 19, 95-99 (1989)
tempering. Overall a good agreement was found 11. Bouchet, R., Siebert, E., Proton conduction
between the measured and simulated results. in acid doped polybenzimidazole, Solid State
Many more investigations are necessary to Ionics, 118, 287-299 (1999)
precisely validate the model set-up. Including 12. Zhang, J., Tang, Y., Song, C., Zhang, J.,
gas solubility at electrochemical active sites is a Polybenzimidazole-membrane-based PEM fuel
task for future work. The model will more cell in the temperature range of 120-200°C, J.
accurately describe the mass limitation effects by Power Sources, 172, 163-171 (2007)
taking the amount of H3PO4 within the RL into
account. 8. Acknowledgements

7. References This work was supported by ‘LE

GOUVERNEMENT DU GRAND-DUCHÉ DE
1. Cheddie, D., Munroe, N., Parametric model of LUXEMBOURG, MCESR – Recherche et
an intermediate temperature PEMFC, J. Power Innovation’, Grant No.: BFR07/007.
Sources, 156, 414-423 (2006)
2. Cheddie, D., Munroe, N., Three dimensional
modeling of high temperature PEM fuel cells, J.
Power Sources, 160, 215-223 (2006)
3. Cheddie, D., Munroe, N., A two-phase model
of an intermediate temperature PEM fuel cell,
Int. J. Hydrogen Energy, 32, 832-841 (2007)
4. Peng, J., Lee, S.J., Numerical simulation of
proton exchange membrane fuel cells at high
operating temperature, J. Power Sources, 162,
1182-1191 (2006)
5. Peng, J., Shin, J.Y., Song, T.W., Transient
response of high temperature PEM fuel cell, J.
Power Sources, 179, 220-231 (2008)
6. Siegel, C., Bandlamudi G., Heinzel, A.,
Proceedings of the European COMSOL
Conference – Vol. 1 Numerical simulation of a
high-temperature PEM (HTPEM) fuel cell, 428-
434, Petit, J.-M., Squalli, O., Grenoble, France
(2007)
7. Xiao, L., Zhang, H., Scanlon, E., Ramanathan,
S., Choe, E.-W., Rogers, D., Apple, T.,
Benicewicz, B.C., High-temperature
polybenzimidazole fuel cell membranes via a
sol-gel process, Chem. Mater., 17, 5328-5333,
(2005)
8. Zhang, J., Xie, Z., Zhang, J., Tang, Y., Song,
C., Navessin, T., Shi, Z., Song, D., Wang, H.,
Wilkinson, D.P., Liu, Z.-S., Holdcroft, S., High
temperature PEM fuel cells, J. Power Sources,
160, 872-891 (2006)
9. Ma, Y.-L., The fundamental studies of
polybenzimidazole/phosphoric acid polymer
electrolyte for fuel cells, Dissertation, Case
Western Reserve University, U.S.A. (2004)