# Stress modelling from land/channel to full PEMFC stack using homogenisation method

Fabien Mons-Quendo, Jean-François Blachot, Jean-Philippe Poirot-Crouvezier, Christophe Carral, Patrice Mele, Sebastien Kawka

## ► To cite this version:

Fabien Mons-Quendo, Jean-François Blachot, Jean-Philippe Poirot-Crouvezier, Christophe Carral, Patrice Mele, et al. Stress modelling from land/channel to full PEMFC stack using homogenisation method. International Journal of Hydrogen Energy, 2026, pp.152592. 10.1016/j.ijhydene.2025.152592. cea-05526169

HAL Id: cea-05526169

https://cea.hal.science/cea-05526169v1

Submitted on 24 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Stress modelling from land/channel to full PEMFC stack using homogenisation method

Fabien MONS-QUENDO¹, Jean François BLACHOT¹, Jean-Philippe POIROT-CROUVEZIER¹, Christophe CARRAL², Patrice MELE², Sébastien KAWKA¹,*

¹Univ. Grenoble Alpes, CEA, Liten DEHT, 38000 Grenoble, France

²Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, Grenoble INP, LEPMI, 38000 Grenoble, France

*Corresponding author: sebastien.kawka@cea.fr, 17 avenue des Martyrs 38054 Grenoble Cedex, France, T +33 4 38 78 01 60

fabien.mons-quendo@cea.fr

sebastien.kawka@cea.fr

jean-francois.blachot@cea.fr

jean-philippe.poirot@cea.fr

christophe.carral@univ-smb.fr

patrice.mele@univ-smb.fr

## 1 Abstract

Proton exchange membrane fuel cell (PEMFC) are composed of stacked cells compressed between two clamping plates. The performance and durability of the PEMFC components within the stack depend on the nominal pressure set by these plates. Too low pressure may result in poor electrical contact between the cells and potential leakage; conversely, too high pressure can lead to over-compression of the membrane electrode assembly, which could impair the transport properties. The nominal pressure can be set through experimental means or by using simple mechanical models at a land/channel scale inferior to 1mm. However, these methods do not take into account the stress heterogeneities that exist within the hundreds of cells within the few decimetres size stack, or do they consider the stress variation during exploitation or during the clamping process. To be able to perform a full 3D calculation of the whole stack calculate the stress distribution it is necessary to decrease the number of finite elements. Indeed, this calculation would be too costly in terms of both time and resources. In this study, a homogenised mechanical model of the stack is presented. It allows the calculation of the stress distribution within the whole stack and simultaneously at each point of each cell into the stack.

## 2 Keywords

Proton Exchange Membrane Fuel Cell, homogenisation method, stack, mechanical simulation

## 3 Introduction

Proton Exchange Membrane Fuel Cells (PEMFC) for diverse applications, such as heavy-duty transport, passenger cars or stationary plants, generally consist of a stack of a few hundred cells compressed by a clamping system. The components of these cells have complex mechanical properties that can have an impact on their functionality. The local pressures affect the local performances at the rib-channel scale inside each cell. Hence the clamping affects the whole PEMFC stack performance [1–3]. A study of the mechanical behaviour at the stack scale is therefore necessary to estimate these stress heterogeneities and determine how they can affect the global performance.

The cells are considered as the sum of two zones: the peripheral area and the active area. The peripheral area is made of a hyperelastic material whose compression ensures the sealing of the device. The active area is composed of numerous components pressed together: BiPolar Plates (BPP), Gas Diffusion Layers (GDL), MicroPorous Layers (MPL), Proton Exchange Membrane (PEM) and Catalyst Layers (CL). The membrane electrode assembly (MEA) comprises two GDLs positioned laterally adjacent to PEM. A detailed understanding of their mechanical properties and of their interfaces is required to anticipate the effects of the hygrothermal loading on the components and to improve the overall clamping process, leading to improved performances and durability of the device.
The GDL is the most affected component, as the clamping pressure leads to permanent changes in its physical properties [4–9]. Zhang et al. show correlation between the contact pressure at the BPP/GDL interfaces and the electrical contact resistance [10]. Mechanical stress also affects the transport properties of the PEMFC [11–13] and the thermal conductivity [14–16].

Several 2D models at a land/channel scale have been developed. Kleemann et al. showed a strong dependency between the compression and the through-plane electrical resistance [11]. García-Salaberri et al. [17] highlighted a strong reduction of porosity under the land during the clamping process. Studies have investigated the effect of the hygrothermal loading on the local stress redistribution in the different components of the cell [18–21] but most of them have only explored the effect of the compression on the performance [22,23]. None of these models simulate the stress heterogeneity that occurs during the clamping process.

To understand the stress repartition inside a stack, an upscaling of the 2D numerical models is necessary. It exists only few 3D models based on the finite element simulation. Moreover, the mechanical properties of the PEMFC components or their geometry are often simplified. Indeed, a complete mechanical calculation is too time consuming, because the characteristic lengths involved are very wide-ranging, from the cell thickness of a few hundred microns to the final object of the order of a few tens of centimetres. Most of these studies assume homogeneous linear elastic and isotropic properties for the components of the fuel cell in order to reduce the degrees of freedom required to solve a 3D model, which greatly decreases the calculation time of these numerical simulations.

Carral and Mele [24] proposed a model of a stack of 2 to 16 cell with an active area of 85 cm² by cell to study the effect of the clamping process on the MEA and the influence of the cell's number inside a stack. They showed that a minimum number of cells is required to limit the stress heterogeneity in the MEA. Firat et al. proposed a model of a 5-cells stack with an active area of 50 cm² by cell, detailing all the components. It can be noted that they used an orthotropic law with a strain dependent Young's modulus for the GDL [25]. Bates et al. proposed a numerical study of a 16-cells stack with an active area of 100 cm² by cell showing over compression at the seal area [26]. Chien et al. investigated the effect of the bolt preload on the performance of the GDL using a 3D model for a 25 cm² active area. They show an optimal value for clamping pressure, at which point performance is maximised. [27]. Ouaidat et al. proposed a land/channel 3D model to investigate the electrical resistance via a design of experiment. A set of parameters is identified as a means of optimising the performance of the PEFMC.[28].

One way to overcome the too numerous degrees of freedom defining the active area is to use a homogenisation method. This method is widely referred in literature, notably for composite materials [29–33]. An important decrease in computational cost maintaining a good accuracy on the stress distribution is reached by replacing the initial heterogeneous media by a homogeneous equivalent medium [34–36]. Only a few articles deal with the application of this method in the mechanical modelling of fuel cells. Charon et al. proposed a mechanical simulation of the stack clamping process through a 3D model using a homogenised model of stamped metal bipolar plates [37]. They obtained a good agreement between the Representative Elementary Volume (REV) with 750,000 degrees of freedom and the equivalent model with 1000 degrees of freedom. Kuts et al. proposed a homogenisation of the linear elastic properties in porous non-woven fibrous materials and show that the porosity and the fiber's intersection ratio are the main factors influencing the mechanical properties of their porous media [38]. These two studies show the interest in using a homogenisation method to simulate complex mechanical behaviour.
This study presents the development of a homogenisation method for the whole active area. The experimental validation of the model is achieved through the implementation of a compression test of the MEA's project, which is conducted between the real land/channel pattern. The model obtained by using this method is then used to simulate the clamping process of a stack with 100 cells of 280 cm². Finally, the strain heterogeneities induced by the clamping process, as determined by the homogenised model, are utilised to return to the land/channel scale.

# 4 Material and methods

## 4.1 Definition of the Representative Elementary Volume

The principle of a homogenisation method is to determine the mechanical properties of a volume containing a complex structure representative of the whole assembly. Due to the periodic pattern of the active area (Figure 1), it is possible to define a REV corresponding to one period of the geometry, with a longitudinal dimension corresponding to the transverse dimension.

It is made of several components: one BPP composed of two stamped metal sheets disposed at the top and bottom of the reference geometry, two GDLs including MPLs, and a membrane equipped with CLs on both sides, forming a Catalyst Coated Membrane (CCM). The geometry of the REV is shown in Figure 1. The geometry of the flow fields corresponds to the ZSW HyFaB generic stack's design [39].

A Representative Elementary Volume (REV) of the rib/channel design is defined, on which a series of loads is numerically performed with a detailed 3D model to obtain some macroscopic mechanical properties. These mechanical properties are then implemented in a homogenized model (Figure 1).

<figure><img src="image_2.png" /><figcaption>Figure 1: Principle of a homogenisation method applied on the active area of a PEMFC: in a cross section of the stack (a), the periodic pattern of the active zone is analysed (b) to define the Representative Elementary Volume (c), whose macroscopic mechanical properties are implemented in the Homogenized model (d)</figcaption></figure>

## 4.2 Materials properties

The mechanical structure of the REV is given by the Bipolar Plates (BPP) made of stamped 316L stainless steel. In between lies the MEA (see flat layers Figure 1). The mechanical parameters of the CCM are taken from the work of Kusoglu et al. [19], neglecting the mechanical behaviour of the catalyst layers. The GDL is the Freudenberg reference H15CX483. Its mechanical behaviour is fitted to experimental data thanks to the model of Carral and Mele [40], which divides the mechanical behaviour of the GDL compression into 3 parts.

- The first part corresponds to the initial behaviour after manufacturing (Eq 1)
$$ \sigma_1(\varepsilon_{zz}) = A_1 \left( \left( \frac{\rho_1}{1 - \varepsilon_{zz}} \right)_1^n - \rho_1^{n_1} \right) + \sigma_c \qquad \text{Eq 1} $$

- The second part corresponds to the original behaviour of the GDL (Eq 2)

$$ \sigma_2(\varepsilon_{zz}) = A_2 \left( \frac{\rho_1}{1 - \varepsilon_{zz}} \right)^{n_2} \quad \text{with} \quad A_2 = \frac{\sigma_p}{\left( \frac{\sigma_p - \sigma_c}{A_1} + \rho_1^{n_1} \right)^{\frac{n_3}{n_1}}} \qquad \text{Eq 2} $$

- The third part corresponds to the mechanical behaviour of the GDL after a first compressive stress (Eq 3)

$$ \sigma_3(\varepsilon_{zz}) = A_3 \left( \frac{\rho_1}{1 - \varepsilon_{zz}} \right)^{n_3} \quad \text{with} \quad A_3 = \begin{cases} \frac{\sigma_p}{\left( \frac{\sigma_c - \sigma_M}{A_1} + \rho_1^{n_1} \right)^{\frac{n_3}{n_1}}} & \text{if } \sigma_M < \sigma_p \\ \frac{\sigma_p}{\left( \frac{\sigma_M}{A_2} \right)^{\frac{n_3}{n_2}}} & \text{if } \sigma_M \ge \sigma_p \end{cases} \qquad \text{Eq 3} $$

With $\sigma_p$ depending of the manufacturing process, $\rho_1$ the initial relative density of the GDL measured at $\sigma_c = 0.1 \text{ MPa}$, $\sigma_M$ the maximum value of stress experienced by the GDL during compression cycles, $A_1, n_1, n_2$ fitted constants. $n_3$ depends on $\sigma_M$, following an exponential law defined in Eq 8, with $n_{30}$ the asymptotic value of $n_3$ when $\sigma_M$ tends to the infinite. B and C are fitted constants.

$$ n_3 = B \exp \left( \frac{\sigma_M}{C} \right) + n_{30} \qquad \text{Eq 4} $$

The values of the parameters fitted to the experimental results of the compression of the GDL take from the MEA are presented in Table 1. The experimental setup, consisting of 5 GDL samples separated by metal plates, is similar to the one used by Le Carre et al. [41].

<table><caption>Table 1: Fitted parameters to model the compression behaviour of the GDL H15CX483</caption><thead><tr><th>Ref.</th><th>σ<sub>p</sub> [MPa]</th><th>ρ<sub>1</sub></th><th>A<sub>1</sub> [MPa]</th><th>n<sub>1</sub></th><th>n<sub>2</sub></th><th>B</th><th>C [MPa]</th><th>n<sub>30</sub></th></tr></thead><tbody><tr><td>precharged H15CX483</td><td>0.8</td><td>0.30</td><td>3290071</td><td>14.6</td><td>12.6</td><td>-43432</td><td>-2932</td><td>43462</td></tr></tbody></table>

The remaining mechanical parameters are derived from the work of Ouerghemmi et al. [42] and are shown in Table 2. It has been observed that the initial thickness of the GDL taken from the MEA does not correspond to the value indicated in the datasheet. This phenomenon is attributable to the assembly process of the two GDL with the CCM by hot pressing. It is evident that the GDL undergoes an initial compression during this process, which results in irreversible strain and a permanent modification of the mechanical behaviour.

The mechanical parameters of the REV components are detailed in Table 2.

<table><caption>Table 2: Physical properties of the REV's components</caption><thead><tr><th>Component</th><th>Material</th><th>Thickness [µm]</th><th>Young Modulus [MPa]</th><th>Poisson's coefficient</th></tr></thead><tbody><tr><td rowspan="2">BPP</td><td>Stainless steel<br>316l</td><td>100</td><td>200E3</td><td>0.3</td></tr><tr><td>Freudenberg<br>H15CX483<br>after assembly</td><td>191 raw[43]<br>164 after<br>assembly</td><td>E<sub>x</sub> = 3000<br>E<sub>y</sub> = 3000<br>E<sub>z</sub>(ε<sub>zz</sub>)</td><td>ν<sub>xz</sub>, ν<sub>yz</sub> = 0<br>ν<sub>xy</sub> = 0.2</td></tr><tr><td>Membrane</td><td>Nafion</td><td>26</td><td>197 [19]</td><td>0.25</td></tr></tbody></table>

To simplify the calculation, the following assumptions are added:
* The mechanical behaviour of the catalyst layers is not considered.

* The BPP and the membrane are assumed isotropic linear materials.

* The GDL's mechanical behaviour is represented by an orthotropic nonlinear elasto-plastic law with strain hardening, similarly to our previous work [45]. The principle is to consider the mechanical behaviour of the discharge as the elastic part of the elasto-plastic law with a low yield limit. The hardening function is then defined to reproduce the behaviour of the charge. This model allows to take into account the irreversible strain of the GDL after compressive stress. The thickness of the GDL of the MEA, i.e. after assembly, has been determined at 0.1 MPa of compression to be 164 ± 2 µm. Due to the porous nature of the GDL, its Poisson's coefficient is neglected for a compression in its through plane/z direction [11]. The contact between the initial land of the BPP and the GDL is modelled using a sharing of topology, as a perfect contact. The contact between the curved part of the BPP and the GDL during the clamping is set as a contact without friction.

* The contact between the various components of the membrane electrode assembly is also modelled using a shared topology.

## 4.3 Stack clamping process

The clamping of the stack is divided in two parts

* The stack is compressed by applying a displacement to the upper end plate to obtain an average value of 1MPa of compressive stress on the active area, named thereafter $P_{imp}$.

* The straps are fastened

* The imposed displacement is released.

The step modelled in this study corresponds to the clamping process that takes place at ambient temperature and relative humidity (25°C and 30% HR). The GDL's physical properties show irreversible changes after compression. This phenomenon has been widely described in the literature [4,5,11], but, it is not taken into account during stack clamping because of the complexity of the involved mechanisms. Due to the stress relaxation that occurs during the second stage, a model that takes into account the irreversible strain of the REV after compression is developed. It corresponds to an elasto-plastic behaviour to model the change in the mechanical behaviour of the REV of the active area induced by the compressive stress. It can model the evolution of the clamping pressure over the life cycles of PEMFCs giving information to improve performance and durability.

### 4.3.1 Homogenized model of the active area

The stress/strain relationship in anisotropic elasticity is given by a Hook's law $\sigma = C\epsilon$, with $C$ the elastic tensor, $\sigma$ the Cauchy stress tensor and $\epsilon$ the strain tensor. A Hook's law can be written as Eq 5 [44].

$$ \begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{zz} \\ \sigma_{xy} \\ \sigma_{xz} \\ \sigma_{yz} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{13} & C_{14} & C_{15} & C_{16} \\ C_{21} & C_{22} & C_{23} & C_{24} & C_{25} & C_{26} \\ C_{31} & C_{32} & C_{33} & C_{34} & C_{35} & C_{36} \\ C_{41} & C_{42} & C_{43} & C_{44} & C_{45} & C_{46} \\ C_{51} & C_{52} & C_{53} & C_{54} & C_{55} & C_{56} \\ C_{61} & C_{62} & C_{63} & C_{64} & C_{65} & C_{66} \end{bmatrix} \begin{bmatrix} \epsilon_{xx} \\ \epsilon_{yy} \\ \epsilon_{zz} \\ \gamma_{xy} \\ \gamma_{xz} \\ \gamma_{yz} \end{bmatrix} \qquad \text{Eq 5} $$

Taking into account the three orthogonal symmetries present in the REV, the elastic tensor $C$ can be simplified to correspond to an orthotropic material (Eq 6).
Eq 6

$$C = \begin{bmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{21} & C_{22} & C_{23} & 0 & 0 & 0 \\ C_{31} & C_{32} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{bmatrix}$$

To determine the nine coefficients of the elastic tensor, the elementary volume is subjected to six loading cases, shown in Figure 2. For the sake of clarity, only the force or the displacement imposed on one side of the representative element are shown; these are represented by blue arrows for displacement and red arrows for force. Cases 4, 5 and 6 are controlled in displacement to simulate a pure shear stress test on the REV. It should be noted that the normal displacement of the opposite side of which the force or displacement is applied is blocked. A release of the boundary conditions is done for each case, providing the mechanical behaviour of the representative element after compression.

<figure><img src="image_3.png" /><figcaption>Figure 2: The six loading cases applied on the active area REV. Blue arrows for displacement and red arrows for the force, the opposite face being a roller condition. Contours in red correspond to the initial state of the REV. The deformation has been amplified by 10 for the cases 2, 3, 4, 5 and 6 for the purpose of visualisation.</figcaption></figure>

The average stress and strain of the elementary volume are computed on the areas bounding the REV for each loading case. Figure 3 shows the calculated evolution of the compressive stress function of the REV during charge and discharge for the loading case 1. The irreversible strain of the GDL under cyclic compressive stresses generates a modification of the mechanical behaviour during the discharge.
<figure><img src="image_4.png" /><figcaption>Figure 3: Computed evolution of the compressive stress function of the thickness of the MEA between the lands for the REV.</figcaption></figure>

For a perfect orthotropic material, the elastic tensor is symmetric. However, even if the assumption of plane symmetries is done to simplify the elastic tensor, it remains non-symmetrical due to the homogenisation method [46]. During the clamping process and the life cycles of PEMFC, the components are mainly constrained in compression: the elastic coefficients determined by the loading case 1 are preferred ($C_{13}$ and $C_{23}$) instead of the symmetric coefficients determined by loading cases 2 and 3 ($C_{31}$ and $C_{32}$). To choose between $C_{12}$ and $C_{21}$, it is assumed that the stack is longer in the x-direction. In this case, the bending of the end plate happening during the clamping process induces more compression in the x-axis, then for this case, $C_{12}$ is more representative than $C_{21}$. Afterwards, the contribution of these terms is negligible: during the clamping process, the influence of the in-plane strain on the stress is negligible compared to the contribution of the out-of-plane strain on the stress.

The coefficients of the elastic tensor (Eq 7 to Eq 12) are then calculated from Eq. 1.

Loading case 1: $\varepsilon_{xx} = \varepsilon_{yy} = 0 \leftrightarrow C_{13} = \frac{\sigma_{xx}}{\varepsilon_{zz}}; C_{23} = \frac{\sigma_{yy}}{\varepsilon_{zz}}; C_{33} = \frac{\sigma_{zz}}{\varepsilon_{zz}}$ Eq 7

Loading case 2: $\varepsilon_{yy} = \varepsilon_{zz} = 0 \leftrightarrow C_{11} = \frac{\sigma_{xx}}{\varepsilon_{xx}}; C_{21} = \frac{\sigma_{yy}}{\varepsilon_{xx}}; C_{31} = \frac{\sigma_{zz}}{\varepsilon_{xx}}$ Eq 8

Loading case 3: $\varepsilon_{xx} = \varepsilon_{zz} = 0 \leftrightarrow C_{12} = \frac{\sigma_{xx}}{\varepsilon_{yy}}; C_{22} = \frac{\sigma_{yy}}{\varepsilon_{yy}}; C_{32} = \frac{\sigma_{zz}}{\varepsilon_{yy}}$ Eq 9

Loading case 4: $\tau_{xz} = \frac{F_z}{l*\hbar}; \gamma_{xz} = \frac{d_z}{l}; C_{66} = \frac{\tau_{xz}}{\gamma_{xz}}$ Eq 10
$$ \text{Loading case 5: } \tau_{yz} = \frac{F_z}{l*h} ; \gamma_{yz} = \frac{d_z}{L} ; C_{55} = \frac{\tau_{yz}}{\gamma_{yz}} \qquad \begin{array}{l} \text{Eq} \\ 11 \end{array} $$

$$ \text{Loading case 6: } \tau_{yx} = \frac{F_z}{l*h} ; \gamma_{yx} = \frac{d_z}{L} ; C_{44} = \frac{\tau_{yx}}{\gamma_{yx}} \qquad \begin{array}{l} \text{Eq} \\ 12 \end{array} $$

The elastic coefficient $C_{33}$ is calculated using a secant modulus and divided in two parts: one for the charge $C_{33charge}$ and one for the discharge $C_{33discharge}$ with a polynomial fit and an exponential fit, respectively. The elastic stress is obtained by Eq 13.

$$ \sigma = C : \varepsilon \qquad \begin{array}{l} \text{Eq} \\ 13 \end{array} $$

It is necessary to define $C_{33discharge}$, $C_{13}$ and $C_{23}$ considering the elastic strain $\varepsilon_{zzel}$ to separate the plastic behaviour from the elastic behaviour. As a reminder, the engineering strain is defined by $\varepsilon_{zz} = \varepsilon_{zzel} + \varepsilon_{zzinel}$. The variations of these coefficients are plotted in Figure 4 as a function of $\varepsilon_{zzel}$.

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: variations of a) $C_{13}$, $C_{23}$ and $C_{3discharge}$ versus the elastic strain $\varepsilon_{zzel}$ and b) $C_{33charge}$, versus $\varepsilon_{zz}$.</figcaption></figure>

Coefficients $C_{13}$, $C_{23}$, $C_{33discharge}$ depend of $\varepsilon_{zzel}$ and $C_{33charge}$ depends of $\varepsilon_{zz}$: this is due to the non-linear GDL's behaviour in compression. This definition of the elastic coefficients allows to describe the non-linear mechanical behaviour of the REV using Hooke's law. A similar method is achieved by Kleemann et al. to model the nonlinear behaviour in compression of the GDL [11].

The discharge part is set to correspond to the elastic part of the law for a homogenized model with a very low yield limit $\sigma_{ys0}$. The best fit is obtained at $\sigma_{ys0}$ equal to 0.12MPa, with a corresponding value of strain $\varepsilon_{ys0}$. A strain hardening function is then created to reproduce the charge behaviour after the criteria is exceeded (Eq 14).

$$ \sigma_h(\varepsilon_{zz}) = C_{33charge}(\varepsilon_{zz}) - C_{33discharge}(\varepsilon_{ys0}) \qquad \text{Eq 14} $$

The value of the remaining elastic coefficients is given in Table 3.
<table><caption>Table 3: Elastic coefficients for the elastic model of the active area</caption><thead><tr><th>Elastic coefficient</th><th>Value [MPa]</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>5.00e4</td></tr><tr><td>C<sub>22</sub></td><td>5.30e3</td></tr><tr><td>C<sub>12</sub></td><td>2.60e3</td></tr><tr><td>C<sub>44</sub></td><td>5.11e3</td></tr><tr><td>C<sub>55</sub></td><td>5.03e-1</td></tr><tr><td>C<sub>66</sub></td><td>4.77</td></tr></tbody></table>

### 4.3.2 Experimental validation of the model for the active area.

An experimental validation has been done to validate the loading case 1 of the REV. Two samples of MEA have been compressed between three representative cuts of bipolar plates (25x45 mm) (Figure 5.a). This stack is compressed to 1 MPa, and the stress is released. The experiment is performed using a 3R Syntax 100kN machine.

The evolution of the MEA's strain $\epsilon_{zz}$ as a function of the pressure imposed on the metallic plates $P_{imp}$ is shown on the Figure 5.b. The tests demonstrated satisfactory repeatability. A comparison between the REV simulation and the experimental measurements demonstrates the model's consistency (Figure 5.b) and validates its use at the stack level.

a)

<figure><img src="image_7.png" /></figure>

b)

<figure><img src="image_8.png" /><figcaption>Figure 5: a) Experimental setup b) comparison between numerical and experimental results of the evolution of $\sigma_{zz}$ as a function of the MEA's thickness</figcaption></figure>

## 4.4 Heterogeneity of the stress induced by the clamping process

### 4.4.1 Geometry

The homogenized model of the active area presented in the precedent section was implemented in a numerical model of the HyFab generic stack designed by ZSW [39] with 100 cells (Figure 6.a and Figure 6.b)
a)

<figure><img src="image_9.png" /></figure>

b)

<figure><img src="image_10.png" /><figcaption>Figure 6: a) Visualization of a cell with the sealing area (red), the active area (green) and the distribution area (grey) c) Schematic of the ZP/ZA pattern. The thickness of the sealing zone corresponds to $h_{ZP}$, and the thickness of the active area corresponds to $h_{ZA}$.</figcaption></figure>

## 4.4.2 Modelling and boundary conditions

The following assumptions are made:

* The region indicated in Figure 6.a as the 'distribution area' is not associated with an increase in compressive rigidity of the cell. Usual experimentation with Fuji paper demonstrated that this latter remains uncompressed during the clamping process. To avoid an effect on the calculated stress distribution, the modulus of elasticity in compression ($C_{33}$) of this part is defined to be one order of magnitude lower than the active area. The other data of the elastic tensor are assumed equivalent to the active area, which facilitates the transmission of the transversal effect of the compression.

* The contact between the cells and the end plates are modelled as a perfect contact i.e. bonded.

* To improve the stress repartition and assure a better homogeneity of the clamping pressure during life loading, two springs are used between the top end plate and the first cell. To simplify the calculation, the springs are modelled with an equivalent stiffness.

* Due to a lack of data about the mechanical behaviour of the sealing area, an equivalent elastic orthotropic material is assumed. Components of the elasticity tensor necessary for the modelling are given in Table 4. To investigate the influence of the compression behaviour of this area, a parametric study is done on the $C_{33}$ coefficient.
<table><caption>Table 4: Elastic coefficients for the elastic model of the sealing area</caption><thead><tr><th>Elastic coefficient</th><th>Value [MPa]</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>4.20E+4</td></tr><tr><td>C<sub>22</sub></td><td>6.10E+3</td></tr><tr><td>C<sub>33</sub></td><td>1e2, 5e2 or 1e3</td></tr><tr><td>C<sub>13</sub></td><td>2.32E+2</td></tr><tr><td>C<sub>23</sub></td><td>7.15E+2</td></tr><tr><td>C<sub>12</sub></td><td>1.91E+3</td></tr><tr><td>C<sub>44</sub></td><td>4.65E+2</td></tr><tr><td>C<sub>55</sub></td><td>12.8</td></tr><tr><td>C<sub>66</sub></td><td>5.74</td></tr></tbody></table>

The materials properties for the other parts are given in Table 5.

<table><caption>Table 5: Materials properties for the components of the stack</caption><thead><tr><th>Component</th><th>Material</th><th>Young Modulus [MPa]</th><th>Poisson's coefficient</th></tr></thead><tbody><tr><td>End plates</td><td>Aluminum</td><td>70 000</td><td>0.33</td></tr><tr><td>Active area</td><td>Homogenized</td><td>-</td><td>-</td></tr><tr><td>Sealing area</td><td>models</td><td>-</td><td>-</td></tr><tr><td>Strap</td><td>HX380LAD</td><td>200 000</td><td>0.29</td></tr><tr><td>Insulation plates</td><td>POM</td><td>2600</td><td>35</td></tr><tr><td>Current Collector</td><td>AlMg3</td><td>68 000</td><td>3</td></tr><tr><td>Spring</td><td>-</td><td>7.72</td><td>0</td></tr></tbody></table>

The thicknesses in the sealing area $h_{ZP}$ (Figure 6.c) is designed in function of the height and rigidity of the seal to control the distribution of stress generating the tightness the compression of the active area ($h_{ZA}$). To highlight this dimensioning, a study of the thickness variation of the seal on the compressive stress distribution has been conducted. Three cases are simulated:

1. ZP's thickness equal to the ZA's thickness: $h_{ZP}/h_{ZA}=1$

2. ZP's thickness equal to 0.98 of the ZA's thickness: $h_{ZP}/h_{ZA}=0.98$

3. ZP's thickness equal to 0.96 of the ZA's thickness: $h_{ZP}/h_{ZA}=0.96$

The clamping process contains three steps that are reproduced numerically. At first the stack without clamping is compressed to obtain about 1MPa in the active area. Then the straps are set up with no force. Finally, the press is released and the elastic part of the cell's deformation tensioned the straps.

# 5 Results and discussion

Firstly, the objective of simulating the clamping process of a stack within a limited timeframe is achieved: for a CAD with about 170000 elements, it takes 40 minutes to run on a workstation equipped with an Intel Xeon Silver 4114 CPU and 64 Go of RAM.

A parametric study indicates that the elastic coefficient of the seal $C_{33}$ influences the stress distribution on the sealing area (Figure 7). On the one hand, a more flexible seal results in a more homogeneous stress on the sealing surface, but reduces the peak value of the compressive stress. On the other hand, a stiffer material, i.e. $C_{33} = 1000$ MPa, results in a low compression or even in a traction behaviour within specific regions of the sealant. This is not compatible with the stack's tightness. A value of 500 MPa for $C_{33}$ seems to be a good compromise between tightness capacity and stress homogeneities.
<figure><img src="image_11.png" /><figcaption>Figure 7: Visualisation of the σ<sub>zz</sub> stress on the sealing zone, depending of the coefficient's value of C<sub>33</sub>. Zoom on the lateral side of the middle cell of the stack. h<sub>zP</sub>/h<sub>ZA</sub>=0.96</figcaption></figure>

The difference between the nominal clamping pressure reached during the first step of the clamping
and the pressure in the active area after the tensioning strap is induced by two parameters: the
$h_{zp}/h_{za}$ ratio and the elastic coefficient $C_{33}$ of the sealing area. Depending on the clamping process,
the final compression varies due to the stress relaxation when the press is released and the strap are
tensioned. The model enables the quantitative setting of the initial compression, ensuring that the
final internal pressure is equal to the nominal value. A stiffer seal leads to a decrease in the stress in
the active area, from 4% up to 17% with the given configurations. As the active area is softer, it is
advantageous to design the sealing area to be thinner and the best scenario is obtained with a ratio
$h_{zP}/h_{ZA} = 0.96$ in this case.
<figure><img src="image_12.png" /><figcaption>Figure 8: Stress $\sigma_{zz}$ at the center of the middle cell as a function of the chosen seal configuration, represented by $C_{33}$, and the initial ratio $h_{ZP}/h_{ZA}$.</figcaption></figure>

The homogenized model gives access to the stress distribution in the cell area and in the casing, as shown in Figure 9. The figure indicates notably that the stress in the two terminal plates remains below 50MPa, which is below the yield stress for aluminium.
<figure><img src="image_13.png" /><figcaption>Figure 9: Stress field in the stack at the cell level: $\sigma_{zz}$ (left color bar) and von Mises stress in the casing (right color bar). A half stack and top casing with a full bottom casing and straps is shown. $h_{ZP}/h_{ZA}=0.96$, $C_{33}=5e8[Pa]$.</figcaption></figure>

The compression stress at the cell scale obtained at the top and middle cell is plotted in Figure 10. In both cases, the active area compression is close to the expected 1MPa, the distribution area compression remains negligible and the seal experiences higher values as expected.

The compression along the seal path is analysed as a value above 2MPa is expected for a reliable sealing. Figure 10 reveals that, in this configuration, the sealing pressure is lower at the extreme parts of the cell in the middle of the stack which constitutes weak points in the seal where the risk of leakage is high.
<figure><img src="image_14.png" /><figcaption>Figure 10: Stress $\sigma_{ZZ}$ at the cell scale for the top cell (top) and the middle cell (bottom), $h_{ZP}/h_{ZA}=0.96$, $C_{33}=5e8$[Pa].</figcaption></figure>

The GDL, being much softer than the bipolar plates, undergoes most part of the deformation and a zoom on this component GDL (Figure 11) highlights the compression heterogeneities at the top or at middle of the stack. The deformation lies around 20% except on the sides where it reaches 30% for the top cell and decrease to 15% for the middle cell. These heterogeneities will affect locally the transport and resistance properties and the associated performance and durability of cells. Hence the model can help conceiving a clamping system that reduce these heterogeneities.

<figure><img src="image_15.png" /><figcaption>Figure 11: Deformation of GDL in the top (a) and middle cell (b), $h_{ZP}/h_{ZA}=0.96$, $C_{33}=5e8$ [Pa].</figcaption></figure>
A notable benefit of the method is its ability to return to REV scale to study the stress heterogeneities. The simulation of the clamping of stacks reveals a significant disparity in the compression ratio of GDL close to the sides ($ε_{zz} ≈ 15 %$) or at the centre of the active area ($ε_{zz} ≈ 30%$) for cells in proximity of the upper clamping plate (Figure 11). For these, Figure 12 shows the difference in the compressive strain distribution between a REV close to the side and one of the central part of the cell. The stress heterogeneities change the value of the intrusion of the GDL within the channel, the contact pressure and the mass transport properties of the porous medium. These data improve multiphysics simulation and enhance comprehension of the impact of stress heterogeneities on performance and especially on durability. It is visible that the lack of symmetry between the anodic and the cathodic bipolar plate induces a bending of the membrane.

The drift in compression occurring during the hygrothermal cycles, could also be taken into account to simulate the performances and seal properties of the stack over time.

<figure><img src="image_16.png" /><figcaption>Figure 12: Compressive strain of the MEA of the REV for a) 15% of deformation and b) 30% of deformation of the GDL under land</figcaption></figure>

# 6 Conclusion

The study of fuel cell stack mechanics is eased by the development of a homogenized model that allows representing the full stack taking into account detailed active area and sealing track. This model was applied to an opensource design of stacks from ZSW. It provides, for each cell, an insight into the compression in both zones, after the clamping process. These results show the influence of the clamping process and of the seal modulus on the compression stress heterogeneities along the stack and all over the active area of each cell. The method presented in this study enables the optimization of a clamping process to minimize these heterogeneities associated to a robust sealing configuration. As a perspective, further developments should be carried out to simulate for example the local current density linked to the stress heterogeneities potentially promoting local aging. The degradation phenomena of PEM cells are complex, but the original data generated by the model could enable to go further in understanding the effect of the stack mechanics.

# 7 Generative AI and AI-assisted technologies

During the preparation of this work the authors used DeepL Write tool to enhance the language and readability of the paper. After using this tool/service, the authors reviewed and edited the content as needed and takes full responsibility for the content of the publication.
## 8 Funding sources

This work was supported by the European Union in the frame of the RealHy-FC project, number 101111904

## 9 CRediT author statement

**Fabien Mons-Quendo:** Writing – original draft, Writing – review & editing, Software, Conceptualization, Methodology, Validation, Formal analysis, Investigation, Data curation, Visualization.

**Jean-François Blachot:** Conceptualization, Methodology, Validation, Investigation, Resources, Writing – Review & Editing, Supervision, Project administration, Funding acquisition

**Jean-Philippe Poirot-Crouvezier:** Conceptualization, Methodology, Writing – Review & Editing, Supervision

**Christophe Carral:** Conceptualization, Methodology, Writing – Review & Editing, Supervision

**Patrice Mele:** Conceptualization, Methodology, Writing – Review & Editing, Supervision

**Sébastien Kawka:** Writing – original draft, Writing – review & editing, Software, Conceptualization, Methodology, Validation, Formal analysis, Investigation, Data curation, Visualization, Supervision

## 10 References

[1] Khetabi EM, Bouziane K, François X, Lachat R, Meyer Y, Candusso D. In-situ experimental investigations to study the impact of mechanical compression on the PEMFC - analysis of the global cell performance. International Journal of Hydrogen Energy 2024;56:1257-72.
https://doi.org/10.1016/j.ijhydene.2023.12.293.

[2] Khetabi EM, Bouziane K, Zamel N, François X, Meyer Y, Candusso D. Effects of mechanical compression on the performance of polymer electrolyte fuel cells and analysis through in-situ characterisation techniques - A review. Journal of Power Sources 2019;424:8-26.
https://doi.org/10.1016/j.jpowsour.2019.03.071.

[3] Wen C-Y, Lin Y-S, Lu C-H. Experimental study of clamping effects on the performances of a single proton exchange membrane fuel cell and a 10-cell stack. Journal of Power Sources 2009;192:475-85.
https://doi.org/10.1016/j.jpowsour.2009.03.058.

[4] Escribano S, Blachot J-F, Ethève J, Morin A, Mosdale R. Characterization of PEMFCs gas diffusion layers properties. Journal of Power Sources 2006;156:8-13.
https://doi.org/10.1016/j.jpowsour.2005.08.013.

[5] Koorata PK, Bhat SD. Compressive cyclic response of PEM fuel cell gas diffusion media. International Journal of Hydrogen Energy 2021;46:5570-9.
https://doi.org/10.1016/j.ijhydene.2020.11.023.

[6] Lee W, Ho C-H, Van Zee JW, Murthy M. The effects of compression and gas diffusion layers on the performance of a PEM fuel cell. Journal of Power Sources 1999;84:45-51.
https://doi.org/10.1016/S0378-7753(99)00298-0.

[7] Mason TJ, Millichamp J, Neville TP, El-kharouf A, Pollet BG, Brett DJL. Effect of clamping pressure on ohmic resistance and compression of gas diffusion layers for polymer electrolyte fuel cells. Journal of Power Sources 2012;219:52-9.
https://doi.org/10.1016/j.jpowsour.2012.07.021.

[8] Millichamp J, Mason TJ, Neville TP, Rajalakshmi N, Jervis R, Shearing PR, et al. Mechanisms and effects of mechanical compression and dimensional change in polymer electrolyte fuel cells – A review. Journal of Power Sources 2015;284:305-20.
https://doi.org/10.1016/j.jpowsour.2015.02.111.
[9] Radhakrishnan V, Haridoss P. Effect of cyclic compression on structure and properties of a Gas Diffusion Layer used in PEM fuel cells. International Journal of Hydrogen Energy 2010;35:11107-18. https://doi.org/10.1016/j.ijhydene.2010.07.009.

[10] Zhang L, Liu Y, Song H, Wang S, Zhou Y, Hu SJ. Estimation of contact resistance in proton exchange membrane fuel cells. Journal of Power Sources 2006;162:1165-71. https://doi.org/10.1016/j.jpowsour.2006.07.070.

[11] Kleemann J, Finsterwalder F, Tillmetz W. Characterisation of mechanical behaviour and coupled electrical properties of polymer electrolyte membrane fuel cell gas diffusion layers. Journal of Power Sources 2009;190:92-102. https://doi.org/10.1016/j.jpowsour.2008.09.026.

[12] Nitta I, Karvonen S, Himanen O, Mikkola M. Modelling the Effect of Inhomogeneous Compression of GDL on Local Transport Phenomena in a PEM Fuel Cell. Fuel Cells 2008;8:410-21. https://doi.org/10.1002/fuce.200700058.

[13] Zhang H, Hu H, Sarker M, Shao X, Zhan Z, Sui P-C, et al. Numerical investigation of effect of mechanical compression on the transport properties of fuel cell microporous layer using a pore-scale model. International Journal of Hydrogen Energy 2024;62:591-600. https://doi.org/10.1016/j.ijhydene.2024.03.102.

[14] Dafalla AM, Jiang F. Stresses and their impacts on proton exchange membrane fuel cells: A review. International Journal of Hydrogen Energy 2018;43:2327-48. https://doi.org/10.1016/j.ijhydene.2017.12.033.

[15] Jia F, Tian X, Zhang G, Chen M, Ye J, Yang C. Thermal and mechanical investigation of proton exchange membrane fuel cells under combined loading conditions. Applied Thermal Engineering 2024;241:122448. https://doi.org/10.1016/j.applthermaleng.2024.122448.

[16] Shi Q, Feng C, Ming P, Tang F, Zhang C. Compressive stress and its impact on the gas diffusion layer: A review. International Journal of Hydrogen Energy 2022;47:3994-4009. https://doi.org/10.1016/j.ijhydene.2021.10.058.

[17] García-Salaberri PA, Vera M, Zaera R. Nonlinear orthotropic model of the inhomogeneous assembly compression of PEM fuel cell gas diffusion layers. International Journal of Hydrogen Energy 2011;36:11856-70. https://doi.org/10.1016/j.ijhydene.2011.05.152.

[18] Al-Baghdadi MARS. A parametric study of assembly pressure, thermal expansion, and membrane swelling in PEM fuel cells. International Journal of Energy and Environment 2016;7:26.

[19] Kusoglu A, Karlsson AM, Santare MH, Cleghorn S, Johnson WB. Mechanical behavior of fuel cell membranes under humidity cycles and effect of swelling anisotropy on the fatigue stresses. Journal of Power Sources 2007;170:345-58. https://doi.org/10.1016/j.jpowsour.2007.03.063.

[20] Silberstein MN, Boyce MC. Hygro-thermal mechanical behavior of Nafion during constrained swelling. Journal of Power Sources 2011;196:3452-60. https://doi.org/10.1016/j.jpowsour.2010.11.116.

[21] Zhou Y, Lin G, Shih AJ, Hu SJ. Assembly pressure and membrane swelling in PEM fuel cells. Journal of Power Sources 2009;192:544-51. https://doi.org/10.1016/j.jpowsour.2009.01.085.

[22] Kulkarni N, Cho JIS, Jervis R, Roberts EPL, Francesco I, Kok MDR, et al. The effect of non-uniform compression on the performance of polymer electrolyte fuel cells. Journal of Power Sources 2022;521:230973. https://doi.org/10.1016/j.jpowsour.2021.230973.

[23] Zhou Y, Lin G, Shih AJ, Hu SJ. Multiphysics Modeling of Assembly Pressure Effects on Proton Exchange Membrane Fuel Cell Performance. Journal of Fuel Cell Science and Technology 2009;6:041005. https://doi.org/10.1115/1.3081426.

[24] Carral C, Mélé P. A numerical analysis of PEMFC stack assembly through a 3D finite element model. International Journal of Hydrogen Energy 2014;39:4516-30. https://doi.org/10.1016/j.ijhydene.2014.01.036.

[25] Firat E, Beckhaus P, Heinzel A. Finite Element Approach for the Analysis of the Fuel Cell Internal Stress Distribution 2011.
[26] Bates A, Mukherjee S, Hwang S, Lee SC, Kwon O, Choi GH, et al. Simulation and experimental analysis of the clamping pressure distribution in a PEM fuel cell stack. International Journal of Hydrogen Energy 2013;38:6481-93. https://doi.org/10.1016/j.ijhydene.2013.03.049.

[27] Chien C-H, Hu Y-L, Su T-H, Liu H-T, Wang C-T, Yang P-F, et al. Effects of bolt pre-loading variations on performance of GDL in a bolted PEMFC by 3-D FEM analysis. Energy 2016;113:1174-87. https://doi.org/10.1016/j.energy.2016.07.075.

[28] Ouaidat G, Cherouat A, Kouta R, Chamoret D. Study of the effect of mechanical uncertainties parameters on performance of PEMFC by coupling a 3D numerical multiphysics model and design of experiment. International Journal of Hydrogen Energy 2022;47:23772-86. https://doi.org/10.1016/j.ijhydene.2022.05.151.

[29] Fedotov AF. Hybrid model of nonlinear homogenisation of anisotropic composites with ellipsoidal inclusions. Composites Part B: Engineering 2021;222:109076. https://doi.org/10.1016/j.compositesb.2021.109076.

[30] Hoang M, Duong PTM, Abbès B, Guo Y-Q. Modèle d'homogénéisation élasto-plastique pour des structures composites n.d.

[31] Mirkhalaf SM, Van Beurden TJH, Ekh M, Larsson F, Fagerström M. An FE-based orientation averaging model for elasto-plastic behavior of short fiber composites. International Journal of Mechanical Sciences 2022;219:107097. https://doi.org/10.1016/j.ijmecsci.2022.107097.

[32] Tian W, Qi L. Unified periodic boundary condition for homogenizing the thermo-mechanical properties of composites. Applied Mathematical Modelling 2023;121:252-69. https://doi.org/10.1016/j.apm.2023.04.024.

[33] Zhao J, Guo C, Zuo X, Román AJ, Nie Y, Su D-X, et al. Effective mechanical properties of injection-molded short fiber reinforced PEEK composites using periodic homogenization. Adv Compos Hybrid Mater 2022;5:2964-76. https://doi.org/10.1007/s42114-022-00518-y.

[34] Nandakumaran AK. An Overview of Homogenization 2007;87.

[35] Devries F, Dumontet H, Duvaut G, Lene F. Homogenization and damage for composite structures. Numerical Meth Engineering 1989;27:285-98. https://doi.org/10.1002/nme.1620270206.

[36] Bensoussan A, Lions J-L, Papanicolaou G. Asymptotic analysis for periodic structures. Amsterdam ; New York : New York: North-Holland Pub. Co. ; sole distributors for the U.S.A. and Canada, Elsevier North-Holland; 1978.

[37] Charon W, Iltchev M-C, Blachot J-F. Mechanical simulation of a Proton Exchange Membrane Fuel Cell stack using representative elementary volumes of stamped metallic bipolar plates. International Journal of Hydrogen Energy 2014;39:13195-205. https://doi.org/10.1016/j.ijhydene.2014.06.125.

[38] Kuts M, Walker J, Newell P. Computational homogenization of linear elastic properties in porous non-woven fibrous materials. Mechanics of Materials 2024;189:104868. https://doi.org/10.1016/j.mechmat.2023.104868.

[39] Häußler F, Wiedemann B, Bergbreiter C, Schneider S, König L, Cerne B, et al. The ZSW HyFaB Generic Stack. Zentrum fuer Sonnenenergie- und Wasserstoff-Forschung Baden-Wuerttemberg (ZSW); 2024. https://doi.org/10.5281/ZENODO.14223363.

[40] Carral C, Mele P. Modeling the original and cyclic compression behavior of non-woven gas diffusion layers for fuel cells. International Journal of Hydrogen Energy 2022;47:23348-59. https://doi.org/10.1016/j.ijhydene.2022.05.121.

[41] Le Carre T, Blachot J-F, Poirot-Crouvezier J-P, Laurencin J. Mechanical response of carbon paper gas diffusion layer under patterned compression. International Journal of Hydrogen Energy 2023;S0360319923040855. https://doi.org/10.1016/j.ijhydene.2023.08.104.

[42] Ouerghemmi M, Carral C, Mele P. Determination of the Orthotropic Mechanical Elastic Properties of Gas Diffusion Layers for PEMFCs. J Electrochem Soc 2023;170:104508. https://doi.org/10.1149/1945-7111/acfc68.

[43] Freudenberg Gas diffusion layers technical data n.d.
[44] Boresi AP, Schmidt RJ, Sidebottom OM. Advanced mechanics of materials. 5. ed. New York: Wiley; 1993.

[45] Mons-Quendo F, Blachot J-F, Poirot-Crouvezier J-P, Kawka S, Carral C, Mele P. Numerical Investigation of the Irreversible Behavior of GDL Under Cyclic Hygrothermal Loading. J Electrochem Soc 2024. https://doi.org/10.1149/1945-7111/ad5a3f.

[46] Jiménez Segura N, Pichler BLA, Hellmich C. Concentration tensors preserving elastic symmetry of multiphase composites. Mechanics of Materials 2023;178:104555. https://doi.org/10.1016/j.mechmat.2023.104555.