COMSOL 6.4 - Fuel Cell Bipolar Plate

This study presents a model that couples the thermal and the structural analysis in a bipolar plate in a proton exchange membrane fuel cell (PEMFC). The fuel cell stack consists of unit cell of anode, membrane, and cathode connected in series through bipolar plates. The bipolar plates also serve as gas distributors for hydrogen and air that is fed to the anode and cathode compartments, respectively. Figure 1 below shows a schematic drawing of a fuel cell stack. The unit cell and the surrounding bipolar plates are shown in detail in the upper right corner of the figure.

Figure 1: Schematic drawing of the fuel cell stack. The unit cell consists of an anode, membrane electrolyte, and a cathode. The unit cell is supported and supplied with reactants through the bipolar plates. The bipolar plates also serve as current collectors and feeders.

The PEMFC is one of the strongest alternatives for automotive applications where it has the potential of delivering power to a vehicle with a higher efficiency, from oil to propulsion, than the internal combustion engine.

The fuel cell operates at temperatures just below 100°C, which means that it has to be heated at startup. The heating process induces thermal stresses in the bipolar plates. This analysis reveals the magnitude and nature of these stresses.

Model Definition

Figure 2 below shows the detailed model geometry. The plate consists of gas slits that form the gas channels in the cell, holes for the tie rods that keep the stack together, and the heating elements, which are positioned in the middle of the gas feeding channel for the electrodes. Due to symmetry, it is possible to reduce the model geometry to 1/8 of the actual size of the cell.

Figure 2: The modeled geometry.

The study consists of a thermal-structural analysis. In the following specification of the problem a number of constants are used; their values are listed in the following table. The constant Q1 represents a power of 200 W distributed over the hole through the 25 bipolar plates in the fuel cell stack.


Constant	Value	Description
R0	5 mm	Radius of tie rod holes and heating element
th	10 mm	Plate thickness
Q1	200/(25·π·R0^2·th) W/m3	Volume power in one heating element
P1	1 MPa	Pressure on the active part
P2	8 MPa	Pressure on the manifold
h1	5 W/(m2·K)	Heat transfer coefficient to the surroundings
h2	50W/(m2·K)	Heat transfer coefficient in the gas channels

Thermal Analysis

The general equation is the steady-state heat equation according to:

where k denotes the thermal conductivity of the different materials, T the temperature, and Q a heat source or heat sink.

The modeled domain consists of three different domains. The first domain corresponds to the active part of the cell, where the electric energy is produced and the current is conducted. The production and conduction of current involve some losses. It is assumed that the cell is heated prior to operation. This means that the model does not account for the heat sources due to the production and conduction of current; see Figure 3.

Figure 3: The active part of the cell bipolar plate is made of titanium.

The second domain corresponds to the heating element in the bipolar plate, which is made of aluminum. The power from this element is assumed to be uniformly distributed in the small cylindrical domain; see Figure 4.

Figure 4: The heating element in the bipolar plate is made of aluminum.

The last domain forms the manifold in the cell and it is made of titanium. In this domain, only heat conduction is present.

Figure 5: The domain that forms the manifold of the cell.

There are different material models present for both titanium and aluminum.

The boundary conditions for the heat transfer analysis are symmetry conditions and convective conditions. The symmetry condition states that the flux is zero through the boundary.

The convective conditions set the heat flux proportional to the temperature difference between the fluid outside and the temperature at the boundary. The heat transfer coefficient is the proportionality constant:

In this equation, h denotes the heat transfer coefficient.

Figure 6 shows the symmetry and insulation boundaries.

Figure 6: Symmetry and insulation boundaries.

Figure 7 shows the convective boundaries. The heat transfer coefficient is different for the different groups of boundaries in the figure.

Figure 7: The convective boundaries.

Structural Analysis

The thermal expansion causes the stresses in the domain. The thermal loads are proportional to the temperature difference between the reference state and the actual temperature.

The constrains are imposed on the symmetry planes, see Figure 8 and Figure 9. There the displacements perpendicular to the boundary are set to zero.

Figure 8: The symmetry boundaries.

Figure 9: Additional symmetry plane.

In addition, the load applied by the tie rods on the stack is applied in the manifold and the active part at the bipolar plate. The pressure is different in these two domains, see Figure 10.

The loads and constraints on all other boundaries are assumed to be negligible. In this context, the only questionable assumption is the possible loads and constraints imposed by the tie rods on the holes in the bipolar plate. In a more detailed analysis, it would be possible to couple a model for the tie rods with a detailed model of the bipolar plate. The probably the pressure on those different parts would show different gradients.

Figure 10: Pressure due to the rod connection.

Results

Figure 11 shows the temperature field at steady state. The slits for the gas present the largest sinks of heat at the boundaries. The temperature decreases almost with the radial distance from the heating element.

Figure 11: Temperature distribution in the plate.

Figure 12 shows that the thermal loads generate stresses that are one order of magnitude larger than those generated by the external pressure loads (which you can compute by turning off the thermal expansion). The loads are far from the critical values, but the displacements can be even more important. The displacements must be small enough to grant that the membrane can fulfill its function in separating hydrogen and oxygen in the anode and cathode compartments, respectively. The z-component of the displacement in shown in Figure 13.

Figure 12: Distribution of the von Mises stress due to temperature and external loads.

Figure 13: Expansion in the z direction.

Structural_Mechanics_Module/Thermal-Structure_Interaction/bipolar_plate

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

The geometry sequence for the model is available in a file. If you want to create it from scratch yourself, you can follow the instructions in the Appendix — Geometry Modeling Instructions section. Otherwise, insert the geometry sequence as follows:

In the toolbar, click and choose .

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

In the toolbar, click

Click the

button in the toolbar.

Click the

button in the toolbar.

In the window, under click .

Definitions

Create an selection of the boundaries for features in both physics.

Symmetry

In the toolbar, click

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

Locate the section. From the list, choose .

Select Boundaries 3, 5, 20, 23, 25, 37, and 43–45 only.

It might be easier to select the correct boundaries by using the window. To open this window, in the toolbar click and choose . (If you are running the cross-platform desktop, you find in the main menu.)

Solid Mechanics (solid)

In the window, under click .

Select Domains 1 and 2 only.

Materials

Titanium

In the window, under right-click and choose .

Select Domains 1 and 2 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	105[GPa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.33	1	Young’s modulus and Poisson’s ratio
Density	rho	4940[kg/m^3]	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	7.5[W/(m*K)]	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	710[J/(kg*K)]	J/(kg·K)	Basic
Coefficient of thermal expansion	alpha_iso ; alphaii = alpha_iso, alphaij = 0	7.06e-6[1/K]	1/K	Basic

In the text field, type Titanium.

Aluminum

Right-click and choose .

Select Domain 3 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	160[W/(m*K)]	W/(m·K)	Basic
Density	rho	2700[kg/m^3]	kg/m³	Basic
Heat capacity at constant pressure	Cp	900[J/(kg*K)]	J/(kg·K)	Basic