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 channels (∂Ωa,inlet) and a current collector (∂Ωa,cc). 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.

(1)

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 you specify the pressure as follows:

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 equations

to 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 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.

The Membrane Transport features of the Hydrogen Fuel Cell interface are used to model the transport of water in the ionomer phase in the membrane domain. The molecular flux of water depends both on chemical potential gradient driven permeation and electroosmotic drag, using experimentally estimated parameters available in the Nafion material in the Fuel Cell and Electrolyzer Material Library.

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.

Figure 5: Reactant mass fractions 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. Water is transported through both diffusion and convection to the membrane on the anode side. 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 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.

Fillet 1 (fil1)

In the toolbar, click

In the window for , locate the section.

Find the subsection. Click to clear the

toggle button.

Click the

button in the toolbar.

On the object , select Points 6 and 7 only.

On the object , select Points 4 and 5 only.

Click to select the

toggle button.

On the object , select Points 6 and 7 only.

On the object , select Points 4 and 5 only.

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

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.

In the window for , type Anode GDL in the text field.

Membrane

In the toolbar, click

Select Domain 2 only.

In the window for , type Membrane in the text field.

Cathode GDL

In the toolbar, click

Select Domain 3 only.

In the window for , type Cathode GDL in the text field.

Anode GDE

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundary 10 only.

In the text field, type Anode GDE.

Cathode GDE

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Boundary 13 only.

In the text field, type Cathode GDE.

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 for , type Anode GDL Variables in the text field.

Locate the section. From the list, choose .

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

In the window for , type Cathode GDL Variables in the text field.

Locate the section. From the list, choose .

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 for , type Anode GDE Variables in the text field.

Locate the section. From the list, choose .

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

In the window for , type Cathode GDE Variables in the text field.

Locate the section. From the list, choose .

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

Materials

This model uses polymer electrolyte material (Nafion 1100, vapor equilibrated) that is available in the material library. Add the material twice, and assign it first to the membrane domain, and next to the membrane boundaries adjacent to the gas phase nodes.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Right-click and choose .

Materials

Nafion, EW 1100, Vapor Equilibrated, Protonated (mat1)

In the window, under click .

In the window for , locate the section.

Click

Select Domain 2 only.

Add Material

Go to the window.

In the tree, select .

Right-click and choose .

In the toolbar, click

to close the window.

Materials

Nafion, EW 1100, Vapor Equilibrated, Protonated 1 (mat2)

In the window for , locate the section.

From the list, choose .

Select Boundaries 10 and 13 only.

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, along with electroosmotic water drag in the membrane. 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, under click .

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. Select the check box.

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

Membrane 1

In the toolbar, 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 .

In the node, the electrolyte conductivity is set to be taken from the node. 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.

The properties for electroosmotic water drag in the node and in the child nodes that added by default are automatically taken from the node.

Initial Values 1

In the window, expand the node, then click .

In the window for , locate the section.

In the T0 text field, type T.

H2 Gas Diffusion Layer 1

Set up the properties of the and nodes.

In the window, under click .

In the window for , locate the section.

In the σs text field, type kappa_s.

Locate the section. In the εg text field, type eps_mac.

In the κg text field, type kappa_p.

O2 Gas Diffusion Layer 1

In the window, click .

In the window for , locate the section.

In the σs text field, type kappa_s.

Locate the section. In the εg text field, type eps_mac.

In the κg text field, type kappa_p.

Thin H2 Gas Diffusion Electrode 1

Set up the properties of the node. The details of electrode kinetics are set in the child node.

In the window, click .

In the window for , locate the section.

In the dgde text field, type l_act.

Thin H2 Gas Diffusion Electrode Reaction 1

In the window, click .

In the window for , locate the section.

From the Eeq list, choose . Due to the user-defined kinetics expression, the equilibrium potential parameter is not in use in this model. Therefore keep the default value of 0 V.

Locate the section. From the iloc,expr list, choose . In the associated text field, type i_a.

Locate the section. In the av text field, type Av.

Thin O2 Gas Diffusion Electrode 1

Similarly, set up the properties of the node. The details of electrode kinetics are set in the child node.

In the window, under click .

In the window for , locate the section.

In the dgde text field, type l_act.

Thin O2 Gas Diffusion Electrode Reaction 1

In the window, click .

In the window for , locate the section.

From the Eeq list, choose . Also for the oxygen reduction reaction, the equilibrium potential parameter is not in use in this model. Therefore keep the default value of 0 V.