Abstract

In this article, the microscopic deformation of carbon fiber papers for proton exchange membrane (PEM) fuel cell gas diffusion layers (GDLs) under compression is investigated numerically. A nonlinear model for mechanical simulations of the fiber structure is described. It is based on beam bending theory applied to a random fiber structure. Additional nonlinear effects are included in the model to ensure realistic mechanical behavior. The effect of the binder material is included in the model via a spring coupling between close fibers. Damage to fibers and binder, as well as additional fiber contacts, are considered. Second-order terms are added to the beam theory. The model can predict the microscopic deformations of the individual fibers and helps to understand the damage accumulating in the fiber structure during the assembly of a fuel cell. The results of the model highlight the importance of including damage terms in the theoretical description of GDLs. The hysteresis in the stress–strain curve seen in compression experiments of GDLs is reproduced qualitatively. It is shown that the hysteresis is based on both binder damage and fiber damage. While damage to the carbon fibers is limited to large compression ratios greater than 0.16, the model predicts damage to the binder throughout the whole compression range. To understand the influence of the individual components of the system, a parameter study is carried out.

1 Introduction

In the quest for sustainable energy solutions, fuel cells and electrolyzers have emerged as key solutions in the transition toward a cleaner and more efficient energy landscape. These technologies harness the principles of electrochemistry to facilitate energy conversion between electrical energy and chemical energy carriers.

The design of the cells must enable the flow of gas and water molecules to and from the electrochemical reaction site as well as the electrical contact of the electrode. This gives rise to the necessity of a porous layer on the outside of the electrode [1]. In proton exchange membrane fuel cells, a gas diffusion layer (GDL) composed of carbon fibers is commonly employed to facilitate this material access and electrical contact. In this study, carbon fiber paper-type GDLs are examined. These are composed of straight horizontal carbon fibers with a random distribution in-plane, which are glued together by a binder material. Positioned on both sides between the flowfield and membrane electrode assembly, the GDL plays a crucial role in ensuring efficient functioning by enabling simultaneous material transport and electrical connectivity to the electrodes.

During the assembly process, the system of layers is pressed together. As the GDL is the most compressible component inside the cell, the majority of deformation and assembly pressure generated is governed by the mechanical properties of the GDL [2]. A high force acting inside the cell improves the contact resistivity between the flowfield and the GDL. However, as the surface of the GDL is not even due to the pores, the GDL also places demands on the mechanical stability of the membrane [3]. As the membrane is compressed between two GDL layers, GDL's design plays a crucial role in preserving the integrity of the membrane within the assembled cell [4]. This becomes increasingly vital as the membrane thickness is reduced to improve cell performance. Although thin membranes are advantageous in reducing ionic resistance, the reduction of the membrane thickness is accompanied by sensitivity to mechanical stress. Careful consideration must be given to the mechanical properties of the GDL design to have a low mechanical impact on the membrane and the entire cell assembly.

In recent years, studies have been conducted to understand the mechanics and the structure of GDLs [5]. To optimize the structure of the GDL, Balakrischnan et al. experimentally investigated the effect of structural changes in the fiber material on the performance of the cell [6]. Several effects of the GDL mechanics on the cell performance have been discussed. A damaged GDL can have a significant effect on the performance of the cell [7]. This has been shown by corrosion experiments [8], as well as by varying the clamping pressure [9,10]. The causes of the effect of the GDL mechanics on the cell performance and the mechanics of the damage inside the GDL are still poorly understood. Local fiber breaking in the GDL can lead to high local stresses on the membrane or poor mechanical and thus electrical contact in some regions [11].

Theoretically, the mechanical behavior of the GDL has been analyzed by Carral and Mélé [12] and Hoppe et al. [13]. These studies have focused on the macroscopic elastic properties of GDLs and have neglected their microstructure or inelastic deformations. Several proposed models, such as those presented by Carral and Mele [14], aim to characterize the elastic properties of the porous fiber structure in GDLs. Norouzifard and Bahrami [15,16] use the beam bending theory to provide an analytical expression of the GDL compression behavior. They derive the mechanical properties from a fiber unit cell. Works of Kleemann et al. [17], Garcia-Salaberri et al. [18], Afrasiab et al. [19], and Le Carre et al. [20] also include nonlinear behavior in the macroscopic model. These models provide analytical expressions for the macroscopic elastic properties and deformations. They take the fiber structure of the GDL into account; however, no prediction for the local deformation of individual fibers is made. Alternatively, the works of Xiao et al. [21] and Zhang et al. [22] seek to reconstruct the actual microstructure of fibers and binders, offering a more detailed perspective and enabling predictions of the microscopic fiber deformations. They perform a full 3D discretization of the fibers and the binder based on real GDL samples. This entails a high computational cost, limiting its application, especially when examining damage or interactions with the membrane. These models are limited to linear effects and small solution domains.

Generated microscopic fiber distribution models with [23,24] and without [25] binder have been proposed for modeling physical quantities. They have been successfully used to model porosities and flow through a membrane [26], as well as electrical conduction [27]. Yi et al. [28] use a unidirectional fiber model including a binder field to model GDL failure. From a mechanical perspective, generated fiber distribution enables an optimized discretization of the fiber and thus reduces the computational cost of the model.

In the prior work, the model of Norouzifard and Bahrami [15] is generalized from a unit cell consisting of two fibers to a randomly generated fiber structure and a binder model is added [29]. The model is able to predict the local deformation on the level of individual fibers. The computational cost of the model is low enough to include damage such as fiber breaking and binder damage to the model. In the following, the model is described and the influence of the GDL parameters on the mechanical properties is discussed. The role of the damage of both GDL and fibers on the mechanical behavior of the GDL is highlighted.

2 Theory

The model consists of a randomly generated fiber structure, a beam theory model for the mechanical description of the fibers, and a spring model for the binder. In the following, the microscopic mechanical models for the fiber and binder in the structure are described. The discretization and solution procedure is then presented.

2.1 Fiber Model

2.1.1 Fiber Structure Creation.

In a carbon paper GDL, the fibers have an orientation pointing in the plane of the GDL sheet. The out-of-plane angle is small [30].

The creation of such a set of fibers for a GDL-like fiber structure was proposed by Thiedmann et al. [23]. In the model, straight fibers are assumed. The orientation of the undeformed fiber structure is uniformly distributed in the x,y-plane or the plane of the GDL sheet with a zero out-of-plane angle. A single fiber in the structure is defined by the angle in the plane, the in-plane distance of the fiber to a fixed origin, and the out-of-plane position.

A cylindrical domain is considered with an in-plane radius RDomain and a height of HDomain. The entire structure is then constructed by a set of random fibers with a uniform distribution of the in-plane orientation. The random in-plane position is uniformly distributed between 0 and RDomain. Deviating from the model proposed by Thiedmann et al. [23], fibers are not arranged into discrete z-planes. In the model, the z-position of the fibers is equidistant between 0 and the domain height HDomain. A sketch of the domain and the randomly distributed fibers is shown in Fig. 1(a).

Fig. 1
(a) Sketch of the simulation domain. The fibers are randomly distributed in a cylindrical domain with radius RDomain and height HDomain. In the undeformed state, all fibers are aligned in plane and (b) Sketch of the binder model. The binder model adds a coupling in the vicinity of fiber crossings up to a maximum distance rMax.
Fig. 1
(a) Sketch of the simulation domain. The fibers are randomly distributed in a cylindrical domain with radius RDomain and height HDomain. In the undeformed state, all fibers are aligned in plane and (b) Sketch of the binder model. The binder model adds a coupling in the vicinity of fiber crossings up to a maximum distance rMax.
Close modal

During the structure's creation, no interaction between fibers is considered. The intersection of the fibers is possible if they are close. Two fibers are considered to be touching if the crossing point of the straight fibers lies within the domain, and the difference in the z-position is smaller than the fiber diameter DF. The total number of fibers is adjusted such that the porosity of the constructed structure matches the real GDL porosity without the binder.

2.1.2 Fiber Deformation.

Norouzifard and Bahrami [15] and Afrasiab et al. [31] proposed using the beam bending theory to model the mechanics of the fibers. This approach is valid if the fibers are oriented in-plane, which is true for the undeformed configuration as well as for small deflections. These studies have used analytical approaches to solve the bending equations and thus were limited to a small number of fibers. The theory can also be applied to large quantities of fibers. In this case, a numerical approach has to be used to solve the resulting equations.

The equation of the deflection of a fiber ω in the z-direction with respect to the position along the beam x is given by the force per length q(x), the elastic modulus of the beam EF, and the second moment of area of the beam IF:
(1)
For high compression ratios, the deflection of the beam is high enough for higher order terms to have a significant contribution. In the case of second order, the beam equation is given by the following equation [32]:
(2)
with the beams' cross-sectional area AF.
At the crossing point xc of two intersecting fibers, the individual equations of the two fibers i and j are coupled. The deflection of these two fibers is set to be equal:
(3)
On the side boundary, a Neumann-type boundary is used for the deflection:
(4)
This corresponds to a fiber slope at the boundary equal to zero. The top and bottom boundaries are given by a Dirichlet-type boundary for the first n fibers on the top and bottom of the structure:
(5)
where ωt,b is the displacement of the top or bottom boundary for the given step. This corresponds to a fixed deflection of the fibers.

This equation only describes the deflection in the out-of-plane direction. For the given case of a compression for a thin GDL layer, the major force is perpendicular to the GDL plane. As the Poisson’s ratio of porous materials such as GDLs is near zero, the compression in this direction will not lead to a large deformation in-plane. Arising forces and displacements in the x-y directions are considered small and are neglected in this model.

2.2 Binder Model.

By adding additional constraints to the deflections of two fibers connected by the binder material, the average length of those free-bending fiber sections is reduced, making the binder a key component when looking at the mechanical properties of GDLs. Different binder materials have been used in GDLs. In the following, a carbon-based binder is assumed. During production of the GDL, the binder accumulates at fiber crossings, depending on the adhesion energy between the binder base material and the fibers [33]. Figure 2 shows a CT image of the GDL in the center of the sheet and one from close to the surface. The center exhibits less binder. This behavior has previously been reported by Odaya et al. [34].

Fig. 2
CT images of the fiber structure (a) at the center of the sample and (b) near the surface. The binder accumulates at fiber crossings and locations with high fiber density. In the center, less binder is present.
Fig. 2
CT images of the fiber structure (a) at the center of the sample and (b) near the surface. The binder accumulates at fiber crossings and locations with high fiber density. In the center, less binder is present.
Close modal
In this section, a method is discussed that includes the binder in a beam model of paper GDLs using a simple spring model. The binder field sets an additional coupling to the equations of two fibers in the vicinity of the fiber crossing. A line force is added between fiber points that are less than a cutoff distance rmax away from each other. Figure 1(b) shows a basic sketch of the binder model. The fiber segments are coupled by a spring force if their distance is less than rmax. If the displacements ωi(x) and ωj(x′) of the fibers differ, the fiber exerts a force Fz,Bi on the segments.
(6)
where ri,j is the distance between the segments in the undeformed configuration, bi,j is the width of the binder between the segments, li is the length of the segment, and Δri,j is the distance variation in the deformed configuration. x and x′ are determined by the condition that x, x′ and the fiber crossing point form an isosceles triangle in the x-y-plane. This takes the fibrous structure of the binder field, as seen in tomography images, into account. Shear in the binder is not considered. Just as in the fiber deformation, only the force and displacements in the z-direction are considered. The force in the z-direction is given by the following equation:
(7)
where θ is the angle between the segment connection and the z-axis.
The Δri,j can be approximated from the difference in the z-displacement ωiωj and the position of the fibers zi and zj:
(8)
The case of zizj < ωiωj, can be neglected in the case where only force in the z-direction is considered. The force of the binder is thus given by:
(9)
As mentioned earlier, the binder is not homogeneously distributed inside the GDL. To account for the inhomogeneous distribution, a bimodal weighing function B(z,r) is used for the thickness distribution depending on the z-coordinate z and the fiber distance r of the binder field:
(10)
where zmin, zmax, and zmid are the boundaries of the domain and the mid-plane.
The binder strength depends on the width of the binder bi,j connecting two fibers. The volume of all binder connections should be equal to the total binder volume VB.
(11)

2.3 Damage and Additional Coupling.

The mechanical properties of the GDL cannot be described as pure elastic and also depend on inelastic deformations such as the breaking of fibers and damage to the connecting binder. Similar to Ref. [28], a maximum stress is assumed for binder and fiber.

The stress in the binder field can be obtained from Eq. (9). If the stress within the binder exceeds the maximum threshold, the corresponding binder is removed from the calculation. The maximum stress inside a bent fiber occurs near the top and bottom sides of the fiber. It is given by the elastic modulus EF, the diameter DF of the fiber, and the second derivative of the displacement:
(12)

If this stress exceeds the breaking point of a fiber, the fiber will be considered broken at this point. During compression, fibers that are not in contact with the undeformed state may become deformed in such a way that additional contact points arise. The number of contact points increases with higher levels of compression.

2.4 Equation Discretization.

Each fiber is discretized into equidistant segments in the fiber direction. The linearization of the differential equations is done by calculating the difference quotient of the derivative at each point. The first-order equation for a free-bending segment at a point i can be discretized as follows:
(13)

This is the lowest-order approximation with equidistant points possible for m = 4. As this includes terms i−2 and i + 2, the boundary consists of two points on each side of a beam. This results in the linear matrix equation. The matrix equation is solved in python using the linalg package of SciPy [35].

To determine a well-defined breaking point, the second derivative for the damage calculation is evaluated in between two discrete points i and i + 1 at x = xi + 0.5Δx. In discretized form, the maximal stress between points i and i + 1 can be calculated as follows:
(14)
Points i and i + 1 are then treated as boundary points at a free-end boundary or a point force F. The point force boundary is expressed as follows:
(15)

2.5 Model Parameters.

The elastic modulus of carbon fibers and of the carbon binder material has been studied extensively in Ref. [36]. It depends on the production process employed. In this study, it is assumed that EF = EB = 230 GPa [37]. However, it was not possible to confirm this value experimentally. The fiber diameter for GDLs has been reported in several studies to be about DF = 7 µm [15]. The cutoff distance for the binder distribution is assumed to be rmax = 17 µm. The maximum strain in a carbon fiber depends heavily on the production process and the precursor material [37]. Most GDLs are based on polyacrylonitrile (PAN) fibers. In this study, the maximum strain in the fibers and binders is assumed to be εB,Fmax=0.015. The porosity of the fiber structure can be determined using CT images. The fiber volume ratio is assumed to be VF = 0.1, and the binder volume ratio is assumed to be VB = 0.02. A cylindrical domain with a radius of RDomain = 600 µm and a height of HDomain = 110 µm is used for the simulation.

3 Results

3.1 Stress–Strain Curve.

The data obtained from the simulation were compared to compression tests of real GDL samples. The experimental compression tests were performed on an ElektroForce DMA 3200. Single samples of 1 cm diameter were compressed, and the force was measured up to a value of P = 1.9 MPa. A second cycle was performed up to P = 3.8 MPa. The simulation curve showed a compression of up to 0.27. Compression cycles were simulated, up to compression ratios of 0.14 and 0.27.

Figure 3 shows the simulation result and deformation of a real GDL as measured in the experiment. The specification of the manufacturer is highlighted in the graph [38]. The pressure deformation curve features a linear region for small deformations and a convex behavior toward higher ones. In the experimental curve, the inelastic behavior of the sample is seen in the second cycle. The inelastic deformations during compression lead to a lower mechanical resistance in the second cycle, up to the point of maximum previously applied pressure. At the point of highest deformation in the first cycle, the curves merge. The curve after deformation exhibits a steep stress–strain behavior. This behavior has been reported in Ref. [39]. The simulation shows a linear behavior in the stress–strain curve if only the first-order terms are used in the beam equation. If higher orders are included, there is a slightly convex stress–strain curve. The simulation predicts a lower quadratic contribution to the compression/pressure curve than the specifications from the manufacturer and the experimental data. Input parameters might not be accurate for the simulation. Especially the binder parameters are not well understood. The binder model used is linear, whereas in the real GDL, the binder is likely to have additional nonlinear effects that are not considered.

Fig. 3
(a) Pressure against compression ratio for the simulation with hysteresis. The specification given by the manufacturer is marked. (b) Experimental result for the pressure compression curve of the GDL with hysteresis.
Fig. 3
(a) Pressure against compression ratio for the simulation with hysteresis. The specification given by the manufacturer is marked. (b) Experimental result for the pressure compression curve of the GDL with hysteresis.
Close modal

The additional cycles with the damaged GDL material show a similar inelastic behavior to the experiment. The inelastic deformation is not as pronounced in the simulation. This might be partially due to the neglected friction between the fibers, which can play a significant role in such structures [40]. The model assumes a perfect fiber geometry without any production defects or imperfect geometries. This is expected to lead to an earlier fiber breaking in the real GDL than predicted in the model. To include this in the model, additional information is necessary concerning these imperfections on a subfiber length scale. The increasing density of fiber crossings leads to an increase in the elastic modulus. This is counteracted by fiber breaking and damage to the binder structure. For small deformations, the damage to the structure is governed by damage to the binder field.

The reaction of the elasticity of the structure to different effects is shown in Fig. 4. The simulation (curve b) is compared to simulations without damage (curve a), without additional fiber crossings during compression (curve c) and without binder field (curve d). The curve without damage overestimates the stiffness of the system. It shows a smooth convex behavior. There are no fluctuations due to breaking points. If the fiber structure is compressed, fibers that were apart from each other now touch and the number of fiber crossings is increased. This effect is omitted in curve c, which means that the fibers will interpenetrate without mechanical hindrance. The stiffness is lower than in the original simulation, and the convex behavior of the simulation and the experiment is not reproduced. The curve d shows the outcome of the simulation if the binder terms are omitted from the equations. The resulting stiffness is more than a magnitude lower than the stiffness observed in the experiment.

Fig. 4
Effects of the binder and damage terms on the compression of the GDL. Curve b shows the pressure compression curve of the original simulation as depicted in Fig. 3(a) without the hysteresis. Curve a shows the predicted curve of the simulation for the case where binder damage and fiber breaking are not considered in the simulation. For curve c, additional fiber contacts during compression are not considered. In the case of curve d, no binder terms are included in the simulation.
Fig. 4
Effects of the binder and damage terms on the compression of the GDL. Curve b shows the pressure compression curve of the original simulation as depicted in Fig. 3(a) without the hysteresis. Curve a shows the predicted curve of the simulation for the case where binder damage and fiber breaking are not considered in the simulation. For curve c, additional fiber contacts during compression are not considered. In the case of curve d, no binder terms are included in the simulation.
Close modal

Figure 5 shows the percentage of broken binder and the density of fiber breaking points for different compressions in the simulation. Parts of the binder field exceed the maximum stress, and inelastic deformations occur at small compressions of ɛ < 0.01. The binder tends to break at locations where the two connected fibers are far apart. Binder connecting fibers horizontally is more apt to break than vertical binder connections. The fiber breaking occurs at higher displacements of ɛ > 0.16.

Fig. 5
Percentage of the binder field broken and density of fiber breaking points for different compressions. At a compression of 0.25, more than one-third of the binder material exceeded the maximum stress and is removed from the simulation. The first broken fiber occurs at a compression of 0.16. The binder damage is the percentage of the binder material that exceeded the maximum stress at the given compression. The fiber damage is given as the number of broken fibers per volume.
Fig. 5
Percentage of the binder field broken and density of fiber breaking points for different compressions. At a compression of 0.25, more than one-third of the binder material exceeded the maximum stress and is removed from the simulation. The first broken fiber occurs at a compression of 0.16. The binder damage is the percentage of the binder material that exceeded the maximum stress at the given compression. The fiber damage is given as the number of broken fibers per volume.
Close modal

3.2 Parameter Study.

Variations in the model parameters, as well as in the discretization and convergence parameters are performed to determine their influence on the outcome of the simulation. The elastic modulus of the fiber EF and the binder EB material, the breaking point of the binder εBmax, and the spatial extent of the binder rmax are varied by 10%. The resulting force, which the fiber structure exerts at a compression of ɛ = 0.1, is compared to the original simulation. Figure 6 shows the effect of parameter variations on the simulation. When looking at low compressions, the mechanical properties of the binder have a larger influence than the parameters of the fibers. The spatial binder distribution has a great impact on the overall mechanics of the GDL. The compressibility of the entire structure decreases with an increase in the spatial extent of the binder. Variations in the elastic modulus of the binder and the binder breaking point lead to a convex response function. Variations in the elastic modulus of the fiber and the spatial extent of the binder lead to a concave response function.

Fig. 6
Effect of variations of the input parameters on the result of the simulation. The value of the respective variables is changed by a factor of 1.1 and 0.9. The pressure at a compression ratio of 0.1 is plotted for the original simulation and the simulations with variations. The parameters are normalized to the values used in the original simulation. Five input variables are investigated: the fiber and binder elastic modulus, the binder breaking point, the maximum binder distance, and the number of fibers in the domain.
Fig. 6
Effect of variations of the input parameters on the result of the simulation. The value of the respective variables is changed by a factor of 1.1 and 0.9. The pressure at a compression ratio of 0.1 is plotted for the original simulation and the simulations with variations. The parameters are normalized to the values used in the original simulation. Five input variables are investigated: the fiber and binder elastic modulus, the binder breaking point, the maximum binder distance, and the number of fibers in the domain.
Close modal

4 Conclusions

Pristine GDL samples show a well-defined stress–strain behavior. In a second compression cycle, the slope of the stress–strain curve drastically increases. If the GDL is deployed inside the cell at a point where the strain has a large influence on the stress, small changes in the cell displacements due to temperature changes during operation can lead to large changes in the contact pressure between the flowfield and the GDL. A combined binder/fiber model was proposed to describe the inelastic deformations of the GDL to get an insight into the mechanics of the structure. The model predicts a dominant influence of the binder on the GDL's mechanics for small deformations. It qualitatively reproduces the inelastic behavior and provides an explanation of the kink in the stress–strain curve of a precompressed GDL at the point of maximum precompression. This hysteresis is explained by damage to the binder occurring in the whole compression range. It predicts that a significant amount of binder material of more than 25% is subject to damage at intermediate compression ratios of 0.16. At large compressions above compression ratios of 0.16, damage to the fibers also adds to the hysteresis. Variations between the experimental and theoretical results still persist. However, the importance of including mechanical effects from the binder material is shown. Failure to do so results in a stress response to compression of more than a magnitude lower than observed in the experiment. The pure fiber structure cannot explain the observed behavior of the GDL. The results show a significant contribution of damage to binder and fibers, as well as additional fiber contact points during compression, to the mechanical behavior of paper-type GDLs.

Acknowledgment

This work was supported by the German Ministry for Digital and Transport as part of the project “GDL-Qualitätssicherung für den Markthochlauf QM-GDL” (No. 03B110016D) and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) (491111487). The author would also like to thank Dieter Froning for his support. The author declares the following financial interests which may be considered as potential competing interests: The experimental samples were provided by the company SGL Carbon in the framework of the project “GDL-Qualitätssicherung für den Markthochlauf QM-GDL.”

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Appendix

Figure 7 shows the numerical parameter study for the number of discrete points and the residuals. The simulation was done with a mesh size of about 80,000 discrete points. The first compression steps were computed with up to four times more spatial accuracy to estimate the numerical error. The numerical error seems to be smaller than 1% for the spatial error. The residual error is small compared to other numerical errors.

Fig. 7
Accuracy of the numerical parameters. The number of discrete points and the value of the maximal residuals are varied, and the simulation is carried out. The predicted pressures at 0.02 compression are plotted for different numbers of discrete points and for different residuals. The number of discrete points and the residuals are normalized to the values of the original simulation.
Fig. 7
Accuracy of the numerical parameters. The number of discrete points and the value of the maximal residuals are varied, and the simulation is carried out. The predicted pressures at 0.02 compression are plotted for different numbers of discrete points and for different residuals. The number of discrete points and the residuals are normalized to the values of the original simulation.
Close modal

References

1.
Barbir
,
F.
,
2012
,
PEM Fuel Cells: Theory and Practice
,
Academic Press
,
Cambridge, MA
, pp.
10
12
.
2.
Yim
,
S.-D.
,
Kim
,
B.-J.
,
Sohn
,
Y.-J.
,
Yoon
,
Y.-G.
,
Park
,
G.-G.
,
Lee
,
W.-Y.
,
Kim
,
C.-S.
, and
Kim
,
Y. C.
,
2010
, “
The Influence of Stack Clamping Pressure on the Performance of PEM Fuel Cell Stack
,”
Curr. Appl Phys.
,
10
(
2
), pp.
S59
S61
.
3.
Ettouhami
,
M. K.
,
Atifi
,
A.
,
Mounir
,
H.
, and
Amadane
,
Y.
,
2019
, “
Numerical Simulation of Effect of Contact Pressure on Mechanical Behavior of Gas Diffusion Layers (GDL) and PFSA Membrane Assembly
,”
Mediterr. J. Chem.
,
8
(
6
), pp.
453
461
.
4.
Bograchev
,
D.
,
Gueguen
,
M.
,
Grandidier
,
J.-C.
, and
Martemianov
,
S.
,
2008
, “
Stress and Plastic Deformation of MEA in Fuel Cells: Stresses Generated During Cell Assembly
,”
J. Power Sources
,
180
(
1
), pp.
393
401
.
5.
Khetabi
,
E. M.
,
Bouziane
,
K.
,
Francois
,
X.
,
Lachat
,
R.
,
Meyer
,
Y.
, and
Candusso
,
D.
,
2024
, “
In-Situ Experimental Investigations to Study the Impact of Mechanical Compression on the PEMFC-Analysis of the Global Cell Performance
,”
Int. J. Hydrogen Energy
,
56
, pp.
1257
1272
.
6.
Balakrishnan
,
M.
,
Shrestha
,
P.
,
Ge
,
N.
,
Lee
,
C.
,
Fahy
,
K. F.
,
Zeis
,
R.
,
Schulz
,
V. P.
,
Hatton
,
B. D.
, and
Bazylak
,
A.
,
2020
, “
Designing Tailored Gas Diffusion Layers With Pore Size Gradients via Electrospinning for Polymer Electrolyte Membrane Fuel Cells
,”
ACS Appl. Energy Mater.
,
3
(
3
), pp.
2695
2707
.
7.
Irmscher
,
P.
,
Qui
,
D.
,
Janßen
,
H.
,
Lehnert
,
W.
, and
Stolten
,
D.
,
2019
, “
Impact of Gas Diffusion Layer Mechanics on PEM Fuel Cell Performance
,”
Int. J. Hydrogen Energy
,
44
(
41
), pp.
23406
23415
.
8.
Ha
,
T.
,
Cho
,
J.
,
Park
,
J.
,
Min
,
K.
,
Kim
,
H.-S.
,
Lee
,
E.
, and
Jyoung
,
J.-Y.
,
2011
, “
Experimental Study on Carbon Corrosion of the Gas Diffusion Layer in Polymer Electrolyte Membrane Fuel Cells
,”
Int. J. Hydrogen Energy
,
36
(
19
), pp.
12436
12443
.
9.
Atkinson
R. W.
, III
,
Garsany
,
Y.
,
Gould
,
B. D.
,
Swider-Lyons
,
K. E.
, and
Zenyuk
,
I. V.
,
2017
, “
The Role of Compressive Stress on Gas Diffusion Media Morphology and Fuel Cell Performance
,”
ACS Appl. Energy Mater.
,
1
(
1
), pp.
191
201
.
10.
El-Kharouf
,
A.
, and
Steinberger-Wilckens
,
R.
,
2015
, “
The Effect of Clamping Pressure on Gas Diffusion Layer Performance in Polymer Electrolyte Fuel Cells
,”
Fuel Cells
,
15
(
6
), pp.
802
812
.
11.
Lai
,
Y.-H.
,
Li
,
Y.
, and
Rock
,
J. A.
,
2010
, “
A Novel Full-Field Experimental Method to Measure the Local Compressibility of Gas Diffusion Media
,”
J. Power Sources
,
195
(
10
), pp.
3215
3223
.
12.
Carral
,
C.
, and
Mélé
,
P.
,
2018
, “
A Constitutive Law to Predict the Compression of Gas Diffusion Layers
,”
Int. J. Hydrogen Energy
,
43
(
42
), pp.
19721
19729
.
13.
Hoppe
,
E.
,
Janßen
,
H.
,
Müller
,
M.
, and
Lehnert
,
W.
,
2021
, “
The Impact of Flow Field Plate Misalignment on thë Gas Diffusion Layer Intrusion and Performance of a High-Temperature Polymer Electrolyte Fuel Cell
,”
J. Power Sources
,
501
, pp.
230036
.
14.
Carral
,
C.
, and
Mele
,
P.
,
2022
, “
Modeling the Original and Cyclic Compression Behavior of Non-Woven Gas Diffusion Layers for Fuel Cells
,”
Int. J. Hydrogen Energy
,
47
(
55
), pp.
23348
23359
.
15.
Norouzifard
,
V.
, and
Bahrami
,
M.
,
2014
, “
Analytical Modeling of PEM Fuel Cell Gas Diffusion Layers Deformation Under Compression: Part 1—Linear Behaviour Region
,”
ECS Trans.
,
61
(
11
), pp.
1
12
.
16.
Norouzifard
,
V.
, and
Bahrami
,
M.
,
2014
, “
Deformation of PEM Fuel Cell Gas Diffusion Layers Under Compressive Loading: An Analytical Approach
,”
J. Power Sources
,
264
, pp.
92
99
.
17.
Kleemann
,
J.
,
Finsterwalder
,
F.
, and
Tillmetz
,
W.
,
2009
, “
Characterisation of Mechanical Behaviour and Coupled Electrical Properties of Polymer Electrolyte Membrane Fuel Cell Gas Diffusion Layers
,”
J. Power Sources
,
190
(
1
), pp.
92
102
.
18.
Garcia-Salaberri
,
P. A.
,
Vera
,
M.
, and
Zaera
,
R.
,
2011
, “
Nonlinear Orthotropic Model of the Inhomogeneous Assembly Compression of PEM Fuel Cell Gas Diffusion Layers
,”
Int. J. Hydrogen Energy
,
36
(
18
), pp.
11856
11870
.
19.
Afrasiab
,
H.
,
Davoodi
,
K. H.
,
Barzegari
,
M. M.
,
Gholami
,
M.
, and
Hassani
,
A.
,
2022
, “
A Novel Constitutive Stress-Strain Law for Compressive Deformation of the Gas Diffusion Layer
,”
Int. J. Hydrogen Energy
,
47
(
75
), pp.
32167
32180
.
20.
Le Carre
,
T.
,
Blachot
,
J.-F.
,
Poirot-Crouvezier
,
J.-P.
, and
Laurencin
,
J.
,
2024
, “
Mechanical Response of Carbon Paper Gas Diffusion Layer Under Patterned Compression
,”
Int. J. Hydrogen Energy
,
50
, pp.
234
247
.
21.
Xiao
,
L.
,
Yin
,
Z.
,
Bian
,
M.
,
Bevilacqua
,
N.
,
Zeis
,
R.
,
Yuan
,
J.
, and
Sui
,
P.-C.
,
2022
, “
Microstructure Reconstruction Using Fiber Tracking Technique and Pore-Scale Simulations of Heterogeneous Gas Diffusion Layer
,”
Int. J. Hydrogen Energy
,
47
(
46
), pp.
20218
20231
.
22.
Zhang
,
Z.
,
He
,
P.
,
Dai
,
Y.-J.
,
Jin
,
P.-H.
, and
Tao
,
W.-Q.
,
2020
, “
Study of the Mechanical Behavior of Papertype GDL in PEMFC Based on Microstructure Morphology
,”
Int. J. Hydrogen Energy
,
45
(
53
), pp.
29379
29394
.
23.
Thiedmann
,
R.
,
Fleischer
,
F.
,
Hartnig
,
C.
,
Lehnert
,
W.
, and
Schmidt
,
V.
,
2008
, “
Stochastic 3D Modeling of the GDL Structure in PEMFCS Based on Thin Section Detection
,”
J. Electrochem. Soc.
,
155
(
4
), p.
B391
.
24.
Daino
,
M. M.
, and
Kandlikar
,
S. G.
,
2012
, “
3D Phase-Differentiated GDL Microstructure Generation With Binder and PTFE Distributions
,”
Int. J. Hydrogen Energy
,
37
(
6
), pp.
5180
5189
.
25.
Liu
,
J.
,
Kerner
,
F.
,
Schlueter
,
N.
, and
Schröeder
,
D.
,
2023
, “
Predicting the Topological and Transport Properties in Porous Transport Layers for Water Electrolyzers
,”
ACS Appl. Mater. Interfaces
,
15
(
46
), pp.
54129
54142
.
26.
Froning
,
D.
,
Gaiselmann
,
G.
,
Reimer
,
U.
,
Brinkmann
,
J.
,
Schmidt
,
V.
, and
Lehnert
,
W.
,
2014
, “
Stochastic Aspects of Mass Transport in Gas Diffusion Layers
,”
Transp. Porous Media
,
103
(
3
), pp.
469
495
.
27.
Ye
,
L.
,
Qiu
,
D.
,
Peng
,
L.
, and
Lai
,
X.
,
2024
, “
Conduction Mechanism Analysis and Modeling of Different Gas Diffusion Layers for PEMFC to Improve Their Bulk Conductivities via Microstructure Design
,”
Appl. Energy
,
362
, pp.
122987
.
28.
Yi
,
P.
,
Peng
,
L.
,
Lai
,
X.
, and
Ni
,
J.
,
2011
, “
A Numerical Model for Predicting Gas Diffusion Layer Failure in Proton Exchange Membrane Fuel Cells
,”
ASME J. Fuel Cell Sci. Technol.
,
8
(
1
), p.
011011
.
29.
Benz
,
F.
,
2024
, “
Mechanical Modeling of Carbon Fiber Paper Structures for PEM Gas Diffusion Layers Including Damage
,”
European PhD Hydrogen Conference (EPHyC) 2024
,
Ghent, Belgium
,
Mar. 20–22
, pp.
405
410
.
30.
Yiotis
,
A. G.
,
Kainourgiakis
,
M. E.
,
Charalambopoulou
,
G. C.
, and
Stubos
,
A. K.
,
2016
, “
Microscale Characterisation of Stochastically Reconstructed Carbon Fiber-Based Gas Diffusion Layers; Effects of Anisotropy and Resin Content
,”
J. Power Sources
,
320
, pp.
153
167
.
31.
Afrasiab
,
H.
,
Gharehhajloo
,
E. E.
, and
Barzegari
,
M. M.
,
2023
, “
Electrical and Mechanical Characterization of the Gas Diffusion Layer During Compression in PEM Fuel Cells
,”
Int. J. Hydrogen Energy
,
48
(
82
), pp.
31996
32010
.
32.
Reddy
,
J.
,
2007
, “
Nonlocal Theories for Bending, Buckling and Vibration of Beams
,”
Int. J. Eng. Sci.
,
45
(
2–8
), pp.
288
307
.
33.
Hinebaugh
,
J.
,
Gostick
,
J.
, and
Bazylak
,
A.
,
2017
, “
Stochastic Modeling of Polymer Electrolyte Membrane Fuel Cell Gas Diffusion Layers—Part 2: A Comprehensive Substrate Model With Pore Size Distribution and Heterogeneity Effects
,”
Int. J. Hydrogen Energy
,
42
(
24
), pp.
15872
15886
.
34.
Odaya
,
S.
,
Phillips
,
R.
,
Sharma
,
Y.
,
Bellerive
,
J.
,
Phillion
,
A.
, and
Hoorfar
,
M.
,
2015
, “
X-Ray Tomographic Analysis of Porosity Distributions in Gas Diffusion Layers of Proton Exchange Membrane Fuel Cells
,”
Electrochim. Acta
,
152
, pp.
464
472
.
35.
Virtanen
,
P.
,
Gommers
,
R.
,
Oliphant
,
T. E.
,
Haberland
,
M.
,
Reddy
,
T.
,
Cournapeau
,
D.
,
Burovski
,
E.
, et al
,
2020
, “
SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python
,”
Nat. Methods
,
17
(
3
), pp.
261
272
.
36.
Naito
,
K.
,
Tanaka
,
Y.
, and
Yang
,
J.-M.
,
2017
, “
Transverse Compressive Properties of Polyacrylonitrile (PAN)-Based and Pitch-Based Single Carbon Fibers
,”
Carbon
,
118
, pp.
168
183
.
37.
Mirdehghan
,
S. A.
,
2021
, “1—Fibrous Polymeric Composites,”
Engineered Polymeric Fibrous Materials
,
M.
Latifi
, ed.,
Woodhead Publishing
,
Sawston, UK
, pp.
1
58
.
38.
Schweiss
,
R.
,
Meiser
,
C.
,
Damjanovic
,
T.
,
Galbiati
,
I.
, and
Haak
,
N.
,
2016
,
SIGRACET Gas Diffusion Layers for PEM Fuel Cells, Electrolyzers and Batteries
,
SGL Carbon Group
,
Wiesbaden, GE
, pp.
6
7
. https://www.sglcarbon.com/data/pdf/SIGRACETWhitepaper.pdf
39.
Gigos
,
P.-A.
,
Faydi
,
Y.
, and
Meyer
,
Y.
,
2015
, “
Mechanical Characterization and Analytical Modeling of Gas Diffusion Layers Under Cyclic Compression
,”
Int. J. Hydrogen Energy
,
40
(
17
), pp.
5958
5965
.
40.
Barbier
,
C.
,
Dendievel
,
R.
, and
Rodney
,
D.
,
2009
, “
Role of Friction in the Mechanics of Nonbonded Fibrous Materials
,”
Phys. Rev. E
,
80
(
1
), p.
016115
.