Improved analytical modeling and mechanical characterization of gas diffusion layers under compression load

Gas diffusion layers (GDLs) are the most rigid layers in a 5‐layer or 7‐layer membrane electrode assembly (MEA) of a proton exchange membrane fuel cell. Therefore, in the fuel cell analysis, the mechanical properties of GDLs have a great impact on the stress distribution of the membrane as well as the performance of the whole cell. However, the mechanical properties of GDLs are not sufficiently studied. Nonlinear behavior of GDLs under cyclic compression is discussed rarely in literature. The existing model takes both constraints into consideration, but due to a geometrical oversimplification of the carbon paper microstructure, it shows some deviation, which cannot be overcome by selecting appropriate parameters. In this paper, the geometry of carbon paper microstructure has been reanalyzed. Based on the improved geometry formula, a modified model of carbon paper GDL is presented and verified by experimental data. Furthermore, both the experimental mechanical characteristics are achieved as support data for the improved model.

conditions, for instance, assemble clamping, membrane swelling, and vibration or impact. The representative models include hygro-thermal-mechanical model 3 and biaxial elastic-viscoplastic model 4 of Nafion the PEM from Boyce and Silberstein's team, poroelastic model 5 and thermo-mechanical model 6 of catalyst layer presented by Cho's team, and so on. They also investigated ex situ temperature and humidity aging effect on the tensile behavior of PEM, 7,8 and proposed a mechanical model of PEMFC electrodes considering effects of temperature and volume agglomerate fraction. 9 Among all these layers, GDL plays a key role in the operation as one of the key components among these layers influencing the whole fuel cell performance and durability. 10 First, to ensure the optimum operation of a PEMFC, reactants should continuously be transported and uniformly distributed by the GDL to the CL, where the electrochemical reactions occur. 11 GDL provides access of reactants. So homogeneous distribution and quick passages of gas should be provided by GDLs on both cathode and anode. The water and heat produced in the reaction should be effectively removed. Moreover, electrons should pass from the CL toward the BPP or in the reverse direction in anode and cathode sides, respectively, to generate an electric current. Then, high thermal and electrical conductivities are also crucial for GDLs. Furthermore, the contact pressure between the bipolar plates and the MEA in the cell assembly is important to the performance and durability of fuel cells. 12 The sufficient mechanical rigidity and strength are essential characteristics of GDLs to support MEA and reduce the effect of inhomogeneous compressive load given by ribs and channels of bipolar plates. 13 However, due to GDL's orthotropic and nonlinear characteristic as well as its complex compositions, it is very difficult to build a parameterized mechanical model of GDLs.
In the literature, there are three main types of GDLs, namely, carbon paper, woven, and felt. And the mechanical, electrical, and thermal properties have also been widely concerned, mainly focusing on carbon paper GDLs. [11][12][13] In a PEMFC, the compressive load on GDLs comes from two sources. One is the clamping force during the assembling process which can produce approximately 3 MPa stresses to maximize overall fuel cell performance by balancing the conductivities (electrical and thermal), mass transport resistance, and potential damage. 14,15 The other one is the swelling of PEM 3,16,17 during fuel cell operation which can induce large strains (up to 40%) of GDLs under the ribs, 18,19 the inhomogeneous pressure can be over 10 MPa. 20 Under compression, the microstructure 21 of GDLs would be changed in terms of pore network 22 and structural porosity. 23 These variations can result in a decrease of gas diffusivity and water transportation 24,25 as well as the electrical 26 and thermal 27 conductivities of GDLs. All the basic characteristics of GDLs, including diffusivity, water transport, thermal conductivity, and electrical contact resistance, are all dependent on compressive deformation. However, the mechanical response of GDLs under compression is nonlinear and nonisotropic. What is worse, the compressive response of GDLs depends on the load history. That is to say, due to the irreversible deformation, the compressive response of GDLs would vary significantly as the compression number of times and load changes. 20 Numerous models have been developed to describe the complex compression behavior of GDL. The nonlinear orthotropic model indicates that the linear isotropic models widely published in the literature may lead to incorrect predictions of contact stress distribution. 28 Another analytical model was built by Norouzifard et al 29 based on a unit-cell approach and statistical analysis. In this model, the strain-stress relationship was described by the geometrical parameters of the microstructure. But only linear response of GDL under compressive load was studied. Afterward, Gigos et al 30 developed a new model on the basis of the Norouzifard's model by converting geometrical parameters into macro-parameters with several assumptions in order to describe the nonlinear response of GDL under cyclic compression. However, in this model, some assumptions are over simplified, and the treatment of cyclic strain-stress relationship is not suitable for further study, neither component's influence of GDL in mechanical load nor entire fuel cell prediction analysis.
Therefore, a more accurate geometrical relationship is developed in the present paper using an analytical approach. And an improved model of the cyclic strain-stress relationship of a carbon paper type commercial GDL is proposed. Finally, the parameters are included as support data for the improved model. The model given in this paper can not only predict a more accurate mechanical response of GDLs under compression but also serve in microstructure study of GDLs to understand the influence of temperature, humidity, aging time, etc 2 | EXPERIMENTAL APPROACH

| Material
As known, a GDL is a thin layer of high porous carbon fibers which is arranged in different ways such as combining nonwoven fibers randomly by binder (carbon fiber papers), weaving fibers orthogonally (carbon fiber cloth), or matting fibers together (carbon fiber felt). The fibers are generally coated with polytetrafluoroethylene (PTFE) to improve hydrophobic characteristic of GDLs. In addition, a carbon microporous layer (CMPL) is often added to minimize the contact resistance between the GDL and catalyst layer, limit the loss of catalyst to the GDL interior, and help to improve water management as it provides effective water transport. In order to get a more accurate mechanical model of GDL fibers | 2801 by avoiding the influence of CMPL, the specimen used in the paper is a commercialized product obtained from Hesen electric co., Ltd., named HCP030N and HCP030P, which consists of only carbon fiber and binder. The detailed characteristics of HCP030N and HCP030P are listed in Table 1.

| Mechanical measurements
Measurements are conducted with an electronic universal testing machine (specific performance parameters are shown in Table 2) manufactured by China Mechanical Testing Equipment Co., Ltd. Since the mechanical behavior of GDLs is sensitive to environmental conditions, we added a temperature control system to ensure the accuracy of test data. Furthermore, due to the thin thickness of GDL, we are equipped with a set of auto-adaptive compression plates to guarantee the parallelism between the pair of loading plate and fixed plate, and the schematic diagram is shown in Figure 1. Ten pieces of Hesen ® GDLs are tested as one specimen to reduce measurement error. In each experiment, since there was no mechanical behavior modification of GDLs after five cycles, a specimen was compressed for five times with a compression amount up to 15 MPa (equivalently about 6 kN in terms of force), which is consistent with the results in literature. 30 The applied load is drawn in Figure 2 for 5 cycles. In the profile, curve-none shows the load applied without specimen, and it represents the deformation of the tester. In addition, curve-1st shows the load applied to a set of specimen for the first time. Similarly, curve-2nd to curve-5th shows the load applied to the same set of specimen for the second to the fifth time.
The load is applied with constant velocity (10 mm/min). Thus, curve-none ascends steeply due to the high rigidness of the test itself. Since the specimen is fluffy before compression, the load of curve-1st increases gently representing the gradual deformation of the specimen. After the first compression, the specimen is compact with the residual strain. So curve-2nd, curve-3rd, curve-4th, and curve-5th are nearly the same shape. GDL samples are cut by a handle die with a cutting size of 20 × 20 mm 2 . After the cut, specimens are measured by camera measure method of 10 μm to obtain size accuracy. Each test was repeated for 3 times and calculated averaged value as the measurement data.
As mentioned in literature, the stiffness moduli of GDL are close to that of the test machine, so the classical assumption of ignoring the deformation of test machine is not suitable for this study. During a compression test, the tested displacement is composed of the deformations of both test machine and specimens. To further reduce deformation test error, we measure the deformation of the test machine in each experiment and subtract the deformation of test machine from the obtained displacement-force curves of specimens. In this method, combining 10 pieces of GDLs as a set of specimen can reduce the test error to less than 5% which is close to the accuracy of test machine itself. Therefore, the experimental data are reliable.

| Data analysis
The displacement of the plate during one experiment of a specimen of GDL is shown in Figure 3 as scatter symbol curves. The symbol of "displacement-none" represents displacement-force relationship without the specimen during compression, which means the deformation of test machine under compression load. The lower scatter symbol curves in the sequence represent the first to the fifth compression of the same specimen in terms of displacement-force curves. By subtracting the machine deformation from 1st to 5th curves, the thickness of the specimen under multiple compressions can be obtained. The process to subtract the deformation of the tester is achieved by a data analysis and graphing software named OriginPro ® . Through a command, Analysis\ Mathematics\Simple Math, the displacement difference between "displacement-none" and "displacement-1st" under the same force can be calculated, where the difference value is the thickness of specimen under the first compression. Similarly, the thickness of specimen under the second, the third, the forth, and the fifth compression can be calculated in the same way. The thickness curves of specimen, for example, HCP030N, under cyclic loads are shown in Figure 3. The essential data for strain-stress curves can be abstract from these thickness curves with specimen area measurement. Figure 4 shows the strain-stress curves of the specimens under multiple compression loads up to 15 MPa. The loading shapes of the curves are similar to responses of GDLs from literature. The strain-stress curves of HCP030P contains 5% of PTFE, as shown in Figure 4B, show a slight soften phenomenon. A stabilization of the response has been observed, and the difference between the curves of fourth and fifth compression is negligible. The irreversible strain can be monitored, and the amount after the first compression is the major part of total irrepressible strain which is consistent with the stabilization phenomenon in this study as well as literature.

| Existing model discussion
In Norouzifard's model, 29 the stress-strain relationship of carbon paper GDL is given in Equation (1). where ε is the strain, σ is the stress, and E is the elastic modulus of carbon fibers in longitudinal direction. And the nondimensional term that represents the geometry structure is given in Equation (2).
where C and μ are the coefficient and random variable mean of variation, l, d, and A pore indicate the unit-cell length, fiber diameter, and pore area mean value, and p is the porosity of the sample, respectively.
Meanwhile, an equation, which links the mechanical strain and the porosity geometrically, 21 is described as Equation (3). where p 0 is the initial porosity of the GDL, and ε is the strain.
Gigos et al 30 represented the relationship between macro property (porosity p) and microstructure dimensions (unitcell length l, fiber diameter d) by making some simplifications and assumptions as follows.

Assumption 1
The GDL fibers distribution was simplified in orthotropic way regularly as shown in Figure 5 (a) and (b). Therefore, the pore area (A pore ), as shown in Figure 5 (c) can be expressed as Equation (4) Assumption 2 The GDL porosity, which is defined as the ratio of void volume (V pore ) over total volume (A total ), was assumed to be equal to the ratio of pore area (A pore ) over total area (A total ) as shown in Figure 5 (e) and Equation (5).
where the subscript YM means PA Gigos' model.

Assumption 3 A residual strain was subtracted from actual
strain to represent irreversible compression. Assumption 4 Along with compressive load increasing, GDL behaved respect to an apparent initial porosity (p app = p 0 ), which only changed when the load was higher than historical maximal load, in first compression.
Based on Equation (2) and the assumptions above, the calibrated GDL model from reference 30 is given as Equation (7).
where μ is a calibration parameter to represent real contact area ratio of carbon fibers in an actual GDL, and λ is the ratio between the apparent initial porosity and the real initial porosity. Assumption 2 is limited by the initial porosity value, only applicable to GDLs with high porosity (above 90%). In fact, as the main component of industrial GDLs, the porosity of carbon paper is generally around 75%. As can be seen from Figure 5D, the exact shapes of fibers and void are not simple cuboids. So there is an error (as shown in Figure 5F) between exact porosity and simplified expression (Equation 6). When the initial porosity of GDL is high, the error can be negligible. However, when the initial porosity is not high enough, which is the ration of l over d is not big enough, the error will be sufficiently large and nonnegligible. Moreover, this error can also affect the calibration of λ, which will be discussed in the following section.
The units of parameters mentioned above are listed as in Table 3.

| Advanced model development
In order to establish an GDL analytical model can be applied in the engineering field, we need to improve the equation linking macro property (porosity p) and microstructure dimensions (unit-cell length l and fiber diameter) while maintaining other assumptions from the existing models.
As assumed that GDL fibers are distributed regularly in an orthotropic manner, a unit cell of the microstructure of carbon paper can be abstracted as a repetitive equivalent element ( Figure 5A-C), and the total volume of the cuboid element can be easily given as Equation (8).
And the fibers in the element are shown in Figure 5D; the fiber volume can be added up by Equation (9).
According to definition, porosity can be expressed as Equation (10).
Considering Assumption 3 and 4 mentioned above along with the calibration method used in Equation (7), the expression of advanced compression model is obtained as Equation (14).
where μ and λ are the same physical significance as in Equation (7). Comparison between existing model and advanced model.

| Model calibration
There are 2 undetermined parameters μ and λ in the advanced model developed in the previous section. They have to be calibrated with the experimental data. The selected conditions are the maximum load point of experiment data and optimal curve fitting. From the averaged maximum load point (ε max = 0.49662, σ max = 15.00373 MPa), the relationship between μ and λ can be obtain as the expressions below.
We calibrated μ and λ in existing model by using Equation (15) and in advanced model by using Equation (16)

| RESULT DISCUSSION
It can be observed in Figure 6 that the advanced model presented in this paper shows a better agreement with experiment data. The curve plotted from existing model showed a higher curvature than experiment data. It leads to a low overlapping range between these two curves. As mentioned previously, the existing model has simplified the expression of porosity. Hence, it is limited to simulate GDLs with high porosity. If the porosity of GDL is not high enough, that is, Hesen ® HCP030N and HCP030P whose porosity is 75% studied in this paper, the error of porosity expression will be nonnegligible. The reasons are discussed as follows: From Equations (6) and (11), we can find that the function variables for simplified porosity p YM and advanced porosity expression p are the same one, namely, length-diameter ratio l d . As we have discussed, for different GDLs whose porosities are different, the average size of pores is different. Accordingly, the value of l d , which changes with pore size, is different. With a simplified model, the error of porosity expression will be too great to be acceptable. By accepting the structure assumptions of existing model, we built an idealized GDL whose pores are uniformly distributed as shown in Figure 5. The real porosity of this idealized GDL can be easily calculated and the relationship between real porosity, p real , and the value of l d can be easily obtained.
We define the error of p YM as Equation (17).
In order to have a visual view, we plot p YM , p real , and error of p YM along l d in Figure 7.
(15) As mentioned in previous section, when l d is big enough, or, equivalently, initial porosity is high enough (around 95% or higher), and the error of p YM is negligible (<5%), which is shown as the magenta dash lines in Figure 7. However, the GDLs used in industrial production are mostly made from carbon papers with lower porosity. What is worse, low porosity will lead to inaccuracy of p YM . For instance, the green dot lines in Figure 7 show that when porosity (p YM ) is 75%, the error is over 15%. Here we need to note that, for these two cases mentioned above, the porosity are calculated by Equation (6). If the advanced porosity expression Equation (11) is utilized, the error will be even bigger. The blue dash lines demonstrated in Figure 7 show that when porosity (p) is 75%, the error of p YM is in excess of 35%.
Although the existing of calibration parameter, λ can reduce the influence of the over simplification of p YM to some extent. It will lead to an inaccuracy of λ itself. What is worse, these inaccuracies will confuse researchers in further study, such as temperature/humidity/aging effect on each parameter of analytical model of GDL as well as mathematic expression of deformation mechanism of GDL under compression load, etc Under compression load as high as 15 MPa, the calibration parameter λ in existed model is only 0.714 as shown in Table 4. λ is the ratio between the apparent initial porosity and the real initial porosity, and it represents the percentage of the pores involved in the compression among total pores of certain specimen. It is obviously abnormal that only 71.4% pores were involved in compression test. In the advanced model presented here, λ equals to 0.949 under the same compression condition. It also proves that the advanced model can represent the actual situation more accurately.
We need to note that the mechanical analytical model of GDL is developed not only to represent the stress-strain relationship of GDLs under compression but also to analyze the mechanism as well as the variation of microstructure in a mathematic way. Hence, an accurate analytical model has a very positive effect in uncovering the mechanism of GDL compression behavior.

| CONCLUSION
In this study, an advanced analytical model of GDL was present based on accurate microstructure equations. The parameters were calibrated by macro compression tests. Compared with the model, the advanced model showed a better agreement with the data when initial porosity is as low as 75%. The compression behavior of certain type GDL was obtained and described accurately by the model represent here. Furthermore, it can also be calibrated to represent other types of carbon paper GDLs.
Besides a better accuracy, the calibrated parameters of the advanced model showed a better rationality as well. This makes it possible to interpret the reality of the compression deformation of a GDL. Moreover, the influence of humidity/ temperature/aging time on compression behavior can be studied quantificationally by defining the parameters as functions of variations mentioned above.
With these advantages, our model can serve in the mechanical simulation of PEMFC cells or stacks in order to get a more description of deformation and stress distribution of GDL as well as other layers in MEA. It can provide theoretical support for PEMFC analysis in durability, reliability, etc F I G U R E 7 Values of real porosity p real , simplified porosity p YM , and error of p YM along l d continuum