# A Mathematical Model of a Direct Propane Fuel Cell

Authors: Hamidreza Khakdaman; Yves Bourgault; Marten Ternan
Journal: Journal of Chemistry
 (2015)

---

## Abstract

A rigorous mathematical model for direct propane fuel cells (DPFCs) was developed. Compared to previous models, it provides better values for the current density and the propane concentration at the exit from the anode. This is the first DPFC model to correctly account for proton transport based on the combination of the chemical potential gradient and the electrical potential gradient. The force per unit charge from the chemical potential gradient (concentration gradient) that pushes protons from the anode to the cathode is greater than that from the electrical potential gradient that pushes them in the opposite direction. By including the chemical potential gradient, we learn that the proton concentration gradient is really much different than that predicted using the previous models that neglected the chemical potential gradient. Also inclusion of the chemical potential gradient made this model the first one having an overpotential gradient (calculated from the electrical potential gradient) with the correct slope. That is important because the overpotential is exponentially related to the reaction rate (current density). The model described here provides a relationship between the conditions inside the fuel cell (proton concentration, overpotential) and its performance as measured externally by current density and propane concentration.


---

## Body

## 1. Introduction

The focus of this study is on the direct propane fuel cell (DPFC), which belongs to the polymer membrane electrolyte fuel cell (PEMFC) family but consumes propane instead of hydrogen as its feedstock. The generation of electrical energy in rural areas is our primary target application for DPFCs. The cost of delivering electrical energy to rural areas is substantially greater than to urban areas, because longer transmission lines are required to serve a comparatively small number of customers. Therefore more costly fuel cells can be justified for use in rural areas compared to urban areas. In addition, the infrastructure for delivering liquefied petroleum gas (LPG) or propane to rural areas already exists. Two major advantages of DPFCs over hydrogen fuel cells are that the expense of hydrogen production plants and of hydrogen transport/storage will be eliminated from the fuel cell energy production cycle. However, a drawback associated with DPFCs is that the propane reaction rate is much slower than that of hydrogen. Liebhafsky and Cairns [1], Bockris and Srinivasan [2], and Cairns [3] reviewed the majority of the DPFC experimental research that had been done in the 1960s. Only a little work has been performed since then. Polarization curves were reported both by Cheng et al. [4] and by Savadogo and Rodriguez Varela [5] for low temperature PEMFCs. Intermediate temperature (100–300°C) proton conducting fuel cells were investigated by Heo et al. [6]. Solid oxide fuel cell (SOFC) studies were performed with propane at 550–650°C by Feng et al. [7] and at 900°C by Yang et al. [8].The approaches to the modeling of fuel cells are summarized below. Weber and Newman [9] have reviewed four groups of fuel cell models that consider transport of water and protons in the electrolyte phase: simple models, diffusive models, hydraulic models, and combination models. The simple models [10–13] describe proton transfer using Ohm’s law with a constant ionic conductivity. These models cannot predict phenomena such as membrane dehydration, in which water content and thus ionic conductivity are variables. For water movement, a numerical value of the net water flux has to be determined as the boundary condition at the interface between the catalyst layers and the membrane.The diffusive models [14–19] predict the movement of dissolved water and protons within the membrane as a result of concentration and electrical potential gradients. They are applicable for the electrolyte systems with low water content (λ
<
14, where λ is moles of water per mole of sulfonic acid sites in the Nafion membrane) where liquid water does not exist. The diffusive models are referred to as single phase models of membranes and can predict proton distribution in the electrolyte phase and membrane dehydration.At high water contents, membrane pores are completely filled with liquid water, and the water content is assumed to be uniform everywhere. Therefore, water diffusion does not occur and the convection mechanism causes proton and water transport. The hydraulic models [20–23] were developed for membranes with high water content. Two phases, liquid water and membrane, are described by the hydraulic models. Water velocity is calculated by Schlögl’s equation [23] which is a function of electrical potential gradient and pressure gradient. Finally, the hydraulic and diffusive models are merged in the combination models [24–26] when calculations covering the whole range of water content are desirable. This approach considers concentration and pressure gradients as driving forces for water and proton transport.There are two possible approaches to dealing with the transport properties in the diffusive models, that is, dilute solution theory and concentrated solution theory [27]. Mass transport in dilute electrolyte systems is usually described by the Nernst-Planck equation [27] in which the flux of a charged species is a function of the concentration gradient of that species as well as the electrical potential gradient. For a noncharged species, the potential gradient term in the Nernst-Planck equation disappears. The membrane transport properties are not required to be constant in this approach.Employing concentrated solution theory leads to rigorous models that consider the interactions between all species. Krishna [19] used Generalized Maxwell-Stefan (GMS) equations to implement this approach for multicomponent electrolyte systems in general. Wöhr et al. [17] also used Maxwell-Stefan (MS) equations to model proton and water transport in PEM fuel cells in which the MS diffusion coefficients are modified as a function of temperature and humidity. Fuller and Newman [28] used the electrochemical potential of each species as the driving force in the MS equations. Fimrite et al. [16] developed a transport model for water and protons based on the binary friction model. The mole fraction and potential gradients were considered in the electrochemical potential gradient expression. Baschuk and Li [15] also used MS equation but they calculated the MS diffusion coefficients based on experimental data available in the literature. Then, they validated those coefficients with experimental data for the electroosmotic drag coefficient.A diffusive model has been developed in the present study to investigate the movement of water and protons in the electrolyte phase of a DPFC where the operation temperature is above the boiling point of water. One possible strategy for increasing the reaction rate in DPFC is to operate at temperatures of 150°C or higher. A membrane that can resist high temperature and show acceptable conductivity (5.0 S m−1) has been developed in our research group [29]. This membrane is composed of porous polytetrafluoroethylene (PTFE) that contains zirconium phosphate (Z
r
H
P
O
4
2
·
H
2
O or ZrP) in its pores. ZrP is a known proton conductor [30]. Concentrated solution theory was used in which the binary interactions between water, protons, and ZrP species were described.We are developing mathematical models of DPFCs in order to understand this phenomenon and hopefully to enhance their performance. The results reported here are major improvements over our previous model [10]. Our previous model, like the vast majority of fuel cell models, only used an electrical potential gradient to describe migration and neglected the proton concentration gradient in accounting for proton transport through the electrolyte layer. As we noted previously [10], neglect of the proton concentration gradient caused the overpotential gradient and the electrical potential gradient in the electrolyte phase to be incorrect. The model being described here, unlike the majority of fuel cell models, includes both a valid electrical potential gradient and a proton concentration gradient to account for proton transport by a combination of migration and diffusion. This model accounts for the influence of the proton concentration in the electrolyte phase and thereby overcomes the deficiencies mentioned above.


## 2. Model Development

This model solved the governing equations for the Membrane Electrode Assembly (MEA), consisting of the membrane layer, anode layer, and cathode layer. A schematic of a typical DPFC is shown in Figure1. The cell is composed of two bipolar plates, two catalyst layers, and a membrane layer. Each bipolar plate has two sets of channels: one for reactants and one for products. The channels are connected to each other through the catalyst layer. Figure 1 shows these channels for the anode bipolar plate. The interdigitated flow fields show a symmetric geometry with repetitive pieces. In order to increase the computational speed, only one of these pieces was considered as the modeling domain. Therefore, the modeling domain can be defined as the part of the MEA that is located between the middle of a feed channel and the middle of its adjacent product channel (cross section in Figure 1). That cross section is shown in Figure 2 as the modeling domain. Its boundaries are shown as a dashed black line.Figure 1
A direct propane fuel cell with interdigitated flow field.Figure 2
Boundaries in the modeling domain.Previously, it was shown that neglecting proton diffusion in the proton conservation equation (an assumption used in many fuel cell models) led to incorrect results for the electrolyte potential and overpotential profiles even though the polarization curve was predicted correctly [10]. The present model includes both proton diffusion and migration.


### 2.1. Governing Equations

Three phases are present in the anode and cathode catalyst layers. They are the “gas phase” containing reactants and products, the “solid catalyst phase” containing the carbon support and platinum, and the “solid electrolyte phase.” The latter consists of a stationary ZrPm
a
t
r
i
x
=
[
Z
r
H
P
O
4
2
·
H
2
O
] containing mobile H2O = [Zr(HPO4)2  ·  2H2O] and mobile H
+
=
[
Z
r
H
P
O
4
2
·
H
3
O
+
] species that can be transported. The membrane layer contains the ZrP electrolyte phase as well as PTFE.Conservation equations for momentum, total mass, and mass of noncharged species were solved for the gas phase in each of the catalyst layers. A list of equations that were used for the gas phase of both anode and cathode catalyst layers is shown as follows:Conservation of mass in gas phase:(1)
∇
·
ε
G
ρ
G
u
→
+
∑
i
n
ν
i
M
W
i
j
z
F
=
0
,
wherei = C3H8, H2O, and CO2 for the anode and O2 and H2O for the cathode.
Conservation of momentum in gas phase:(2)
-
∇
P
=
150
μ
G
1
-
ε
G
2
D
p
2
ε
G
3
u
→
.
Conservation of noncharged species in gas phase:(3)
∇
·
ε
G
c
G
u
→
y
i
-
∇
·
ε
G
c
G
D
i
∇
y
i
+
ν
i
j
z
F
=
0
,
wherei = C3H8 and CO2 for the anode and O2 and H2O for the cathode.
Conservation of species in the electrolyte phase:
for water,(4)
-
∇
·
c
E
L
Y
B
H
2
O
–
H
2
O
′
-
B
H
2
O
–
H
+
′
∇
x
H
+
+
∇
·
c
E
L
Y
B
H
2
O
–
H
+
′
F
x
H
+
R
T
∇
ϕ
E
L
Y
-
j
z
F
=
0
;
for proton,(5)
∇
·
c
E
L
Y
B
H
+
–
H
+
′
-
B
H
+
–
H
2
O
′
∇
x
H
+
+
∇
·
c
E
L
Y
B
H
+
–
H
+
′
F
x
H
+
R
T
∇
ϕ
E
L
Y
+
j
z
F
=
0
.
Butler-Volmer equation in the anode:(6)
j
A
=
j
A
0
A
P
t
exp
⁡
α
A
F
η
A
R
T
-
exp
⁡
-
α
C
F
η
A
R
T
,
where(7)
j
A
0
=
j
C
3
O
x
0
r
e
f
p
C
3
p
C
3
r
e
f
exp
⁡
Δ
G
C
3
O
x
‡
R
1
T
r
e
f
-
1
T
,
(8)
η
A
=
Δ
ϕ
A
-
Δ
ϕ
A
E
Q
=
ϕ
P
t
A
-
ϕ
E
L
Y
A
-
ϕ
P
t
A
E
Q
-
ϕ
E
L
Y
E
Q
.
Butler-Volmer equation in the cathode:(9)
j
C
=
j
C
0
A
P
t
exp
⁡
α
A
F
η
C
R
T
-
exp
⁡
-
α
C
F
η
C
R
T
,
where(10)
j
C
0
=
j
O
2
R
d
0
r
e
f
p
O
2
p
O
2
r
e
f
exp
⁡
Δ
G
O
2
R
d
‡
R
1
T
r
e
f
-
1
T
,
(11)
η
C
=
Δ
ϕ
C
-
Δ
ϕ
C
E
Q
=
ϕ
P
t
C
-
ϕ
E
L
Y
C
-
ϕ
P
t
C
E
Q
-
ϕ
E
L
Y
E
Q
.Equation (1) describes the total mass conservation in the gas phase of the catalyst layers. The second term in this equation is the sink or source term describing the mass consumption or production in the gas phase caused by electrochemical reactions. Equation (2) is the linear form of the Ergun equation. It was used to calculate the pressure profiles in the gas phase of the catalyst layers because they are packed beds. At the conditions used in this study, the magnitude of the quadratic velocity term in the Ergun equation was much smaller than the linear term. Hence only the linear term in velocity was used in (2). Equations (1) and (2) were solved together to calculate the velocity and pressure profiles in the gas phase of the catalyst layers. Mass balances for each of the individual gas phase species account for convection, diffusion, and reaction, as shown in (3).Equations (4) and (5) describe, respectively, water and proton conservation in the electrolyte phase of the membrane and catalyst layers. Diffusion was described by concentrated solution theory through the use of the GMS equations. The following paragraphs illustrate the derivation of (4) and (5).A general procedure for the calculation of mass fluxes in multicomponent electrolyte systems was presented by Krishna [19]. It has been proven that the Nernst-Planck equation is a limiting case of the GMS equations. The GMS equations can be written as follows:(12)
d
→
i
=
∑
j
=
1
j
≠
i
n
x
i
J
→
j
-
x
j
J
→
i
c
E
L
Y
Ð
i
j
i
=
1,2
,
…
,
n
-
1
,where d
→
i is a generalized driving force for mass transport of species i. Because the summation of the n driving forces is equal to zero due to the Gibbs-Duhem limitation [31], only n
-
1 driving forces are independent. The equation to calculate the generalized driving force has been derived based on nonequilibrium thermodynamics [31]. A simplified expression for a solid stationary electrolyte (no convection term) [19] can be written as(13)
d
→
i
=
∇
x
i
+
x
i
z
i
F
R
T
∇
ϕ
E
L
Y
.For a noncharged species such as water z
i is equal to zero, and, according to (13), the concentration gradient will be the only driving force.The migration term in (13) was obtained by representing ion mobility by the Nernst-Einstein relation (D
i
=
R
T
u
i). This equation is applicable only at infinite dilution. However, it can be used in concentrated solutions if additional composition-dependent transport parameters, such as the B
′ parameters in (19), are used to calculate the flux of ions [27]. It will be shown in the following paragraphs that (18) represent the composition-dependent parameters.Equation (12) results in (n
-
1) independent equations that can be written in matrix form for convenience:(14)
c
E
L
Y
d
→
1
⋮
d
→
n
-
1
=
-
B
1
1
⋯
B
1
,
n
-
1
⋮
⋱
⋮
B
n
-
1
,
1
⋯
B
n
-
1
,
n
-
1
J
→
1
⋮
J
→
n
-
1
,where the elements of the matrix of inverted diffusion coefficients [B] are given by(15)
B
i
i
=
∑
j
=
1
j
≠
i
n
x
i
Ð
i
j
i
=
1,2
,
…
,
n
-
1
,
B
i
j
=
-
x
i
Ð
i
j
i
=
1,2
,
…
,
n
-
1
;
i
≠
j
.The fluxes of species, J
→
i, can be calculated from (16) which is the inversion of (14):(16)
J
→
1
⋮
J
→
n
-
1
=
-
c
E
L
Y
B
11
⋯
B
1
,
n
-
1
⋮
⋱
⋮
B
n
-
1,1
⋯
B
n
-
1
,
n
-
1
-
1
d
→
1
⋮
d
→
n
-
1
.For the present electrolyte system containing three species, mobile H2O and H+ plus immobile solid ZrP, (16) may be written as(17)
J
→
H
2
O
J
→
H
+
=
-
c
E
L
Y
B
H
2
O
–
H
2
O
′
B
H
2
O
–
H
+
′
B
H
+
–
H
2
O
′
B
H
+
–
H
+
′
d
→
H
2
O
d
→
H
+
,where [B
′] is the inverse of the matrix of inverted diffusion coefficients. Because Ð
H
2
O
–
H
+
=
Ð
H
+
–
H
2
O, the elements of [B
′] are calculated using (18) which are functions of the GMS diffusivities and the species mole fractions in the electrolyte phase:(18)
B
H
2
O
–
H
2
O
′
=
x
H
2
O
Ð
H
+
–ZrP
+
Ð
H
2
O
–
H
+
x
H
+
+
Ð
H
+
–ZrP
/
Ð
H
2
O
–ZrP
x
H
2
O
+
Ð
H
2
O
–
H
+
/
Ð
H
2
O
–ZrP
,
B
H
2
O
–
H
+
′
=
x
H
2
O
Ð
H
+
–ZrP
x
H
+
+
Ð
H
+
–ZrP
/
Ð
H
2
O
–ZrP
x
H
2
O
+
Ð
H
2
O
–
H
+
/
Ð
H
2
O
–ZrP
,
B
H
+
–
H
2
O
′
=
x
H
+
Ð
H
2
O
–ZrP
x
H
2
O
+
Ð
H
2
O
–ZrP
/
Ð
H
+
–ZrP
x
H
+
+
Ð
H
2
O
–
H
+
/
Ð
H
+
–ZrP
,
B
H
+
–
H
+
′
=
x
H
+
Ð
H
2
O
–ZrP
+
Ð
H
2
O
–
H
+
x
H
2
O
+
Ð
H
2
O
–ZrP
/
Ð
H
+
–ZrP
x
H
+
+
Ð
H
2
O
–
H
+
/
Ð
H
+
–ZrP
.Combining sets of (17) and (13) results in two independent equations that can be used to calculate the fluxes of mobile species (J
→
H
2
O and J
→
H
↑
+
↓
↑) within the electrolyte phase:(19)
J
→
H
2
O
=
-
c
E
L
Y
B
H
2
O
–
H
2
O
′
∇
x
H
2
O
-
c
E
L
Y
B
H
2
O
–
H
+
′
∇
x
H
+
+
F
x
H
+
R
T
∇
ϕ
E
L
Y
,
(20)
J
→
H
+
=
-
c
E
L
Y
B
H
+
–
H
2
O
′
∇
x
H
2
O
-
c
E
L
Y
B
H
+
–
H
+
′
∇
x
H
+
+
F
x
H
+
R
T
∇
ϕ
E
L
Y
.Equations (19) and (20) show that diffusion flux of each species is a function of the concentration gradient of all species as well as of the potential gradient. There are five unknowns in (19) and (20): J
→
H
2
O, J
→
H
+, x
H
2
O, x
H
+, and ϕ
E
L
Y. Therefore, three more equations are required.ZrP is immobile. As a result, the diffusion phenomenon will effectively be the interchange of H+ and H2O species. Therefore, for diffusion purposes we will only consider the domain of the mobile species, H+ and H2O, and will ignore the immobile species, ZrP. On that basis, (21) can be used as a third equation. Nevertheless, the presence of ZrP is important because of its interaction with the mobile species. Specifically, the values of the B
′ coefficients for H+ and H2O were influenced by the presence of ZrP:(21)
x
H
2
O
+
x
H
+
=
1.0
.The differential equations for H2O and H+ mass conservation in the electrolyte phase can be expressed in molar units as(22)
∇
·
J
→
H
2
O
=
-
j
z
F
,
∇
·
J
→
H
+
=
j
z
F
,where j is the volumetric current production. This quantity which appears in (1), (3) to (5) and (22) is the rate of production of protons in the anode. Therefore, it is positive in the anode, j
A, and negative in the cathode, j
C. It was calculated using the Butler-Volmer equation for the anode and cathode, (6) and (9), respectively. Exchange current densities at the anode and cathode are a function of the reactants’ partial pressure and the operating temperature as shown in (7) and (10). The Butler-Volmer equation and its parameters for both propane oxidation and oxygen reduction were described in our previous communication [10]. Complete conversion of C3H8 to CO2 was reported in experiments by Grubb and Michalske [34]. Equations (19) to (22) were combined and are shown as (4) and (5).


### 2.2. Numerical Procedure

The numerical solution procedure is illustrated in Figure3. Equations (1)–(11) define the problem at steady state. However, a time derivative was appended to each partial differential equation and a backward Euler time stepping method was used to increase stability while converging to the steady-state solution. The Finite Element Method was used to discretize the partial differential equations in space, with all dependent variables discretized by a linear finite element except for the pressure that is taken as a quadratic.Figure 3
Modeling procedure.FreeFEM++ software has been successfully used to solve two-dimensional partial differential equations (1)–(11). It is open-source software and is based on the Finite Element Method developed by Hecht et al. [32]. The calculated results from FreeFEM++  were exported to ParaView visualization software [35] for postprocessing. ParaView is also open-source software.There is no proton loss through the exterior boundaries of the domain (Figure2). Therefore, the total rate of proton production in the anode, ∫
A
n
o
d
e
j
d
V, has to be equal to the total rate of proton consumption in the cathode, ∫
C
a
t
h
o
d
e
(
-
j
)
d
V. In each case, the electrical potential of the catalyst phase of the anode, ϕ
P
t
A, and that of the cathode, ϕ
P
t
C, had individual constant values. Then all the variables in the whole domain were calculated. However, having fixed electrical potentials of the anode and cathode catalyst phases does not guarantee that the proton production at the anode will equal the proton consumption at the cathode. The difference between the rate of proton production and consumption can be minimized by shifting ϕ
E
L
Y by a constant value because the production and consumption rates are functions of the electrical potential in both of their respective catalyst phases, ϕ
P
t
A and ϕ
P
t
C, and in the electrolyte phase, ϕ
E
L
Y. Therefore, the Newton method was used to force equal proton production and consumption. In other words, balancing ∫
A
n
o
d
e
j
d
V and ∫
C
a
t
h
o
d
e
(
-
j
)
d
V acts as a constraint for the conservation of protons in the electrolyte phase.The equations for the conservation of momentum, total mass, and individual species in the gas phase of the anode and cathode were solved by assuming there was no species crossover through the membrane. Electrical potential, proton, and water concentrations in the electrolyte phase of the anode, cathode, and membrane layers were coupled to each other. These variables were calculated by solving (4), (5), and (21) iteratively in each layer. Then, the Robin method [10] was used to couple the solutions between layers. In the Robin method, both of the following transfer conditions are progressively satisfied on the anode catalyst/membrane interface and the membrane/cathode catalyst interface through iterations of (a) the continuity of the variable (e.g., potential) and (b) the continuity of the flux (e.g., electrical current).Figure2 shows four types of boundary conditions for the modeling domain, that is, inlet, outlet, wall of the land, and the midchannel symmetry boundaries. The flux of species in the gas phase is zero at the walls because there is no transfer through walls. The zero flux condition is also true at the midchannel symmetry boundaries. The compositions of the gaseous species are known at the inlet of the anode and cathode catalyst layers. It was assumed that no change in the composition of gas mixture occurred after leaving the catalyst bed. Therefore, the composition gradients are zero in the direction normal to the catalyst layer at the outlet boundaries. The zero flux condition is applied at all exterior boundaries for the species in the electrolyte phase.


### 2.3. Input Parameters

The parameters used for the simulations are shown in Table1. The GMS diffusivities, Ð
i
j, which are used in (18) have to be calculated from the Fickian diffusion coefficients, D
i
j. For ideal solutions, the Fickian diffusion, D
i
j, can be used as Ð
i
j in the Stefan-Maxwell equations [26] because the concentration dependence of Fickian diffusion coefficients is ignored. Experimental values for D
H
+
–ZrP and D
H
2
O
–
H
+ are given in Table 1. Note that the diffusivity of protons in ZrP is approximately two orders of magnitude smaller than the diffusivity of protons in water. The movement of protons causes the electroosmotic flow of water [9]. It was assumed that one water molecule is dragged by each proton, H3O+, that travels from anode to cathode. Therefore, the diffusivity of water in ZrP was set equal to the diffusivity of protons in ZrP [36], the smaller of the two proton diffusivities in Table 1. Proton diffusivity and proton mobility are different quantities. The three diffusivities in Table 1 were the ones used to calculate the B
′ parameters in (18).Table 1
Operational, electrochemical, and design parameters for simulations.
Property
Value
Temperature,T
423–503 K
Pressure,P
101.3 k Pa
Proton–ZrP diffusivity,D
H
+
–ZrP
3.1 × 10−12 m2 s−1 [29]
Proton–water diffusivity,D
H
2
O–
H
+
2.9 × 10−10 m2 s−1 [12]
Ionic conductivity in membrane,σ
ZrP
/
PTFE
5.0 S m−1 [24]
Electrical resistivity in membrane,R
PTFE
1.0 × 1016 Ω m
Charge transfer coefficients,α
A and α
C
1.0 [30]
Equilibrium potential of catalyst phase at the anode,ϕ
Pt
A
EQ
0.136 V [1]
Equilibrium potential of catalyst phase at the cathode,ϕ
Pt
C
EQ
1.229 V
Equilibrium potential of electrolyte phase,ϕ
ELY
EQ
0.136 V
Apparent bulk density of carbon catalyst support,ρ
CAT
⁡
0.259 gcatalyst
mL
catalyst
-
1
Specific surface area of carbon catalyst support in the anode and cathode,A
CAT
⁡
255m
catalyst
2 g
catalyst
-
1
Gas phase volume fraction in anode and cathode,ε
G
0.5
Electrolyte phase volume fraction in anode and cathode,ε
ELY
0.4
Effective particle diameter in anode and cathode,D
p
5μm
Land width,L
W
2–8 mm
Anode and cathode thickness, ThA, ThC
200–400μm
Membrane thickness,T
h
M
100–200μm
Fluid channels width in bipolar plates
0.4 mm


### 2.4. Model Validation

The model predicts the performance of a DPFC that (i) has interdigitated flow fields, (ii) has zirconium phosphate as the electrolyte, and (iii) operates over a temperature range of 150–230°C. As there are no experimental data for DPFCs having zirconium phosphate electrolytes and interdigitated flow fields, the model results have been compared to published results for DPFCs with other types of electrolytes and flow fields.Figure4 compares the modeling results for zirconium phosphate electrolyte with the experimental data for other types of electrolytes [34, 37]. The figure shows that the polarization curve for ZrP-PTFE electrolyte is somewhat comparable to that for the other electrolytes. The difference between the polarization curves can be partially explained by the difference between conductivities of the electrolytes. The proton conductivity of a nonmodified Nafion 117 approaches 10 S m−1 at 80°C [38]. The conductivity of the 95% H3PO4 electrolyte is 35 S m−1 at 200°C [39]. However, the proton conductivity for the best ZrP-PTFE that has been developed in our laboratory is about 5 S m−1 at 150°C.Figure 4
Polarization curves of direct propane/oxygen fuel cell using Pt anode and cathode. (a) Experimental results [31] using Nafion 117 at 95°C. (b) Experimental results [32] using 95% H3PO4 at 200°C. (c) The present proton migration and diffusion model results for a solid ZrP-PTFE electrolyte at 150°C.


## 2.1. Governing Equations

Three phases are present in the anode and cathode catalyst layers. They are the “gas phase” containing reactants and products, the “solid catalyst phase” containing the carbon support and platinum, and the “solid electrolyte phase.” The latter consists of a stationary ZrPm
a
t
r
i
x
=
[
Z
r
H
P
O
4
2
·
H
2
O
] containing mobile H2O = [Zr(HPO4)2  ·  2H2O] and mobile H
+
=
[
Z
r
H
P
O
4
2
·
H
3
O
+
] species that can be transported. The membrane layer contains the ZrP electrolyte phase as well as PTFE.Conservation equations for momentum, total mass, and mass of noncharged species were solved for the gas phase in each of the catalyst layers. A list of equations that were used for the gas phase of both anode and cathode catalyst layers is shown as follows:Conservation of mass in gas phase:(1)
∇
·
ε
G
ρ
G
u
→
+
∑
i
n
ν
i
M
W
i
j
z
F
=
0
,
wherei = C3H8, H2O, and CO2 for the anode and O2 and H2O for the cathode.
Conservation of momentum in gas phase:(2)
-
∇
P
=
150
μ
G
1
-
ε
G
2
D
p
2
ε
G
3
u
→
.
Conservation of noncharged species in gas phase:(3)
∇
·
ε
G
c
G
u
→
y
i
-
∇
·
ε
G
c
G
D
i
∇
y
i
+
ν
i
j
z
F
=
0
,
wherei = C3H8 and CO2 for the anode and O2 and H2O for the cathode.
Conservation of species in the electrolyte phase:
for water,(4)
-
∇
·
c
E
L
Y
B
H
2
O
–
H
2
O
′
-
B
H
2
O
–
H
+
′
∇
x
H
+
+
∇
·
c
E
L
Y
B
H
2
O
–
H
+
′
F
x
H
+
R
T
∇
ϕ
E
L
Y
-
j
z
F
=
0
;
for proton,(5)
∇
·
c
E
L
Y
B
H
+
–
H
+
′
-
B
H
+
–
H
2
O
′
∇
x
H
+
+
∇
·
c
E
L
Y
B
H
+
–
H
+
′
F
x
H
+
R
T
∇
ϕ
E
L
Y
+
j
z
F
=
0
.
Butler-Volmer equation in the anode:(6)
j
A
=
j
A
0
A
P
t
exp
⁡
α
A
F
η
A
R
T
-
exp
⁡
-
α
C
F
η
A
R
T
,
where(7)
j
A
0
=
j
C
3
O
x
0
r
e
f
p
C
3
p
C
3
r
e
f
exp
⁡
Δ
G
C
3
O
x
‡
R
1
T
r
e
f
-
1
T
,
(8)
η
A
=
Δ
ϕ
A
-
Δ
ϕ
A
E
Q
=
ϕ
P
t
A
-
ϕ
E
L
Y
A
-
ϕ
P
t
A
E
Q
-
ϕ
E
L
Y
E
Q
.
Butler-Volmer equation in the cathode:(9)
j
C
=
j
C
0
A
P
t
exp
⁡
α
A
F
η
C
R
T
-
exp
⁡
-
α
C
F
η
C
R
T
,
where(10)
j
C
0
=
j
O
2
R
d
0
r
e
f
p
O
2
p
O
2
r
e
f
exp
⁡
Δ
G
O
2
R
d
‡
R
1
T
r
e
f
-
1
T
,
(11)
η
C
=
Δ
ϕ
C
-
Δ
ϕ
C
E
Q
=
ϕ
P
t
C
-
ϕ
E
L
Y
C
-
ϕ
P
t
C
E
Q
-
ϕ
E
L
Y
E
Q
.Equation (1) describes the total mass conservation in the gas phase of the catalyst layers. The second term in this equation is the sink or source term describing the mass consumption or production in the gas phase caused by electrochemical reactions. Equation (2) is the linear form of the Ergun equation. It was used to calculate the pressure profiles in the gas phase of the catalyst layers because they are packed beds. At the conditions used in this study, the magnitude of the quadratic velocity term in the Ergun equation was much smaller than the linear term. Hence only the linear term in velocity was used in (2). Equations (1) and (2) were solved together to calculate the velocity and pressure profiles in the gas phase of the catalyst layers. Mass balances for each of the individual gas phase species account for convection, diffusion, and reaction, as shown in (3).Equations (4) and (5) describe, respectively, water and proton conservation in the electrolyte phase of the membrane and catalyst layers. Diffusion was described by concentrated solution theory through the use of the GMS equations. The following paragraphs illustrate the derivation of (4) and (5).A general procedure for the calculation of mass fluxes in multicomponent electrolyte systems was presented by Krishna [19]. It has been proven that the Nernst-Planck equation is a limiting case of the GMS equations. The GMS equations can be written as follows:(12)
d
→
i
=
∑
j
=
1
j
≠
i
n
x
i
J
→
j
-
x
j
J
→
i
c
E
L
Y
Ð
i
j
i
=
1,2
,
…
,
n
-
1
,where d
→
i is a generalized driving force for mass transport of species i. Because the summation of the n driving forces is equal to zero due to the Gibbs-Duhem limitation [31], only n
-
1 driving forces are independent. The equation to calculate the generalized driving force has been derived based on nonequilibrium thermodynamics [31]. A simplified expression for a solid stationary electrolyte (no convection term) [19] can be written as(13)
d
→
i
=
∇
x
i
+
x
i
z
i
F
R
T
∇
ϕ
E
L
Y
.For a noncharged species such as water z
i is equal to zero, and, according to (13), the concentration gradient will be the only driving force.The migration term in (13) was obtained by representing ion mobility by the Nernst-Einstein relation (D
i
=
R
T
u
i). This equation is applicable only at infinite dilution. However, it can be used in concentrated solutions if additional composition-dependent transport parameters, such as the B
′ parameters in (19), are used to calculate the flux of ions [27]. It will be shown in the following paragraphs that (18) represent the composition-dependent parameters.Equation (12) results in (n
-
1) independent equations that can be written in matrix form for convenience:(14)
c
E
L
Y
d
→
1
⋮
d
→
n
-
1
=
-
B
1
1
⋯
B
1
,
n
-
1
⋮
⋱
⋮
B
n
-
1
,
1
⋯
B
n
-
1
,
n
-
1
J
→
1
⋮
J
→
n
-
1
,where the elements of the matrix of inverted diffusion coefficients [B] are given by(15)
B
i
i
=
∑
j
=
1
j
≠
i
n
x
i
Ð
i
j
i
=
1,2
,
…
,
n
-
1
,
B
i
j
=
-
x
i
Ð
i
j
i
=
1,2
,
…
,
n
-
1
;
i
≠
j
.The fluxes of species, J
→
i, can be calculated from (16) which is the inversion of (14):(16)
J
→
1
⋮
J
→
n
-
1
=
-
c
E
L
Y
B
11
⋯
B
1
,
n
-
1
⋮
⋱
⋮
B
n
-
1,1
⋯
B
n
-
1
,
n
-
1
-
1
d
→
1
⋮
d
→
n
-
1
.For the present electrolyte system containing three species, mobile H2O and H+ plus immobile solid ZrP, (16) may be written as(17)
J
→
H
2
O
J
→
H
+
=
-
c
E
L
Y
B
H
2
O
–
H
2
O
′
B
H
2
O
–
H
+
′
B
H
+
–
H
2
O
′
B
H
+
–
H
+
′
d
→
H
2
O
d
→
H
+
,where [B
′] is the inverse of the matrix of inverted diffusion coefficients. Because Ð
H
2
O
–
H
+
=
Ð
H
+
–
H
2
O, the elements of [B
′] are calculated using (18) which are functions of the GMS diffusivities and the species mole fractions in the electrolyte phase:(18)
B
H
2
O
–
H
2
O
′
=
x
H
2
O
Ð
H
+
–ZrP
+
Ð
H
2
O
–
H
+
x
H
+
+
Ð
H
+
–ZrP
/
Ð
H
2
O
–ZrP
x
H
2
O
+
Ð
H
2
O
–
H
+
/
Ð
H
2
O
–ZrP
,
B
H
2
O
–
H
+
′
=
x
H
2
O
Ð
H
+
–ZrP
x
H
+
+
Ð
H
+
–ZrP
/
Ð
H
2
O
–ZrP
x
H
2
O
+
Ð
H
2
O
–
H
+
/
Ð
H
2
O
–ZrP
,
B
H
+
–
H
2
O
′
=
x
H
+
Ð
H
2
O
–ZrP
x
H
2
O
+
Ð
H
2
O
–ZrP
/
Ð
H
+
–ZrP
x
H
+
+
Ð
H
2
O
–
H
+
/
Ð
H
+
–ZrP
,
B
H
+
–
H
+
′
=
x
H
+
Ð
H
2
O
–ZrP
+
Ð
H
2
O
–
H
+
x
H
2
O
+
Ð
H
2
O
–ZrP
/
Ð
H
+
–ZrP
x
H
+
+
Ð
H
2
O
–
H
+
/
Ð
H
+
–ZrP
.Combining sets of (17) and (13) results in two independent equations that can be used to calculate the fluxes of mobile species (J
→
H
2
O and J
→
H
↑
+
↓
↑) within the electrolyte phase:(19)
J
→
H
2
O
=
-
c
E
L
Y
B
H
2
O
–
H
2
O
′
∇
x
H
2
O
-
c
E
L
Y
B
H
2
O
–
H
+
′
∇
x
H
+
+
F
x
H
+
R
T
∇
ϕ
E
L
Y
,
(20)
J
→
H
+
=
-
c
E
L
Y
B
H
+
–
H
2
O
′
∇
x
H
2
O
-
c
E
L
Y
B
H
+
–
H
+
′
∇
x
H
+
+
F
x
H
+
R
T
∇
ϕ
E
L
Y
.Equations (19) and (20) show that diffusion flux of each species is a function of the concentration gradient of all species as well as of the potential gradient. There are five unknowns in (19) and (20): J
→
H
2
O, J
→
H
+, x
H
2
O, x
H
+, and ϕ
E
L
Y. Therefore, three more equations are required.ZrP is immobile. As a result, the diffusion phenomenon will effectively be the interchange of H+ and H2O species. Therefore, for diffusion purposes we will only consider the domain of the mobile species, H+ and H2O, and will ignore the immobile species, ZrP. On that basis, (21) can be used as a third equation. Nevertheless, the presence of ZrP is important because of its interaction with the mobile species. Specifically, the values of the B
′ coefficients for H+ and H2O were influenced by the presence of ZrP:(21)
x
H
2
O
+
x
H
+
=
1.0
.The differential equations for H2O and H+ mass conservation in the electrolyte phase can be expressed in molar units as(22)
∇
·
J
→
H
2
O
=
-
j
z
F
,
∇
·
J
→
H
+
=
j
z
F
,where j is the volumetric current production. This quantity which appears in (1), (3) to (5) and (22) is the rate of production of protons in the anode. Therefore, it is positive in the anode, j
A, and negative in the cathode, j
C. It was calculated using the Butler-Volmer equation for the anode and cathode, (6) and (9), respectively. Exchange current densities at the anode and cathode are a function of the reactants’ partial pressure and the operating temperature as shown in (7) and (10). The Butler-Volmer equation and its parameters for both propane oxidation and oxygen reduction were described in our previous communication [10]. Complete conversion of C3H8 to CO2 was reported in experiments by Grubb and Michalske [34]. Equations (19) to (22) were combined and are shown as (4) and (5).


## 2.2. Numerical Procedure

The numerical solution procedure is illustrated in Figure3. Equations (1)–(11) define the problem at steady state. However, a time derivative was appended to each partial differential equation and a backward Euler time stepping method was used to increase stability while converging to the steady-state solution. The Finite Element Method was used to discretize the partial differential equations in space, with all dependent variables discretized by a linear finite element except for the pressure that is taken as a quadratic.Figure 3
Modeling procedure.FreeFEM++ software has been successfully used to solve two-dimensional partial differential equations (1)–(11). It is open-source software and is based on the Finite Element Method developed by Hecht et al. [32]. The calculated results from FreeFEM++  were exported to ParaView visualization software [35] for postprocessing. ParaView is also open-source software.There is no proton loss through the exterior boundaries of the domain (Figure2). Therefore, the total rate of proton production in the anode, ∫
A
n
o
d
e
j
d
V, has to be equal to the total rate of proton consumption in the cathode, ∫
C
a
t
h
o
d
e
(
-
j
)
d
V. In each case, the electrical potential of the catalyst phase of the anode, ϕ
P
t
A, and that of the cathode, ϕ
P
t
C, had individual constant values. Then all the variables in the whole domain were calculated. However, having fixed electrical potentials of the anode and cathode catalyst phases does not guarantee that the proton production at the anode will equal the proton consumption at the cathode. The difference between the rate of proton production and consumption can be minimized by shifting ϕ
E
L
Y by a constant value because the production and consumption rates are functions of the electrical potential in both of their respective catalyst phases, ϕ
P
t
A and ϕ
P
t
C, and in the electrolyte phase, ϕ
E
L
Y. Therefore, the Newton method was used to force equal proton production and consumption. In other words, balancing ∫
A
n
o
d
e
j
d
V and ∫
C
a
t
h
o
d
e
(
-
j
)
d
V acts as a constraint for the conservation of protons in the electrolyte phase.The equations for the conservation of momentum, total mass, and individual species in the gas phase of the anode and cathode were solved by assuming there was no species crossover through the membrane. Electrical potential, proton, and water concentrations in the electrolyte phase of the anode, cathode, and membrane layers were coupled to each other. These variables were calculated by solving (4), (5), and (21) iteratively in each layer. Then, the Robin method [10] was used to couple the solutions between layers. In the Robin method, both of the following transfer conditions are progressively satisfied on the anode catalyst/membrane interface and the membrane/cathode catalyst interface through iterations of (a) the continuity of the variable (e.g., potential) and (b) the continuity of the flux (e.g., electrical current).Figure2 shows four types of boundary conditions for the modeling domain, that is, inlet, outlet, wall of the land, and the midchannel symmetry boundaries. The flux of species in the gas phase is zero at the walls because there is no transfer through walls. The zero flux condition is also true at the midchannel symmetry boundaries. The compositions of the gaseous species are known at the inlet of the anode and cathode catalyst layers. It was assumed that no change in the composition of gas mixture occurred after leaving the catalyst bed. Therefore, the composition gradients are zero in the direction normal to the catalyst layer at the outlet boundaries. The zero flux condition is applied at all exterior boundaries for the species in the electrolyte phase.


## 2.3. Input Parameters

The parameters used for the simulations are shown in Table1. The GMS diffusivities, Ð
i
j, which are used in (18) have to be calculated from the Fickian diffusion coefficients, D
i
j. For ideal solutions, the Fickian diffusion, D
i
j, can be used as Ð
i
j in the Stefan-Maxwell equations [26] because the concentration dependence of Fickian diffusion coefficients is ignored. Experimental values for D
H
+
–ZrP and D
H
2
O
–
H
+ are given in Table 1. Note that the diffusivity of protons in ZrP is approximately two orders of magnitude smaller than the diffusivity of protons in water. The movement of protons causes the electroosmotic flow of water [9]. It was assumed that one water molecule is dragged by each proton, H3O+, that travels from anode to cathode. Therefore, the diffusivity of water in ZrP was set equal to the diffusivity of protons in ZrP [36], the smaller of the two proton diffusivities in Table 1. Proton diffusivity and proton mobility are different quantities. The three diffusivities in Table 1 were the ones used to calculate the B
′ parameters in (18).Table 1
Operational, electrochemical, and design parameters for simulations.
Property
Value
Temperature,T
423–503 K
Pressure,P
101.3 k Pa
Proton–ZrP diffusivity,D
H
+
–ZrP
3.1 × 10−12 m2 s−1 [29]
Proton–water diffusivity,D
H
2
O–
H
+
2.9 × 10−10 m2 s−1 [12]
Ionic conductivity in membrane,σ
ZrP
/
PTFE
5.0 S m−1 [24]
Electrical resistivity in membrane,R
PTFE
1.0 × 1016 Ω m
Charge transfer coefficients,α
A and α
C
1.0 [30]
Equilibrium potential of catalyst phase at the anode,ϕ
Pt
A
EQ
0.136 V [1]
Equilibrium potential of catalyst phase at the cathode,ϕ
Pt
C
EQ
1.229 V
Equilibrium potential of electrolyte phase,ϕ
ELY
EQ
0.136 V
Apparent bulk density of carbon catalyst support,ρ
CAT
⁡
0.259 gcatalyst
mL
catalyst
-
1
Specific surface area of carbon catalyst support in the anode and cathode,A
CAT
⁡
255m
catalyst
2 g
catalyst
-
1
Gas phase volume fraction in anode and cathode,ε
G
0.5
Electrolyte phase volume fraction in anode and cathode,ε
ELY
0.4
Effective particle diameter in anode and cathode,D
p
5μm
Land width,L
W
2–8 mm
Anode and cathode thickness, ThA, ThC
200–400μm
Membrane thickness,T
h
M
100–200μm
Fluid channels width in bipolar plates
0.4 mm


## 2.4. Model Validation

The model predicts the performance of a DPFC that (i) has interdigitated flow fields, (ii) has zirconium phosphate as the electrolyte, and (iii) operates over a temperature range of 150–230°C. As there are no experimental data for DPFCs having zirconium phosphate electrolytes and interdigitated flow fields, the model results have been compared to published results for DPFCs with other types of electrolytes and flow fields.Figure4 compares the modeling results for zirconium phosphate electrolyte with the experimental data for other types of electrolytes [34, 37]. The figure shows that the polarization curve for ZrP-PTFE electrolyte is somewhat comparable to that for the other electrolytes. The difference between the polarization curves can be partially explained by the difference between conductivities of the electrolytes. The proton conductivity of a nonmodified Nafion 117 approaches 10 S m−1 at 80°C [38]. The conductivity of the 95% H3PO4 electrolyte is 35 S m−1 at 200°C [39]. However, the proton conductivity for the best ZrP-PTFE that has been developed in our laboratory is about 5 S m−1 at 150°C.Figure 4
Polarization curves of direct propane/oxygen fuel cell using Pt anode and cathode. (a) Experimental results [31] using Nafion 117 at 95°C. (b) Experimental results [32] using 95% H3PO4 at 200°C. (c) The present proton migration and diffusion model results for a solid ZrP-PTFE electrolyte at 150°C.


## 3. Results and Discussion

Figure5(a) shows the two-dimensional variation of the proton concentration in the electrolyte phase of the entire domain, that is, the anode catalyst layer (AN), the membrane layer (ML), and the cathode catalyst layer (CA). The proton concentration at the anode inlet close to the feed gas channel has the highest value. This would be expected because the propane’s partial pressure is higher at the anode inlet and that causes a higher propane oxidation reaction rate, according to Butler-Volmer equation (6). Because protons are produced in the anode catalyst layer and consumed in the cathode catalyst layer, the proton concentration is greater at the anode than the cathode. The resulting proton concentration gradient is the driving force for protons to diffuse from the anode to the cathode.Figure 5
(a) Proton concentration in the electrolyte phase of the anode, membrane, and cathode layers. (b) Electrical potential profile for the electrolyte phase of the anode, membrane, and cathode layers. (c) Protonic flux from anode to cathode in the electrolyte phase. The vectors lengths indicate the flux magnitude which varies from 0 to 17 mA cm−2 in this case.
(a)
(b)
(c)The electrical potential variation in the electrolyte phase of the catalyst layers and membrane is shown in Figure5(b). As the reaction rate in the catalyst layers is not uniform, current density and electrical potential will be variable. Figure 5(b) shows that the electrical potential is higher at the cathode electrolyte phase than at the anode electrolyte phase. That electrical potential gradient is the driving force for protons to migrate from the cathode to the anode. This proton migration (caused by the electrical potential gradient) is in the opposite direction to the proton diffusion (caused by the proton concentration gradient) that was discussed above. In reality, protons are known to be transported from the anode to the cathode. Therefore the dominant driving force is the proton concentration gradient. Furthermore it can be concluded that the electrical potential gradient is not the dominant driving force for proton transport.Figure5(c) shows the magnitude and direction of protonic flux in the electrolyte phase of the anode, cathode, and membrane layers. Protons are produced in the anode and travel from the anode, through the membrane layer, and to cathode, where they are consumed. As discussed above, in Figure 5(a), the concentration driving force for proton flux was from anode to cathode and in Figure 5(b) the electrical potential driving force for protons was in the opposite direction, from cathode to anode. Finally, Figure 5(c) demonstrates that the net flux of protons is from the anode toward the cathode. As the net flux is the summation of two driving forces that are in opposite directions, again one can conclude that proton diffusion is dominant over proton migration. For the fuel cell to operate, the net transport of protons must be from the anode to the cathode. Therefore, the rate of proton diffusion must exceed the rate of proton migration. Figure 5(c) also shows that the arrows’ lengths are becoming longer (indicating that the proton flux increases) in the y-direction from the anode land/anode catalyst interface to the anode catalyst/membrane interface as more protons are produced throughout the anode catalyst layer. Similarly, the arrows’ length becomes shorter (as the proton flux decreases) in the y-direction from membrane/cathode catalyst interface to the cathode catalyst/cathode land interface.There are two routes by which electrons can flow from the anode to the cathode. The electron flux through the electrolyte is shown in Figure6. The electron flow rate through the electrolyte will be many orders of magnitude smaller than the electron flow rate through the external circuit. Although the vast majority of electrons flow through the external circuit, the production and consumption of the miniscule number of electrons that flow through the electrolyte have a distribution (Figure 6) that is similar to the distribution of protons (Figure 5(c)).Figure 6
Electronic flux from anode to cathode in electrolyte phase. The vectors lengths indicate the flux magnitude which varies from 0 to1
e
-
11 mA cm−2 in the same case as in Figure 5(c).It is constructive to compare this model (migration plus diffusion) with a migration-only model [10]. A cross section of Figure 5(b) along the y-direction at the middle of the domain (x
=
L
W
/
2) is shown in Figure 7(a), where the electrical potential for the migration plus diffusion model in the electrolyte phase (the left axis in Figure 7(a)—solid line) is compared with that in the two solid catalyst phases (the right axis in Figure 7(a)—dashed lines). The electrical potentials in each of the two solid catalyst phases (dashed line) are almost constant throughout their layers because these phases have high electrical conductivities. The greater electrical potential at the cathode than at the anode (both in the catalyst phases and the electrolyte phase) provides a driving force that (a) pushes positively charged protons from the cathode to the anode, via the electrolyte and (b) pushes negatively charged electrons from the anode to the cathode via both the external circuit (almost all the electrons) and the electrolyte (a miniscule quantity of electrons). The flow rate of negatively charged electrons through the electrolyte phase from the anode to the cathode will be miniscule.Figure 7
Electrical potential profiles in they-direction for the electrolyte and catalyst phases located at the middle of the domain x-direction for the cathode and anode catalyst layers and membrane layer. The arrows point in the direction of the ordinate scale that applies to each of the three curves. (a) Proton migration plus diffusion within the electrolyte phase (the present model) (b) Proton migration only within the electrolyte phase [5].
(a)
(b)The results of the migration plus diffusion model shown in Figure7(a) correctly describe these phenomena. In contrast, the results from the migration only model [10] are seen in Figure 7(b). Those calculations showed that the migration-only model produced incorrect results. Specifically, the electrical potential gradient in the electrolyte has the wrong slope. The slope (gradient) predicted by the migration-only model incorrectly drives the positively charged protons in the electrolyte from cathode to anode. In reality, they move from the anode to the cathode in the electrolyte.Figure8 compares the anodic and cathodic overpotential for two cases. The solid lines in Figure 8 are the results from the migration plus diffusion model. The dashed lines are the results from a migration only model. The dashed lines (migration-only) have a negative slope whereas the solid lines (migration plus diffusion) have a positive slope. Since the overpotential is the electrochemical driving force for the reaction (see (6) and (9)), it will always have its largest value adjacent to the anode land and decrease toward the membrane. In summary, the migration plus diffusion model predicted the correct behaviour, while the migration-only model predictions were incorrect.Figure 8
Overpotential profile in the anode and cathode alongy-axis at the middle of the modeling domain. Solid lines (migration plus diffusion). Dashed lines (migration only) [5].Figure9 shows the propane mole fraction in the gas phase of the anode catalyst layer along the x-direction. For similar operating conditions, the migration plus diffusion model predicted different propane concentrations than the migration-only model. This difference is caused by the different overpotential profiles predicted by the two models. The difference in overpotentials for migration plus diffusion compared to migration-only model is shown in Figure 8. Those differences are small. However, those small differences are in exponential terms, as shown in (6) and (9). It is the exponential terms that cause the large differences in concentration shown in Figure 9. If proton diffusion in the electrolyte phase is ignored, the prediction of species distribution within the gas phase of the catalyst layers becomes incorrect. In other words, the migration-only model can not correctly calculate either the proton concentration in the electrolyte phase or the propane concentration in the gas phase.Figure 9
Propane mole fraction in the gas phase of the anode catalyst layer along thex-direction at the middle of the anode catalyst layer. (a) Proton migration plus diffusion within the electrolyte phase (the present model). (b) Proton migration only within the electrolyte phase [5].In Figure10, the polarization curves for the migration plus diffusion model are compared with the migration-only model. At a specific cell potential, the cell current density predicted by the migration plus diffusion model is lower than that of the migration-only model. That is because the steady-state value for concentration occurs in the equation for the exchange current density, (7) and (9). This deviation may appear to be small at some conditions. In Figure 10, at a cell potential of 0.4 V, the migration plus diffusion model predicts a current density near 50 mA cm−2. In contrast, the migration-only model predicts nearly 70 mA cm−2. That is, one cannot conclude that a reasonable prediction of the fuel cell overall performance can be obtained using simple models that ignore the proton diffusion phenomenon in the electrolyte. In addition, there are other phenomena for which the migration-only model predicts results that are completely erroneous.Figure 10
Modeling results for polarization curves of direct propane/oxygen fuel cells using a solid ZrP-PTFE electrolyte at 150°C. (a) Proton migration and diffusion within the electrolyte phase (the present model). (b) Proton migration only within the electrolyte phase [5].It would be desirable to expand the range of the polarization curve in Figure10 to greater current densities and to smaller cell potentials. Many attempts to obtain such a wider range of values were made. Unfortunately they were all unsuccessful. As the current density increased, convergence to an acceptable numerical solution of the equations became progressively more difficult. Convergence was not obtained at values of current densities greater than those shown in Figure 10. The difficulty was caused by the exponential nature of the Butler-Volmer equation in combination with the complex Generalized Maxwell-Stefan equations. Small changes in cell potential cause the current density calculated from the Butler-Volmer equation to vary enormously. The search for superior convergence techniques is a topic that is being actively pursued in our laboratory.Activation overpotential and ohmic polarization are the major sources of potential drop in a direct propane fuel cell. Any change in the operating conditions or cell design that results in a decrease in activation overpotential and ohmic polarization will improve the cell performance. Figure11 shows the performance of a DPFC predicted by the model at different operating temperatures. It also shows the performance of a hydrogen PEM fuel cell at 80°C [40] and that of a DPFC at 200°C having a phosphoric acid electrolyte [34]. As temperature is increased from 150°C to 230°C, the rate of reaction increases according to (7) and (10). This leads to a decrease in the overpotential term in the Butler-Volmer equation and a major improvement in the cell performance. It can be concluded that the predicted performance of a DPFC operating at 230°C can approach that of a hydrogen PEMFC at 80°C when both operate at current densities less than 40 mA cm−2.Figure 11
(a), (b), and (c) Predicted polarization curves for a direct propane/oxygen fuel cell at different operating temperatures, (d) experimental data for a typical hydrogen/oxygen PEMFC [33], and (e) experimental data for the best performed DPFC at 200°C [32].


## 4. Conclusions

The migration plus diffusion model, described in this work, was shown to be superior to the migration-only model that is used in many fuel cell modeling studies. Specifically, the migration-only model predicted values of electrical potential in the electrolyte that are erroneous. The gradient of the electrolyte electrical potential predicted by the migration-only model was in the wrong direction. The incorrect values of the electrical potential in the electrolyte caused the values for the overpotential to be incorrect. Incorrect overpotential values caused the values calculated for the propane concentration to be incorrect. This work has shown that the predicted values for steady-state current density and steady-state propane concentration become substantially different when the effect of proton diffusion in the electrolyte is included in the model. The migration plus diffusion model described here has been shown to be a major improvement over the migration-only model that was used in earlier studies.Many important phenomena that occur in fuel cells are not described by polarization curves. Meaningful values for variables internal to the fuel cell, for example, overpotential and reactant concentration, are essential for the understanding of fuel cell performance. At some operating conditions, variables external to the fuel cell, for example, current density and the exit concentration of propane, are substantially different when proton diffusion in the electrolyte is included in the model. The insight obtained using the migration plus diffusion model is far more useful than that obtained from the migration-only model.



---

*Source: 102313-2015-12-21.xml*