A Mechanical Model for Carbon Paper Gas Diffusion Layers for Proton Exchange Membrane Fuel Cells Including Fiber and Binder Failure

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.

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.

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 R_Domain and a height of H_Domain. 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 R_Domain. 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 H_Domain. A sketch of the domain and the randomly distributed fibers is shown in Fig. 1(a).

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 D_F. The total number of fibers is adjusted such that the porosity of the constructed structure matches the real GDL porosity without the binder.

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.

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.

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].

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.

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.

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].

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 E_F = E_B = 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 D_F = 7 µm [15]. The cutoff distance for the binder distribution is assumed to be r_max = 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 V_F = 0.1, and the binder volume ratio is assumed to be V_B = 0.02. A cylindrical domain with a radius of R_Domain = 600 µm and a height of H_Domain = 110 µm is used for the simulation.

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.

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.

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.

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 E_F and the binder E_B material, the breaking point of the binder $εBmax$⁠, and the spatial extent of the binder r_max 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.

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.

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.”

There are no conflicts of interest.

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

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.