# Continuum Mathematical Model of Fuel Cells

Mathematical modeling plays an important role in the design of proton ex­change membrane fuel cells. In particular, there is a wide variety of continuum macroscopic models that are considered successful as they agree well with the

Constants used in the model are declared at the top.

 A В C E F G H 1 Electroclierni мі Model of =EM Fuel Cel 3 7 2 3 Eo 1.229 V T 333.15 4 R 8.314 P H2 0 26 5 F 96487 P 02 1 6 Rs 0.0035 7 xl 1 -0.948 C 02 8.7767E-007 3 xi 2 3.10E-003 n c 1 9 xl 3 7.60E-005 10 xl 4 -1.9ЭЕ-004 11 в 0.2 12 I L 42 13 14 15 E N V activ V cone V ohms V cell 16 1 1.20966506 0.268339344 0.00481951 0.0035 0.933006226 17 2 1.20966508 0 312907267 0.009758033 0.007 0.87999976 13 0.338977862 0.014821594 0.0105 0.845365624 19 20966508 0.357475229 0.020016692 0.014 0.818173159 20 1.20966508 0.371822902 0.025350341 0.0175 0.794991836

The current is in column B, the open circuit voltage (calculated using the Nernst expression) is in column C, the activation voltage drop is in column D, the effects of reduced concentrations at the electrodes are in column E, and the effects of the membrane and electrode resistance are in column F. Combining all of the these effects results in the cell’s voltage in column G.

FIGURE 7.6

Spreadsheet implementation of a simple electrochemical model of a fuel cell.   experimental data. (Although it should be noted that when a model or sim­ulation has a number of free parameters, then agreement with experimental data is often either fortuitous or obtained through fitting the free parameters.) In general, the model or simulation is considered valid when the experimental and modeled polarization curves are in good agreement. However, it is often found that different models and simulations, that are significantly different from one another, are found to agree with the same experimental data. In­creasingly complex models are continually being developed, incorporating full dynamic three-dimensional fuel cell models, multiple phase fluid flows, heat transfer and other processes, that capture the important physics of fuel cells. Rather than explore these increasingly complex models and simulations that are emerging within this field, however, we will look at an earlier model which is more suited to spreadsheet modeling.

The model we will consider was originally developed by Springer et al (J. Electrochemical Soc., 1991). The model is a simple one-dimensional model which is capable of capturing the fuel cell current-voltage performance, and was found to agree reasonably well with experimental data. This model is important, not only because it is the cornerstone from which later models were descended, but also because it emphasized the importance of keeping the membrane well hydrated, and the importance of water management. In particular, the model is steady-state, isothermal and one-dimensional, with the fuel cell consisting of five layers (membrane, diffusion media, and catalyst layers), and importantly this model considered the water content in the mem­brane (allowing for variable properties in the membrane such as conductivity and the water diffusion coefficient).

In the one-dimensional model all quantities are considered to vary only in the direction perpendicular to the anode and cathode interfaces. Recall that at the anode the hydrogen reacts to yield electrons to the platinum and protons to the membrane material

H2 ^ 2H + + 2e – (7.23)

At the cathode, the oxygen reacts with the protons coming through the mem­brane to form molecules of water

O – + 2H + ^ H2O (7.24)

The crux of this model, therefore, is to capture the water transport throughout the cell. This includes transport through the porous electrodes, and transport through the membrane electrolyte, to obtain the steady-state water distribu­tion throughout the cell.

A schematic of how the model works is shown in Figure 7.8. The interfaces are labeled 1 through 4, as we go from the anode to the cathode, and the water concentrations are obtained at these outer interfaces first. Then the water content is calculated using differential equations (that we will cover shortly) starting from interface 1 and working our way through to interface 3,  FIGURE 7.8

Schematic of the continuum mathematical model of a fuel cell.

while at the same time the water content is calculated starting from interface 4 and working our way through to interface 3 through the cathode. In this manner we obtain two separate solutions for the water content at interface 3. The ratio of water flux to the molar flux of H2 through the anode is a free parameter, which can be varied until the two solutions agree, giving us a consistent distribution of water content throughout the fuel cell. Note that not only does water enter the fuel cell as humidified fuel and air streams, but it is also generated at the cathode by the oxygen and hydrogen reaction.

The cathode inlet stream is assumed to be saturated with water vapor, and the mole fraction of water is given by

P SAT

xwc _ p (7.25)

where the subscripts W and C represent water and cathode, respectively. PSAT is the saturated pressure at the cathode and Pc is the pressure at the cathode. Essentially this sets up a Dirichlet boundary condition (where we assume that we know the water content at this boundary).

The mole fraction of water and oxygen at interface 4 depend on the mole fraction of water in the cathode inlet stream. In particular, the mole fraction of water at interface 4 is set up as a ratio of the water flow and the total flow. _ xWcvO + 2(1 + a) (l

vo + (2a + 1) 1 – xWc xon

and, similarly for oxygen

xo4 _ (vo – 1)(1 – xWe) x°N (7.27)

vO + (2a + 1) 1 — xwc xON

where vo is a stochiometric coefficient and the ratio of inlet flux to the total flux across interface 4, and xoN is the inlet dry gas mole fraction of oxygen (and is assumed to be known). The above equations, in a relatively straight­forward manner, give us the water and oxygen mole fractions at interface 4.

 Note, however, that this depends on the free parameter, a, which is the ratio of water flux to the molar flux of H2. Once we have the mole fractions of oxygen and water at interface 4 we can obtain the mole fractions at interface 3 by integrating dxwc dz

 RTI ~PC

 (7.28)

 and dxo dz

 RTI ~PPC

 (7.29)

 where xwc and xo are the mole fractions of water and oxygen in the cathode, respectively. z is the distance through the cathode, R is the gas constant, T is temperature, I is current density, Pc is the pressure at the cathode, and a is again the ratio of water flux to the molar flux of H2. DWN is the binary diffusion coefficient between water and nitrogen, DoN is the binary diffusion coefficient between oxygen and nitrogen, and Dwo is the binary diffusion coefficient between water and oxygen. We can integrate the above equations using a simple finite difference approximation. In other words, we can write

 xi+i xwc

 xwc+ RTI   (7.33)

 xVVA a (1 xw^ + VH 1

where vH is a stochiometric coefficient and the ratio of inlet hydrogen flux to the total flux across interface 1.

The mole fraction of water at interface 2 can be obtained through inte­grating the following equation across the anode.

 dxwA dz

 RTI — —————— [xwa(1 + a) – a] PA DWH

 (7.34)

where xWa is the water mole fraction in the anode and DWHis the binary diffusion coefficient between water and hydrogen. This differential equation is much easier to solve than Equations 7.28 and 7.29, and a simple analytical equation can be obtained to give us the mole fraction of water at interface 2

 exp ^ + _a_ 1 + a) pV Pa Dwh) + 1 + a

 (7.35)

 XW 2

 XW 1

where tA is the anode thickness.

Now we have obtained the water content at interface 2, we must obtain the water content at interface 3, which we will then be able to compare with the solution from earlier which involved the numerical integration of water content from interface 4 to interface 3. When the parameter a is chosen such that these two solutions agree then the simulation is solved.

To get from interface 2 to interface 3 we have to go through the membrane, and for that we must consider A, the membrane water content or the local ratio H2O/SO3 in the membrane. This is calculated directly from the mole fraction of water content. First we calculate the water vapor activity using xW P

PSAT

where xW is the mole fraction of water, P is pressure and PSaT is the satura­tion pressure. Next we calculate the water content using a phenomenological equation of the form ( 0.043 + 17.81a – 39.85a2 + 36a3 if a < 1;

14 + 1.4(a – 1) if 1 < a < 3.

It is this water content which we will integrate across the membrane using the differential equation

dA л IM,

dz = [2ndrag – a] PdryD(A) (7.38)

where ndrag is the number of molecules per proton, M is the molecular weight of the membrane, Pdry is the density of dry membrane, and the diffusion coefficient (which depends on water content) is given by 1

(7.39)

From the above system of equations we can obtain the water content at interface 3, A3, both through the integration of water content from the anode and the cathode. We can vary a, the ratio of water flux to the molar flux of H2, until these two solutions converge and we have a complete calculation of the water content throughout the fuel cell. From here we can then calculate the voltage across the fuel cell, and hence obtain the polarization curve.       The conductivity of the membrane as a function of water content, A, has been experimentally determined to be of the form   where T is temperature. From this conductivity, we can calculate the resis­tance of the membrane using    where the integration is over the entire membrane. To obtain the polarization curve (voltage-current curve) the model accounts for kinetic losses and ohmic losses in the membrane. To calculate kinetic losses, the Tafel equation is used, which relates the current density to the overpotential, q.

where jo is the exchange current density (referenced to pure oxygen at 1 atm at the open cell potential, Voc). Pa is the pressure at the cathode, xliq is the mole fraction of liquid which is generally very small (and often ignored), F is Faraday’s constant, and n is given by n Voc Vcell JR

where Voc is the open cell voltage, and vcell is the voltage across the cell, which we are trying to obtain in order to generate the polarization curve.

To summarize, the simulation involves calculating the water content at interface 3, A3, from both the anode and cathode sides of the fuel cell such that the two solutions agree. Then we can take the water content throughout the cell and calculate the resistance and voltage across the cell. One of the parameters used throughout these calculations is the current density. There­fore, we can obtain the voltage across the cell as a function of current density, and generate the polarization curve.

Let’s turn our attention to the spreadsheet implementation of the model described above. Figure 7.9 depicts the top of the spreadsheet, where the many constants are contained. Most of these numbers are obtained directly from the literature, but the current density in cell B8 is changed to obtain the polarization (or voltage-current) curve. To the right of Figure 7.9 the The model uses a large number of constants, most of which are obtained experimentally. The current density in cell В8 is the variable changed when obtaining the voltage-current curves.

from the cathode. When an a is found which results in these two values being close together then the simulation is considered solved.

FIGURE 7.9

Spreadsheet implementation of the continuum mathematical model of a fuel cell.

variable a is determined (a is stored in cell H3). As a function of this variable, the water content at interface 3, A3, is obtained twice (as described in the description of the model above). These two A’s are copied from elsewhere in the spreadsheet model (which we will come back to later) in to the cells H5 and H6. The difference between these two solutions is calculated in H8. The objective, therefore, is to find a value of a which would minimize this difference. While this could be done numerically, using for example the bisection method, in the current implementation this is done by hand (although it doesn’t take too many iterations to converge).

The calculation of the two A3’s, and ultimately the voltage across the cell, is depicted in Figure 7.10. In cells A20 and A23 the mole fractions of wa­ter and oxygen, respectively, are calculated at interface 4. In terms of these concentrations, the mole fractions of water and oxygen are calculated by inte­grating over distance. These integrations (the equations of which are covered above) occur in cells D20:BY20 and D23:BY23, for the mole fraction of water and oxygen, respectively. For example, the code contained in cell E20 is of the form

The water and oxygen levels are calculated through the cathode to obtain the membrane water content at the cathode-membrane interface (interface 3)

 A В C D E F G H I J 17 Г 7 18 о 1 2 3 4 5 6 19 x w4 x у 0 0.000005 0.00001 0.000015 0.00002 0.000025 0.00003 20 0.231848637 X wc 0.231848637 0.231980124 0.232111596 0.232243054 0.232374498 0.232505927 0.232637343 21 22 x o4 23 0.139302061 X о 0.139302061 0.139214072 0.139126065 0.139038102 0.138950121 0.138862144 0.136774169 24 25 0 1 26 z 0 0.000005 27 x w1 x w2 a Lambda 2 D(lambda) 5.4463E-006 5.4144E-006 28 0 066523777 0.065967365 0.423410556 3.172445763 lambda 3.172445763 3.229143285 29 30 / sigma 0.023594587 0.02412166 31 /. dR 0.000211913 0.000207283 32 33 R 0.004354066 34 V cell 0.759394462 35

The water content is calculated through the anode and through the membrane to again obtain the membrane content at the cathode-membrane interface.

FIGURE 7.10

Spreadsheet implementation of the continuum mathematical model of a fuel cell (continued).

( (1-D20-D23)*(1+\$H\$3)/\$B\$12 + (0.5*D20+D23*(1+\$H\$3))/\$B\$13)

and corresponds with the calculation of the mole fraction of water in Equation 7.28 above, while for the mole fraction of oxygen we discretize Equation 7.29 and the code in cell E23 is of the form

=D23-\$E\$15*(\$B\$6*\$B\$7*\$B\$10/\$B\$11) *

((D23*(1-\$H\$3)+0.5*D20)/\$B\$14 + (1-D20-D23)/\$B\$13)

In cell BY20, at the end of the spreadsheet to the right (not shown), the mole fraction of water at interface 3 is obtained. From this it is easy to obtain the water content, and this is then displayed in cell H5 (Figure 7.9).

In cell A28 the mole fraction of water at interface 1 is calculated, and the mole fraction of water at interface 2 is determined in cell C28. The code in cell C28 is

=(A28-(H3)/(1+H3))*EXP(B6*B7*B10*E7/(E5*E6)) + (H3)/(1+H3)

and the integration of mole fraction of water is obtained analytically. To ob­tain the water content at interface 3, we first calculate the water content at interface 2 and integrate this through the membrane to interface 3. The water vapor activity and water content at interface 2 are calculated in cells E28 and

F28, respectively. The diffusion coefficient of water is calculated from the wa­ter content in cells I27:AR27. For example, the diffusion coefficient of water calculated in cell J27 uses the code

=0.00000309381428211097*

(2.563-0.33*J28+0.0264*J28*J28-0.000671*J28*J28*J28)

The water content throughout the membrane is now calculated in cells I28:AR28. For example, the code in cell J28 is

=I28+\$E\$15*(2*2.5*128/(22) – \$H\$3)*(\$B\$10*\$E\$10/(\$E\$11*I27))

which is a simple finite difference approximation of the differential equation for water content (Equation 7.38). This ultimately yields the water content at interface 3, in cell AR28 (not shown) and this is then displayed in cell H6 (Figure 7.9). By comparing these two A’s (in cells H5 and H6) we can find a value of a such that the two solutions reasonably agree, and the water content throughout the membrane.

Once we have the water content throughout the membrane, we can deter­mine the resistance of the membrane and finally the voltage across the cell. Cells I30:AR30 calculate the conductivity of the membrane (as a function of the water content). This is then integrated in cells I31:AR31 to obtain the resistance of the membrane in cell I33. In other words, the numerical inte­gration (essentially a summation) of the resistance of small portions of the membrane are added together to give the total resistance of the membrane. From this the voltage across the cell is obtained in cell I34. By systematically varying the current density in cell B8, and solving the above simulation to obtain the voltage across the cell, the polarization curve can be obtained. The polarization cell, obtained from the spreadsheet model, is depicted in Figure 7.11. For high currents, the concentration losses play a more important role and these are not taken into consideration here. Therefore the model is more accurate for regions of the plot with low to intermediate current. The curve, however, is just as we might expect, and this model has been successfully fitted to experimental data, validating the model to some degree.

Updated: September 25, 2015 — 8:58 pm