Pore-Scale Modeling of Air–Water Two Phase Flow and Oxygen Transport in Gas Diffusion Layer of Proton Exchange Membrane Fuel Cell

Abstract

Understanding multiphase flow and gas transport occurring in electrodes is crucial for improving the performance of proton exchange membrane fuel cells. In the present study, a pore-scale model using the lattice Boltzmann method (LBM) was proposed to study the coupled processes of air–water two-phase flow and oxygen reactive transport processes in porous structures of the gas diffusion layer (GDL) and in fractures of the microscopic porous layer (MPL). Three-dimensional pore-scale numerical results show that the liquid water generation rate is gradually reduced as the oxygen consumption reaction proceeds, and the liquid water saturation in the GDL increases, thus the constant velocity inlet or pressure inlet condition cannot be maintained while the results showed that at t = 1,200,000 iterations after 2900 h running time, the local saturation at the GDL/MPL was about 0.7, and the maximum value was about 0.83, while the total saturation was 0.35. The current density reduced from 2.39 to 0.46 A cm⁻². Effects of fracture number were also investigated, and the results showed that for the fracture numbers of 8, 12, 16, and 24, the breakthrough point number was 4, 3, 3, and 2, respectively. As the fracture number increased, the number of the water breakthrough points at the GDL/GC interface decreased, the liquid water saturation inside the GDL increased, the GDL/MPL interface was more seriously covered, and the current density decreased. The pore-scale model for the coupled multiphase reactive transport processes is helpful for understanding the mechanisms inside the porous electrodes of PEMFC.

1. Introduction

During the past few decades, the proton exchange membrane fuel cell (PEMFC) has received much attention due to several advantages including high energy efficiency, room-temperature operation and low pollution. Improving fuel cell performance highly depends on effective thermal and water management of PEMFC, which are based on aa comprehensive understanding of water transport mechanisms inside the PEMFC [1]. The water transport process in GDL plays an important role in the mass transport and fuel cell performance [2]. Too much liquid water will cause flooding and hinder the oxygen transport. Therefore, it is desirable to accelerate the liquid water migration from the GDL and to reduce the liquid water saturation in GDL. In practice, the GDL is usually treated with hydrophobic polytetrafluoroethylene (PTFE) to facilitate the liquid water transport [3]. An microporous layer (MPL) is added to the GDL to help improve the transport of liquid water through the GDL [4].

Advanced visualization techniques have been adopted to conduct in situ measurement of multiphase flow inside GDL such as x-ray computed tomography (XCT), nuclear magnetic resonance, and neutron imaging. With great improvement in spatial and temporal resolution, XCT has been adopted to reconstruct the porous structures of GDL and to investigate liquid water dynamic behaviors [5,6,7,8]. Supplementary to the experiments, numerical simulation also plays important roles in investigating the effects of structural and operating parameters on liquid water transport in GDL. Models for the multiphase flow in GDL can be divided into continuum-scale models and pore-scale models [9]. At the continuum-scale, models based on the concept of representative elementary volume (REV) have been proposed for transport processes in porous media, in which the extension of the Darcy law involving relative permeability has been adopted for multiphase flow. A typical example of such continuum-scale model is the multiphase mixture model proposed by Wang et al. [10,11]. Such a model is computationally more efficient and can be used for modeling the entire cells. In this model, the capillary pressure-saturation relationship such as the Leverett–Udell function is adopted for closing the governing equations [12,13]. It has been recognized the Leverett–Udell function obtained from porous media in geology or petrology is not suitable for describing multiphase flow in GDLs. A more accurate pressure–saturation relationship has also been proposed in the literature [14,15,16,17].

The mesh size in the continuum-scale cell model is too large to consider the microscopic details of the GDL structures. To gain a deep understanding of the interactions between porous structures, multiphase flow processes, and GDL performance, pore-scale modeling has been developed and conducted. Among the many numerical methods for pore-scale modeling of multiphase flow, the lattice Boltzmann method (LBM) is preferred due to its power capacity of handling the complicated structures and of capturing the dynamic evolution of the phase interface. Using the LBM, the realistic porous structures of the GDL can be easily taken into account. This method has been most widely adopted for simulating multiphase flow in GDL [18,19]. Hao and Cheng [20] found that by adding a hydrophilic columnar through the GDL, liquid water will migrate through the GDL by the hydrophilic columnar, leaving the remaining hydrophobic part available for the reactant gas. Yu et al. [21] performed three-dimensional (3D) pore-scale simulations to study effects of distributions of PTFE on the liquid water transport inside and breakthrough the GDL. With PTFE distributed heterogeneously along the through-plane direction, the liquid water saturation curves show different shapes. The pore-scale results suggest that placing PTFE near the GDL/MPL interface facilitates the liquid water transport in GDL [21]. Effects of the PTFE distribution on the removal of one single droplet were also studied by Kakaee et al. [22]. Jinuntuya et al. [23] numerically studied at the pore-scale the liquid water transport in three types of GDL structures obtained from XCT and the effects of the GDL structure, the wettability, and the pressure difference applied on the liquid water flow pattern, breakthrough time and saturation inside the GDLs were investigated. Pore-scale two-phase flow in compressed GDL was also studied [24]. It was found that as compression ratio increases, the location of water breaking through the GDL moves from the channel to the channel/rib interface. Under a relatively high compression ratio, liquid water tends to form aa water film inside GDL, which greatly hinders the reactant mass transport [24]. Deng et al. [25] adopted the LBM to simulate liquid water transport inside and at the interface of GDL and MPL. It was found that with the addition of MPL, liquid water flooding in GDL is greatly reduced.

In PEMFCs, oxygen transport through the GDL and MPL arrives at the CL and participates in the electrochemical reaction that generates water. The water is then transported in the opposite direction from the CL to the GDL, and finally removed out of the fuel cell through the GC. However, only Zhang et al. conducted a pore-scale study of such coupled processes between multiphase flow and reactive transport [26]. As the liquid water saturation increases in the GDL, the oxygen transport is gradually hindered, leading to a lower generation rate of liquid water. Therefore, uniform liquid inlet velocity or constant pressure applied at the bottom surface of GDL, which are widely adopted in the current pore-scale simulations, is not sustainable.

From the above review, it can be found that pore-scale modeling considering the coupling processes of air–water two-phase flow and oxygen reactive transport is really rare in the literature. Understanding such coupling mechanisms is of great importance for enhancing mass transport, improving water management, and increasing the cell performance. Therefore in this study, a 3D multiphase reactive transport model was developed to study the coupled multiphase flow and oxygen transport in the GDL and MPL fractures. Furthermore, a 1D model is proposed at the GDL bottom surface to consider the transport resistance in MPL and CL. Effects of the number of the fractures on the multiphase flow are investigated in detail.

2. The Computational Domain and the Physicochemical Processes

In the present study, first, a three-dimensional domain composed of GDL and a fracture in the MPL was constructed to simulate oxygen transport as well as liquid water dynamic behaviors. The computational domain was rectangular with dimensions of 100 µm× 100 µm× 100 μm. Here, the resolution of one lattice was 1 μm. As can be seen from Figure 1a, the GDL showed layered structure characteristics. Therefore, the GDL was reconstructed by inserting the straight carbon fibers layer by layer in the through-plane direction, with the constraining parameters as the GDL porosity and the diameter of the carbon fiber. Details of the reconstruction method can be found in [27]. From the SEM image of MPL, as shown in Figure 1b, there were many fractures inside the MPL. There, fractures can provide preferred pathways for liquid water from the CL and can reduce the water flooding [5]. Therefore, at the bottom of the reconstructed GDL, an elongated rectangle was added to represent one fracture in the MPL. The final computational domain is shown in Figure 2.

The physicochemical processes taking place inside the GDL can be described as follows. Oxygen from the top boundary diffuses into the GDL and arrives at the GDL bottom surface for the electrochemical reaction. In practice, the electrochemical reaction of oxygen actually takes place inside the CL. After diffusing through the GDL, the oxygen transports through the MPL, and arrives at the CL. Inside the CL, the oxygen first diffuses in the macropores between the carbon particles within the CL, and then into the local structures around a carbon particle, before finally being consumed at the triple-phase boundary (TPB). During the above processes, the oxygen flux is the same and can be expressed as follows:

$$D_{\mathrm{GDL}}\frac{\partial C_{\mathrm{oxygen}}}{\partial n}=D_{\mathrm{MPL}}\frac{\partial C_{\mathrm{oxygen}}}{\partial n}=D_{\mathrm{CL},\mathrm{macropore}}\frac{\partial C_{\mathrm{oxygen}}}{\partial n}=D_{\mathrm{CL},\mathrm{local}}\frac{\partial C_{\mathrm{oxygen}}}{\partial n}=k_{\mathrm{elec}}{C}_{\mathrm{oxygen}}$$

(1)

In Equation (1), D is the diffusivity and C is the oxygen concentration. The five terms refer to the flux at the GDL/MPL interface, through the MPL, the macropores in CL, the local structures in CL, and that at the TPB due to oxygen consumption, respectively. Equation (1) can be rearranged as follows:

$$J_{\mathrm{oxygen}}=2J_{\mathrm{water}}=D_{\mathrm{GDL}}\frac{\partial C_{\mathrm{oxygen}}}{\partial n}=\frac{1}{(\frac{1}{k_{\mathrm{MPL}}}+\frac{1}{k_{\mathrm{CL},\mathrm{macropore}}}+\frac{1}{k_{\mathrm{CL},\mathrm{local}}}+\frac{1}{k_{\mathrm{elec}}})}{C}_{\mathrm{oxygen}}$$

(2)

where the four terms at the denominator represent the resistance through MPL and CL. J refers to the flux.

The typical pore size GDL is about 10 μm, that of MPL is about 50 nm, while that inside the CL is from 1~10 nm. Simultaneously resolving the pores in the GDL as well as that in the MPL and CL requires an enormous number of grids, which cannot be affordable by current computational resources. Therefore, in the present study, only GDL and fractures in MPL were taken into account, while the solid matrix with nanoscale pores of MPL as well as CL structures were not resolved. Thus during the simulation, Equation (1) was adopted as the boundary condition at the GDL bottom surface. This boundary condition takes into account the transport resistance during the pathway from the GDL/MPL interface to the TPB. In short, the transport processes numerically resolved in the 3D structures of GDL were combined with the 1D boundary condition described by Equation (1), and such a model was called the 3D + 1D model in the present study.

The electrochemical reaction generates water, which is assumed to be in liquid phase. The liquid water generated in the CL migrates from the fracture in MPL into the GDL [5]. Therefore, in the present study during each simulation step, the liquid water generated by Equation (2) is summarized and then averaged into every node in the fracture. Liquid water then grows from the fracture, and enters the GDL. For the air–liquid two phase flow in GDL, capillary number Ca is about 10⁻⁵~10⁻⁸ [19,28]. With such low Ca and the dynamic viscosity ratio between water and air at about 18, the two-phase flow in the GDL is dominated by the capillary force and effects of other forces such as gravitational force, viscous force, and inertial force can be neglected. Correspondingly, the flow is located in the capillary fingering region [29]. Therefore, neglecting the density ratio between air and liquid and setting the density ratio as unity is acceptable for simulating two-phase flow in GDL [30].

As shown in Figure 2., periodic boundary conditions were employed for the x–y plane and x–z plane. For the top and bottom boundaries at the y–z plane, no slip boundary condition was imposed. The electrochemical reaction takes place at the surface of GDL/MPL, as described above. As liquid water migrates, it will cover the GDL/MPL interface, and under such a circumstance, the local interface is not available to the electrochemical reaction. Besides, the pores occupied by the liquid water are also not allowed for oxygen diffusion. For the inner solid fiber surface, no-slip and non-flux boundary conditions were adopted for fluid flow and mass transport, respectively. The contact angle was set as 140° for the multiphase flow simulation.

3. The Lattice Boltzmann Method

The multiphase flow and oxygen reactive transport processes in the porous structures described in Section 2 was solved using the lattice Boltzmann method (LBM). Details of the LB model are introduced as follows.

3.1. Multiphase Flow Model

In the present study, the pseudopotential multiphase LB model was employed to study air–water two-phase flow. A source term is added to consider the water generation due to electrochemical reaction. For the kth component, the evolution equation for the density distribution function with multiple-relaxation-time (MRT) collision term is defined as [31]:

$$\begin{array}{l}{f}_i^k(\mathbf{x}+c{\mathbf{e}}_i\Delta t,t+\Delta t)-f_i^k(\mathbf{x},t)=\\ -\left({\mathbf{T}}^{-1}\cdot {\overline{\mathbf{D}}}^k\cdot \mathbf{T}\right)\left[f_i^k(\mathbf{x},t)-f_i^{\mathrm{eq},k}(\mathbf{x},t)\right]+\left({\mathbf{T}}^{-1}\cdot (\mathbf{I}-\frac{{\overline{\mathbf{D}}}^k}{2})\cdot \mathbf{T}\right)(F_i^k+S_i^k)\end{array}$$

(3)

where $f_i^k$ is the kth density distribution function at the lattice site x and time t.c = Δx/Δt is the lattice speed with Δx and Δt as the lattice length and time step, respectively. k equals 1 for gas and equals 2 for liquid water. The equilibrium distribution functions f^eq is as follows

$$f_i^{\mathrm{eq},k}={\omega }_i{\rho }^k\left[1+\frac{3}{c^2}\left({\mathbf{e}}_i\cdot {\mathbf{u}}^{\mathrm{eq}}\right)+\frac{9}{2c^4}{\left({\mathbf{e}}_i\cdot {\mathbf{u}}^{\mathrm{eq}}\right)}^2-\frac{3}{2c^2}{({\mathbf{u}}^{\mathrm{eq}})}^2\right]$$

(4)

For the D3Q18 lattice model used in this study, the values of the weight coefficient w_i were w_i = 1/3, i = 0; w_i = 1/18, i = 1, 2, …, 6; w_i = 1/36, i = 7, 8, …, 18. The transformation matrix T in Equation (3) is a (N + 1) × (N + 1) matrix, with T⁻¹ as the inverse matrix of T [31]. I is the unit matrix. The diagonal relaxation matrix ${\overline{\mathbf{D}}}^k$ is defined as [32].

$${\overline{\mathbf{D}}}^k={\mathrm{diag}(\mathrm{s}}_0{,\mathrm{s}}_1,\dots {,\mathrm{s}}_{17}{,\mathrm{s}}_{18})$$

(5a)

$${\mathrm{s}}_{0~8,10,12,16~18}=1, {\mathrm{s}}_9={\mathrm{s}}_{11}={\mathrm{s}}_{13-15}=\frac{1}{{\tau }^k}$$

(5b)

More details of T and ${\overline{\mathbf{D}}}^k$ can be found in Ref. [31]. The density, velocity, and viscosity are calculated by [30]

$${\rho }^k={\displaystyle \sum f_i^k}, {\rho }^k{\mathbf{u}}^k={\displaystyle \sum f_i^k}{\mathbf{e}}_i+\frac{\delta t}{2}{\mathbf{F}}^k,{\upsilon }^k=\frac{1}{3}({\tau }^k-0.5)\frac{\Delta x^2}{\Delta t}$$

(6)

The effective velocity ${\mathbf{u}}^{\mathrm{eq}}$ in Equation (4) is as follows

$${\mathbf{u}}^{\mathrm{eq}}={\displaystyle \sum s_0^k{\rho }^k}{\mathbf{u}}^k/{\displaystyle \sum s_0^k{\rho }^k}$$

(7)

In Equation (3), $F_i^k$ is the force term and $S_i^k$ is the source term in the LB framework. $F_i^k$ is calculated by

$$F_i^k=\frac{{\mathbf{F}}^k\cdot ({\mathbf{e}}_i-{\mathbf{u}}^{\mathrm{eq}})}{{\rho }^k{c}_s{}^2}{f}_i^{eq,k}\Delta t$$

(8)

where $F^k$ (also in Equation (6)) includes fluid–solid interaction force, fluid–fluid interaction force, and external body force

$${\mathbf{F}}^k={\mathbf{F}}_f^k+{\mathbf{F}}_s^k+{\mathbf{F}}_e^k$$

(9)

For the particles of the kth component at lattice site x, the total fluid–fluid surface tension force is expressed as

$${\mathbf{F}}_f^k=-{\psi }^k\left({\rho }^k(\mathbf{x})\right)g_{kk\prime }{\displaystyle \sum _{\mathbf{x}\prime } }{\displaystyle \sum _{k\prime }^{s}w(\mathbf{x}\prime )}{\psi }^k\left({\rho }^{k\prime }(\mathbf{x}\prime )\right)(\mathbf{x}\prime -\mathbf{x})$$

(10)

where ψ is the effective mass or interparticle potential and is defined as ψ_k(ρ_k) = 1-exp(-ρ_k). $g_{kk\prime }$ is the interaction strength and iteffectively controls the immiscibility of the binary mixture and the surface tension force. The value of $g_{kk\prime }$ was taken as 1.95 in this study. If the interaction force of four nearest neighbors is considered, the weight factor $w(\mathbf{x}\prime )$ is 1/3 and 1/12 for $|\mathbf{x}\prime -\mathbf{x}{|}^2=1$ and $|\mathbf{x}\prime -\mathbf{x}{|}^2=2$, respectively. ${\mathbf{F}}_s^k$ is the fluid–solid interaction force

$${\mathbf{F}}_s^k=-{\psi }^k\left({\rho }^k(\mathbf{x})\right)g_s{\displaystyle \sum _{\mathbf{x}\prime }w(\mathbf{x}\prime )s(\mathbf{x}\prime )(\mathbf{x}\prime -\mathbf{x})}$$

(11)

where s is an indicator function, with 0 and 1 for pore and solid, respectively. g_s controls the fluid-solid strength, By adjusting $g_s$, different wettability (contact angle) can be obtained.

The LB source term in Equation (3) is calculated by [33]

$$S_i^k={\omega }_i{S}^k\left[1+\frac{3}{c^2}\left({\mathbf{e}}_i\cdot {\mathbf{u}}^{\mathrm{e}\mathrm{q}}\right)+\frac{9}{2c^4}{\left({\mathbf{e}}_i\cdot {\mathbf{u}}^{\mathrm{e}\mathrm{q}}\right)}^2-\frac{3}{2c^2}{({\mathbf{u}}^{\mathrm{e}\mathrm{q}})}^2\right]\Delta t$$

(12)

where $S^k$ is the actual source term considering theelectrochemical reaction at the Pt–ionomer interface. This source term will be discussed in Section 3.3. Incorporating $S^k$ into the evolution equation (Equation (3)) through the form of Equation (12), the Galilean invariance can be guaranteed according to previous study [33].

3.2. Mass Transport Model

Oxygen transports from the GC into the GDL and arrives at the bottom surface of GDL for the electrochemical reaction. The mass transport is solved using the LB mass transport model as follows

$$g_i(\mathbf{x}+c{\mathbf{e}}_i\Delta t,t+\Delta t)-g_i(\mathrm{x},t)=Q^{-1}\mathsf{\Lambda }Q(g_i(\mathbf{x},t)-g_i^{\mathrm{eq}}(\mathbf{x},t))$$

(13)

where g_i is the concentration distribution function at the lattice site x and time t in the ith direction. $g_i^{\mathrm{eq}}$ is the equilibrium distribution function. The D3Q7 lattice model is adopted in this study, and $g_i^{\mathrm{eq}}$ was calculated by $g_i^{\mathrm{eq}}=w_{i}C$, with $w_0=1/4\mathrm{and}{w}_{1-6}=1/8$. The transformation matrix Q in Equation (13) transfers the distribution function in velocity space into momentum space. Q⁻¹ is the inverse matrix of Q, and for moredetails of Q and Q⁻¹, one can found in Ref. [34]. Λ is the relaxation coefficient matrix.

$$\mathsf{\Lambda }=\left[\begin{array}{ccccccc}{\tau }_0,& 0,& 0,& 0,& 0,& 0,& 0\\ 0,& {\tau }_{xx},& {\tau }_{xy},& {\tau }_{xz},& 0,& 0,& 0\\ 0,& {\tau }_{yx},& {\tau }_{yy},& {\tau }_{yz},& 0,& 0,& 0\\ 0,& {\tau }_{zx},& {\tau }_{zy},& {\tau }_{zz},& 0,& 0,& 0\\ 0,& 0,& 0,& 0,& {\tau }_4,& 0,& 0\\ 0,& 0,& 0,& 0,& 0,& {\tau }_5,& 0\\ 0,& 0,& 0,& 0,& 0,& 0,& {\tau }_6\end{array}\right]$$

(14)

where the relaxation coefficients ${\tau }_{\alpha \beta }$ is related to the diffusivity

$$D_{\alpha \beta }=\zeta ({\tau }_{\alpha \beta }-\frac{1}{2}{\delta }_{\alpha \beta })\frac{\mathsf{\Delta }{x}^2}{\mathsf{\Delta }t}$$

(15)

where ${\delta }_{\alpha \beta }$ is the kronecker symbol. For isotropic mass transport, ${\tau }_{xx}={\tau }_{yy}={\tau }_{zz}$ and ${\tau }_{\alpha \beta }(\alpha \ne \beta )=0$. In the present study, the values of ${\tau }_0$, ${\tau }_4$, ${\tau }_5$, and ${\tau }_6$ in Equation (16) were set as unity. $\zeta$ in Equation (16) was 1/4.

3.3. Source/Sink Term for Multiphase Flow

Now, our attention turned to determining the source/sink term in Equation (12) for multiphase flow. As discussed above, liquid water enters the GDL from the factures inside the MPL. During each time step, the total amount of liquid water generated can be determined by the summation of all the flux described by Equation (2). Then, all water generated is relocated into the fractures.

$$S^1=-\frac{1}{V_{\mathrm{facture}}}{\displaystyle \sum M_{\mathrm{water}}{J}_{\mathrm{water}}A}$$

(16)

$$S^2=\frac{1}{V_{\mathrm{facture}}}{\displaystyle \sum M_{\mathrm{water}}{J}_{\mathrm{water}}A}$$

(17)

where M_water is the molar mass; A is the surface area of each computational node; and V_fracture is the total volume of the fracture.

The general modeling procedure of each time step during the pore-scale modeling are as follows: (1) The pore-scale two-phase flow modeling is conducted, by which the liquid water distribution is determined; (2) The pore-scale oxygen mass transport modeling is conducted, and such modeling is not allowed in the liquid phase and only conducted in the gas phase; and (3) Based on the results in Step 2, the liquid water generation rate is determined, which is then adopted in the two-phase flow modeling in Step 1. Steps 1–3 are repeated during the simulation until the liquid water touches the top wall of the GC or the liquid water generation rate is extremely low due to the coverage of the bottom surface by liquid water or the oxygen transport is seriously hindered by the liquid water. It is worth mentioning that the pore-scale two-phase model and the oxygen reactive transport model have been well validated in our previous work, and thus the validation is not repeated here for brevity [35,36].

4. Results and Discussion

4.1. Multiphase Reactive Transport in the 3D Domain

The 3D simulation of oxygen reactive transport and air–liquid water two-phase flow was performed by a self-developed LBM code parallelized based on domain decomposition scheme using message passing interface (MPI). The domain was cut into 72 sub-domains and 72 CPU cores were adopted.

Figure 3 shows the air–liquid water interface distribution in the GDL at different times. As can be seen in Figure 3, as the oxygen reduction reaction proceeds, liquid water gradually accumulates inside the MPL fracture and enters the GDL as tiny droplets. As time goes on, these tiny droplets merge to big droplets. Affected by the complex porous structures, some droplets advance deep into the GDL, while the growth of other droplets is suppressed. The liquid water front moves in both in-plane and through-plane direction, although irregularly, but following the capillary fingering mechanism, which means that the liquid water front always advances into the pores with the largest pore size. Finally at t = 1,200,000 iterations, there was excessive liquid water in the GDL, and the GDL was severely flooded at the bottom half. Figure 4 further shows the bottom view of the computational domain, which provides information of coverage by the liquid water. It can be found that during the migration of the liquid water inside the GDL, the bottom surface was gradually covered by the liquid water, leading to a reduced reactive surface area. For the liquid water distribution at t = 1,200,000 iterations, it can be found that the transport pathway for the oxygen was seriously blocked and the reactive surface area was greatly reduced.

Note that in the present simulation, the multiphase flow and the oxygen reactive transport were coupled. Therefore, the more the oxygen consumed, the more liquid water generated. As the liquid water saturation increases in the GDL, the pores inside the GDL are gradually blocked, and thus less oxygen is available to the reaction at the bottom surface. Therefore, the generation rate of the liquid water gradually decreases as time proceeds. In the literature, the majority of the pore-scale study of liquid water transport in the GDL adopts either a velocity inlet or constant pressure boundary condition for the liquid water entering GDL. Such a boundary condition leads to a continuous and steady supply of liquid water into the GDL even when the GDL is seriously flooded, which in fact neglects the coupling mechanisms between air–water two-phase flow and the oxygen reactive transport. In fact, in the present study, due to the extremely slow generation rate of the liquid water in the latter stage due to the inhibition of the oxygen reactive transport process, the running time for the results shown in Figure 3 is as long as four months, leading to a total of about 2900 h for each CPU core, which is really time consuming. Even with such a long time simulation, the liquid water did not break through the GDL at t = 1,200,000 with liquid saturation in the entire GDL at about 0.35. In the pore-scale study of the coupled multiphase flow and oxygen transport by Zhang et al., the liquid water mainly transports in the through-plane direction, and breakthrough of the liquid water was observed with a saturation of about 0.23 [26]. The pore-scale results in the present study as well as that of Zhang et al. suggest that liquid water transport and distribution should be controlled effectively and that the liquid water transports in the through-plane direction should be facilitated while that in the in-plane direction should be suppressed. Perforated GDL is one such scheme in which liquid water penetrates the GDL by the perforated pores [37].

The through-plane distributions of the liquid water saturation at different times are shown in Figure 5. The saturation is defined as the ratio of volume occupied by the liquid water to the entirepore volume in a given cross-section. It can be found that liquid water saturation generally decreases from the GDL/MPL interface to the GDL/GC interface. Behind the liquid water front, the liquid water saturation is below unity, indicating that the void space is partially filled by the liquid water, confirming the capillary fingering mechanism. In the literature, for pore-scale modeling adopting the uniform velocity inlet as the boundary condition, usually the local water saturation at the GDL/MPL approaches unity. The results in Figure 5 show that the local saturation when t = 1,200,000 was only about 0.7. The maximum value of local saturation located around x* = 0.25 was about 0.83. However, liquid water saturation in the remaining parts, especially x* > 0.6, was very low. Such high local high saturation leads to a small space for oxygen transport, greatly reducing the effective diffusivity of oxygen through the entire GDL, even though the averaged saturation in the entire GDL was only about 0.35 at t = 1,200,000. Figure 6 displays the time evolution of the current density obtained, which reduced from the initial value of 2.39 to 0.46 A cm⁻² at t = 1,200,000, indicating the negative effects of liquid water on the reactive transport processes.

4.2. Effects of the Fracture Number

The above pore-scale 3D simulation clearly show the coupling mechanisms of multiphase flow and oxygen transport processes inside the GDL. Since a 3D simulation is really time consuming, in this section, 2D simulations were conducted to investigate the effects of facture number on the coupled processes. The computational domain is shown in Figure 7 with a size of 2048 × 456 lattices with the resolution of one lattice as 1 μm. The domain contained a GDL, a MPL with different numbers of fractures, and a free region on the top of the GDL representing the GC. The height of the MPL, the GDL, and the free region was 7, 156, and 290, respectively. The porosity of the GDL was 0.8004. The contact angle of the solid particles in GDL, the GDL/MPL interface as well as the top wall of the GC was the same as 140°. Different number of fractures in the MPL was investigated including 8, 12, 16, and 24 fractures. The width and depth of each fracture was 10 and 4, respectively. The gap between two neighboring fractures was 227, 157, 120, and 81 for the four different number of fractures studied. The dynamic behaviors of the liquid water, the evolution of the saturation, the saturation distribution along the y direction, the coverage area of the GDL/MPL interface, and the current density will be discussed in detail.

Figure 7 shows the dynamic behaviors of liquid water in the domain with the fracture number of 8. As the reaction proceeded, liquid water is gradually generated and eight liquid clusters were observed connected to the eight fractures (t = 40,000). Each cluster advances inside the GDL following the mechanism of capillary fingering, namely always searching for the pore with the largest pore size. Affected by the local complicated porous structures, some clusters move mainly in the through-plane direction, while others were mainly in the in-plane direction. At t = 200,000, clusters A1–A3 merged into one cluster; cluster A4 found its pathway along the in-plane direction; clusters A5, A6, and A8 were trapped by local structures with relative low pore size, and cluster A7 advanced the most along the through-plane direction. At t = 550,000, cluster A1–A4 merged into one cluster, A5–A6 merged into one cluster, cluster A8 still advanced separately, while cluster A7 broke through the GDL and formed a tiny droplet in the GC. At t = 750,000, clusters A1–A4 merged and supported the growth of two droplets. At t = 1,500,000, clusters A5 and A6 merged and supported the growth of one droplet. Clusters A7 and A8 also merged and supported the growth of the first droplet forming in the GC. After the breakthrough, continuous pathways were created for the liquid water inside the GDL, while the dynamic change of the local liquid water distribution inside the GDL still could be observed. Finally, at t = 3,500,000, there were four droplets in the GC connected to different clusters. The growth rate of the four droplets in the GC were different, and it was found that the earliest breakthrough did not necessarily mean the highest growth rate. The growth rate of the droplet in the GC highly depends on the flow rate of the water cluster connected to it. Finally, it can be seen that liquid water distribution was not uniform inside the GDL. Locally, the liquid water amount at the top of GDL was low while the bottom of the GDL was seriously flooded. The GC was also severely blocked by the four droplets.

Figure 8 shows the corresponding oxygen concentration field at two t = 45,000 and t = 350,000. It was found that due to the consumption at the GDL bottom, the oxygen concentration gradually decreased from the top to the bottom of the domain. The liquid water hinders the transport of oxygen, and causes the reduction of the reactive surface area.

Figure 9, Figure 10 and Figure 11 further show the dynamic behaviors of liquid water in the domain with fracture numbers of 12, 16, and 24. The dynamic behaviors of growth, coalescence, and breakthrough discussed in Figure 7 can also be observed and are not repeated here. It can be found that the fractures inside the MPL significantly affects the liquid water dynamic behaviors in the GDL and GC. In particular, the breakthrough location, breakthrough time, and breakthrough point number (or droplet number in the GC) are closely related to the fracture numbers. For the breakthrough time, it is expected that the breakthrough time for the case with the most fractures is the longest. For the breakthrough point number, the results are interesting. Intuitively, more fractures in the MPL will lead to more clusters in the GDL, and thus finally more droplets appearing in the GC. However, it can be found in Figure 7, Figure 9, Figure 10 and Figure 11 that the number of the breakthrough points and the thus the corresponding droplets in the GC decreased as the fracture number in the MPL increased. For the fracture numbers of 8, 12, 16, and 24, the breakthrough point number was 4, 3, 3, and 2, respectively. Note that although there were two droplets in the GC for the case with fracture numbers of 12, the left droplet was actually nourished by two breakthrough points, as can be seen at t = 1,200,000 and 2,100,000. This is because as more fractures were added in the MPL, the gap between the neighboring fractures decreased, and the neighboring clusters were more easily to merge, thus resulting in fewer clusters in the GDL and finally few breakthrough points at the GDL/GC interface.

Figure 12 shows the evolution of the liquid water saturation along the y direction for the four cases with different fracture numbers. It can be found that the liquid water saturation in the GDL increases as time proceeds. After the breakthrough, the saturation distribution inside the GDL still varies. These results are in consistent with the dynamic behaviors of liquid water in Figure 7, Figure 9, Figure 10 and Figure 11. For all the four images in Figure 12, it can be found that the saturation decreased and then increased across the GDL/GC interface, with a local minimum value obtained. The minimum value was usually less than 0.1. However, at the GDL/MPL interface, the local flooding was serious, and the local saturation value was higher than 0.8. Comparing the four images, it can be found that as the fracture number increases, the local liquid water saturation at the GDL/MPL increases.

Figure 13a shows variation of the total liquid water saturation inside the GDL and the coverage of the GDL/MPL interface. Here, the total saturation is defined as the ratio between the total liquid water volume in the GDL and total volume of the void space in GDL. It can be found that the coverage ratio increases as more fractures are considered, in consistent with the local saturation at the GDL/MPL interface as discussed in Figure 12. For the total saturation, in all cases, the total saturation in the GDL gradually increased and reached a constant value. As the fracture number increased, the final value of the total saturation generally decreased. In particular, the saturation of the case with 24 fractures was much higher than the others. As discussed previously in Figure 7, Figure 9, Figure 10 and Figure 11, more fractures led to fewer droplets in the GC. Therefore, it seems that for the range of fracture numbers studied in the present study, as the fracture number increased, the flooding in GDL was more severe while that in the GC was alleviated, with the former one reducing the cell performance while the latter was desirable. Finally, Figure 13b shows the time evolutions of the current density. Corresponding to the above discussion, the current density decreases as the time proceeded due to the increased liquid water saturation. As the fracture number increased, the final value of the current density decreased.

5. Conclusions

Understanding multiphase reactive transport processes inside porous electrodes of PEMPC is of great importance for improving performance of PEMFC. In this study, a pore-scale model adopting the LBM was proposed for the coupled air–water two-phase flow and the oxygen transport processes in GDL and fractures of MPL. In the model, liquid water generated by the electrochemical reactions enters the GDL from the MPL fractures. The coupled mechanisms between multiphase and oxygen reactive transport are discussed in both 3D and 2D simulations. The simulation results show that as the liquid water saturation in the GDL increases, the oxygen transport is hindered, and thus the generation rate of liquid water decreases. The liquid water cannot maintain constant velocity inlet or pressure inlet boundary conditions that are widely adopted in the literature. For the 3D domain studied in the present study, the migration of the liquid water in the in-plane direction caused serious flooding, leading to breakthrough failure of the liquid water from the GDL. Through-plane migration of liquid water is highly desirable. Furthermore, the effects of fracture number on the coupled processes and some important variables such as reactive surface coverage, saturation, and current density were discussed. It was found that liquid water from the different fractures of MPL shows complex growth, coalescence, and breakthrough behaviors in the GDL. Generally, as the fracture number in the GDL increased, the number of liquid water breakthrough points decreased, the total saturation in the GDL as well as the local saturation at the GDL/MPL interface increased, and the current density decreased. The optimization of the GDL structure and wettability based on the pore-scale modeling results is undergoing in our group.