COMSOL 6.3 - Low-Temperature PEM Fuel Cell with Serpentine Flow Field

Low-Temperature PEM Fuel Cell with Serpentine Flow Field

Water management is crucial for low-temperature polymer electrolyte fuel cell (PEMFC) operation. Water is produced on the air (cathode) side in the oxygen reduction reaction, but may permeate through the membrane to the hydrogen (anode side). Running the cell under too wet conditions may result in mass transport limitations of gases due to flooding of liquid water in the pores, whereas running the cell under too dry conditions may result in poor performance due to a low ohmic conductivity in the ionomer (polymer electrolyte) used in the membrane and catalytic layers.

This tutorial explores the current distribution in a low-temperature PEMFC when using serpentine flow field patterns, in combination with operating the cell in counter-flow mode so that the oxygen and hydrogen inlet flow streams are located at opposite sides (in the in-plane direction of the membrane) of the cell. Relatively dry inlet gas compositions are used so that the cell relies on self-humidification for achieving good performance.

For a more detailed analysis and discussion of the local transport phenomena of the fuel cell membrane–electrode assembly (MEA), see also the Transport Phenomena in a Polymer Electrolyte Fuel Cell Membrane Electrode Assembly tutorial.

Model Definition

Figure 1 shows the model geometry. The MEA is sandwiched between the anode and cathode gas diffusion layers, and the flow field channels, with the location of the gas stream inlets as indicated in the figure.

Figure 1: Model geometry. The fuel cell MEA is sandwiched between two gas diffusion layers, and the hydrogen and oxygen serpentine flow channels. The air side is located above the MEA, the hydrogen side is located below the MEA. The gas inlet positions are indicated in the figure.

The model solves for the charge balance (the electrode and electrolyte phase potentials) in the gas diffusion layers and the membrane, as well as the mass transfer (the molar fractions) and momentum transport (pressure and velocity) in the gas phase on each side of the membrane. The membrane transport of water, due to both diffusion (permeation) and migration (electroosmotic drag) is also included in the model.

The cell temperature is 70°C, with the relative humidities or the hydrogen and air inlet streams humidified to 25 and 75%, respectively. The molar flow rates of hydrogen and oxygen are set to be proportional to the total current, with a 20% excess of hydrogen and a 150% excess of oxygen (that is, using a hydrogen and oxygen flow stoichiometry of 1.2 and 2.5, respectively).

The model is defined using the Hydrogen Fuel Cell interface and solved using an Auxiliary sweep, ramping up the average cell current density from 0.01 to 1 A/cm2. See the Notes About the COMSOL Implementation and the Modeling Instructions below, and the Transport Phenomena in a Polymer Electrolyte Fuel Cell Membrane Electrode Assembly tutorial for more details on the model setup.

Results and Discussion

Figure 2 shows the cell polarization curve.

Figure 2: Polarization plot.

Figure 3 and Figure 4 show the streamlines of hydrogen and oxygen, and the corresponding molar fractions for current density of 1 A/cm2. The molar fractions decrease toward the outlets.

Figure 3: Hydrogen total flux streamlines and molar fraction.

Figure 4: Oxygen total flux streamlines and molar fractions.

Figure 5 shows the water activity (which equals the relative humidity in the gas phase) in the channels and in the membrane at 1 A/cm2. For both gas streams, the water activity increases toward the outlet. On the oxygen side the increased water activity is a direct effect of the water being produced in the cell. On the hydrogen side, the water activity increase is related mainly to the depletion of hydrogen — that is, the water fraction (at 25% relative humidity at the inlet) of the gas stream increases as hydrogen is consumed — but also to water membrane transport between the two gas compartments.

Figure 5: Left: Relative humidity in the channels. Right: Water activity in the membrane.

Figure 6 shows the electrolyte conductivity and the electrolyte current density in the z direction in the membrane. The conductivity depends on the water activity, which has an effect on the current density distribution in the cell, with generally lower current densities in the low conductivity regions. However, the current density distribution also relates strongly to the oxygen levels, lowering the current densities in areas less accessible to oxygen.

Figure 6: Left: Electrolyte conductivity in the membrane. Right: Electrolyte current density in the z direction.

Finally, Figure 7 shows the membrane water flux in the z direction. Close to the oxygen inlet/hydrogen outlet, where the hydrogen side is more humidified than the oxygen side, the flux is positive, indicating that water is transported from the hydrogen to the oxygen side. Close to the oxygen inlet/hydrogen outlet, the water is transported in the negative z direction through the membrane.

Figure 7: Water flux in the z direction of the membrane.

Notes About the COMSOL Implementation

The geometry is defined using an Assembly in the final Union node of the geometry sequence. This creates a geometry consisting of two parts, with a common boundary pair in the middle of the membrane. In this way, nonmatching meshes may be used on each side of the membrane. An automatically defined Identity Pair and a Continuity boundary node set up the needed boundary condition at the boundary pair between the assembly parts.

The gas diffusion electrodes of the MEA are not included explicitly in the geometry, and are instead defined as Thin Gas Diffusion Electrodes boundary nodes when defining the physics. This saves memory since the gas diffusion electrodes need not to be meshed.

The gas diffusion layers use anisotropic electronic conductivities, featuring about an order of magnitude higher conductivities in the in-plane (x and y) directions compared to the through-plane (z) direction).

To save computational time, Darcy’s law is used to solve for the flow in both the flow channels as well as in the gas diffusion layers. In the channels, considering the constant rectangular cross section, the permeability was derived based on an expression by Boussinesq (Ref. 1) originally derived for straight channels. See the Solid Oxide Electrolyzer tutorial for how define the flow using the Navier–Stokes and Brinkman equations instead.

The study sequence is set up in three steps: The first step solves for the potential variables for the initial gas composition, the second step solves for velocity distribution by computing the pressure variables, the third step solves the full problem, ramping up the current using an Auxiliary sweep.

The default Direct solver is disabled by enabling an iterative multigrid solver in study step 3. This reduces the memory required for solving the model.

The geometry is fully parameterized. In the modeling instructions a smaller version of the geometry is solved for first before solving for the full geometry.

Reference

1. J. Boussinesq, “Mémoire sur l'influence des Frottements dans les Mouvements Réguliers des Fluids,” J. Math. Pures Appl., vol. 13, no. 2, pp. 377–424, 1868.

Fuel_Cell_and_Electrolyzer_Module/Fuel_Cells/pemfc_serpentine_flow_field

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

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file pemfc_serpentine_flow_field_geom_sequence.mph.

In the toolbar, click

In the window, under click .

Global Definitions

Geometry Parameters

The geometry is fully parameterized. Create a smaller version of the geometry to use while setting up the physics.

In the window, under click .

In the window for , type Geometry Parameters in the text field.

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


Name	Expression	Value	Description
W_plate_min	10[mm]	0.01 m	Minimum plate width
N_ch	1	1	Number of channels
N_repeat	1	1	Number of repeating units

Geometry 1

In the toolbar, click

In the window, under click .

In the window, collapse the node.

Click the

button in the toolbar.

Global Definitions

Physics Parameters

In the toolbar, click

and choose > .

In the window for , type Physics Parameters in the text field.

Locate the section. Click

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

Change the hydrogen stoichiometry (inlet flow rate) and inlet relative humidity. This will make it easier to solve for the smaller geometry. You will go back and revert these changes later before solving for the full model.

In the table, enter the following settings:


Name	Expression	Value	Description
stoich_H2	3	3	Hydrogen flow stoichiometry
RH_an	50[%]	0.5	Inlet relative humidity, anode side

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > ®, EW 1100, Vapor Equilibrated, Protonated.

Right-click and choose .

In the toolbar, click

to close the window.

Hydrogen Fuel Cell (fc)

In the window for , locate the section.

Find the subsection. Select the checkbox.

Locate the section. Select the checkbox.

Membrane 1

In the toolbar, click

and choose .

Hydrogen Fuel Cell (fc)

In the window, expand the > node, then click > > .

In the window for , locate the section.

From the list, choose .

Geometry 1

In the window, collapse the > node.

Hydrogen Fuel Cell (fc)

In the window, under click .

In the window for , click to expand the section.

Select the checkbox.

Initial Values 1

In the window, under > > click .

In the window for , locate the section.

In the aw,0 text field, type (RH_cath+RH_an)/2.

In the T0 text field, type T_hum.

Water Absorption-Desorption, H2 Side 1

In the window, click .

In the window for , locate the section.

From the list, choose ®.

Water Absorption-Desorption, O2 Side 1

In the window, click .

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 .

Locate the section. From the list, choose .

The GDLs feature anisotropy with regards to the electron conductivity. The first, second and third values on the diagonal refer to the conductivities in the x, y, and z directions, respectively.

Specify the σs matrix as


sigmas_GDL_IP	0	0
0	sigmas_GDL_IP	0
0	0	sigmas_GDL_TP

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

In the κg text field, type kappag_GDL.

O2 Gas Diffusion Layer 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Specify the σs matrix as


sigmas_GDL_TP	0	0
0	sigmas_GDL_TP	0
0	0	sigmas_GDL_IP

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

In the κg text field, type kappag_GDL.

H2 Gas Flow Channel 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the H text field, type H_ch.

In the W text field, type W_ch.

O2 Gas Flow Channel 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the H text field, type H_ch.

In the W text field, type W_ch.

Thin H2 Gas Diffusion Electrode 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the dgde text field, type L_CL.

Thin H2 Gas Diffusion Electrode Reaction 1

In the window, click .

In the window for , locate the section.

In the i0,ref(T) text field, type i0_H2_ref.

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

Thin O2 Gas Diffusion Electrode 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the dgde text field, type L_CL.

Thin O2 Gas Diffusion Electrode Reaction 1

In the window, click .

In the window for , locate the section.

In the i0,ref(T) text field, type i0_O2_ref.

In the αa text field, type 3.

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

Electronic Conducting Phase 1

In the window, under > click .

Initial Values, O2 Domains 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Electronic Conducting Phase 1

In the window, click .

Electric Ground 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Select the checkbox.

In the Rc text field, type Rc.

Electronic Conducting Phase 1

In the window, click .

Electrode Current 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Is,total text field, type -I_tot.

Locate the section. Select the checkbox.

In the Rc text field, type Rc.

Initial Values 1

In the window, under > > click .

In the window for , locate the section.

In the x0,H2O text field, type x_H2O_an.

H2 Gas Phase 1

In the window, click .

H2 Inlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Specify the mass flow rates to make the inlet gas flow proportional to the cell current.

Locate the section. From the list, choose .

In the J0,H2O text field, type m_H2O_an.

In the ω0,bnd,H2O text field, type w_H2O_an.

Locate the section. From the list, choose .

In the J0 text field, type m_an.

H2 Gas Phase 1

In the window, click .

H2 Outlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Initial Values 1

In the window, under > > click .

In the window for , locate the section.

In the x0,H2O text field, type x_H2O_cath.

In the x0,N2 text field, type x_N2_cath.

O2 Gas Phase 1

In the window, click .

O2 Inlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the J0,H2O text field, type m_H2O_cath.

In the J0,N2 text field, type m_N2.

In the ω0,bnd,H2O text field, type w_H2O_cath.

In the ω0,bnd,N2 text field, type w_N2.

Locate the section. From the list, choose .

In the J0 text field, type m_cath.

O2 Gas Phase 1

In the window, click .

O2 Outlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Global Definitions

Default Model Inputs

In the window, under click .

In the window for , locate the section.

In the tree, select > .

Find the subsection. In the text field, type T.

Mesh 1

A user-defined mesh is required for this model. Use mapped and swept meshed in order to reduce the total number of mesh element.

With the last node in the geometry sequence set to , the mesh nodes do not need to match along the pair boundaries between the two parts of the assembly at the center of the membrane. In this way different swept meshes can be used on each side of the membrane.

Mapped 1

In the toolbar, click

and choose .

Size

In the window, click .

In the window for , click to expand the section.

In the text field, type W_ch/2.

Mapped 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Distribution 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the text field, type 4.

In the text field, type 2.

Select the checkbox.

Swept 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Size 1

Right-click and choose .

In the window for , locate the section.

Click the button.

Locate the section.

Select the checkbox. In the associated text field, type W_ch/1.1.

Click

Mapped 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click

Free Triangular 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click

Boundary Layers 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Boundary Layer Properties

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 2.

From the list, choose .

In the text field, type W_rib/10.

Click

Free Tetrahedral 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Size 1

Right-click and choose .

In the window for , locate the section.

Click the button.

Locate the section.

Select the checkbox. In the associated text field, type W_ch/2.1.

Click

Swept 2

In the toolbar, click

Distribution 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 4.

Distribution 2

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the text field, type 4.

Click

In the window, collapse the node.

Study 1

Step 1: Current Distribution Initialization

Change the current distribution type to also include activation overpotentials (secondary). This will make the initial values for the subsequent study steps more accurate.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Stationary - Pressures Only

In the window, click .

In the window for , type Stationary - Pressures Only in the text field.

Create a third study step that will solve for the whole problem, for a range of currents densities.

Stationary - All Physics

In the toolbar, click

and choose > .

In the window for , type Stationary - All Physics in the text field.

Click to expand the section. Select the checkbox.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
I_avg (Average cell current density)	I_avg_init range(I_avg_final/10, I_avg_final/10, I_avg_final)	A/cm^2

Solution 1 (sol1)

Generate the default solver and make some manual changes.

In the toolbar, click

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

The second study step should only solve for pressure variables.

In the window, expand the > > > node, then click .

In the window for , locate the section.

Clear the checkbox.

In the window, under > > > click .

In the window for , locate the section.

Clear the checkbox.

In the window, under > > > click .

In the window for , locate the section.

Clear the checkbox.

In the window, under > > > click .

In the window for , locate the section.

Clear the checkbox.

In the window, under > > > click .

In the window for , locate the section.

Clear the checkbox.

In the window, under > > > click .

In the window for , locate the section.

Clear the checkbox.

In the window, under > > > click .

In the window for , locate the section.

Clear the checkbox.

In the window, under > > > click .

In the window for , locate the section.

Clear the checkbox.

In the window, under > > > click .

In the window for , locate the section.

Clear the checkbox.

In the window, under > > > click .

In the window for , locate the section.

Clear the checkbox.

All the modified nodes are marked with a small * in the corresponding icon in the model tree. Inspect all subnodes under the node and make sure that all nonpressure variables are marked with a (*), and that all 4 pressure variables are left unmodified. (Two of the pressure variables are the pressures on each side of the membrane. The other two pressures are auxiliary variables used to implement the total mass flow rate inlet conditions.)

In the window, collapse the node.

In the toolbar, click