Pore-Scale Investigation of Mass Transport in Compressed Cathode Gas Diffusion Layer of Proton Exchange Membrane Fuel Cells

Abstract

Proton exchange membrane fuel cells (PEMFCs) are considered a promising energy source in the field of transport and distributed power generation. Fundamental research into their key components is needed to improve PEMFC performance and accelerate commercialization. Binder addition and compression induced by assembly pressure can significantly change the microstructure of the gas diffusion layer and affect mass transport. A two-dimensional multicomponent lattice Boltzmann (LB) model considering the cathode electrochemical reaction was developed, and a GDL was reconstructed numerically and considering a binder structure. The effects of the binder and compression on mass transport and electrochemical performance within the GDL were investigated. The results showed that an increase in binder volume fraction led to more chain-like structures and closed pores that were unfavorable for mass transport. Compression increased the mass transfer resistance of the GDL in the region under the rib, leading to a decrease in oxygen concentration and local current density.

1. Introduction

A proton exchange membrane fuel cell (PEMFC) is a device that directly converts the chemical energy of fuel into electrical energy. With the advantages of a high fuel energy conversion rate, low noise, and zero emissions, PEMFCs can be widely used in vehicles, such as cars, airplanes, and trains, and stationary power stations [1]. The development of core components in PEMFCs limits the process of scaling up production [2]. The gas diffusion layer (GDL) is between the bipolar plate and the catalyst layer, which is a porous media layer made of carbon fiber as the backbone and carbonized resin material as the binder, with a thickness of about 100–300 μm [3]. GDL is the place for gas diffusion, water penetration, electron conduction, and heat transfer, and it also needs to have mechanical strength, corrosion resistance, and a long service life [4]. The electrochemical reaction process is an important link governing fuel cell performance, and optimizing the GDL structure to improve the reactive gas transport process is an effective way to improve the performance of PEMFCs [5].

The manufacturing process of carbon paper GDL includes (1) chopped carbon fiber, (2) manufacturing carbon paper, (3) impregnation with carbonized resin, (4) graphitization, (5) water-repellent treatment, and (6) GDL lamination to membrane electrode assembly (MEA). By impregnating carbon paper with a thermosetting resin, i.e., the binder, the material strength and porosity required for a GDL can be obtained [6]. Several algorithms have been developed to generate digitized GDL structures containing binders. The morphologically closed spherical structure element method was proposed by Daino et al. [7] to generate a binder, where equal-diameter spheres were filled into the pore space of the fiber skeleton, and the unfilled area was defined as the binder. Burganos et al. [8] suggested that a binder is present in the form of a web-like structure between fiber skeletons. Zamel et al. [9] considered that a binder fills the smallest pores preferentially. Equal circles tangent to a fiber skeleton were placed in the pore space, and the closed area was filled with a binder by Kakaee et al. [10]. Thiedmann et al. [11,12] divided a plane into blocks according to the distribution of carbon fibers, and random blocks were filled with a binder using either the filling method or the dilation method. Altogether, to present the wetting characteristics of a binder, a binder generated by the algorithms usually collected at the intersections of carbon fibers and filled small pores preferentially.

However, the presence of binders has been ignored in most studies, which may lead to a misunderstanding of mass transport processes within a GDL, and few scholars have carried out studies related to binders in GDLs. Nabovati et al. [13] found that an increase in binder and PTFE volume fractions at constant porosity resulted in GDLs having higher IP and TP permeability. While the binder addition increased the TP effective thermal conductivity of a GDL, the IP effective thermal conductivity was not affected by the binder [9]. A study of the water removal process within a GDL was carried out by Kakaee et al. [10], who found that a binder hindered water droplet removal, but it is worth noting that the authors did not state a constant porosity. The effective thermal conductivity and permeability of a binder to anisotropy were evaluated by Simaafrookhteh et al. [14]. The process of liquid water penetration into a GDL was simulated by Jinuntuyaa et al. [15]. The results showed that a change in the structure of the GDL by a binder made a serious impact on liquid water breakthrough characteristics. Besides that, scholars have focused more on the effect of a binder on GDL transport parameters such as permeability, effective diffusivity, effective thermal conductivity, and tortuosity [16,17,18].

As mentioned earlier, keeping the porosity of the GDL constant, an increase in binder content makes the number of carbon fibers decrease and, consequently, has a significant effect on the GDL microstructure. Evaluating the mechanism of the change in the microstructural properties of a GDL on mass transport at the cathode of fuel cells at the pore scale has become a popular research topic. The lattice Boltzmann method has become the dominant approach for pore-scale investigations due to its advantages in code realization, parallelization, treatment of complex boundaries, and freedom from fluid continuity assumptions. Therefore, LBM has been used by many scholars to study mass transport in porous media inside fuel cells. The degree to which different structural parameters affect the breakthrough velocity of liquid water within a GDL was evaluated by Lai et al. [19]. Oxygen transport in solid oxide fuel cells was studied by Wu et al. [20] and Shimura et al. [21], respectively. Yin et al. [22] developed a three-dimensional multiphase LB model to study water transport in a perforated GDL. The effect of a wavy cathode surface on mass transport and electrochemical reaction processes in a PEMFC cathode was investigated by Ashorynejad et al. [23,24]. The effects of activation overpotential, inlet and outlet pressure difference, rib–channel width ratio, and channel width on cathode reactive gas distribution and current density distribution were evaluated by Molaeimanesh et al. [25]. Molaeimanesh et al. [26] also simulated single multicomponent gas transport within a carbon cloth GDL and obtained the current density distribution on the surface of the catalyst layer. Chen et al. [27] explored the effect of a porous structure on material transport and electrochemical reactions and simulated the removal of liquid water in GDL. Structural variations due to parameters such as the fiber diameter, thickness, porosity, fiber orientation, etc. on the electrochemical reaction process within a GDL were also evaluated [28,29,30]. The changes in the microstructure of a GDL by a binder are indisputable, but the mechanism by which a binder affects a cathode electrochemical reaction and mass transport has not been reported.

The preload applied during PEMFC assembly will cause a GDL below the rib to deform in compression. The deformation of a GDL due to compression changes the mass transfer characteristics, so many scholars have focused on the transport parameters [31,32,33,34,35] and water management [36,37,38,39,40] of a compressed GDL. Anisotropic tortuosity and permeability of a compressed GDL were calculated by Froning et al. [30]. The effective diffusivity and effective thermal conductivity of a compressed GDL were investigated by Zhu et al. [35]. It was found that the effective gas diffusion coefficient decreased with increasing compression ratio, and the effective electrical and thermal conductivity coefficients increased with increasing compression ratio. The dynamic transport behavior of liquid water within a compressed GDL was investigated by Lee et al. [39], and the rib structure was found to have a significant effect on the water breakthrough point location. Ira et al. [36] simulated the effects of wettability heterogeneity and compression on liquid water transport within a GDL, with compression leading to the appearance of more discontinuous water clusters. Liquid water removal behavior within a compressed GDL was evaluated by Molaeimanesh et al. [37]. It is worth mentioning that there are few studies on multicomponent transport and electrochemical reactions within a compressed GDL.

Both a binder and compression significantly change the microstructure of a GDL, yet the effects of both on mass transport within a GDL are extremely poorly studied. In this study, GDL structures with binder volume fractions of 0, 10%, 20%, and 30% were generated using the stochastic reconstruction method, and a compression model of GDLs was established to generate a compressed GDL considering compression ratios of 0, 10%, 20%, 30%, and 40%. A two-dimensional, single-phase, multicomponent LB model that takes into account electrochemical reactions was developed to investigate the effects of a binder and compression on cathode mass transport. The distribution of oxygen and water vapor within the GDL and the current density distribution are exhibited and discussed. The present study contributes to a deeper understanding of the effects of a binder and compression on mass transport mechanisms within GDLs and contributes to the structural design of GDLs.

2. Methodology

2.1. Two-Dimensional Single-Phase Multicomponent Fluid Flow LB Model

In this paper, the multicomponent fluid flow phenomena of a PEMFC cathode assembled with an interdigitated flow field were investigated, and the geometrical structure is schematically shown in Figure 1. The PEMFC cathode involves multiple component flow and reaction processes: air enters the GDL from the inlet and is driven by the pressure difference between the inlet and outlet, the oxygen component reacts in the catalyst layer to form water vapor, and the water vapor is taken out of the outlet by being coerced by the gas mixture. The present study set the inlet air to consist of 21% oxygen and 79% nitrogen, thus involving the transport of a total of three components: oxygen, nitrogen, and water vapor.

The D2Q9 LB model of the single relaxation time collision operator was applied, with nine possible velocity directions for the particles at each lattice. Macroscopic quantities of the fluid such as density and velocity can be expressed by the distribution function, and the Boltzmann equation for iteration of the distribution function is simplified as follows [41]:

$$f_{k,i}(x+c_i\Delta t,t+\Delta t)=f_{k,i}(x,t)-\frac{\Delta t}{{\tau }_k}[f_{k,i}(x,t)-f_{k,i}^{eq}(x,t)]$$

(1)

where f and f^eq are the distribution function and equilibrium distribution function, respectively, subscript k refers to the different components, and c refers to the velocity vector. The equilibrium distribution function can be calculated as follows:

$$f_{k,i}^{eq}(x,t)=w_i{\rho }_k\left[1+\frac{c_{i}u}{{c_s}^2}+\frac{{(c_{i}u)}^2}{2{c_s}^2}-\frac{u^2}{2{c_s}^2}\right]$$

(2)

The c_s is the lattice sound velocity, and w represents the weighting factor for each velocity direction, where the values of c and w are given as follows:

$$\left[\begin{array}{ccc}\begin{array}{ccc}\begin{array}{cccc}\begin{array}{cc}{c}_1& c_2\end{array}& c_3& c_4& c_5\end{array}& c_6& c_7\end{array}& c_8& c_9\end{array}\right]=\left[\begin{array}{c}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{ccc}1& 0& -1\end{array}& 0& 1& -1\end{array}& -1& 1& 0\end{array}\\ \begin{array}{cccc}\begin{array}{cccc}\begin{array}{ccc}0& 1& 0\end{array}& -1& 1& 1\end{array}& -1& -1& 0\end{array}\end{array}\right]$$

(3)

$$w_i=\left\{\begin{array}{cc}1/9\hfill & i=1,2,3,4\hfill \\ 1/36\hfill & i=5,6,7,8\hfill \\ 4/9\hfill & i=9\hfill \end{array}\right.$$

(4)

ρ and u are the mass and velocity of the fluid, respectively, and can be calculated from the distribution function:

$${\rho }_k={\displaystyle \sum _{i=1}^9{f}_{k,i}(x,t)}$$

(5)

$$u=\frac{{\displaystyle \sum _k\frac{1}{{\tau }_k}}{\displaystyle \sum _{i=1}^9{f}_{k,i}(x,t)c_i}}{{\displaystyle \sum _k\frac{1}{{\tau }_k}}{\rho }_k}$$

(6)

τ is the relaxation time and is related to the LB viscosity of the fluid (ν_LB) as follows:

$${\tau }_k=0.5+3{\nu }_{k,LB}$$

(7)

The conversion between LB viscosity and physical viscosity (ν_Ph) is as follows:

$${\nu }_{k,LB}=\frac{{\nu }_{k,P\mathrm{h}}\Delta t}{\Delta x^2}$$

(8)

2.2. Boundary Conditions

Dry air (21% oxygen and 79% nitrogen) was set to flow in through the inlet, while a mixture of oxygen, nitrogen, and water vapor flowed out through the outlet, with Zou and He’s [42] constant pressure boundary conditions applied at both the inlet and outlet. At the inlet and outlet boundaries, f₄, f₇, and f₈ are unknown and can be calculated specifically as follows:

$$u_k=\frac{f_{k,1}+f_{k,3}+f_{k,9}+2(f_{k,2}+f_{k,5}+f_{k,6})}{{\rho }_k}-1$$

(9)

$$f_{k,4}=f_{k,2}-\frac{2}{3}{\rho }_k{u}_k$$

(10)

$$f_{k,7}=f_{k,5}+\frac{1}{2}(f_{k,1}-f_{k,3})-\frac{1}{6}{\rho }_k{u}_k$$

(11)

$$f_{k,8}=f_{k,6}-\frac{1}{2}(f_{k,1}-f_{k,3})-\frac{1}{6}{\rho }_k{u}_k$$

(12)

At the outlet, the total pressure was constant, but the densities of the components were unknown, so the partial pressures of the components were assumed to be the same as at the previous lattice of the outlet. Since the representative computational element was chosen, symmetrical boundary conditions were set at the left and right boundaries, which were implemented as follows:

Left boundary:

$$f_{k,1}=f_{k,3}, f_{k,5}=f_{k,6}, f_{k,8}=f_{k,7}$$

(13)

Right boundary:

$$f_{k,3}=f_{k,1}, f_{k,6}=f_{k,5}, f_{k,7}=f_{k,8}$$

(14)

No-slip bounce-back boundary conditions were employed on all solid surfaces. The boundary at y = 0 is the simplified catalyst layer and, therefore, needs to be set as a reactive interface where oxygen is consumed to produce water vapor. As the oxygen particles move to the surface of the catalyst layer, a fraction $K_{sr}^{LB}$ of oxygen is thought to be consumed, while the remaining oxygen rebounds into the pore space. $K_{sr}^{LB}$ is considered the electrochemical reaction rate in lattice units, and the conversion relation from the electrochemical reaction rate in physical units is as follows [43]:

$$K_{sr}^{LB}=\left(\frac{6k_{sr}\Delta t}{\Delta x}\right)/\left(1+\frac{k_{sr}\Delta x}{2D}\right)$$

(15)

where k_sr is the electrochemical reaction rate constant, and D is the diffusivity of oxygen. Based on the current density generated by the electrochemical reaction, k_sr can be calculated by the following equation:

$$k_{sr}=\frac{J}{4F{\rho }_{O2}}$$

(16)

where J is the current density at the surface of the catalyst layer, while F is the Faraday constant, and ρ_O2 is the oxygen density at the surface of the catalyst layer. According to the Butler–Volmer equation, the current density can be calculated as follows:

$$J=a_r{J}_{ref}\frac{{\rho }_{O2}}{{\rho }_{O,ref}}\left[\exp \left(\frac{{\alpha }_{f}F\eta }{R_{u}T}\right)-\exp \left(-\frac{{\alpha }_{r}F\eta }{R_{u}T}\right)\right]$$

(17)

where a_r is the roughness factor of the catalyst layer, J_ref is the reference current density, and α_f and α_r are the forward and reverse oxygen reduction transfer coefficients, respectively.

Thus, at the surface of the catalyst layer, the unknown distribution function is calculated as follows:

Oxygen:

$$f_{O_2,2}=(1-K_{sr}^{LB})f_{O_2,4}$$

(18)

$$f_{O_2,5}=(1-K_{sr}^{LB})f_{O_2,7}$$

(19)

$$f_{O_2,6}=(1-K_{sr}^{LB})f_{O_2,8}$$

(20)

Water vapor:

$$f_{H_{2}O,2}=\frac{2M_{H_{2}O}}{M_{O_2}}{K}_{sr}^{LB}{f}_{O_2,4}+f_{H_{2}O,4}$$

(21)

$$f_{H_{2}O,5}=\frac{2M_{H_{2}O}}{M_{O_2}}{K}_{sr}^{LB}{f}_{O_2,7}+f_{H_{2}O,7}$$

(22)

$$f_{H_{2}O,6}=\frac{2M_{H_{2}O}}{M_{O_2}}{K}_{sr}^{LB}{f}_{O_2,8}+f_{H_{2}O,8}$$

(23)

2.3. Reconstruction of GDL

2.3.1. Stochastic Reconstruction of Fiber Skeleton

The fiber backbone of a GDL typically consists of randomly stacked carbon fibers with diameters of 7–10 µm. The 2D GDL is usually reconstructed in two ways: either by obtaining slices from the 3D GDL [25,27,37] or by simplifying the carbon fibers into circles that do not overlap [23,24,36]. In the present study, the GDL has a width of 2000 µm and a thickness of 196 µm. Considering the need to add the binder into the GDL in the subsequent steps, the second reconstruction method of the 2D GDL was adopted, with all carbon fibers having a diameter of 7 µm. Keeping the porosity of all uncompressed GDLs (ε_GDL) to be consistently 0.78, the relationship between the porosity of the fiber skeleton (ε_f) and the binder volume fraction (BVF) can be expressed as follows:

$${\epsilon }_f=BVF(1-{\epsilon }_{GDL})+{\epsilon }_{GDL}$$

(24)

2.3.2. Addition of Binder

The material of the binder is usually carbonized resin; therefore, in this study, the density of the binder was considered the same as the carbon fibers, and the binder volume fraction was defined as the ratio of the volume occupied by the binder to the volume of solids within the GDL [7,31]. The morphologically closed spherical structure element method proposed by Daino et al. [7] was adopted to generate the binder in the following steps:

• Calculate the shortest distance d_c from each lattice to the carbon fiber;
• Select an initial value r_c as the circle radius, draw circles with radius r_c at the lattices where d_c ≥ r_c, and mark the area covered by the circles;
• The area not covered by the circles is marked as binder, and the binder volume is counted;
• Decrease r_c if the counted binder volume is greater than the preset value, increase r_c if it is less, and repeat steps (2) (3) until the preset binder volume is reached.

2.3.3. Compression Model

The preload force causes the GDL below the rib to deform in compression, and based on the experimental observation and simulation study by Jeon et al. [44], the deformation of the GDL in the region not below the rib was considered negligible. The compression ratio (CR) was defined as the ratio of the reduced GDL thickness due to compression to the initial GDL thickness. Compression ratios of 0, 10%, 20%, 30%, and 40% were considered in this paper. Further details on the compression model can be found in our previous work [45].

2.4. Computational Domain

The computational domain for this study is shown in Figure 1, with a size of 2000 × 220 µm² and 1 lattice unit representing 1 µm. The computational domain was taken from the interdigitated flow field, taking half of the inlet and outlet, the rib, the gas diffusion layer, and the catalyst layer, which was simplified to a thin layer at the bottom of the computational domain. Dry air flows in through the inlet, and water is difficult to condense in dry conditions, so it was assumed that all the water in the GDL is present as water vapor [25,27,28,30,43]. Carbon fibers were simplified as circles with the same radius of 3.5 µm. The overall porosity of the uncompressed GDL was kept constant at 0.78 and is shown in Figure 1 for GDL with a binder volume fraction of 20%; the binder is marked pink, and the carbon fibers are marked black. The distribution of the binder corresponds to the characteristics of preferentially filling small pores and presenting wettability. The calculation ends when the change in current density is less than 1 × 10⁻⁴ after an adjacent 10,000 step. All codes were written in MATLAB and run on a server equipped with two 24-core Intel Xeon Gold 6248R 3.0–4.0 GHz CPUs and 128 GB of RAM.

3. Results and Discussion

Based on the GDL stochastic reconstruction method of Section 2.3, computational domains were reconstructed for binder volume fractions of 0, 10%, 20%, and 30% and compression ratios of 0, 10%, 20%, 30%, and 40%. The two-dimensional LB multicomponent fluid flow and electrochemical reaction models developed in Section 2.1 and Section 2.2 were applied to simulate the flow of oxygen, nitrogen, and water vapor inside the cathode, and the concentration distributions of oxygen and water vapor, as well as the current density distributions, were obtained. The two-dimensional multicomponent LB model developed in this study has been validated in our previous studies [43,46]. The physical parameters used in the model are listed in

Table 1. Structural and operational parameters used in this paper.

Parameter	Value
Width of the GC	1000 µm
Width of the rib	1000 µm
GDL thickness	200 µm
GDL porosity	0.78
Operating temperature	353 K
Operating pressure	1.5 atm
Differential pressure between inlet and outlet	0.01 atm
Oxygen molar fraction at the inlet	0.21
Nitrogen molar fraction at the inlet	0.79
Dynamic viscosity of oxygen	2.34 × 10−5 kg/(m·s)
Dynamic viscosity of nitrogen	2.01 × 10−5 kg/(m·s)
Dynamic viscosity of water vapor	1.20 × 10−5 kg/(m·s)
Diffusivity of oxygen in the mixture	1.891 × 10−5 m2/s
Roughness coefficient	2000
Reference current density	1.3874 × 10−2 A/m2
Forward oxygen reduction transfer coefficients	0.5
Reverse oxygen reduction transfer coefficients	1
Binder volume fraction	0, 10%, 20%, 30%
Compression ratio	0, 10%, 20%, 30%, 40%

. The results of the investigation are discussed next.

3.1. Effect of Binder on Mass Transport within Uncompressed GDL

The distribution of oxygen concentration within the GDL at different binder volume fractions is demonstrated in Figure 2a. Observing the structure of GDLs at different binder volume fractions, the solids in the GDLs were more uniformly distributed when the BVF was low, such as when the BVF was 0 or 10%. When the BVF was high, such as 20 or 30%, on the one hand, the solids became more easily distributed in chunks or chains, and several lumps of solids tended to form a chain and led to the blockage of the gas transport throat. On the other hand, more binder filled the small pores, and connected binder appeared between multiple adjacent fibers, leading to the emergence of closed pores, which is very unfavorable for gas transport within GDLs. The oxygen concentration distribution showed a gradual decrease from the inlet to the outlet because oxygen was consumed at the catalyst layer. Oxygen-poor regions occurred within the GDL close to the outlet side, where oxygen was almost depleted, which can lead to oxygen starvation and the intensified polarization loss of the PEMFC being exacerbated. The water vapor distribution at different binder volume fractions is given in Figure 2b. Water vapor was generated in the catalyst layer and was carried by a mixture of oxygen and nitrogen toward the outlet so that there was an accumulation of water vapor in the area close to the outlet. In addition, at binder volume fractions of 20% and 30%, in the region below the inlet, a chain-like structure due to the increase in binder volume fraction was observed to be a clear boundary of the water vapor concentration distribution, indicating that the water vapor concentration distribution was affected by the GDL structure.

The distribution of oxygen and water vapor concentrations along the x-direction and the current density distribution on the surface of the catalyst layer is depicted in Figure 3. For ease of comprehension, the inlet, rib, and outlet are illustrated in the figure by dashed lines, indicating that the area illustrated is the GDL below that location. From Figure 3a,b, the distribution of oxygen concentration and water vapor concentration below the exit region shows an oscillating pattern. Oxygen was supplied in time below the inlet so that from x = 0 to x = 500 µm, the local current density was high, and the changes in oxygen and water vapor concentrations were not drastic. Due to the obstruction of the rib, the lateral transport of oxygen was weakened so that the supply of oxygen below the rib was gradually weakened, and a rapid decline in oxygen concentration and current density occurred. Unusual and sudden fluctuations in oxygen and water vapor concentrations in the area below the rib and outlet were due to the formation of closed pores. In Figure 3c, the obtained current density is compared to the results of Chen et al. [27] and Molaeimanesh et al. [25]. The reason for not being able to match exactly is due to differences in the microstructure of the GDL, but the trend and the range of the distribution of the current density are very close to each other, corroborating the applicability of the LB model developed in this study. For the binder volume fractions of 0, 10%, 20%, and 30% calculations, the average oxygen concentrations were 2.4070, 2.3984, 2.3184, and 3.0770 mol/m³, respectively, and the average water vapor concentrations were 15.0747, 15.0948, 15.2399, and 15.7054 mol/m³, and the average current densities were 1.6402, 1.4819, 1.3432, and 1.1062 A/cm². Compared with the no binder case, the current density decreased by 8.43%, 18.1%, and 32.55% in the cases with BVF = 10%, 20%, and 30%, respectively, and the higher the binder volume fraction, the more significant the decrease in current density. Cathode electrochemical performance decreased with increasing binder volume fraction, and the presence of block and chain structures resulted in a weaker mass transfer. However, this study is an analysis of electrochemical performance, and a more comprehensive judgment would require studies of other aspects, such as strength and resistance.

3.2. Effect of Binder and Compression on Oxygen Distribution

Oxygen concentration distributions for GDLs with binder volume fractions of 0, 10%, 20%, and 30% at different levels of compression are shown in Figure 4. Commonly, as the compression ratio increased, the oxygen-poor region represented by the blue moved toward the inlet, representing an expansion of the oxygen-poor region. At lower binder volume fractions (BVF = 0, 10%), there were not many closed pores due to compression. The pore connectivity decreased from 100% to 999.9939% and 99.8826%, respectively, and for the samples with BVF = 10%, closed pores were significant only at a compression ratio of 40%. For samples with BVF = 20% and 30%, closed pores were already present in uncompressed samples, and pore connectivity was 99.6933% and 99.8301%, respectively. After compression, the closed pores were further increased, and pore connectivity was reduced, which greatly reduced mass transport efficiency.

The average oxygen concentration along the x-direction within the GDLs for different binder volume fractions and compression ratios is given in Figure 5. Since no deformation of the GDL occurred below the inlet area, compression had little effect on the oxygen concentration in this area. As the rib intrusion into the GDL resulted in a blockage of the lateral oxygen transport, a portion of the oxygen was forced to flow downward, and the increase in the compression ratio resulted in a slight increase in the oxygen concentration below the inlet. A sharp decrease in the oxygen concentration in the region under the rib occurred. The increase in compression ratio decreased the porosity of the GDL under the rib, and the mass transfer resistance increased, which led to the difficulty in supplying the oxygen being consumed, and therefore, the oxygen concentration decreased with the increase in compression ratio. On the other hand, the increase in compression ratio led to the emergence of more closed pores, as shown in Figure 5b at CR = 40%, and the apparently closed pores at x = 630 μm gave a sudden increase in local oxygen concentration, but this is actually detrimental to oxygen transport. Sudden fluctuations in oxygen concentration occurred more frequently in samples with BVF = 30% and 40%, implying more closed pores.

3.3. Effect of Binder and Compression on Water Vapor Distribution

The average water vapor concentration along the x-direction within the GDLs for different binder volume fractions and compression ratios is given in Figure 6. Similarly, below the ribs and outlets were areas where water vapor was enriched. In the region below the inlet, for samples with different binder volume fractions, compression resulted in a slight decrease in water vapor concentration, which may be due to the fact that compression reduced the porosity of the GDL below the ribs, the gas flow rate below the inlet region was accelerated, and the water vapor was carried away faster after it was generated. As the electrochemical reaction rate remained constant, indicating that water vapor was being generated at a flat rate, compression had no effect on the concentration of water vapor in the GDL below the outlet. For the GDL below the rib, compression led to an increase in the concentration of water vapor, and the porosity of the GDL below the rib was enriched by the increase in gas transport resistance due to the compression being reduced, making it difficult for water vapor to be removed. The increase in compression ratio also resulted in a leftward shift in the location of the maximum water vapor concentration; for example, the locations of maximum water vapor concentration within the GDL with compression ratios of 0, 10%, 20%, 30%, and 40% were x = 1962 µm, 1465 µm, 1313 µm, 1172 µm, and 1017 µm, respectively, as an example of cases for BVF = 0. Untimely removal of water vapor may block the pathway for oxygen supply and have a negative impact on the performance of the PEMFC. Compression changes to the location of the maximum water vapor concentration should be taken into account to design a GDL structure to improve water–gas management. The sharp decrease in water vapor concentration in Figure 6b–d does not imply that water vapor is favorably excluded but is a phenomenon resulting from the closure of pores.

3.4. Effect of Binder and Compression on the Local Current Density Distribution

The surface current density distribution of the catalyst layer at different compression ratios is demonstrated in Figure 7. In the region of 0–500 µm, the local current density almost did not change due to compression. In the 1500 μm to 2000 μm region, on the other hand, a region of near-zero current density occurred due to the basic depletion of oxygen. In the region of 1500–2000 μm, a phenomenon of near-zero current density occurred because oxygen was almost exhausted. In the region below the rib, for samples with different binder volume fractions, it was observed that the higher the compression ratio, the lower the local current density. This suggests that the GDL structure under the rib needs to be paid attention to for a proper structural design to improve mass transfer capability.

The average current density for all samples is given in Figure 8. The current density decreased close to linearly with increasing compression ratio, where slight fluctuations may be the effect of structural differences. For the samples with BVF = 0, 10%, 20%, 30%, and 40%, the average current density decreased by 32.55%, 31.48%, 31.59%, 30.70%, and 30.85% at a compression ratio of 40% compared with uncompressed, which was at a similar level.

4. Conclusions

In this paper, a GDL was generated based on the stochastic reconstruction method, and a binder was added with the morphologically closed spherical structure element method. A two-dimensional single-phase multicomponent LB model was developed to investigate the effect of binder volume fraction and compression ratio on mass transport within the GDL, and more specific conclusions are summarized below:

• The increase in binder volume fraction led to more block and chain structures that blocked the gas transport throat and closed pores that are not favorable for gas transport;
• For uncompressed GDLs, the increase in binder volume fraction led to an increase in oxygen-poor and water-vapor-enriched regions and a decrease in average current density;
• An increase in compression ratio decreased the porosity of the GDL under the rib, and mass transfer resistance increased, which led to the difficulty in supplying the oxygen being consumed. Therefore, the oxygen concentration decreased with an increase in compression ratio;
• For samples with different binder volume fractions, the higher the compression ratio, the lower the local current density.