Pore-Scale Investigation of Coupled Two-Phase and Reactive Transport in the Cathode Electrode of Proton Exchange Membrane Fuel Cells

A three-dimensional multicomponent multiphase lattice Boltzmann model (LBM) is established to model the coupled two-phase and reactive transport phenomena in the cathode electrode of proton exchange membrane fuel cells. The gas diffusion layer (GDL) and microporous layer (MPL) are stochastically reconstructed with the inside dynamic distribution of oxygen and liquid water resolved, and the catalyst layer is simplified as a superthin layer to address the electrochemical reaction, which provides a clear description of the flooding effect on mass transport and performance. Different kinds of electrodes are reconstructed to determine the optimum porosity and structure design of the GDL and MPL by comparing the transport resistance and performance under the flooding condition. The simulation results show that gradient porosity GDL helps to increase the reactive area and average concentration under flooding. The presence of the MPL ensures the oxygen transport space and reaction area because liquid water cannot transport through micropores. Moreover, the MPL helps in the uniform distribution of oxygen for an efficient in-plane transport capacity. Crack and perforation structures can accelerate the water transport in the assembly. The systematic perforation design yields the best performance under flooding by separating the transport of liquid water and oxygen.


Introduction
Hydrogen energy is recognized as one of the most promising approaches for achieving carbon neutrality [1][2][3]. Proton exchange membrane (PEM) fuel cell can convert hydrogen energy into electricity in an efficient and clean manner [4], and it has been widely applied in fuel cell vehicles. To ensure a smooth operation of PEM fuel cells, liquid water inside an electrode, which significantly affects the performance, should be managed properly, particularly under a high-current condition [5,6]. In particular, water produced from the cathode catalyst layer (CL) flows through the microporous layer (MPL) and gas diffusion layer (GDL) to the flow channel. A certain amount of water in the CL can help the membrane to stay hydrated and ensure proton conductivity, but excessive water could cause flooding in the electrode, which greatly increases the transport resistance and leads to severe performance degradation [7,8]. Therefore, the microstructure of the GDL/ MPL assembly should be delicately optimized for a good water management strategy, which is based on in-depth insight into the dynamic behavior of liquid water in the cathode electrode and its influence on oxygen transport and performance.
To resolve the water distribution inside the electrode, X-ray tomography [9,10], neutron imaging [11], and nuclear magnetic resonance [12] are widely used in experiments. Ko et al. [13] investigated the effect of the GDL with gradient porosity on the water transport in PEM fuel cells via X-ray visualization and found that the GDL with medium gradient porosity yields a good performance due to the low transport resistance. Lee et al. [14] also employed X-ray tomography to observe the influence of the MPL on water transport, and the results reveal that the presence of the MPL leads to less liquid water stuck at the MPL/CL interface. Markötter et al. [15] studied the transport path of liquid water in the GDL and MPL via in situ imaging technology and found that liquid water mainly flows through the large pore in the GDL and cracks in the MPL. Gerteisen et al. [16] investigated the effect of vertical laser perforations in the GDL on water transport and found that these perforations can alleviate flooding by separating the transport of the reactant and 1 3 liquid water. The mentioned experimental techniques are very useful in revealing the dynamic behavior of liquid water in the cathode electrode. However, conducting an in situ observation of the two-phase flow in a porous microstructure is still difficult. Furthermore, the spatial and temporal resolutions still need to be further developed to resolve the detailed dynamics of water transport, particularly in the MPL and CL.
To complement the experimental issues, numerical modeling has been widely employed to study the two-phase flow in the electrodes of PEM fuel cells [5,17]. The volume-of-fluid (VOF) method [18], pore network (PN) model [19], and lattice Boltzmann model (LBM) [20] are currently the most popular numerical methods. The VOF method is widely applied to study the two-phase flow in GDLs [21,22]. However, the VOF method is limited in modeling the nanoscale flow where the assumption of fluid continuity is not validated. As the macroscopic fluid dynamic results from the interaction among various fluid particles, a mesoscale model is thus expected to describe the two-phase flow in a fundamental view [23]. Recently, Cetinbas et al. [19] proposed a PN model [24] to investigate the water evolution process in GDL/MPL assemblies, with cracks in the MPL addressed. The results agreed with those in the work by Markötter et al. [15] and showed that liquid water tends to transport through the cracks in the MPL. However, the PN model as a rule-based model can only give a general prediction of the transport phenomena [25]. LBM, rooted in kinetic theory, can resolve the mesoscale phenomena with an acceptable computation cost [26,27] and has been widely applied to study the two-phase flow in porous media. Kim et al. [28] developed an LBM to resolve the water distribution in GDL/MPL assemblies and showed that a thick MPL leads to less water in the GDL/MPL. Hou et al. [29], for the first time, studied the effect of the GDL microstructure on the droplet dynamics in the flow channel of PEM fuel cells and found that the presence of the GDL microstructure affects the droplet motion direction and hinders the droplet movement. In their later work, Deng et al. [30] investigated the hydrophobicity of the GDL and MPL on water distribution. The results show that liquid water tends to flow through macrocracks, leading to less water in the electrode, which agrees with the results of the work by Cetinbas et al. [19].
However, these works only captured the two-phase flow with its effect on the reactive transport and performance neglected. Recently, Zhang et al. [31] employed the LBM to simulate the coupled two-phase and mass transport phenomena in GDL/MPL assemblies. The results show that the oxygen transport benefits from the presence of the MPL for less water in the GDL. However, the MPL is simplified as two dimensional, which cannot give a realistic description of the structural effect. In this work, a numerical stochastic algorithm is adopted to reconstruct the realistic three-dimensional (3D) microstructure of GDL/MPL assemblies, and the structural characteristics can be flexibly adjusted. To achieve a sufficiently low spurious velocity at the gas-liquid interface for model stability and accuracy, a multiple-relaxation-time (MRT) LBM is employed to simulate the coupled transport of liquid water and reactant gas in GDL/MPL assemblies. The electrochemical reaction in the CL is also addressed by applying the source term at the MPL/CL interface. The effects of assembly porosity, macrocracks in the MPL, and perforation design are investigated in terms of dynamic water transport, reactant distribution, and output performance. Further consideration of reactive transport provides deep insights into interlinks among flooding, transport resistance, and performance.
The remainder of the paper is structured as follows: the stochastic reconstruction model of GDL/MPL assemblies is introduced in Sect. 2. The developed numerical tool, LBM, and conducted validation tests are described in Sect. 3. The simulation results are presented and discussed in Sect. 4. Finally, conclusions are summarized in Sect. 5.

Stochastic Reconstruction Model
Based on the literature, reconstruction methods can be generally categorized into experimental approaches, such as X-ray scanning [32] and focused ion beam/scanning electron microscopy [33], and numerical approaches [30]. Considering the advantages of cost saving and flexibility, the stochastic reconstruction method is adopted in this work.
The GDL and MPL are constructed by micron-scale carbon fiber and nanoscale carbon black, respectively. Polytetrafluoroethylene (PTFE) is usually added for bonding and surface hydrophobicity adjustment. The pore size in GDL is in micrometer scale, whereas the MPL has nanometer scale pores built by carbon particles and micrometer scale pores Fig. 1 Reconstructed GDL, crack-free MPL, and cracked MPL built by cracks. In this work, the carbon fiber is simplified as a long-straight cylinder, the carbon black is simplified as a sphere, and the generation of PTFE is based on the local pore size to serve as a binder. The specific reconstruction algorithms can be found in our previous works [30,34]. In this work, all carbon fibers in the GDL have the same diameter and are set as 10 μm. The MPL is 40-μm thick with the width of vertical cracks ranging from 8 to 10 μm. The reconstructed GDL, crack-free MPL, and cracked MPL are shown in Fig. 1.

Two-Phase Model
Considering the outstanding advantages of the implemented boundary condition and parallel computing, the LBM is employed as a numerical tool. The single relaxation time (SRT) collision operator is mostly used in previous works [28,31,35]. However, in the MPL, where the local pore is much smaller than that in the GDL, the SRT operator could lead to a high spurious velocity at the two-phase interface and thus model instability. Therefore, the MRT LBM is employed to ensure the model stability, and the evolution equation is expressed as where f , is the density distribution function of component ; f eq is its equilibrium form; represents the discrete direction in the model; x is the lattice position; and t is the time. In the D3Q19 model, the discrete velocity e is given as Δt is the time step and Δx is the lattice spacing. M is a transformation matrix between the velocity and momentum space, and M −1 is the inverse form of M . In this work, M is expressed as I is the identity matrix, and S is the diagonal relaxation matrix. For the D3Q19 model, S is descried as The relaxation matrix can be adjusted to decrease the spurious velocity, and in this work, −1 The bulk viscosity and kinematic viscosity v are given by and F represents the force term in the model, which is given as where u eq is the barycentric velocity of the fluid, which indicates the velocities at the equilibrium state; F SC(s) is the component force, which is the combination of the fluid-fluid interaction force F ,int (x) and fluid-solid interaction force F ,ads (x) . F ,int (x) and F ,ads (x) are calculated by where G , is the interaction strength parameter between water and air and G ,w is the interaction strength between the water and solid phases. The hydrophobicity of the solid wall can be controlled by adjusting the G ,w value. is a pseudopotential, which equals the local density. w( | | e | | 2 ) is the weighted parameter, shown as

Mass Transport Model
Regarding reactant transport, the D3Q7 model is selected for its simplicity [36]. The evolution equation is described as where c (x, t) represents the concentration distribution function and c eq (x, t) is its equilibrium distribution. The transformation matrix M c and corresponding relaxation matrix S c are expressed as The corresponding relaxation coefficients are shown as where D is the macroscopic diffusion coefficient, and the local concentration C can be obtained by In this study, the CL is simplified as an infinitely thin layer attached to the MPL [31,37], and the electrochemical reaction occurring at the thin interface is implemented by the boundary condition as follows: where j is the current density; r ′′ is the oxygen consumption rate; and F is the Faraday constant.
According to the Butler-Volmer equation [38], the reaction constant k sr of the activated surface in the CL can be described as [39] where C O 2 is the local oxygen concentration; a s represents the roughness factor of the CL; j ref is the reference current density; C ref is the reference oxygen concentration; f and r are the transfer coefficients for the forward reaction and reverse reaction, respectively; and is the overpotential.
The modified bounce-back conditions proposed by Kamali et al. [39] are adopted to the surface reaction. Then, the real physical electrochemical reaction rate k sr is converted into the nondimensional lattice unit k LB sr by Then, the boundary condition can thus be implemented.

Water-Oxygen Coupling Scheme
Because the solubility and diffusivity of oxygen in liquid water are orders of magnitude lower than those in the air, this work ignores the presence of oxygen in the water. When liquid water moves into the oxygen node, the local oxygen concentration will be cleared to zero, and the original local distribution will stream to the neighbor node to ensure mass conservation. The averaged extrapolation-refilling scheme proposed by Peng et al. [40] is applied in this study, which is given by where N is the number of discrete directions.

Validation of Laplace's Law
The Young-Laplace equation shows the relationship among the pressure, surface tension, and droplet raids. The internal and external pressure differences of droplets with different radii are compared to verify the model. For the 3D case, Laplace's law is written as where P in − P out represents the pressure difference on both sides of the droplet; is the surface tension; and R is the radius. The domain size is 150 × 150 × 150 lu (lattice unit), and the periodic boundary is applied in each direction. Droplets with different radii are initialized in the center of the domain. When the droplet is stable, the pressure difference across the liquid-gas interface is measured. As shown in Fig. 2a, the simulation results agree well with the analytical results, which validate the two-phase model.

Validation of the Mass Transport Model
To verify the mass transport model and reactive boundary condition, the one-dimensional diffusion benchmark is conducted, schematically shown in Fig. 2b. The distance between the two infinite reaction planes is 2l , and the initial oxygen concentration in the computation domain is set as C 0 . Due to the symmetry structure, the gradient of oxygen concentration in the middle point should be 0, and the governing equations can be presented as The analytical results can be calculated as where C * (x, t) represents the nondimensional oxygen concentration at time t and position x and Fo is the Fourier number obtained from the equation Fo = Dt∕l 2 ; n is the Eigen parameter and can be determined by n tan( n ) = Da , Da is Damköhler number and can be calculated by Da = lk sr ∕D . The comparisons between analytical and simulation results under different Fo values are shown in Fig. 2c, and a good agreement is achieved, which validates the mass transport model.

Computational Domain and Boundary Conditions
The computation domain is shown in Fig. 3a, including the buffer space, GDL, MPL, and CL. The thicknesses of the buffer space A, GDL, MPL, and buffer space B are set to 160, 180, 40, and 20 μm, respectively. As mentioned above, the CL is simplified as a thin interface attached to the bottom of the MPL for implementing a reactive boundary condition (labeled as "reaction surface" in Fig. 3a). The buffer spaces A and B are added to the inlet of the gas and liquid water, respectively, to ensure model stability [41]. The computational domain is discretized into 200 × 200 × 200 grids, and the grid resolution is 2 μm. The unit conversion method can be found in Ref. [30].
To reduce the scale effect, periodic boundary conditions are applied in the y and z directions, and the internal solid structure (carbon lattices) is defined as the half-bounce-back boundary [26]. For the oxygen diffusion process, constant concentration is set at the top surface. The reactive boundary scheme proposed by Kamali et al. [39] is employed to describe the electrochemical reaction on the MPL/CL interface. Regarding the two-phase flow, liquid water flows into the domain from the inlet holes in buffer space B with a constant velocity of 0.0024 m/s, as shown in Fig. 3b. The boundary implementations are based on the Zou-He method [42].

Results and Discussion
In this work, 12 cases are utilized to examine the effects of GDL/MPL assemblies on the water and oxygen transport in terms of the structural porosity, MPL cracks, and perforation design, and the specific quantitative summary is shown in Table 1.

Effect of GDL Porosity
In this study, GDLs with porosities from 0.55 to 0.75 are reconstructed, and the inside dynamic transport processes of liquid water and oxygen are illustrated in Fig. 4. In this work, the solid node is hydrophobic with a contact angle of 120°.
First, oxygen diffuses into the GDL, which is consumed in the CL, and the oxygen distribution in the computation domain gets stable at 90000 lu time. Then, liquid water flows into the computation domain from the CL side. As presented in Fig. 4a, liquid water breaks through the GDL and occupies the original pore space. For porosities ranging from 0.55 to 0.75, the breaking through of liquid water to the GDL with a low porosity takes less time, which can be explained by the relatively small average pore size. Moreover, the water distribution inside the GDL will become stable after the breaking through, as the liquid water prefers to transport through large pores, which agrees well with the results in Ref. [30]. The flooding effect on the oxygen distribution is shown in Fig. 4b. Before the flooding (90000 lu time), the GDL with a large porosity has a high average oxygen concentration in the CL layer for a low transport resistance. However, after the flooding, much pore space is occupied by liquid water, which increases the transport resistance, and the average oxygen concentration in the CL layer is significantly decreased, as shown in Fig. 4c. The results show that the concentration drop is the highest in the 0.75 porosity case, indicating severe flooding. The gradient porosity design yields the highest average oxygen concentration and the lowest gap after the flooding, which proves its advantage on water management.  Figure 5a shows the evolution of the current density distribution under flooding, and Fig. 5b quantitatively presents the change in the activated reaction area ratio in CL. The findings show that liquid water can easily form a water film and block the CL surface under a large porosity near the CL side, resulting in a high reactive area loss, which agrees with the results of Ko et al. [13]. Figure 5c shows the variation of the average current density against the flooding and proves that the flooding leads to serious performance degradation. After the flooding, the current density of the GDL with 0.65 porosity is slightly higher than the GDLs with 0.55 and 0.75 porosities. Moreover, the GDL with gradient porosity yields the best performance against flooding, and thus, it is suggested to improve water management. On the one hand, as shown in Fig. 4c, the gradient porosity GDL exhibits good liquid water distribution, resulting in a high average oxygen concentration in the CL after flooding. On the other hand, as shown in Fig. 5b, the gradient porosity GDL has a low porosity near the CL side, and less liquid water accumulates at the GDL/CL interface, resulting in a large activated area. Therefore, it yields the best performance.

Effect of MPL Porosity
In this section, the GDL porosity is set as 0.65, and three crack-free MPLs with porosities of 0.5, 0.6, and 0.7 are reconstructed and attached to the GDL.
As shown in Fig. 6a, liquid water can hardly transport through the microscopic pores in the MPL, and the presence of the MPL effectively ensures the oxygen transport space under the flooding, leading to a high average oxygen concentration in the CL layer and low reactive area loss. Such results agree well with the conclusions by Lee et al. [14]. Figure 6b shows that the concentration drop is low with a high porosity MPL due to the presence of many water transport paths. For the same reason, liquid water can easily be stuck in the MPL/CL interface under the 0.5 porosity MPL, hence resulting in a high activation area loss, as shown in Fig. 6c, d. Moreover, the MPL makes the oxygen concentration distribution and current density uniform due to its strong in-plane transport capacity. Figure 6e shows that the MPL helps to ensure the performance under flooding, and the crack-free MPL with high porosity is preferred.

Effect of the GDL/MPL Structure
MPLs with cracks and perforations are reconstructed to investigate the effect of the GDL/MPL structure. Micronscale cracks are usually formed inside the MPL microstructure during its production [43], and studies have proven that laser perforation can optimize water management in PEM fuel cells [21]. The crack morphology is based on Ref. [44]. In Cases 9 and 10, perforations are only constructed in the GDL, whereas in Case 11, they penetrate through the GDL and MPL. In Case 10, the perforation is located at the water breakthrough point. By contrast, in Case 9, it is randomly located. The diameter of the perforation is 20 μm.
As shown in Fig. 7a, liquid water mainly flows through the GDL/MPL assembly by cracks and perforations, and the perforation design can lead the liquid water transport through the assembly effectively. Moreover, the perforation design significantly decreases the concentration drop caused by the flooding, as shown in Fig. 7b, as less liquid water blocks the pore space in the assembly and the transport resistance is small. Such optimization is evident when the perforation is located at the liquid water breakthrough point (Case 10). Among the various designs, the systematic perforation design yields the lowest concentration drop by separating the transport of liquid water and reactant gas.
As shown in Fig. 8a, b, the crack and perforation design can increase the activated area for a low water transport resistance. The reaction area is the highest under the systematical perforation design (Case 11) for an organized water transport path. Moreover, the systematical perforation design ensures a complete uniform distribution of the reaction rate, as shown in Fig. 8a, and yields the best performance for simultaneously increasing the activated reaction area and oxygen concentration near CL, as shown in Fig. 8b, c. As mentioned above, the gradient porosity GDL and high porosity MPL, with the presence of cracks and systematic perforations, are favorable for the performance under flooding. Case 12 simulates all the above optimal conditions. The MPL with large porosity and the GDL with gradient porosity will retain more pores in the oxygen transmission path after the liquid water breakthrough, resulting in a great average oxygen concentration in the CL, as shown in Fig. 7b. However, the MPL with a large porosity also leads to a small activation

Conclusions
In this work, a 3D multicomponent multiphase LBM is established to examine the coupled two-phase and reactive transport in the stochastically reconstructed GDL/MPL assembly. The MRT collision operator is used to ensure the model stability and accuracy. The flooding effect on the mass transport and performance is studied with the inside water, oxygen, and reaction rate distributions resolved. The simulation results show that water quickly breaks through the GDL in a low porosity structure and the water distribution becomes stable after breaking through the GDL, as water prefers to transport through large pores. The gradient porosity GDL can increase the average reactant concentration under the flooding condition. The presence of the MPL effectively ensures the oxygen transport space under the flooding and achieves a low reactive area loss. Moreover, the MPL promotes the uniform distribution of oxygen concentration and current density for a good in-plane transport capacity. The crack and perforation designs can accelerate the water transport in the assembly. The systematic perforation design ensures the separated transport of water and gas and yields the best performance under flooding.