A Mathematical Model of a Direct Propane Fuel Cell

1Department of Chemical and Biological Engineering, University of Ottawa, 161 Louis-Pasteur, Ottawa, ON, Canada K1N 6N5 2Center for Catalysis Research and Innovation, University of Ottawa, 30 Marie-Curie Street, Ottawa, ON, Canada K1N 6N5 3Department of Mathematics and Statistics, University of Ottawa, 585 King Edward Avenue, Ottawa, ON, Canada K1N 6N5 4EnPross Incorporated, 147 Banning Road, Ottawa, ON, Canada K2L 1C5


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][11][12][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][15][16][17][18][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][21][22][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][25][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 (Zr(HPO 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.

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 Figure 1.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.
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.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: where  = C 3 H 8 , H 2 O, and CO 2 for the anode and O 2 and H 2 O for the cathode.
Conservation of momentum in gas phase: Conservation of noncharged species in gas phase: where  = C 3 H 8 and CO 2 for the anode and O 2 and H 2 O for the cathode.
Conservation of species in the electrolyte phase: for water, for proton, Butler-Volmer equation in the anode: where Butler-Volmer equation in the cathode: where 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: where ⃗   is a generalized driving force for mass transport of species .Because the summation of the  driving forces is equal to zero due to the Gibbs-Duhem limitation [31], only  − 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 For a noncharged species such as water   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  = u  ).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   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 (−1) independent equations that can be written in matrix form for convenience: ) , (14) where the elements of the matrix of inverted diffusion coefficients [] are given by The fluxes of species, ⃗   , can be calculated from ( 16) which is the inversion of ( 14): ) .

(16)
For the present electrolyte system containing three species, mobile H 2 O and H + plus immobile solid ZrP, ( 16) may be written as where [  ] is the inverse of the matrix of inverted diffusion coefficients.Because Đ , the elements of [  ] are calculated using (18) which are functions of the GMS diffusivities and the species mole fractions in the electrolyte phase: Combining sets of ( 17) and ( 13) results in two independent equations that can be used to calculate the fluxes of mobile ) within the electrolyte phase: 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) ,  H + , and  ELY .Therefore, three more equations are required.
ZrP is immobile.As a result, the diffusion phenomenon will effectively be the interchange of H + and H 2 O species.Therefore, for diffusion purposes we will only consider the domain of the mobile species, H + and H 2 O, 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   coefficients for H + and H 2 O were influenced by the presence of ZrP: The differential equations for H 2 O and H + mass conservation in the electrolyte phase can be expressed in molar units as where  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,  A , and negative in the cathode,  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 C 3 H 8 to CO 2 was reported in experiments by Grubb and Michalske [34].Equations ( 19) to (22) were combined and are shown as ( 4) and (5).

Numerical
Procedure.The numerical solution procedure is illustrated in Figure 3. 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.
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 opensource software.
There is no proton loss through the exterior boundaries of the domain (Figure 2).Therefore, the total rate of proton production in the anode, ∫ Anode , has to be equal to the total rate of proton consumption in the cathode, ∫ Cathode (−).In each case, the electrical potential of the catalyst phase of the anode,  Pt A , and that of the cathode,  Pt 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  ELY by a constant value because the production and consumption rates are functions of the electrical potential in both of their respective catalyst phases,  Pt A and  Pt C , and in the electrolyte phase,  ELY .Therefore, the Newton method was used to force equal proton production and consumption.In other words, balancing ∫ Anode  and ∫ Cathode (−) acts as a constraint for the conservation of protons in the electrolyte phase.

Journal of Chemistry
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).
Figure 2 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.

Input Parameters.
The parameters used for the simulations are shown in Table 1.The GMS diffusivities, Đ  , which are used in (18) have to be calculated from the Fickian diffusion coefficients,   .For ideal solutions, the Fickian diffusion,   , can be used as Đ  in the Stefan-Maxwell equations [26] because the concentration dependence of Fickian diffusion coefficients is ignored.Experimental values for  H + -ZrP and  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, H 3 O + , 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   parameters in (18).

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.
Figure 4 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% H 3 PO 4   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.

Results and Discussion
Figure 5(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.
The electrical potential variation in the electrolyte phase of the catalyst layers and membrane is shown in Figure 5(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.Figure 5(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 -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 -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 Figure 6.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)).
It is constructive to compare this model (migration plus diffusion) with a migration-only model [10].A cross section of Figure 7: Electrical potential profiles in the -direction for the electrolyte and catalyst phases located at the middle of the domain -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].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.
The results of the migration plus diffusion model shown in Figure 7(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 migrationonly 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.
Figure 8 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 9 shows the propane mole fraction in the gas phase of the anode catalyst layer along the -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.
In Figure 10, 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.
It would be desirable to expand the range of the polarization curve in Figure 10 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.Figure 11 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 .

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.Total pressure (kPa) PTFE: Polytetrafluoroethylene :

Greek Letters
A and  C : Anodic and cathodic charge transfer coefficients : Volume fraction : Overpotential (V) : Moles of water per mole of sulfonic acid sites : Dynamic viscosity (kg m −1 s −1 ) ]  : Stoichiometric coefficient of species , positive for reactants and negative for products : Mass density (kg m −3 )  CAT : Apparent bulk density of catalyst support (kg catalyst m −3 catalyst )  ZrP/PTFE : Ionic conductivity in membrane (S m −1 ) : Electrical potential (V) EQ Pt : Equilibrium potential of catalyst phase (V)  EQ ELY : Equilibrium potential of electrolyte phase (V).

Figure 4 :
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% H 3 PO 4 at 200 ∘ C. (c) The present proton migration and diffusion model results for a solid ZrP-PTFE electrolyte at 150 ∘ C.

Figure 5 :
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.

Figure 6 :
Figure 6: Electronic flux from anode to cathode in electrolyte phase.The vectors lengths indicate the flux magnitude which varies from 0 to 1 − 11 mA cm −2 in the same case as in Figure 5(c).

Figure 5 (
b) along the -direction at the middle of the domain ( =   /2) is shown in Figure7(a), where the electrical potential for the migration plus diffusion model in the electrolyte phase (the left axis in Figure7(a)-solid line) is compared with that in the two solid catalyst phases (the right axis in Figure7(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

Figure 8 :
Figure 8: Overpotential profile in the anode and cathode along axis at the middle of the modeling domain.Solid lines (migration plus diffusion).Dashed lines (migration only) [5].

Figure 9 :Figure 10 :
Figure 9: Propane mole fraction in the gas phase of the anode catalyst layer along the -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].

2 Figure 11 :
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].

Table 1 :
Operational, electrochemical, and design parameters for simulations.
Ox: Propane oxidation reaction on Pt catalyst ELY: Electrolyte phase in the membrane, anode, and cathode catalyst layers containing solid ZrP and mobile H 2 O and H + EQ: Equilibrium state G: Gas mixture : Species in gas or solid phase: propane, water, CO 2 , O 2 , H + , and ZrP ML: Membrane layer O 2 Rd: Oxygen reduction reaction on platinum