Investigation on Mechanical Properties of a Carbon Paper Gas Diffusion Layer through a 3 ‐ D Nonlinear and Orthotropic Constitutive Model

: The mechanical loads that gas diffusion layers (GDLs) withstand in polymer electrolyte membrane fuel cell (PEMFC) stacks are sensitive to the assembly and working conditions. The me ‐ chanical properties of GDLs mostly depend on their composition materials, microstructural charac ‐ teristics, operation conditions, etc. An accurate and comprehensive understanding of the mechani ‐ cal performance of GDLs is significant for predicting the stress distribution and improving the as ‐ sembly technology of PEMFC stacks. This study presented a novel 3 ‐ D nonlinear and orthotropic constitutive model of a carbon paper GDL to represent the material stiffness matrix with its com ‐ pressive, tensile, and shear properties. Numerical simulations were performed based on the 3 ‐ D constitutive model, and the proposed 3 ‐ D model was validated against the experimental data re ‐ ported previously. It is found that the simulation results of the 3 ‐ D constitutive model show a good agreement with the experimental results. Besides, the novel 3 ‐ D nonlinear and orthotropic model was applied in the overall stress simulation of a simplified PEMFC unit cell, compared to a conven ‐ tional 3 ‐ D linear and isotropic model, and the simulation results of the two models show a signifi ‐ cant


Introduction
As one of the most promising fuel cells, polymer electrolyte membrane fuel cells (PEMFCs) have attracted increasingly more attention owing to their excellent properties, such as low emission, noiselessness, high efficiency, high specific power, etc. A PEMFC unit [1] typically consists of a polymer electrolyte membrane or proton exchange membrane (PEM) in the middle, a catalyst layer (CL), a gas diffusion layer (GDL), and a bipolar plate (BPP) in the anode and cathode regions, respectively. Most commercial GDLs [2][3][4][5][6][7] consist of a substrate (made of carbon paper, carbon cloth, metallic foams, etc.) and a microporous layer (MPL) coated with carbon powder and polytetrafluoroethylene (PTFE). Located between a BPP and a CL, the GDL takes crucial responsibilities for PEMFCs, such as providing mechanical supports, collecting the current, reducing contact resistance, transporting gas and water, etc. Especially for carbon paper GDLs, they are widely applied in PEMFCs due to their outstanding electrical conductivity, corrosion resistance, and water management. Significantly, carbon paper GDLs have high porosity and permeability because their porous structures result in efficient passageways for water and gas across PEMFCs [8]. However, excessively high porosity can lead to a decline in the mechanical performance of carbon paper GDLs [9]. Once mechanical fractures or breakages occur, the internal carbon fiber connections change correspondingly, resulting in changes in the thermal and electrical conductivity [10][11][12], porosity and permeability [13,14], and the current density of fuel cells [15]. Besides, it has been studied that the existence of cracks in the MPL is beneficial to the transport properties [16] but detrimental to the mechanical durability of PEMFCs [17]. In general, mechanical failures in GDLs, such as fractures or breakages of carbon fibers in the substrate and cracks in the MPL, greatly influence the overall performance and mechanical durability of fuel cells. Studying the mechanical performance of GDLs is fundamental, which can guide the performance improvement of PEMFC stacks.
GDLs offer mechanical supports for fuel cells, and the mechanical properties of GDLs play a key role in influencing the internal stress distribution. According to the typical PEMFC assembly procedures, GDLs are highly subject to clamping loads in the throughplane (or thickness) direction. The compressive properties of GDLs have been highlighted. It has been reported that the fiber volume fraction, resin volume fraction, and porosity vary with the compressive loads applied on GDLs, which eventually affects their effective mechanical properties, such as longitudinal elastic moduli and shear elastic moduli [18]. Some investigations considered that GDLs were linear and isotropic materials without indepth studies, and the constant Young's modulus [19,20] was employed in the fuel cell stress simulation. In practice, GDLs have been experimentally characterized with nonlinear compressive properties [21][22][23]. Additionally, some researchers have made considerable efforts to develop reliable numerical models [24][25][26][27] to explain the nonlinearity in the compressive performance of GDLs.
In a real PEMFC device, GDLs under the gas channel ribs of BPPs suffer more complex and inhomogeneous stress besides compression [28,29]. Rapid changes of stress occur at the edges of ribs, leading to bending deflection, tensile breakage, and shear failure even with irreversible damage of GDLs. In addition to the compressive properties, study of the tensile and shear characteristics of GDLs is needed and important, as well. It is presented that the mechanical performance of GDLs under in-plane tensile and flexural loads is linear and isotropic [30] or anisotropic [31]. The shear behavior of GDLs has been investigated with either isotropic [32] or anisotropic [30] characteristics. Due to a lack of a comprehensive 3-D mechanical model, GDLs were generally assumed to be isotropic without the consideration of the real and intrinsic material properties during PEMFC mechanical simulation [33,34].
This study proposed a sophisticated 3-D nonlinear and orthotropic constitutive model of a commercial carbon paper GDL for the PEMFCs stress simulation and prediction. The novel 3-D constitutive model was comprehensively expressed by the material stiffness matrix, where the compressive, tensile, and shear numerical models were included, especially considering the microstructural characteristics, material properties, and assembly configurations in PEMFC stacks. The established 3-D constitutive model was validated against the correspondingly experimental data published elsewhere [31]. Significantly, the 3-D nonlinear and orthotropic constitutive model was embedded into a simplified PEMFC unit for the mechanical stress simulation, compared to the existing elastic and isotropic model of GDLs.

3-D Nonlinear and Orthotropic Constitutive Model of a Carbon Paper GDL
In the mechanics of 3-D materials, the generalized Hook's law is representatively used to state the deformation or strain caused by an arbitrary combination of stress. The mechanical stress vs. strain relationship of a 3-D material body is generally written in matrix form, where 36 stiffness matrix components are used to describe the material properties, as given here: 11 where σ is the stress vector, C is the stiffness matrix of a 3-D material, and ε is the strain vector. Owing to the strain energy density function and the symmetric stiffness and compliance matrices of materials, the material stiffness matrix can be expressed by 21 independent stiffness components, especially for anisotropic materials. As representative engineering materials, fiber-reinforced composites normally behave in an orthotropic manner [35,36]. In general, at least two symmetric orthogonal planes exist in an orthotropic material, which makes its material properties be directionally determined within each plane. Therefore, only nine independent variables are required to represent the stiffness matrix of orthotropic materials. In compliance with the generalized Hook's law, the 3-D constitutive model of an orthotropic composite material can be defined as: where Ex, Ey, and Ez are Young's moduli in the X, Y, and Z directions, respectively; Gyz, Gzx, and Gxy are the shear moduli in the YZ, ZX, and XY planes, respectively; and νxy, νyx, νyz, νzy, νzx, and νxz are Poisson's ratios in the corresponding directions. Carbon paper GDLs, as thin composite materials, are typically made of random carbon fibers in the substrate and coated with PTFE and carbon powder in the MPL. On the whole, the arrangement and architecture of carbon fibers and uniform coating can result in a quasi-planar and symmetric structure, which makes GDLs have orthotropic characteristics. In addition, the mechanical performance of GDLs has been experimentally and theoretically investigated with orthotropic characteristics [24,37,38]. Considering the actual assembly conditions of GDLs in PEMFCs, a GDL is located between a BPP and a CL, so GDLs have to withstand assembly loads in the through-plane direction. However, the real mechanical loads that GDLs suffer are more complex and inhomogeneous, including tensile and shear stress, besides compressive loads. In consequence, a 3-D constitutive model with nine independent stiffness components to describe the orthotropic characteristics of a carbon paper GDL is required for the accurate and reliable stress distribution of PEMFCs. Due to the porous structures of GDLs, Poisson's ratios are generally considered as zero [39,40]. Mathematically, Poisson's ratios of GDLs can be written as: Thus, a 3-D constitutive model of a GDL could be described as: Given the configurations of GDLs in fuel cells and the original size provided by the supplier, the direction of a commercial carbon paper GDL is defined as shown in Figure  1.Error! Reference source not found. In this study, σzz, the stress in the Z (thickness or through-plane) direction is expressed by its compressive characteristic; σxx and σyy, the stress values in the X and Y (or in-plane) directions, are represented through its tensile properties; and τyz, τzx, and τxy are stated with its shear performance in the YZ, ZX, and XY planes, respectively.

Compressive Model
In PEMFC stacks, GDLs mainly experience compressive loads in the through-plane direction due to the assembly and clamping process. To date, the compressive behavior of GDLs has been experimentally investigated with nonlinearity [21,[41][42][43]. Besides experimental investigations, theoretical approaches were also employed, such as curve fitting through a polynomial function [30] or piece-wise functions [24,44,45], to describe the nonlinear compressive stress vs. strain relationships of GDLs. Even though data fitting can be used for representation of the nonlinear compressive performance, it heavily depends on the experimental results without any explanations of the structural characteristics of materials. It has been presented that the microstructural parameters of paper-type GDLs, such as porosity and fiber diameter, play a significant role in the determination of the compressive stress-strain curves [46]. Considering the connections of carbon fibers in the substrate, a numerical model including geometrical parameters (like carbon fiber diameter, elastic modulus of the fiber, pore size distribution, etc.) was created based on the beam theory to state the nonlinear compression behavior of GDLs [25] up to a compressive pressure of 2 MPa, which is expressed by: where εc is the compressive strain, σc is the compressive stress, E is Young's modulus in the longitudinal direction of carbon fibers, and eff l d is the effective ratio of the length over the diameter of a carbon fiber model. In view of carbon fiber arrangements, the porous characteristics, and the relationship between porosity and compressive strain of GDLs, P. A. Gigos et al. [26] developed the analytical model expressed in Equation (6), and applied the modified model in a wider pressure with a maximum cyclic compressive load up to 12 MPa. The developed compressive model was stated as: where p0 is the initial porosity of a pristine GDL, and μ and λ are calibration parameters to present the structural differences between a real GDL and the unit cell model, and the ratio between the apparent initial porosity and the real initial porosity, respectively. Subsequently, Yang Xiao et al. [27] optimized the relationship between porosity and microstructural parameters, such as the length (l) and diameter (d) of carbon fibers in the unit cell model listed in Equation (6): by the description of: Thus, with a more precise expression of the porosity in Equation (9) instead of Equation (8), the analytical compressive model of a GDL can be further deduced by: However, in reality, it can be clearly observed in Figure 1b that carbon fibers are bonded by the binder. Even though the binder occupies a lesser part among the total volume of a GDL unit cell, compared to carbon fibers, it still influences the effective porosity. Considering the effects of the binder on porosity, the porosity can be further improved with higher accuracy by the expression of: where ζ (ζ<1) contributes the influence of the binder to the effective porosity of a GDL unit cell structure. Correspondingly, the compressive model can be improved by: In both two representative compressive models of GDLs in Equations (7) and (10), λ was defined as a constant to represent the apparent initial porosity and the real initial porosity of GDLs to be filled the available porosity volume during compression. However, the porosity is geometrically proportional to the mechanical strain [47] when GDLs are exposed to compressive loads. Thus, λ can be described in a more accurate way: where k and c are fitting constants. Putting Equation (13) into Equation (12), the improved nonlinear compressive model of a carbon paper GDL can be eventually expressed by:

Tensile Model
Different from the nonlinear compressive behavior of GDLs, it has been proved that the mechanical performance of GDLs in the in-plane directions (X and Y) behaves in a linear and isotropic manner [30]. However, due to the preferential alignments of carbon fibers resulting from the original fabrication process [48], the in-plane mechanical properties have been theoretically and experimentally confirmed to be elastic and anisotropic [31,32,49]. In summary, the in-plane tensile strain of GDLs was proportional to their tensile stress regardless of the isotropy or anisotropy. In this research, the in-plane tensile model of carbon paper GDLs with the consideration of anisotropy is given as: where σxx and σyy are the tensile stress in the X and Y directions, εxx and εyy are the tensile strain in the X and Y directions, and Ex and Ey are the tensile moduli in the X and Y directions.

Shear Model
Compared to the compressive properties, the shear characteristics of GDLs seem to take less responsibility for providing the mechanical support for fuel cells, and the studies of the shear performance are very limited. However, in PEMFC stacks, the shear performance of GDLs takes a vital role in influencing the local stress distribution. Especially under the gas channel ribs, the shear failure of GDLs easily occurs with irreversible cracks of carbon fibers, resulting in the large-scale path of transporting gas and water, which would heavily affect the performance of fuel cells. To the authors' knowledge, for carbon paper GDLs, such as thin and brittle composites, it is extremely difficult to achieve their net shear characteristics due to the lack of reliable measurements. Great efforts have been made to predict and measure the shear performance of GDLs. The shear characteristics of GDLs have been theoretically studied with isotropy [32] or anisotropy [30], and the shear moduli are constants regardless of isotropy or anisotropy. In addition to theoretical methods to predict the shear performance of GDLs, it can be experimentally characterized by referring to the punch-and-die method [50][51][52]. The illustration of the shear property measurement of a GDL is illustrated in Figure 2. Thus, the shear strain of a GDL specimen can be defined as: where γ is the shear strain, Δl is the gap, and Δd is the deformation of a GDL in the vertical direction. It has been experimentally measured that the deformation of a GDL sample (Δd) was proportional to the punch load (F) in characterizing the shear performance of GDLs [31], which is given as: The shear stress can be calculated as: where τ is shear stress and A is the cross-section area of a specimen. Combining Equations (16)-(18) together, the shear stress can be calculated as: where C is a constant related to the specimen size, material properties, and the setting gap during the shear test based on the improved punch-and-die measurement. Hence, the shear moduli Gyz and Gxz (in YZ and XZ planes) can be derived as: Due to the limitation of the experimental measurement, the shear modulus Gxy in the XY plane was employed from [53] for the 3-D constitutive model simulation.

Validation of the 3-D Nonlinear and Orthotropic Constitutive Model of a GDL
In Equation (5), it is found that a 3-D material matrix of an orthotropic material, including its compression, tension, and shear descriptions, is typically employed to state the overall mechanical performance. In this research, a comprehensive 3-D nonlinear and orthotropic constitutive model of a carbon paper GDL was established by combining the compression, tension, and shear numerical models. Next, we present the validation of the proposed 3-D constitutive model by comparing the model results with experimental data published in [31]. A commercial carbon paper GDL (JNT-30-A1, Hwasung, South Korea) was employed in this study. Detailed material specifications and property parameters of JNT-30-A1 GDL and the novel 3-D constitutive model are listed in Table 1.  Figure 3 shows the compressive performance of a JNT-30-A1 GDL in the throughplane direction by experimental and numerical results. It can be seen that the GDL exhibits distinctly nonlinear characteristics. Thus, the simplification and assumption that the compressive properties of a carbon paper GDL were constant [19,20] were not accurate or acceptable. The improved compressive numerical model of the carbon paper GDL was verified by the accordingly experimental data, which confirms that the numerical data match the experimental results very well, as given in Figure 3Error! Reference source not found.. In conclusion, instead of a linear and isotropic model, the nonlinear compressive numerical model can better present the actual compressive characteristics of GDL, and be reliably employed to predict the stress distribution of PEMFCs with higher accuracy.  The GDL has different tensile moduli in different directions, which proves that it behaves in an anisotropic manner. Further, that the tensile performance of carbon paper GDLs was directly treated as isotropic [30] is not reliable. Importantly, even though the in-plane mechanical characteristics of the GDL are anisotropic, it is linear and elastic (about 720 MPa in the X direction and 420 MPa in the Y direction, [31]), and entirely different from its compressive properties in the through-plane direction. In general, the numerical results of the in-plane tensile numerical model show a good agreement with its corresponding experimental results, as shown in Figure 4. The in-plane linear and anisotropic tensile model can be implemented in the 3-D mechanical stress simulation of PEMFCs.  In the current, the shear properties of GDLs were weakly investigated due to the lack of reliable experimental measurements. The shear performance of GDLs was theoretically characterized as isotropic and linear [32], or anisotropic and linear [30]. This study established the shear numerical model of a carbon paper GDL based on an improved punchand-die measurement. The shear characteristics of a JNT-30-A1 GDL in the ZX and YZ planes with numerical data and experimental results are given in Figure 5. Different from the achievements in the publications [30,32], the shear performance of the GDL is confirmed, which displays strong nonlinearity, as well as slightly anisotropic characteristics. Besides, the shear numerical results fit the experimental results pretty well. As one part of the 3-D nonlinear constitutive model of a GDL, the nonlinear and anisotropic shear numerical model can be used for the accurate 3-D stress prediction of PEMFCs. In summary, the comprehensive 3-D nonlinear and orthotropic constitutive model of a carbon paper GDL, including the nonlinear compressive model, in-plane linear and anisotropic tensile model, and nonlinear and anisotropic shear model, was validated by the comparison with the related experimental data. It is concluded that the simulation results of the presented 3-D model show quite good agreement with the corresponding experimental results with higher coefficients of determination as listed in Table 2. Compared to a traditional linear and isotropic model, the proposed 3-D nonlinear and orthotropic model is more closed to its intrinsic and real material properties where the microstructural characteristics and stress-bearing states are considered.

Application of the 3-D Nonlinear and Orthotropic Constitutive Model of a GDL
In addition, the proposed 3-D nonlinear and orthotropic constitutive model of a carbon paper GDL, compared to a traditional linear and isotropic model, was applied in a simplified PEMFC unit cell, which was subject to an assembly load for the mechanical stress distribution analysis. The material specifications of the components used in the simplified PEMFC unit are summarized in Error! Reference source not found. 3. All the components (such as BPP, GDL, CL, and membrane) were assumed to be isotropic and independent of temperatures in this study. According to the symmetry of the PEMFC unit structure, a simplified PEMFC unit cell with detailed size and direction descriptions was established as given in Figure 6a, which was employed for the mechanical stress simulation. In order to simplify the simulation process, the symmetry boundary condition was defined as shown in Figure 6b, where a uniformed load (1 MPa in this study) was applied on the top of the BPP to stimulate the fuel cell clamping load. The material properties of all the remaining components (such as the BPP, CL, and membrane), as well as the boundary conditions, had the same values as listed in Table 3 under the simulation process for the comparison with the proposed 3-D nonlinear and orthotropic model and a traditional linear and isotropic model of GDLs. Subsequently, the linear and isotropic and the 3-D nonlinear and orthotropic constitutive models of GDLs were respectively implemented into a simplified PEMFC unit cell for the mechanical stress simulation, the corresponding simulation results of which are shown in Figure 7. As given in Figure 7Error! Reference source not found., the maximum stress concentration of the two models (the linear and isotropic and the 3-D nonlinear and orthotropic constitutive models of GDLs) occurs at the same location where the gas channel ribs connect to the GDL, similar to the findings in [24]. Significantly, the magnitude of the maximum stress obtained by the two models of GDLs is quite different from each other. Even though the maximum stress concentration happens at the same place, the value of the maximum stress in the 3-D nonlinear and orthotropic model is nearly 10 times larger than that of the linear and isotropic model. In summary, the numerical models closer to their actual material properties should be characterized comprehensively and thoroughly, and the nonlinearity and orthotropy in the mechanical performance of GDLs must be taken into consideration for a reliable stress simulation and prediction of fuel cells.

Conclusions
In this work, a novel 3-D nonlinear and orthotropic constitutive model of a carbon paper GDL, including the compression, tensile, and shear numerical models, was presented by a detailed and comprehensive description of its material stiffness matrix. Three types of numerical models were elaborately utilized and developed, considering the microstructure of carbon paper GDLs, the intrinsic properties of composite materials, and the actual configurations of PEMFC stacks. In the compressive numerical model, nonlinearity was described in the through-plane direction. In the tensile numerical model, the anisotropic and linear properties were stated in the in-plane (X and Y) directions. In the shear numerical model, anisotropy and nonlinearity were included in the ZX and YZ planes. All the developed numerical models were verified through the correspondingly experimental data with high coefficients of determination.
Additionally, compared to a conventional linear and isotropic model, the proposed 3-D nonlinear and orthotropic constitutive model of a carbon paper GDL was implemented in a simplified PEMFC unit cell for the mechanical stress simulation. It is concluded that the maximum stress magnitude of the 3-D nonlinear and orthotropic constitutive model is almost 10 times larger than the magnitude of the maximum stress of the linear and isotropic model, even though the maximum stress occurs at the same location. Overall, the stress distribution in PEMFCs is strongly affected by the internal components and material models. For reliable stress simulation and prediction of fuel cells, the nonlinearity and orthogonality in the mechanical performance of GDLs must be considered.