Species Transport in the Gas Diffusion Layers of a PEM

This example focuses on the species transport within the gas diffusion layers (GDLs) of a proton exchange membrane (PEM) fuel cell. The geometry models a cell with two adjacent flow channels of different pressures, a situation that may occur in a cell with serpentine flow channels, or in a cell using a interdigitated flow field design. The model uses current balances, mass transport equations (Maxwell-Stefan diffusion for reactants, water and nitrogen gas), and momentum transport (Darcy’s law for the gas flows) to simulate a PEM fuel cell’s behavior.

Model Definition

The modeled section of the fuel cell consists of three domains: an anode (Ωa), a proton exchange membrane (Ωm), and a cathode (Ωc) as indicated in Figure 1.

Figure 1: Model geometry with domain and boundary labels.

Each of the electrodes (gas diffusion layers) is in contact with an interdigitated gas distributor, which has an inlet channel (∂Ωa,inlet), a current collector (∂Ωa,cc), and an outlet channel (∂Ωa,outlet). The same notation is used for the cathode side.

Humidified hydrogen and air are supplied to the inlet channels of the anode and cathode, respectively. Hydrogen reacts and is consumed at the anodic active layer to form protons that carry the ionic current to the cathode. Each proton is assumed to drag three molecules of water from the anode to the cathode, resulting in the following reaction formula for the anode:

(1)

Where λH2O is the number of water molecules that are transported per proton through the membrane by water-proton interaction.

At the cathode, oxygen reacts together with the protons to form water at the active layer according to:

(2)

Both feed gases (humidified hydrogen and humidified air) are treated as ideal and are transported through diffusion and convection. The electrodes are treated as homogeneous porous media with uniform morphological properties such as porosity and permeability. The gas within each of the electrodes exists as a continuous phase so Darcy’s law applies.

An agglomerate model describes the electrode reactions in the active catalyst layers. The agglomerates consist of catalyst and carbon particles embedded in polymer electrolyte. The equations for the agglomerate model originate from the analytical solution of a diffusion-reaction problem in a spherical porous particle (Ref. 5 and Ref. 6). At the anodic active catalyst layer, hydrogen is the diffusing and reacting species in the agglomerates, while oxygen is the diffusing and reacting species in the agglomerates at the cathode. An agglomerate model of the cathodic active catalyst layer of a PEM fuel cell has been presented by Broka and others (Ref. 7 and Ref. 8).

Charge Balances

The Hydrogen Fuel Cell interface is used for modeling the potential distributions in the three domains, with the charge-transfer current density of the active layers described as interior boundary conditions, using Thin H2/O2 Gas Diffusion Electrode nodes. The local current density expressions for the anode and cathode, ia and ic, are specified according to the equation

where the index e stands for “a” (anode) or “c” (cathode), Ragg is the agglomerate radius (SI unit: m) and jagg,a and jagg,c (SI unit: A/m3) are the current densities given by the agglomerate model. Further, Lact is the active layer (thin gas diffusion electrode) thickness (m) and εmac its porosity (the macroscopic porosity). The specific surface area of the thin gas diffusion electrode (SI unit: 1/m) is specified according to the following equation as

Agglomerate Model for Anode and Cathode

The agglomerate model describes the current density in an active layer consisting of agglomerates of ionic conductor material and electrically conducting particles covered partially with catalyst. The local current density can be expressed analytically by solving a combination of the diffusion equation and the Butler-Volmer electrode kinetic equation for an agglomerate with constant electric and ionic potentials. The resulting equations for the current density in the anode and cathode are (Ref. 7)

where, again, the index e stands for “a” (anode) or “c” (cathode), and

In these equations, Dagg is the agglomerate gas diffusivity (SI unit: m2/s), ne is a “charge transfer” number (1 for the anode and –2 for the cathode), S is the specific area of the catalyst inside the agglomerate (SI unit: 1/m), and F is Faraday’s constant (SI unit: C/mol). Furthermore, ci,ref are the reference concentrations of the species (SI unit: mol/m3), ci,agg are the corresponding concentrations in the agglomerate surface (SI unit: mol/m3), i0a and i0c are the exchange current densities (SI unit: A/m2), R is the gas constant, T is the temperature (SI unit: K), and the overvoltages at the anode and the cathode are given by

where Eeq (SI unit: V) denotes the equilibrium voltage.

You set the anodic and cathodic reference states equal to the molar fractions at the inlet channels of the anode and cathode, respectively, at 1 atm. The dissolved hydrogen and oxygen concentrations at the surface of the agglomerates are related to the molar fractions of the respective species in the gas phase through Henry’s law

where K is Henry’s constant (SI unit: Pa·m3/mol).

Charge Balances, Continued

The potential difference between the cathode and anode current collectors corresponds to the total cell voltage. Choose the potential at the anode current collector as the reference level by setting it to zero. Then the total cell voltage serves as the boundary condition at the cathode current collector:

For the other boundaries you have electric insulation boundary conditions.

Porous Media Fluid Flow

To model the gas flows in the gas backings, this example uses the Darcy’s law in the Hydrogen Fuel Cell interface. The gas velocity is given by the continuity equation according to

where ρ is the mixture density of the gas phase (SI unit: kg/m3) and u denotes the gas velocity (SI unit: m/s). Darcy’s law for porous media states that the gradient of pressure, the viscosity of the fluid, and the structure of the porous media determine the velocity:

Here kp denotes the electrode’s permeability (SI unit: m2), η represents the gas viscosity (SI unit: Pa·s), and p is the pressure (SI unit: Pa). The ideal gas law gives the gas phase’s mixture density, ρ:

(3)

In this equation, R denotes the gas constant (SI unit: J/(mol·K)), T is the temperature (SI unit: K), M is the molar mass (SI unit: kg/mol), and x is the mole fraction.

At the inlets and outlets you specify the pressure:

At the thin gas diffusion electrode boundary for the anode and cathode, the gas velocity is calculated automatically by the Hydrogen Fuel Cell interface, from the total mass flow given by the electrochemical reaction rate and the stoichiometric coefficients of Equation 1 and Equation 2, using Faraday’s law.

Maxwell-Stefan Mass Transport

The model takes into account two species in the anode — H2 and H2O — and three at the cathode — O2, H2O, and N2. The Hydrogen Fuel Cell interface uses Maxwell-Stefan multicomponent diffusion governed by the following equation

o solve for the mass fractions, wi. This particular PEM fuel cell model assumes that the temperature-driven diffusion is insignificant and sets the source term, R, to zero. The Maxwell-Stefan binary diffusion coefficients, Dij (SI unit: m2/s), are calculated automatically by the interface.

The feed-gas mole fractions are specified at the inlets. At the outlets, convective flux boundary conditions are applied, meaning that the flux is convection dominated. At the thin gas diffusion electrode boundaries, the mass fluxes of the species are automatically determined from the electrochemical reaction rate and stoichiometric coefficients, using Faraday’s law. The Stefan velocity contributions are also automatically calculated by the interface.

Results and Discussion

Figure 2 shows the current distribution in the PEM fuel cell. There are significant current spikes present at the corners of the current collectors.

Figure 2: Current density (surface plot) and current vector field (arrow plot) in the fuel cell operating at 0.7 V. The anode is on the left and the cathode is on the right.

To further analyze the cell’s behavior, plot the current density at the active layer as a function of cell height (y) as in Figure 3.

Figure 3: Current-density distribution at the active layer at the anode.

The current density is uneven with the highest density in the cell’s upper region. This means that the oxygen-reduction reaction rate in the cathode determines the current-density distribution. The maximum current density arises close to the air inlet.

The convective fluxes generally dominate mass transport in the cell. To study the convective effects, plot the velocity field as in Figure 4.

Figure 4: Gas velocity field in the anode and cathode compartments.

The flow-velocity magnitude attains its highest values at the current collector corners.

Figure 5 shows the reactant (oxygen and hydrogen) weight fractions in the cathode and anode gases. Surprisingly, the hydrogen fraction increases as the anode gas flows from the inlet (at the bottom) to the outlet (at the top). This is the result of the electro-osmotic drag of water through the membrane, which results in a higher flux than the consumption of hydrogen. This means that the resulting convective flux of anode gas toward the membrane causes the weight fraction of hydrogen to go up. In the cathode gas, there is an expected decrease in oxygen content along the flow direction.

Figure 5: Reactant mass fractions, normalized by their inlet values, on the anode side (left) and cathode side (right). The reactant in the anode is hydrogen and that in the cathode is oxygen.

Although oxygen consumption is small, the concentration overvoltage in the agglomerates gives a substantial contribution to the concentration overvoltage. A small change in the oxygen flow gives a substantial change in cell polarization.

Figure 6 depicts the water mass fraction in the anode and cathode gases as well as the diffusive flux of water in the anode. It is apparent that water is transported through both diffusion and convection to the membrane on the anode side. The results show a minimum occurring in the upper corner of the membrane on the anode side. This is known to limit fuel cell performance. If the anode gas becomes too dry, the membrane dries out, resulting in decreasing ionic conductivity and the cell subsequently fails.

On the other hand, on the cathode side water levels increase with the direction of flow, and a local maximum in water current occurs in the lower corner to the membrane. This might also be critical because water droplets can clog the pores and effectively hinder gas transport to the active layer.

Figure 6: Water mass fraction in the anode (left, surface plot) and the cathode (right, contour plot).

References

1. W. He, J.S. Yi, and T.V. Nguyen, “Two-Phase Flow Model of the Cathode of PEM Fuel Cells Using Interdigitated Flow Fields,” AIChe J., vol. 46, pp. 2053–2063, 2000.

2. C. Marr and X. Li, “Composition and Performance Modelling of Catalyst Layer in a Proton Exchange Membrane Fuel Cell,” J. Power Sources, vol. 77, pp. 17–27, 1999.

3. P. Futerko and I.-M. Hsing, “Two-Dimensional Finite Element Method Study of the Resistance of Membranes in Polymer Electrolyte Fuel Cells,” Electrochimica Acta, vol. 45, pp. 1741–1751, 2000.

4. D.M. Bernardi and M.W. Verbrugge, “Mathematical Model of a Gas Diffusion Electrode Bonded to a Polymer Electrolyte,” AIChe J., vol. 37, pp. 1151–1163, 1991.

5. H. Scott Fogler, Elements of Chemical Reaction Engineering, 3rd ed., Prentice Hall, 1999.

6. R.B. Bird, W.E. Stewart, and E.N. Lightfoot, Transport Phenomena, John Wiley & Sons, 1960.

7. K. Broka and P. Ekdunge, “Modelling the PEM fuel cell cathode,” J. Appl. Electrochem., vol. 27, pp. 281–289, 1997.

8. K. Dannenberg, P. Ekdunge, and G. Lindbergh, “Mathematical model of the PEMFC”, J. Appl. Electrochem., vol. 30, pp. 1377–1387, 2000.

Fuel_Cell_and_Electrolyzer_Module/Fuel_Cells/pem_gdl_species_transport_2d

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

Create the geometry using rectangles.

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.5e-4.

In the text field, type 2e-3.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 1e-4.

In the text field, type 2e-3.

Locate the section. In the text field, type 2.5e-4.

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2.5e-4.

In the text field, type 2e-3.

Locate the section. In the text field, type 3.5e-4.

Rectangle 4 (r4)

In the toolbar, click

In the window for , locate the section.

In the text field, type 1e-4.

In the text field, type 5e-4.

Locate the section. In the text field, type -1e-4.

Array 1 (arr1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

In the text field, type 2.

Locate the section. In the text field, type 7e-4.

In the text field, type 1.5e-3.

Union 1 (uni1)

In the toolbar, click

and choose .

Select the objects , , and only.

In the window for , locate the section.

Clear the check box.

Union 2 (uni2)

In the toolbar, click

and choose .

Select the objects , , and only.

In the window for , locate the section.

Clear the check box.

In the toolbar, click

Global Definitions

Proceed to load a set of global model parameters from a text file provided with the Application Library.

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file pem_gdl_species_transport_2d_parameters.txt.

Definitions

Set up a number of selections on the geometry. These will be used later when setting up the physics.

Anode GDL

In the toolbar, click

Select Domain 1 only.

Right-click and choose .

In the dialog box, type Anode GDL in the text field.

Click .

Membrane

In the toolbar, click

Select Domain 2 only.

Right-click and choose .

In the dialog box, type Membrane in the text field.

Click .

Cathode GDL

In the toolbar, click

Select Domain 3 only.

Right-click and choose .

In the dialog box, type Cathode GDL in the text field.

Click .

Anode GDE

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundary 10 only.

Right-click and choose .

In the dialog box, type Anode GDE in the text field.

Click .

Cathode GDE

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundary 13 only.

Right-click and choose .

In the dialog box, type Cathode GDE in the text field.

Click .

Anode GDL Variables

Define a number of domain- and boundary-specific variables. These will be used both for setting up the physics and in postprocessing.

In the toolbar, click

In the window, right-click and choose .

In the dialog box, type Anode GDL Variables in the text field.

Click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
wReact	fc.wH2		Mass fraction, reactant

Cathode GDL Variables

In the toolbar, click

Right-click and choose .

In the dialog box, type Cathode GDL Variables in the text field.

Click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
wReact	fc.wO2		Mass fraction, reactant

Anode GDE Variables

In the toolbar, click

In the window, right-click and choose .

In the dialog box, type Anode GDE Variables in the text field.

Click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
cH2_agg	fc.p*fc.xH2/KH2		Henry's law hydrogen agglomerate concentration
eta_a	fc.phis-fc.phil-E_eq_a	V	Anodic overpotential
beta_a	cH2_agg-cH2_refexp(-2F_consteta_a/(R_constT))
lda_a	sqrt(i0_aSR_agg^2/(2F_constcH2_ref*D_agg))		Anodic current density subexpression
i_a	K(1-lda_acoth(lda_a))beta_a(R_agg/3)		Anode current density

Cathode GDE Variables

In the toolbar, click

Right-click and choose .

In the dialog box, type Cathode GDE Variables in the text field.

Click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
cO2_agg	fc.p*fc.xO2/KO2		Henry's law oxygen agglomerate concentration
eta_c	fc.phis-fc.phil-E_eq_c	V	Cathodic overvoltage
lda_c	sqrt(i0_cSR_agg^2exp(-F_consteta_c/(2R_constT))/(4F_constcO2_ref*D_agg))		Cathodic current density subexpression
i_c	-2K(1-lda_ccoth(lda_c))cO2_agg*(R_agg/3)		Cathode current density

Hydrogen Fuel Cell (fc)

Set up the current distribution and transport model. Include mass transport using Maxwell-Stefan diffusion and momentum transport using Darcy’s Law in both the anode and cathode gas mixtures. Also, define the reference pressure level in the interface properties. Note that the default gas species are hydrogen and water on the anode side, and oxygen, nitrogen and water on the cathode side. Start with adding the relevant domain nodes.

In the window for , locate the section.

Find the subsection. Select the check box.

Locate the section. Select the check box.

Click to expand the section. In the pref text field, type 0.

Membrane 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

H2 Gas Diffusion Layer 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

O2 Gas Diffusion Layer 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Thin H2 Gas Diffusion Electrode 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Thin O2 Gas Diffusion Electrode 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Electrolyte Phase 1

Set up the electrolyte conductivity in the node.

In the window, click .

In the window for , locate the section.

In the σl text field, type kappa_m.

Inspect the settings in the and nodes. Note that the density and viscosity of the gas mixture, and the binary diffusion coefficients are calculated automatically when the respective default settings are used.

H2 Gas Diffusion Layer 1

Set up the properties of the and nodes.