Abstract
Porous carbon fiber papers, which are formed by randomly distributed carbon fibers, are commonly used as gas diffusion layers of proton exchange membrane fuel cells. Several key properties of the gas diffusion layers, for example, air permeability, thermal conductivity, and electrical resistance, are inevitably controlled by the strain it experienced, which is always a nonlinear function of the applied compressive stress. In the present study, a logarithmic-type constitutive model is proposed to characterize the nonlinear stress–strain relation of carbon fiber papers under compressive loading. A hexahedral configuration of carbon fibers is adopted as a typical microstructure of carbon fiber papers. Contact and bending deformation of carbon fibers are accounted for with use of Hertz’s contact model and beam bending model, respectively. The logarithmic-type constitutive model is eventually obtained by solving a differential equation considering the variation of structural porosity. The stress–strain curves of the present model and experimental data are compared for a number of commercial gas diffusion layer samples. The comparisons indicate that the present model is capable of predicting the mechanical performance of carbon fiber papers, and may be useful for analytical and numerical analysis of gas diffusion layers and fuel cells.
Introduction
Proton exchange membrane fuel cells (PEMFCs) are considered as promising power sources for motor vehicle, stationary and portable equipment due to their high power density and environmentally friendly properties. 1 Gas diffusion layer (GDL), which is responsible for current collection, mechanical supporting, and gas and thermal diffusion, is one of the main components of PEMFCs. 2 Understanding the functionality of GDLs is essential for its design and optimization to meet various service conditions.
Typical GDLs are highly porous papers or cloth constituted by randomly distributed carbon fibers. 3 Numerous experimental and analytical studies4–11 indicated that several key properties of the carbon fiber papers are inevitably influenced by its deformation. For instance, the air permeability is significantly influenced by the porosity of carbon fiber papers,4,5 which decreases when the applied compressive strain is increased. The compressive stress can also change the contact electrical resistance of carbon fiber papers,6–8 and therefore influence the efficiency of energy output. Last but not least, the thermal conductivity of carbon fiber papers is also highly dependent on the magnitude of compressive stress.9–11 Consequently, it is imperative to take mechanical deformation into consideration to determine the electrical, thermal, and fluid performance of carbon fiber papers.
The mechanical behaviors of carbon fiber papers are usually analyzed by finite element (FE) models.12–16 One type of the commonly used FE models considers the carbon fiber papers and bipolar plates as elastic solids, which follow a linear elastic stress–strain or a data-based nonlinear stress–strain relationship. For instance, Zhou et al. 12 analyzed the water content in carbon fiber papers with use of a linear constitutive model. Using a data-based nonlinear constitutive model, Zhou and Wu13,14 investigated the deformation, porosity, and contact resistance of one piece of carbon fiber paper under the compression of two bipolar plates with smooth surfaces numerically. It is found that FE simulations with use of linear constitutive model of carbon fiber papers can only provide rough results 12 while more refined constitutive model of carbon fiber papers is indispensable if more accurate results are required.13,14 Other types of the FE models considered the GDLs and bipolar plates as elastic solids with non-smooth surface morphologies. Yoo and Jang 15 built a finite element model with micro geometric details to examine the deformation and contact resistance between GDLs and bipolar plates. By considering the carbon fiber papers as linear elastic layers, Lee et al. 16 studied the contact behavior between GDLs and bipolar plates with different surface morphologies. In general, the FE models with a data-based nonlinear stress–strain relation or detailed surface morphologies,13–16 however, can only be constructed by resorting to a costly and time-consuming testing process.17,18 Moreover, the transferability of the FE models13–16 is also a problem that cannot be addressed easily. An analytic constitutive model, which can predict the nonlinear stress–strain relationship of the carbon fiber papers accurately, is highly desirable to overcome the aforementioned difficulties.
With the aim of developing an analytic constitutive model for carbon fiber papers, Garcia-Salaberri et al. 19 proposed an empirical nonlinear stress–strain relationship by fitting the experimental data. However, Garcia-Salaberri’s model still requires considerable experimental data for curve fitting, which is very time-consuming in nature. Recently, Bahrami used representative volume element method with cubical unit cell to analyze the material properties (e.g. permeability, 20 thermal conductivity, 21 and elastic modulus) 2 of carbon fiber papers. Although Bahrami’s model is successful for taking the microstructures and deformations of carbon fiber papers into consideration, the Bahrami’s model, which is linear in nature, can only fit the experimental data within a small range. Therefore, a nonlinear analytical constitutive model, which can accurately predict the behavior of carbon fiber papers, is highly required. There is no doubt that a nonlinear analytic constitutive model will also be very helpful for analyzing the mechanical deformation, contact resistance, permeability, and conductivity of carbon fiber papers as well as the membrane electrode assembly.
In the present study, a logarithmic-type constitutive model of carbon fiber papers is proposed to describe the nonlinear stress–strain relation of carbon fiber papers under compressive loading. In our analysis, a hexahedral configuration of carbon fibers is adopted as a typical microstructure of carbon fiber papers. The effect of relevant geometric parameters, for example, the fiber diameter, the edge length, and the unit cell area of the microstructures, on the material properties of carbon fiber papers are analyzed in detail. The contact and bending deformation of carbon fibers are also accounted for with use of Hertz’s contact model and beam bending model. A logarithmic-type stress–strain relation is finally constructed by solving a differential equation of compressive stress and strain. For a number of commercial GDL samples, it is found that the proposed stress–strain curves are in good agreement with the experimental data. The strengthening behavior of structural modulus can be well explained by the effect of strain-induced porosity of carbon fiber papers.
Analytical model
Unit cell description
The carbon fiber papers concerned in the present study are highly porous papers constituted by randomly distributed carbon fibers. By observing the carbon fiber papers closely, it can be found that the carbon fiber papers are usually of quadrilateral shapes and formed by the carbon fibers with random edge lengths and areas (see Figure 1). Consequently, a typical unit cell of the carbon fiber papers is identified and depicted in Figure 2. Figure 2(a) shows a sketch of several randomly distributed carbon fibers at different levels of height. As a tailored part of Figure 2(a), Figure 2(b) shows a typical hexahedral unit cell of the carbon fiber papers in an isometric view, which consists of two horizontal fibers and two vertical fibers at different heights. The total height of the hexahedral unit cell is four times of the diameter of carbon fiber. Meanwhile, schematic illustrations of deformation modes of the carbon fibers are presented in Figure 2(c).

Top views of commercial GDLs, which reveal carbon fibers with random edge lengths: (a) SEM images, (b) detailed unit cell. 3

Schematics of the proposed geometrical model for the GDL: (a) Randomly distributed carbon fibers at different levels of height; (b) a typical hexahedral unit cell of carbon fiber papers; (c) contact deformation of a carbon fiber and bending deformation of a carbon fiber.
Contact deformation analysis
As shown in Figure 2(c), the horizontal carbon fibers will have an intimate contact with the neighboring vertical fibers when compressive stress is applied on the hexahedral unit cell. The nonlinear contact deformation is one of the main deformation modes of carbon fiber papers, which is an important source of the nonlinear mechanical responses of carbon fiber papers.
Table 1 presents the mechanical properties and geometric parameters of several commercial GDL samples. Based on these data, it can be found that the average edge lengths (84–136 μm) are much larger than the average fiber diameters (6.1–7.0 μm). Therefore, the contacting fibers within the unit cell can be modeled as infinitely long cylinders. Based on the Hertz’s contact theory, 22 the compressive load-displacement relation of the unit cell of carbon fiber papers can be calculated as
Mechanical properties and geometric parameters of commercial GDL samples.
GDL: gas diffusion layer.
where P is the applied compressive loading on the unit cell, u is the compressive displacement of the unit cell, and E, D, and v denote the elastic modulus, diameter, and Poisson’s ratio of carbon fiber, respectively. Then the average compressive stress and average contact-deformation-induced compressive strain of the unit cell can be formulated as
where σ is the average compressive stress of the unit cell of carbon fiber papers, εc is the average compressive strain induced by contact deformation, and Apore denotes the average area of the top surface of the hexahedral unit cell. Therefore, the relation between the compressive stress and contact-deformation-induced strain of the unit cell can be expressed in the form of
where
Bending deformation analysis
As shown in Figure 2(c), the carbon fibers crisscross at different heights. Consequently, each section of the fiber also acts as a three-point bending beam where the above and below contact points act as the support. The bending deformation is another main deformation mode of the carbon fiber papers.
Norouzifard and Bahrami 2 successfully obtained an elastic constitutive model of cubical unit cell of carbon fiber papers with use of the Euler-Bernoulli beam theory. Since the applied moment and deformation of bended fiber are independent on intersection angle between fibers, the Euler-Bernoulli beam theory can also be applied to crossed fibers in a hexahedral unit cell. Therefore, the relation between the compressive stress and bending-deformation-induced strain 2 of unit cell can be formulated as
where εb is the average compressive strain induced by bending deformation, l denotes the average edge length of unit cell, µ is the porosity of carbon fiber papers, and Cl and CD are the coefficient of variation of l and D, respectively. Let
Consequently, by taking stress as independent variable, a differential relation between the incremental compressive stress and bending-deformation-induced strain (i.e. dσ and dεb, respectively) can be described as
A logarithmic-type constitutive model
The total deformation of the aforementioned hexahedral unit cell is a superposition of the aforementioned contact deformation and bending deformation. Therefore, by combining equations (5) and (8), the following differential relation between dε and dσ of the unit cell can be established
However, the mechanical behavior of commercial GDL samples is usually different from the unit cells. Since the commercial GDL samples ordinarily contain large voids and pores, which may cause higher porosity of GDL samples than unit cell. As an example, the porosities of GDL unit cells and GDL samples are compared as follows. The porosity of the unit cell µu can be calculated approximately as
From the geometric parameters in Table 1, it can be estimated that the porosities of unit cell of SGL 24AA and SGL 25 AA samples are 83% and 86%, while the porosities of commercial samples are 88% and 92%, respectively. The difference in porosities verifies that there are large voids in commercial GDL samples, which may appear as long hanging fibers (much longer than the average edge length of unit cell). Under the assumption that the large voids are uniformly distributed, an efficiency factor η can be defined as
where the efficiency factor η is expressed as the ratio of GDL modulus to unit cell modulus.
Furthermore, the effect of strain-induced porosity on GDL modulus will be considered as well. The total volume of carbon fiber papers, which have a direct influence on the porosity and therefore the GDL modulus, is controlled by the applied compressive strain. Specifically, for the soft GDL samples, the nominal compressive strain can reach the value of 0.32 under the applied compressive stress of 1.5 MPa, thus the total volume and porosity of carbon fiber papers may change by an amount of 30% and 43%, respectively. By assuming that the Poisson’s ratio of GDL samples can be neglected, the strain-dependent efficiency factor η* can then be evaluated by the following equation
where ε is the nominal compressive strain of GDL. Furthermore, equation (12) indicates that the compressive modulus of carbon fiber papers will be increased when compressive strain is applied since under this circumstance ε < 0. Substituting equation (12) into equation (9), the strain-dependent differential relation of carbon fiber papers can be described as
Separating of variables in equation (13) yields
By integrating both sides of equation (14), the following logarithmic-type stress–strain relation is eventually obtained
Results and discussion
Verification of the logarithmic-type model
In this section, comparisons of different stress–strain relations with experimental results will be made. In order to verify the accuracy of the present logarithmic-type constitutive model, experimental data are used as benchmarks. The stress–strain data of five types of commercial GDL samples, for example, SEGRACET SGL 24AA, 25AA, 10 BA, TORAY TGP-H-60, and TGP-H-120, are obtained from the literature.2,23,24 Considering the fact that the SGL 10 BA, TGP-H-60, and TGP-H-120 GDLs contain PTFE or thermal-plastic resin, which will inevitably affect the mechanical behavior, these three types of GDLs are excluded, thus only SGL 24AA and 25AA GDLs are used for comparison in the following discussions.
According to the parameters listed in Table 1, several stress–strain curves of SGL 24AA and 25AA GDLs, including the experimental results, 2 bending model 2 and the present logarithmic-type model, are shown in Figure 3. In Figure 3, the black scatters are experimental data of SGL 24AA and 25AA GDLs, 2 where the symbols and error bars represent the average values and ranges of the compressive strain, respectively. The red and blue curves with symbols are corresponding to the present model and bending model, respectively. Notably, compressive stress and strain are used in Figure 3.

Comparison of stress–strain curves between experimental results, bending model, and the logarithmic-type model: (a) SGL 24AA GDLs; (b) SGL 25AA GDLs.
Figure 3 indicates that both the experimental data and that obtained from the suggested logarithmic-type model expose a similar strengthening trend. This can be explained by the fact that in both models the elastic modulus of carbon fiber papers all keep increasing as the applied strain increases. The bending model, however, demonstrates a linear stress–strain relation with constant modulus. Meanwhile, it is clear that the present logarithmic-type model, which could consider the strain-induced strengthening effect on modulus, has smaller deviation from the experimental data than the bending model. It is suggested that the strengthening effect is essential to predict the nonlinear mechanical behavior of carbon fiber papers. The above preliminary comparison indicates relatively lower error of the logarithmic-type model.
For clear verification, the relative errors of the bending model and the logarithmic-type model at various stress values are shown in Figure 4. Since 0.77 MPa is suggested to be an optimum clamping pressure, 25 typical stress values are selected as 0.4, 0.8, and 1.2 MPa to cover a relatively wide stress range. In Figure 4, the black scatters are the experimental average strain data of SGL 24AA and 25AA GDLs, 2 where the error bars represent the experimental strain range. The blue dotted line and dashed line are corresponding to the absolute value of relative error to the average strain value from the present model and bending model, respectively.

Relative errors of the bending model and the logarithmic-type model: (a) SGL 24AA GDLs; (b) SGL 25AA GDLs.
The relative errors in Figure 4 illustrate that the logarithmic-type model is more accurate in various conditions. For the SGL 24AA GDL samples, the relative error to the average experimental strain value of the bending model are −71% to −48% when the applied compressive stress is 0.4 to 1.2 MPa; at the same time, the relative errors of the logarithmic-type model are −32% to 8%, which is much lower than the bending model. Positive relative errors are used when the predicted value is higher than experimental results, while the negative relative errors suggest the opposite. Moreover, the relative error to the closest value within the strain range of the bending model are −69% to −44% in the same condition, while the corresponding relative errors of the logarithmic-type model are −26% to 1%, which is much lower as well. Especially, the relative errors to the experimental strain limit value of the logarithmic-type model are approximately equal to zero, when the applied compressive stress is 0.8 to 1.2 MPa. Therefore, the logarithmic-type model is proper for simulation of fuel cells with the optimum clamping pressure around 0.8 MPa. The accuracy of the logarithmic-type model can be explained by the fact that in this model the contact deformation and strain-induced strengthening effect under compressive loading have been considered appropriately.
For SGL 25AA GDL samples, similar global trend is revealed that the logarithmic-type model is more accurate. A better fitting is conducted for the logarithmic-type model to the experimental results, which the maximum absolute value of relative error reduces from 32% for SGL 24AA GDLs to 23% for SGL 25AA GDLs. However, the relative errors to the average strain value of the bending model is higher for SGL 25AA GDLs than that for SGL 24AA GDLs. The aforementioned difference between the bending model and the logarithmic-type model indicates that the logarithmic-type model based on a hexahedral unit cell is more accurate than the cubic-unit-cell-based bending model, when the porosity of carbon fiber papers is relatively high, for example, 92% for SGL 25AA GDLs. An explanation is that the carbon fiber GDLs with high porosity possess scattered edge length and varied fiber included angle which is close to hexahedral configuration rather than cubic configuration.
In addition, the logarithmic-type model predicts higher strain results than the experimental data, when the compressive stress is relatively high. The intense strengthening effect of GDL samples may be caused by the decreasing of edge length, as the unit cells can be squeezed and folded when the applied compressive stress is high. The abatement in edge length will increase stiffness of the carbon fiber papers, according to equation (6). The effect of edge length variation will be considered in the future study.
To summarize, the accuracy of the present constitutive model has been verified in this section. The comparison between the experimental data and the present model suggests that the logarithmic-type constitutive model does have the capability of predicting the overall nonlinear mechanical behavior of carbon fiber papers with a reasonable accuracy.
Effect of contact deformation
In this section, by neglecting the contact deformation, a simplified constitutive model of carbon fiber papers is first proposed to simplify the relatively complex constitutive model as shown in equation (15). Then the effect of contact deformation is examined by comparing the logarithmic-type model in the “Analytical model” section and the simplified model.
By ignoring the contact component, the differential strain–stress relation of GDL unit cell can be formulated as
Solving equation (16), the simplified constitutive model of carbon fiber papers can be expressed in the form of
For the average compressive strain, it can be eventually obtained
The simplified constitutive model in equation (18) can obtain the stress–strain relation of carbon fiber papers only after a small quantity of tests, since the whole stress–strain curve can be accounted for by one stress–strain pair. Consequently, the simplified model is preferable for engineering use once its accurate is verified.
The stress–strain pairs of the logarithmic-type model and the simplified model for SGL 24AA and 25AA GDLs are shown and compared in Figure 5, using the mechanical properties and geometric parameters in Table 1. For detailed comparison, more stress values are selected from 0.2 to 1.2 MPa to cover a large interval of loading conditions. The blue dashed line and dotted line in Figure 5 are the absolute value of relative errors of the simplified model and the logarithmic-type model to the experimental results, respectively.

Comparison between compressive strain values of the logarithmic-type model and the simplified model: (a) SGL 24AA GDLs; (b) SGL 25AA GDLs.
The results in Figure 5 illustrate that the simplified model may produce higher relative errors under low applied stress, yet lower relative errors under high applied stress. The relative errors of logarithmic-type model are in closed intervals [–32.6%, 4.1%] and [–31.6%, 4.4%], respectively, for SGL 24AA and 25AA GDL samples. As for the simplified model those figures are [–42.9%, 0.5%] and [–40.1%, 3.1%]. The lower strain values obtained by the simplified model than the logarithmic-type model can be easily explained, as the simplified model neglected the effect of contact deformation.
Specifically, the relative errors of the simplified model to the average experimental strain are −42.9% to −14.7% for SGL 24AA samples when the applied compressive stress is 0.2 to 0.8 MPa, while the relative errors are −3.9% to 0.5% when the applied compressive stress is 1.0 to 1.2 MPa, based on the experimental results in the “Results and discussion” section. The similar results are illustrated for SGL 24AA samples as well. Furthermore, the simplified model is more accurate than the logarithmic-type model when the stress is higher than 0.8 MPa.
It is suggested that the simplified model can be used to predict the strain of carbon fiber papers when the compressive stress is relatively high (i.e. higher than the suggested optimum clamping pressure). An explanation for the agreement between the experimental results and the simplified model is that the neglecting of contact deformation, which decreased the predicted strain, coincidentally fit the strengthening effect of carbon fiber papers as mentioned in the “Results and discussion” section.
The results in this section demonstrated clearly that the contact deformation contributes significantly to the total deformation of carbon fiber papers when the applied stress is relatively low. In addition, the contribution of contact deformation to total deformation generally decreases with the increasing of compressive stress. Specifically, the contribution of contact deformation may reduce to about 10% when the compressive stress is over 1.0 MPa. In summary, the simplified model is able to characterize the stress–strain relation of carbon fiber papers under relatively high stress.
Refinement of the logarithmic-type model
As mentioned in the “Results and discussion” section, actually there are extra PTFE or thermal-plastic resin in the SGL 10 BA, TGP-H-60, and TGP-H-120 GDLs. The existence of these ingredients will inevitably affect the mechanical behavior of these GDLs. For instance, the SEGRACET SGL 10 BA GDL samples contain 95% of carbon fiber and a 5% PTFE loading as hydrophobized substrate. Therefore, a modifying factor is necessary to demonstrate the influence of PTFE or thermal-plastic resin. Under the assumption that the PTFE, which makes scarce contribution to the modulus, is uniformly distributed in GDL samples, the modified logarithmic-type model can be described as
where
The experimental23,24 and predicted stress–strain curves of SGL 10 BA, TGP-H-60, and TGP-H-120 GDLs are shown in Figure 6, on the basis of the parameters in Table 1. The red and blue curves are the present logarithmic-type model and the bending model, respectively, 2 and the black scatters correspond to the experimental data of average compressive strain. Similarly, compressive stress and strain are used in Figure 6.

Comparison of stress–strain curves between experimental results and the modified logarithmic-type model: (a) SGL 10 BA GDLs; (b) TGP-H-60 GDLs; (c) TGP-H-120 GDLs.
As shown in Figure 6, strengthening behaviors are exhibited for both the experimental stress–strain curves and logarithmic-type model, which repeatedly indicates the importance of considering porosity-depend modulus. The relative error of the logarithmic-type model is relatively lower for TGP-H-120 GDLs, yet higher for SGL 10 BA and TGP-H-60 GDLs. Particularly, the ranges of relative errors for SGL 10 BA, TGP-H-60, and TGP-H-120 GDLs are in closed intervals [–55%, −26%], [–52%, −32%], and [–31%, −5%], respectively. The relative higher errors could be explained mainly by non-uniform distribution of the extra inclusions and inclusion-caused geometry change of unit cell. For example, the soft inclusions may locate at the contact points of carbon fibers, and therefore their existence will increase the compressive strain of the GDL samples. This well explains the negative relative errors associated with the modified model.
In summary, the results in this section illustrate that the modified logarithmic-type model is still capable of predicting the nonlinear mechanical behavior of carbon fiber papers with extra PTFE or resin, yet with some relative errors. A more refined model that can consider the effect of extra inclusions will be addressed in our future study.
Conclusion
In the present study, a logarithmic-type constitutive model is proposed to obtain the nonlinear stress–strain relation of carbon fiber papers under compressive loading. A hexahedral configuration is used as a typical microstructure of carbon fiber papers. Considering both contacting and bending deformations of carbon fiber papers, the logarithmic-type constitutive model is obtained by solving a differential equation. First, the comparison between experimental data and the present logarithmic-type model indicates that the present model is accurate in predicting the nonlinear mechanical performance of carbon fiber papers, and may be helpful for both analytical and numerical study of GDLs and fuel cells. The strengthening behavior of structural modulus is well explained by the effect of strain-induced porosity of carbon fiber papers.
Moreover, the effect of inclusions and contact deformation are also examined. It is found that the modified model is still capable of predicting the overall mechanical behavior of carbon fiber papers with inclusions, for instance, extra PTFE or resin, although with higher relative error. By neglecting the contact deformation, the simplified model can also be used to establish the stress–strain relation of carbon fiber papers under relatively high stress levels.
Footnotes
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors are grateful to the financial supports of National Natural Science Foundation (11702048, 11872138, 51875073 and 11772076), Program for Changjiang Scholars, Innovative Research Team in University (PCSIRT), 111 Project (B14013), Liaoning Natural Science Foundation Guidance Plan (20170520293) and Fundamental Research Funds for the Central Universities, China.
