Proton exchange membrane fuel cells (PEMFCs) require mechanical compression to ensure structural integrity, prevent leakage, and to minimize the electrical contact resistance. The mechanical properties and dimensions of the fuel cell vary during assembly due to manufacturing tolerances and during operation due to both temperature and humidity. Variation in stack compression affects the interfacial contact pressures between components and hence fuel cell performance. This paper presents a one-dimensional equivalent stiffness model of a PEMFC stack capable of predicting independent membrane and gasket contact pressures for an applied external load. The model accounts for nonlinear component compression behavior, thickness variation due to manufacturing tolerances, thermal expansion, membrane expansion due to water uptake, and stack dimensional change due to clamping mechanism stiffness. The equivalent stiffness model is compared to a three-dimensional (3D) finite element model, showing good agreement for multicell stacks. Results demonstrate that the correct specification of gasket thickness and stiffness is essential in ensuring a predictable membrane contact pressure, adequate sealing, and avoiding excessive stresses in the bi-polar plate (BPP). Increase in membrane contact pressure due to membrane water uptake is shown to be significantly greater than the increase due to component thermal expansion in the PEMFC operating range. The predicted increase in membrane contact pressure due to thermal and hydration effects is 18% for a stack containing fully hydrated Nafion® 117 membranes at 80 °C, 90% relative humidity (RH) using an eight bolt clamping design and a nominal 1.2 MPa assembly pressure.

## Introduction

Proton exchange membrane fuel cells (PEMFCs) have seen increased use in recent years as an alternative clean energy source in multiple applications ranging from consumer electronics to automotive powertrains. To obtain a usable output power and voltage, multiple cells are combined in series to produce a fuel cell stack. The design of the fuel cell stack is critical in maximizing performance and efficiency; the method of clamping the individual stack components together is one such area where careful design consideration is required. The purpose of the clamping load in a fuel cell stack is to ensure adequate gas tight sealing, reduce electrical and thermal contact resistance between the different components, and provide some structural integrity [1,2]. However, overcompression can lead to reduced performance and damage [3].

The vast majority of fuel cell stacks are clamped together using multiple bolts to apply a compressive load between two stiff endplates. Wen et al. [4] and Lee et al. [3] used polarization curves to demonstrate how different clamping pressures, applied through varying bolt torque, influenced the cell performance. Wen et al. [4] compared two, four, and six bolt stack designs at different torques. The six-bolt design at the maximum tested torque of 16 N·m gave the best overall performance in a ten-cell stack. Lee et al. [3] used a 10 cm2 single cell design to observe the cell performance at different clamping loads with different types of gas diffusion layer (GDL). An optimum tightening torque was observed with performance decreasing at high bolt torques, most likely due to overcompression of the GDL inhibiting the transport of reactant gases and product water.

The main cause of the improved performance with increased bolt torque seen by Wen et al. [4] and Lee et al. [3] is due to additional clamping pressure reducing the interfacial contact resistance between the fibers of the gas diffusion layer and bi-polar plate (BPP). Mason et al. [5] showed that increasing the clamping pressure on a single cell with Toray H120 carbon paper GDL from 0.25 to 2.50 MPa reduced the electrical resistance from 29 to 10 mΩ/cm2. Lin et al. [6] determined that the optimum compression for two types of carbon cloth GDL was between 59% and 64%. Meanwhile, Zhou et al. [7] also characterized the relationship between interfacial contact resistance and contact pressure using a microscale numerical model, showing good agreement with experimental results.

In PEMFC stacks with a large cell active area ($≥$ 100 cm2), it is important to consider the distribution of contact pressure across the surface of the cell created by the clamping method used. Multiple studies have investigated how uneven clamping pressure from endplate deflection influences cell performance. Bates et al. [8] used a combination of finite element analysis (FEA) simulations and experimental tests with pressure sensitive films, characterizing the pressure distribution in a 16-cell stack with 100 cm2 active area cells. Results demonstrated reduced clamping pressure at the center of the cell caused by deflection of the stainless steel endplates. Carral and coworkers [9,10] also considered how the number of cells in a stack influences the contact pressure distribution, demonstrating both numerically [9] and experimentally [10] how increasing the number of cells from 1 to 16 reduces the variation in contact pressure. Strain gauges and optical methods were used by Carral et al. [10] to measure endplate deflection and compressive load in the tie rods. Montanini et al. [11] used a piezoresistive sensor array to perform real-time measurements of the contact pressure distribution, demonstrating pressure distribution between the membrane electrolyte assembly (MEA) and gasket at different bolt torques. Gatto et al. [12] used the same method as Montanini et al. [11], along with electrochemical testing, to investigate the bolt torque for optimum performance of four gaskets across a range of stiffness and thickness.

When determining an appropriate clamping load, it is important to consider both the assembly and operational conditions that the fuel cell stack will experience. Increased temperature and relative humidity (RH) during operation will cause different rates of hygrothermal expansion and a change in material properties, influencing the contact pressure between components. For the fuel cell membrane, these effects can be significant. Tang et al. [13] demonstrated between 25 °C 30% RH and 85 °C 90% RH, a 10% swelling and 77% drop in the Young's modulus of a Nafion® 112 membrane. Zhou et al. [14] conducted one of the few studies to consider the influence of both temperature and hydration on clamping. Their work used a two-dimensional finite element model of a single BPP channel and land to investigate membrane and GDL compression under different assembly pressures and operating conditions. The deformed model geometry under assembly conditions was then used in a computational fluid dynamics package to obtain cell polarization curves.

Current computational limitations mean that it is not feasible to investigate the clamping performance of large multicell fuel cell stacks (>20 cells) using commercial FEA packages without making simplifying assumptions. One simplifying method, put forward by Lin et al. [15], is to represent each component of the fuel cell stack as a single spring, with a constant stiffness equivalent to the solid bulk material. Using this concept, a one-dimensional model of the stack can be produced to determine stack and component compression for a chosen clamping load, as well as the force distribution between the gasket and MEA. By applying thermal constraints postassembly, Lin et al. [16] calculated the optimum bolt torque range to give an ideal cell clamping pressure at both assembly and operational temperatures. The influence of nonlinear compressive behavior and humidity was not considered.

The equivalent stiffness concept has since been used by several other authors: Liu et al. [17] produced an equivalent stiffness model of a fuel cell stack using a steel belt clamping method to predict behavior under external accelerations. Qiu et al. [18] represented channel height variation of a stamped metallic bi-polar plate using an equivalent stiffness model. This showed good agreement with experimental results obtained using pressure sensitive film. The equivalent stiffness concept was also used by Ahmad et al. [19] to determine the required stack clamping force for a desired GDL compression when developing an assembly method for a 60-cell stack.

The literature has shown that clamping load and associated contact pressure have a significant effect on the performance of PEMFC stacks. Correct specification of the gasket and GDL properties are important in ensuring appropriate distribution of clamping pressures during both assembly and operation. While the influence of operational temperature on clamping pressure has been explored previously, the influence of membrane water content, nonlinear spring stiffness, and manufacturing tolerances is not yet fully understood at a stack level. In this paper, a new equivalent stiffness model of a fuel cell stack is presented. The model is the first stack level model to include the effects of nonlinear component compression, membrane hydration, and manufacturing tolerances. Sections 2 and 3 detail the model and its comparison to three-dimensional (3D) FEA simulations and literature data, respectively, Sec. 4 discusses the results of the model.

## Equivalent Stiffness Model

The equivalent stiffness model used is this work is based on Hooke's law (Eq. (1)). The force Fi required to compress each of the separate components within the fuel cell stack i, by a displacement xi, can be found using an individual stiffness value $ki(x)$ based on its mechanical properties and dimensions. Using this method, a simple model can be produced, capable of determining the compression and contact pressures of the different components within the fuel cell stack for a known overall clamping force
$Fi=ki(x)xi$
(1)

The compressive behaviors of five components for each cell in the fuel cell stack are considered in the model. These are the proton exchange membrane (PEM), GDL, gasket, BPP, and clamping mechanism. Each of the stack components is represented by a single spring, which is then assembled to construct a representative fuel cell stack model. The relationship between the material properties and equivalent spring stiffness $ki(x)$ is shown in Eq. (2) and illustrated in Fig. 1. Where $Ei(x)$ is Young's Modulus, Ai is material cross-sectional area, and $ti(x)$ is material thickness. The model makes the following assumptions:

Fig. 1
Fig. 1
Close modal
1. (1)

All components are homogeneous.

2. (2)

Clamping load is applied evenly (infinitely stiff endplates).

3. (3)

Influence of coolant channels on bi-polar plate stiffness is ignored.

4. (4)
The PEM under the gasket area is independent of the PEM under the gas diffusion layer
$ki(x)=Ei(x)Aiti(x)$
(2)

At low levels of compression (0.01–1.00 MPa), the stress–strain relationship of the gasket and GDL is nonlinear. The Young's modulus increases with strain, approaching a linear relationship in the 1–3 MPa regions [20]. The different relative stiffnesses between components in the fuel cell stack mean that the initial nonlinear region can have a significant influence on the contact pressure distribution, particularly when small dimensional tolerances are considered. The model accounts for nonlinear stiffness in the gasket and GDL using a stiffness look-up table as a function of component displacement x from experimental data. This is applied using a scaling factor from 0 to 1 multiplied by the maximum stiffness. Experimental compression data for an SGL10BA GDL and Samco silicone gasket taken from Ismail et al. [20], shown in Fig. 2, are used in this work. All other components are assumed to be linear elastic.

Fig. 2
Fig. 2
Close modal

Figure 3 shows the arrangement of the different fuel cell components considered in the model. The central cell unit is repeated for the number of cells in the fuel cell stack. By representing each component, as a single spring with stiffness $k(x)$, the structure can be simplified to that of Fig. 4 where the repeating cell unit is enlarged for clarity.

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal

Between each BPP, the clamping force must travel through either the gasket assembly or the MEA, which are in parallel. The gasket assembly consists of two gaskets and a PEM in series (Gasket/PEM/Gasket). The MEA consists of two gas diffusion layers and a PEM in series (GDL/PEM/GDL).

The stiffness of the MEA and gasket assembly are calculated as a function of GDL displacement (xgdl) and gasket displacement (xg), while accounting for nonlinear stiffness behavior. The MEA stiffness (kmea) is
$kmea(xgdl)=(kpem,ckgdl(xgdl)2kpem,c+kgdl(xgdl))$
(3)
The gasket assembly stiffness (kgass) is
$kg,ass(xg)=(kpem,gkg(xg)2kpem,g+kg(xg))$
(4)
where $kpem,c$ is the PEM stiffness in the MEA area and $kpem,g$ is the PEM stiffness in the gasket area. The force through the MEA area (Fmea) is determined by
$Fmea=kmea(xgdl)xmea=kpem,ckmea(xgdl)2xgdlkpem,c−kmea(xgdl)$
(5)
where GDL displacement is obtained through
$xmea=2xgdl+xpem,c$
(6)
Using the same method as Eq. (5), force through the gasket assembly ($Fg,ass$) in terms of gasket displacement is obtained as
$Fg,ass=kg,ass(xg)xg,ass=kpem,gkg,ass(xg)2xgkpem,g−kg,ass(xg)$
(7)
As illustrated in Fig. 4, the total equivalent stiffness of a single cell ($kcell,n$) can then be found as
$kcell,n(xgdl,xg)=(kpem,ckgdl(xgdl)2kpem,c+kgdl(xgdl))+(kpem,gkg(xg)2kpem,g+kg(xg))$
(8)

The area of the PEM under the gasket is different from the area of the PEM under the GDL, meaning two different PEM stiffness values are required in Eqs. (3)(8). $kpem,g$ is used for the area of PEM under the gasket and $kpem,c$ for the area of PEM under the GDL. This method implies that the displacement of the PEM between the gasket assembly ($xg,ass$) is independent of the PEM displacement between the GDLs (xmea), and does not consider that the PEM is connected. This assumption is deemed reasonable since the PEM exhibits good flexibility and the relative displacement in each cell will be minimal, due to the high stiffness of the bi-polar plates. The calculation of all individual component stiffness values used in the model is shown in Table 1.

Table 1

Component area and stiffness calculations

ParameterEquationValue
$kbpp,main$$EbppAbpp,maintbpp,main$$6.05×1010$ N/m
$kbpp,edge$$EbppAbpp,edgedrib$$2.10×1010$ N/m
$kbpp,ribs$$EbppAbpp,ribdrib$$4.60×1010$ N/m
$kpem,cell(T,λ)$$Epem(T,λ)Abpp,ribtpem$$2.15×1010$ N/m (25 °C 30% RH)
$kpem,g(T,λ)$$Epem(T,λ)Agtpem$$6.69×109$ N/m (25 °C 30% RH)
$kgdl(xgdl)$$Egdl(xgdl)Acelltgdl$$1.88×108$ N/m (38% strain)
kclamp$EclampAclampLstack$$2.08×108$ (8 × 8 mm dia steel bolts)
$kg(xg)$$Eg(xg)Agtg$$4.16×108$ N/m (17% strain)
$Abpp,main$$hbppwbpp$0.0121 m2
$Abpp,edge$$hbppwbpp−hcellwcell$0.0021 m2
$Abpp,rib$$wribhribnrib$0.0046 m2
Aclamp$πrclamp2nclamp$$4.02×10−4$ m2 (8 × 8 mm dia)
ParameterEquationValue
$kbpp,main$$EbppAbpp,maintbpp,main$$6.05×1010$ N/m
$kbpp,edge$$EbppAbpp,edgedrib$$2.10×1010$ N/m
$kbpp,ribs$$EbppAbpp,ribdrib$$4.60×1010$ N/m
$kpem,cell(T,λ)$$Epem(T,λ)Abpp,ribtpem$$2.15×1010$ N/m (25 °C 30% RH)
$kpem,g(T,λ)$$Epem(T,λ)Agtpem$$6.69×109$ N/m (25 °C 30% RH)
$kgdl(xgdl)$$Egdl(xgdl)Acelltgdl$$1.88×108$ N/m (38% strain)
kclamp$EclampAclampLstack$$2.08×108$ (8 × 8 mm dia steel bolts)
$kg(xg)$$Eg(xg)Agtg$$4.16×108$ N/m (17% strain)
$Abpp,main$$hbppwbpp$0.0121 m2
$Abpp,edge$$hbppwbpp−hcellwcell$0.0021 m2
$Abpp,rib$$wribhribnrib$0.0046 m2
Aclamp$πrclamp2nclamp$$4.02×10−4$ m2 (8 × 8 mm dia)

Both the gasket and MEA areas contact the bi-polar plates between cells. The flow channels in a bi-polar plate mean it is not a homogeneous block of material and therefore the method of Lin et al. [15] is used to determine the appropriate stiffness. The bi-polar plate is separated into three sections, the main bulk section, the rib section, and the edge section (which meets the gasket), shown in Fig. 5. Each section has a different area and therefore stiffness value. The edge ($kbpp,edge$) and rib area ($kbpp,ribs$) act in parallel with the bulk section ($kbpp,main$) in series; overall stiffness of the bi-polar plate is therefore found using Eq. (9). The method for calculating individual component stiffness is shown in Table 1.

Fig. 5
Fig. 5
Close modal
The presence of flow channels on the bi-polar plate reduces the MEA contact area, which influences the load distribution between the gasket and MEA area. The contact area between the BPP and the GDL is equal to the surface area of the ribs, determined in Table 1. It is assumed that the contact area between the GDL and PEM is equal to the contact area between the BPP and GDL. This assumption is reasonable since both the thickness and stiffness of the gas diffusion layers are significantly lower than the BPP. For the geometry used in this study, the cell active area is 100 cm2 and the contact area is 46 cm2 after considering the bi-polar plate flow channels
$kbpp=(1kbpp,ribs+kbpp,edge+1kbpp,main)−1$
(9)
Using the above methods, the overall force (Fstack) required to displace the fuel cell stack by a distance xstack can be determined using Eq. (10), where ncells is the number of cells and nbpp is the number of bi-polar plates in the stack
$Fstack=[∑n=1ncells1kcell,n(xgdl,xg)+∑n=1nbpp1kbpp,n]−1xstack$
(10)

### Unequal Thickness.

It is unlikely that the gasket and GDL will have the same thickness, especially when manufacturing tolerances are taken into account. During assembly, the gasket and GDL will contact with the BPP at different compression levels, meaning one component is precompressed before the other comes into contact. The equivalent stiffness model accounts for precompression using Eq. (11) to establish the stiffness of each component as a function of displacement $k(x)$ where k is the standard equivalent stiffness found using Eqs. (2), (3), or (4), and δ is the difference between the initial thickness of the MEA and gasket assembly. When not in contact with the BPP, the stiffness value is set to 0, meaning no force can be transmitted through the MEA until the gasket is compressed to the point where contact occurs. The opposite is used for the case when the MEA touches before the gasket
$kmea(x)={kmeax−δ≥00x−δ<0$
(11)

### Hygrothermal Effects.

During operation, the temperature and humidity of the fuel cell stack will increase, changing the mechanical properties of the components within the stack. The equivalent stiffness model considers thermal expansion of the PEM, GDL, gasket, and BPP and expansion of the PEM due to water uptake. Thermal expansion in the model is considered using Eq. (12) to determine the material thickness under operating conditions, where α is the thermal expansion coefficient and T temperature. The new thickness is used to determine the change in overall stack compression
$toperation=tassembly(1+α(Toperation−Tassembly))$
(12)
Expansion of the membrane due to water uptake has been modeled using the swelling coefficient β put forward by Zhou et al. [14]. A constant value of $β=0.0115 λ−1$ has been used in this study based on experimental data for the dimensional change of a Nafion® 112 membrane from Ref. [13]. Where λ is the water content of the membrane expressed in number of water molecules per charge site $SO3−$$H+$ (see Eq. (14)). Expansion of the membrane due to swelling is calculated using the below equation:
$tpem,operation=tpem,assembly(1+β(λpem,operation−λassembly))$
(13)
Water uptake and temperature increase during operation lead to a significant decrease in the strength of the membrane. A surface fit using data from Ref. [14] is used to calculate the membrane Young's modulus (Epem) from water content and temperature, shown in Eq. (15) (R2 = 0.97). The relationship of Springer et al. [21] is used to calculate water content as a function of water activity (aw). The Young's modulus of Eq. (15) is used to evaluate the membrane stiffness kpem at both assembly and operational conditions
$λ=0.0043+17.81aw−39.85aw2+36aw3$
(14)
$Epem=363.8−46.08λ−2.60T+2.61λ2+0.03λT+0.01T2$
(15)

### Manufacturing Tolerances.

Variation in component thickness due to manufacturing tolerances will influence the load distribution between the gasket and MEA of each cell, partly contributing to cell-to-cell performance variation seen in large fuel cell stacks. The model accounts for manufacturing tolerances of the gasket, GDL, and PEM thickness by assigning a normally distributed random number for each thickness. For each component, the mean distribution value is equal to the nominal thickness and the standard deviation specified to give the desired tolerance range. In this work, the standard deviation is set to 0.5% of the nominal thickness value. For a 200 μm thick GDL, this gives a 1 μm standard deviation. A single cell assembly consists of two gaskets, two GDLs, and one PEM. Each of the five components in each cell of the stack is assigned a different thickness, meaning every individual component within the stack has a unique stiffness value and each cell has a different contact spacing δ used in Eq. (11).

### Clamping Method.

When hygrothermal effects are applied to the stack, the distribution of pressure and overall forces will vary depending on the method used in keeping the stack clamped together. In this work, three different clamping methods are considered, these are:

1. (1)

Fixed displacement (stack length is fixed and cannot expand)

2. (2)

Fixed clamping force (stack is free to change in length but clamping force is fixed)

3. (3)

Variable displacement and force (stack can change both length and force depending on stiffness of clamping method)

For the fixed displacement case, any expansion caused by operational conditions is taken up in the form of additional compression of the fuel cell stack components. Additional compression and hence change in force distribution is calculated for each separate component within the stack. For the fixed clamping force case, overall clamping pressure remains constant and the stack is free to expand or contract. The distribution of clamping force between the gasket and GDL of each cell is recalculated due to differing expansion coefficients. For the variable displacement and force clamping case, the equivalent stiffness of the clamping mechanism must be taken into account. As with the fuel cell stack, this is achieved using Eqs. (1) and (2). Equation (16) gives the equivalent stiffness for a clamping design with nbolts bolts of diameter d and stack length Lstack
$kclamp=Eclampπd2nbolt4Lstack$
(16)

The change in stack length for clamping stiffness kclamp is found by calculating the force required to displace the stack by a distance x in both assembly and operational conditions. The point x = 0 corresponds to the initial displacement of the stack before compression under assembly conditions. For assembly compression xa with force Fa, the compression under operational conditions is the intersection of the operational compression vector and a vector with gradient $−kclamp$ starting at point (xa, Fa). This method, illustrated in Fig. 6, accounts for changes in component dimensions, spacing, and nonlinear stiffness between the two conditions.

Fig. 6
Fig. 6
Close modal

The model described in Sec. 2 was implemented using matlab and populated with the geometry and material properties of a typical fuel cell stack with 100 cells each with a 100 cm2 cell active area. Values used are shown in Table 2 unless otherwise stated.

Table 2

Fuel cell stack parameters

ParameterValue
Stack
Number of cells100
Assembly temperature$25 °C$
Assembly humidity30%
Bi-polar plate
Width0.11 m
Height0.11 m
Rib area height0.1 m
Rib area width0.1 m
Rib width1 mm
Number ribs46
Thickness3 mm
Rib depth1 mm
Young's modulus10 GPa [14]
Poisson's ratio0.25 [14]
Thermal expansion7.9 × 10−6 K−1 [15]
GDL
Width0.1 m
Height0.1 m
Thickness200 μm
Young's modulusVariable [20] (15 MPa at 38% strain)
Poisson's ratio0.25 [14]
Thermal expansion$7.9×10−6$ K−1 [15]
PEM
Width0.1 m
Height0.1 m
Thickness50 μm
Young's modulusVariable [14]
Poisson's ratio0.4 [22]
Thermal expansion$90×10−6$ K−1 [15]
Hydration expansion0.0115 $λ−1$
Width (outer)0.108 m
Width (inner)0.1 m
Height (outer)0.108 m
Height (inner)0.1 m
Thickness200 μm
Young's modulusVariable [20] (50 MPa at 17% strain)
Poisson's ratio0.3 [22]
Thermal expansion$77×10−6$ K−1 [15]
Clamp
Bolt diameter8 mm
Number of bolts8
Young's modulus180 GPa
Poisson's ratio0.3
Thermal expansion$12×10−6$ K−1 [15]
Endplate (FEA only)
Width0.144 m
Width0.144 m
Thickness0.02 m
Young's modulus69 GPa
Poisson's ratio0.3
Thermal expansion$77×10−6$ K−1 [15]
ParameterValue
Stack
Number of cells100
Assembly temperature$25 °C$
Assembly humidity30%
Bi-polar plate
Width0.11 m
Height0.11 m
Rib area height0.1 m
Rib area width0.1 m
Rib width1 mm
Number ribs46
Thickness3 mm
Rib depth1 mm
Young's modulus10 GPa [14]
Poisson's ratio0.25 [14]
Thermal expansion7.9 × 10−6 K−1 [15]
GDL
Width0.1 m
Height0.1 m
Thickness200 μm
Young's modulusVariable [20] (15 MPa at 38% strain)
Poisson's ratio0.25 [14]
Thermal expansion$7.9×10−6$ K−1 [15]
PEM
Width0.1 m
Height0.1 m
Thickness50 μm
Young's modulusVariable [14]
Poisson's ratio0.4 [22]
Thermal expansion$90×10−6$ K−1 [15]
Hydration expansion0.0115 $λ−1$
Width (outer)0.108 m
Width (inner)0.1 m
Height (outer)0.108 m
Height (inner)0.1 m
Thickness200 μm
Young's modulusVariable [20] (50 MPa at 17% strain)
Poisson's ratio0.3 [22]
Thermal expansion$77×10−6$ K−1 [15]
Clamp
Bolt diameter8 mm
Number of bolts8
Young's modulus180 GPa
Poisson's ratio0.3
Thermal expansion$12×10−6$ K−1 [15]
Endplate (FEA only)
Width0.144 m
Width0.144 m
Thickness0.02 m
Young's modulus69 GPa
Poisson's ratio0.3
Thermal expansion$77×10−6$ K−1 [15]

## Comparison to Three-Dimensional Finite Element Model

The equivalent stiffness model was compared to a 3D finite element analysis (3D FEA) model to (a) determine if the assumptions made with respect to even clamping load are reasonable, and (b) demonstrate how the equivalent stiffness model can be used as an alternative to the 3D FEA with a low computational cost.

### Three-Dimensional Finite Element Model.

A 3D model of the fuel cell stack detailed in Table 2 was produced using the commercial FEA package abaqus with the addition of 20 mm thick aluminum endplates. The model considers the influence of the endplates, clamping bolts, gasket, GDL, PEM, and BPP including flow channels. The current collector, reactant manifolds, or coolant channels in the BPP are not accounted for, as these will not significantly affect the contact pressure distribution of the MEA. By applying symmetry in the x, y, and z planes, it is possible to perform calculations on 1/8th of the geometry only, significantly reducing the computational time. The z plane symmetry is applied to the middle cell at half of the PEM thickness. Figure 7 shows the simulation geometry for a five-cell stack. The GDL is fixed to the PEM to represent the MEA structure after hot pressing. Contact between all other components is accounted for using the general contact model with a 0.3 tangential coefficient of friction and exponential overclosure contact in the normal direction. Variation of the friction coefficient was seen to have minimal (<1%) effect on the cell contact pressure distribution. Nonlinear behavior of the GDL and gasket was modeled using the hyperelastic material model based on experimental uniaxial compression data [20]. Different mesh seeds were used to account for the variation in geometry between components. The coarsest mesh of 3.0 mm was used for the endplate and the finest mesh seed of 0.5 mm was applied to the PEM and GDL. All components were meshed using hexahedral elements with reduced integration (C3D8R). The number of cells in the fuel cell stack was varied from one to nine, the five-cell stack containing 136,157 elements. Eight bolts of 8 mm diameter were used to apply a compressive load of 1.25 kN per bolt. Figure 8 shows the MEA contact pressure distribution for 1/4 of the end cell in a five-cell stack, and the influence of BPP channel geometry on localized contact pressure can be clearly seen.

Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal

### Comparison.

Contact pressures of the cell and gasket in the 3D FEA model were compared to those predicted by the equivalent stiffness model with the same overall clamping force applied. Figure 9 compares the maximum and minimum contact pressures of the MEA and gasket from the end cell of the 3D model, to those predicted by the equivalent stiffness model. For a stack with only a few cells, the variation in contact pressure, especially in the gasket, is high and the assumption of uniform clamping pressure is not valid. However, as the number of cells in the stack increases, the contact pressure becomes more uniform as the additional cells act to distribute the load. For multiple cell stacks, the equivalent stiffness model shows good agreement with the 3D FEA model, demonstrating that the model is suitable for use as a predictive design tool for fuel cell stack development with a low computational cost.

Fig. 9
Fig. 9
Close modal

### Comparison to Literature Data.

The 3D FEA model of Caral and Mélé [9] also compared the contact pressure distribution as a function of the number of cells in a stack using a validated endplate model and nonlinear MEA stress–strain curve. Their model predicted contact pressure variation of 98%, 87%, and 72% across the surface of the MEA for the end cell of a 2-cell, 5-cell, and 16-cell stacks, respectively. The FEA model in the current work predicted a lower contact pressure variation of 80%, 40%, and 24%, respectively, for a one-cell, five-cell, and nine-cell stacks with different geometry. The work of Caral and Mélé [9] demonstrates a higher level of variation due to the use of a two-bolt endplate design compared to an eight-bolt design in the current work. The authors also did not consider the influence of the sealing gasket in their study, which will help to disperse load concentrations from the edge of the MEA caused by bending of the endplate under load.

Mason et al. [23] used a controlled compression rig to measure the change in clamping pressure of a Nafion® 117 membrane when transitioning from a dry to wet state by exposing the cell to a humidified nitrogen gas feed at 80 °C. After an initial compression assembly of 0.2 MPa, the displacement was fixed while humidified gas was supplied. MEA compression pressure was then seen to increase by 0.55 MPa, from 0.55 to 1.10 MPa due to hydration, stabilizing after 150 s. Replicating the same initial compression, clamping method, and membrane thickness in the equivalent stiffness model, the MEA contact pressure was seen to increase by 0.45 MPa when 80 °C, 95% RH conditions were applied, showing good agreement to the experimental work of Mason et al. [23].

## Results

### Assembly.

Fig. 10
Fig. 10
Close modal

Selection of a thicker and softer gasket can be used to ensure a gas tight seal is formed prior to contact with the MEA. Figures 10(c) and 10(d) show the compression and contact pressure for a 250 μm gasket with a 50% lower Young's modulus. Contact between the bi-polar plate and MEA does not occur until a stack clamping pressure of 0.9 MPa has been applied and the gasket has been compressed by 20%. Similarly, Figs. 10(e) and 10(f) represent the case where the gasket is thinner and stiffer than the GDL ($t=150 μm$, 200% stiffness). This method can be used to achieve a desired compression of the GDL or achieve equal contact pressure of the MEA and gasket at a specific stack clamping pressure.

To demonstrate how the inclusion of nonlinear stiffness influences stack behavior, compression and clamping pressure were also calculated using constant stiffness values for the GDL and gasket. Constant stiffness values at 38% and 13% strain were used for the gasket and GDL, respectively. Results for linear compression are shown as dashed lines in Fig. 10, denoted by (L). The linear model is seen to underpredict GDL and gasket compression in all cases compared to the nonlinear model. This occurs since the variable stiffness is lower than the constant stiffness value at low levels of compression. The increased stiffness at low compression also leads to the linear model overestimating the clamping pressure required to make contact with both the GDL and gasket, compared to the nonlinear model. Disregarding nonlinear compression effects may therefore lead to the underprediction of component compression and overprediction of clamping pressures required for good MEA and gasket contact. Sections 4.2 and 4.4 all use the nonlinear model.

### Operational Conditions.

Changes in temperature and humidity under operational conditions will influence the clamping pressure experienced by the MEA due to thermal expansion of the stack components, water uptake, and stiffness changes of the membrane. Figure 11 shows how different operational conditions influence MEA contact pressure using the equivalent stiffness model populated with parameters from Table 2. Assembly conditions are 1.2 MPa MEA pressure at 20 °C 30% RH using 8 × 8 mm diameter bolts. The MEA contact pressure can be seen to increase in a near linear relationship with increasing temperature. This indicates that the contribution from thermal expansion of the stack components is greater than reduction in PEM stiffness with temperature. The change in contact pressure due to humidity variation is seen to be much greater than the change due to temperature, the membrane swelling effects being greater than the thermal expansion. At 80 °C and 90% RH, the MEA contact pressure increase was 0.21 MPa, 17.9% relative to the assembly condition.

Fig. 11
Fig. 11
Close modal

The MEA contact pressure change under operational conditions will vary with stack compression during assembly. This is shown in Fig. 12 for three different operating conditions. At low initial compression, the change in MEA contact pressure under operational conditions is low. This is because the reduced stiffness of the GDL, and gasket at low strain accommodates component expansion. A peak in MEA pressure change under operational conditions is seen around 2.7% compression. At this point, the GDL and gasket stiffness are in the linear region and the thermal/hydration expansion of the stack components offsets the reduction in PEM stiffness. At high levels of stack compression, the change in MEA pressure is seen to reduce. This occurs as the higher initial clamping force means the PEM is further compressed as its stiffness reduces under operational conditions, reducing bolt stretching and offsetting the increase in load through thermal expansion.

Fig. 12
Fig. 12
Close modal

### Clamping Method.

The previous results have used 8 × 8 mm diameter bolts. Clamping mechanism stiffness will directly influence the MEA contact pressure under assembly conditions as it restricts the stacks ability to expand. The change in MEA pressure when operating at 90% RH using four different clamping mechanisms is shown in Fig. 13. Assembly conditions are 2.3% original compression (1.2 MPa MEA contact pressure) at 25 °C 30% RH. The case where the stack is fixed and unable to expand gives the largest increase in MEA clamping pressure below 80 °C since all hygrothermal expansions must lead to a change in contact pressure. The fixed case is very similar to a conventional bolted design using 4 × 4 mm diameter bolts. At low operating temperatures, the change in MEA clamping pressure due to operational conditions can be reduced by using a clamping method with a low stiffness kclamp. Potential methods include replacing the bolts with multiple thinner wires, as in the 8 × 2 mm diameter case, or by placing constant stiffness “tie rod springs” between the stack endplate and nut used to clamp the stack together [2]. Above 80 °C operating temperature, a reduction in clamping stiffness is seen to lead to an increased change in MEA clamping pressure. At this temperature, the reduction in stiffness of the PEM contributes to stack shrinkage and additional compression under operating conditions. A lower clamping stiffness means further compression of the stack is required to balance the clamping forces, influencing the load distribution between the gasket and MEA contact areas. At the point where the stack transitions from expansion to contraction (80 °C for the parameters in this study), the change in MEA contact pressure is independent of the clamping stiffness.

Fig. 13
Fig. 13
Close modal

### Manufacturing Tolerances.

By applying a normalized random variation to the thickness of each component, as described in Sec. 2.3, the influence of manufacturing tolerances can be modeled and methods of mitigating their impact investigated. Variation in cell MEA contact pressure for a 100-cell stack with properties of Table 2 is shown in Fig. 14(a) for a mean contact pressure of 1.2 MPa. Contact pressure between the cells is seen to vary by ±10%, despite the component thickness standard deviation being only 0.5%. The high contact pressure variation between cells relative to the thickness tolerance is caused by the similarity in nominal thickness and Young's modulus of the GDL and gasket. This means small variations can lead to large changes in force distribution between the two load paths. In comparison, Fig. 14(b) shows the variation in MEA contact pressure if a thicker gasket is used with a lower Young's modulus under the same assembly and tolerance conditions. Cell-to-cell variation is seen to reduce by 50% compared to Fig. 14(a). The thickness difference between the gasket and GDL means that the relative error has a reduced effect on the overall contact pressure leading to a more uniform pressure distribution throughout the stack.

Fig. 14
Fig. 14
Close modal

The influence of material stiffness tolerance on MEA contact pressure variation can also be explored by applying the same normalized distribution to the material Young's modulus. Using this method, mean cell-to-cell MEA contact pressure variation for a Young's modulus standard deviation of 0.5% is ten times less than the same thickness standard deviation. A small variation in Young's modulus creates a proportional variation in component stiffness (Eq. (2)) and a corresponding small change in force for the same displacement. However, a small variation in component thickness can lead to a large variation in displacement, corresponding to a large change in force for the same stiffness. This is because the thickness error (1 μm standard deviation for a 200 μm GDL) is large relative to the mean component compression (37.3 μm for the GDL in Fig. 14(a)). Therefore, ensuring good control over thickness tolerance is more important than material stiffness tolerance in stack designs utilizing similar GDL and gasket thickness.

## Conclusions

An equivalent stiffness model of a PEM fuel cell stack has been presented, which acts as a useful tool for the development of contact pressure distribution in fuel cell stack designs at a low computational cost. Good agreement was seen when compared to 3D FEA simulations for multicell stacks and experimental hydration data from the literature. The model is the first of its kind to consider nonlinear component compression, thermal, and hydration effects under operational conditions and manufacturing tolerances at a stack level. Results show that change in MEA contact pressure due to humidity is greater than the change due to temperature across the operating region of a PEMFC. Reduction in membrane stiffness at operational conditions also has a strong influence on the change in contact pressure. In some cases, this can lead to a reduction in MEA contact pressure and stack shrinkage at high clamping loads. Thickness tolerances with small standard deviations ($±1 μm$) can lead to MEA contact pressure variations of up to 10% between cells. The influence of manufacturing tolerances on contact pressure variation can be minimized through selecting gaskets with different thickness and stiffness to that of the GDL.

## Acknowledgment

The underlying research used in this publication can be found online.2

## Funding Data

• Engineering and Physical Sciences Research Council (Grant No. EP/M023508/1).

## Nomenclature

• A =

area (m2)

•
• aw =

water vapor activity

•
• d =

diameter (m)

•
• E =

Young's modulus (N/m2)

•
• F =

force (N)

•
• h =

height (m)

•
• k =

stiffness (N/m)

•
• Lstack =

stack length (m)

•
• nbolt =

number of bolts

•
• t =

thickness (m)

•
• T =

temperature (°C)

•
• w =

width (m)

•
• x =

displacement (m)

•
• α =

thermal expansion coefficient (K−1)

•
• β =

hydration expansion coefficient ($λ−1$)

•
• δ =

thickness difference (m)

•
• λ =

membrane water content

### Subscripts

Subscripts

• bpp =

bi polar plate

•
• c =

cell

•
• g =

•
• gdl =

gas diffusion layer

•
• g, ass =

•
• mea =

membrane electrolyte assembly

•
• pem =

proton exchange membrane

## References

1.
Zhang
,
W.
, and
Wu
,
C.-W.
,
2014
, “
Effect of Clamping Load on the Performance of Proton Exchange Membrane Fuel Cell Stack and Its Optimization Design: A Review of Modeling and Experimental Research
,”
ASME J. Fuel Cell Sci. Technol.
,
11
(2), p.
020801
.
2.
Millichamp
,
J.
,
Mason
,
T. J.
,
Neville
,
T. P.
,
Rajalakshmi
,
N.
,
Jervis
,
R.
,
Shearing
,
P. R.
, and
Brett
,
D. J.
,
2015
, “
Mechanisms and Effects of Mechanical Compression and Dimensional Change in Polymer Electrolyte Fuel Cells a Review
,”
J. Power Sources
,
284
, pp.
305
320
.
3.
Lee
,
W.-K.
,
Ho
,
C.-H.
,
Van Zee
,
J.
, and
Murthy
,
M.
,
1999
, “
The Effects of Compression and Gas Diffusion Layers on the Performance of a PEM Fuel Cell
,”
J. Power Sources
,
84
(
1
), pp.
45
51
.
4.
Wen
,
C.-Y.
,
Lin
,
Y.-S.
, and
Lu
,
C.-H.
,
2009
, “
Experimental Study of Clamping Effects on the Performances of a Single Proton Exchange Membrane Fuel Cell and a 10-Cell Stack
,”
J. Power Sources
,
192
(
2
), pp.
475
485
.
5.
Mason
,
T. J.
,
Millichamp
,
J.
,
Neville
,
T. P.
,
El-kharouf
,
A.
,
Pollet
,
B. G.
, and
Brett
,
D. J.
,
2012
, “
Effect of Clamping Pressure on Ohmic Resistance and Compression of Gas Diffusion Layers for Polymer Electrolyte Fuel Cells
,”
J. Power Sources
,
219
, pp.
52
59
.
6.
Lin
,
J.-H.
,
Chen
,
W.-H.
,
Su
,
Y.-J.
, and
Ko
,
T.-H.
,
2008
, “
Effect of Gas Diffusion Layer Compression on the Performance in a Proton Exchange Membrane Fuel Cell
,”
Fuel
,
87
(
12
), pp.
2420
2424
.
7.
Zhou
,
Y.
,
Lin
,
G.
,
Shih
,
A.
, and
Hu
,
S.
,
2007
, “
A Micro-Scale Model for Predicting Contact Resistance Between Bipolar Plate and Gas Diffusion Layer in PEM Fuel Cells
,”
J. Power Sources
,
163
(
2
), pp.
777
783
.
8.
Bates
,
A.
,
Mukherjee
,
S.
,
Hwang
,
S.
,
Lee
,
S. C.
,
Kwon
,
O.
,
Choi
,
G. H.
, and
Park
,
S.
,
2013
, “
Simulation and Experimental Analysis of the Clamping Pressure Distribution in a PEM Fuel Cell Stack
,”
Int. J. Hydrogen Energy
,
38
(
20
), pp.
6481
6493
.
9.
Carral
,
C.
, and
Mélé
,
P.
,
2014
, “
A Numerical Analysis of PEMFC Stack Assembly Through a 3D Finite Element Model
,”
Int. J. Hydrogen Energy
,
39
(
9
), pp.
4516
4530
.
10.
Carral
,
C.
,
Charvin
,
N.
,
Trouvé
,
H.
, and
Mélé
,
P.
,
2014
, “
An Experimental Analysis of PEMFC Stack Assembly Using Strain Gage Sensors
,”
Int. J. Hydrogen Energy
,
39
(
9
), pp.
4493
4501
.
11.
Montanini
,
R.
,
,
G.
, and
Giacoppo
,
G.
,
2011
, “
Measurement of the Clamping Pressure Distribution in Polymer Electrolyte Fuel Cells Using Piezoresistive Sensor Arrays and Digital Image Correlation Techniques
,”
J. Power Sources
,
196
(
20
), pp.
8484
8493
.
12.
Gatto
,
I.
,
Urbani
,
F.
,
Giacoppo
,
G.
,
Barbera
,
O.
, and
Passalacqua
,
E.
,
2011
, “
Influence of the Bolt Torque on PEFC Performance With Different Gasket Materials
,”
Int. J. Hydrogen Energy
,
36
(
20
), pp.
13043
13050
.
13.
Tang
,
Y.
,
Karlsson
,
A. M.
,
Santare
,
M. H.
,
Gilbert
,
M.
,
Cleghorn
,
S.
, and
Johnson
,
W. B.
,
2006
, “
An Experimental Investigation of Humidity and Temperature Effects on the Mechanical Properties of Perfluorosulfonic Acid Membrane
,”
Mater. Sci. Eng., A
,
425
(
1–2
), pp.
297
304
.
14.
Zhou
,
Y.
,
Lin
,
G.
,
Shih
,
A.
, and
Hu
,
S.
,
2009
, “
Assembly Pressure and Membrane Swelling in PEM Fuel Cells
,”
J. Power Sources
,
192
(
2
), pp.
544
551
.
15.
Lin
,
P.
,
Zhou
,
P.
, and
Wu
,
C.
,
2009
, “
A High Efficient Assembly Technique for Large PEMFC Stacks
,”
J. Power Sources
,
194
(
1
), pp.
381
390
.
16.
Lin
,
P.
,
Zhou
,
P.
, and
Wu
,
C.
,
2010
, “
A High Efficient Assembly Technique for Large Proton Exchange Membrane Fuel Cell Stacks—Part II: Applications
,”
J. Power Sources
,
195
(
5
), pp.
1383
1392
.
17.
Liu
,
B.
,
Wei
,
M.
,
Zhang
,
W.
, and
Wu
,
C.
,
2016
, “
Effect of Impact Acceleration on Clamping Force Design of Fuel Cell Stack
,”
J. Power Sources
,
303
, pp.
118
125
.
18.
Qiu
,
D.
,
Yi
,
P.
,
Peng
,
L.
, and
Lai
,
X.
,
2015
, “
Assembly Design of Proton Exchange Membrane Fuel Cell Stack With Stamped Metallic Bipolar Plates
,”
Int. J. Hydrogen Energy
,
40
(
35
), pp.
11559
11568
.
19.
,
M.
,
Harrison
,
R.
,
Meredith
,
J.
,
Bindel
,
A.
, and
Todd
,
B.
,
2015
, “
Analysis of the Compression Characteristics of a PEM Stack, Development of an Equivalent Spring Model and Recommendations for Compression Process Improvements
,”
Sixth International Renewable Energy Congress
(
IREC
), Sousse, Tunisia, Mar. 24–26, pp.
1
6
.
20.
Ismail
,
M. S.
,
Hassanpour
,
A.
,
Ingham
,
D. B.
,
Ma
,
L.
, and
Pourkashanian
,
M.
,
2012
, “
On the Compressibility of Gas Diffusion Layers in Proton Exchange Membrane Fuel Cells
,”
Fuel Cells
,
12
(
3
), pp.
391
397
.
21.
Springer
,
T. E.
,
Zawodzinski
,
T. A.
, and
Gottesfeld
,
S.
,
1991
, “
Polymer Electrolyte Fuel Cell Model
,”
J. Electrochem. Soc.
,
138
(8), pp.
2334
2342
.
22.
Wu
,
C.-W.
,
Liu
,
B.
,
Wei
,
M.-Y.
, and
Zhang
,
W.
,
2015
, “
Mechanical Response of a Large Fuel Cell Stack to Impact: A Numerical Analysis
,”
Fuel Cells
,
15
(
2
), pp.
344
351
.
23.
Mason
,
T. J.
,
Millichamp
,
J.
,
Neville
,
T. P.
,
Shearing
,
P. R.
,
Simons
,
S.
, and
Brett
,
D. J.
,
2013
, “
A Study of the Effect of Water Management and Electrode Flooding on the Dimensional Change of Polymer Electrolyte Fuel Cells
,”
J. Power Sources
,
242
, pp.
70
77
.