Microstructure Reconstruction and Gas-Liquid Two-Phase Transport Mechanism within Porous Electrodes of PEM Fuel Cells

Abstract

The structure of porous media is composed of skeleton particles and pores. Its micro-pores and solid skeleton characteristics lead to the capillary fingering movement of fluid in its porous media driven by capillary pressure. Currently, the methods of constructing porous media are mainly random construction and multi-scale imaging construction. The porous structure constructed by these two methods can show the real microstructure characteristics. The research on multiphase flow in microporous structure mainly includes VOF, MC, LBM, and other methods. In this chapter, taking the classic porous structure of polymer electrolyte membrane (PEM) fuel cell gas diffusion layer (GDL) as an example, GDL porous microstructure is constructed through random algorithm, and multiphase LBM is used to study two-phase flow in porous media to explore the relationship between porous structure characteristics and multiphase flow transport.

1. Introduction

In the field of energy and environment, the flow, heat, and mass transfer processes occurring in porous media generally exist [1, 2], such as the migration and distribution of light non-aqueous phase liquid (LNAPL) [3] and heavy non-aqueous phase liquid (DNAPL) [4] in porous media soil, the heat and catalytic reaction process of porous media catalyst particles in chemical granular bed chemical reactors [5], gas-liquid electric thermal multiphysics simulation heat, and mass transfer process in fuel cell porous structure assembly [6, 7, 8]. The heat and mass transfer process in the above porous media is closely related to the energy and environment on which human beings rely for survival, and involves all aspects of human life and industrial generation. Therefore, the study of the transport process in porous media is of great significance for rational energy exploitation and utilization, energy conservation and emission reduction, and environmental protection [9, 10]. To gain a deeper understanding of the influence of porous media structure on its internal transport, this chapter will take the porous media structure of proton exchange membrane (PEM) fuel cell GDL as an example to explore the gas-liquid two-phase transport process inside it.

PEM fuel cells are a type of power generation equipment that differs from traditional fossil fuels. It has the advantages of environmental protection, high efficiency, and fast response [11], and can be widely applied in various fields [12]. It is considered one of the most promising energy technologies. Despite the rapid development of PEM fuel cell technology in recent years, its large-scale commercialization still faces great challenges. Improving the water management capability of PEM fuel cell and strengthening the drainage performance and anti-aging performance of each component are important methods to improve the overall performance of PEM fuel cell [13, 14]. As an important component of PEM fuel cell, the GDL has the following functions: removing liquid water from the fuel cell, providing a pathway for the transport of reactant gas to the catalyst layer, providing structural support for the membrane, and promoting the transfer of electrons and heat [15, 16]. GDL is usually composed of a certain amount of hydrophilic carbon fibers (with a contact angle of 50°). To obtain more hydrophobic surface properties (with a contact angle of 130°), hydrophobic treatment with polytetrafluoroethylene (PTFE) is usually used to obtain high drainage performance [17, 18]. Optimizing GDL structural parameters with reasonable PTFE content and distribution is of great significance for improving the overall performance of fuel cells [19].

Various experimental works have been conducted to investigate the effects of PTFE content and distribution on the drainage and performance of fuel cells [20, 21]. However, the complex structure makes the reaction transport process in porous media much more complex than in a single media. To detect the gas-liquid two-phase transport process under the microstructure of GDL and obtain the influence mechanism of microstructure on drainage performance, microscopic visualization techniques such as neutron imaging, scanning electron microscopy (SEM), and X-ray technology were introduced [22, 23, 24]. Meyer et al. [24] studied GDL with different PTFE contents through neutron imaging and X-ray computed tomography. X-ray CT experimental results showed that PTFE forms like fibers attached to GDL carbon fibers, and high PTFE content will become an obstacle to water transport in GDL. Yu et al. [25] degraded GDL through accelerated stress testing and observed it through SEM. They report that the loss of carbon fibers and PTFE after GDL degradation leads to a decrease in hydrophobicity, seriously affecting the mass transfer characteristics of GDL. However, due to the limited spatial or temporal resolution of the existing experimental techniques, it is still difficult to study the transport process in microscopic porous media, especially the dynamic process of multiphase flow.

The science of studying the fluid transport process in porous media is called mechanics of flow through porous media [26, 27]. At present, the study on meso- and micro-level by using pore scale numerical simulation method has become the frontier topic in the field of microporous media study in the world, such as MC [28], PNM [29] and LBM [30, 31, 32] methods. LBM is suitable for simulating gas-liquid two-phase flow in microporous media structure, especially for simulating the heat and mass transfer process of gas-liquid electricity thermal coupling multiphysics simulation in porous media structure components of PEM fuel cell. Ira et al. [33] used pseudo-potential LBM to simulate the transport of liquid water in the microporous layer and diffusion layer of PEM fuel cell. The results showed that the content of hydrophobic solid skeleton had a great impact on the transport of gas-liquid two-phase in its interior and pointed out that GDL porous media with a certain degree of hydrophilic hydrophobic mixed structure had better liquid water transport performance. Han et al. [34] simulated the liquid water transport in the porous layer of PEM fuel cell through three-dimensional two-phase LBM. The results show that the porosity, pore structure, and liquid water saturation level in the porous layer have a great impact on the liquid water transport, which has guiding significance for the design of porous layer porous media structure optimization. Jithin et al. [35] used multi-relaxation LBM to simulate the fluid, heat, and material transport and oxygen reduction reaction in the cathode porous material assembly of PEM fuel cell. The numerical simulation considered the permeability, porosity, effective diffusion coefficient, and other properties of porous media, and the results showed that the solid skeleton material properties of the porous layer had an important impact on its internal heat and mass transport. Zhou et al. [36] built the three-dimensional porous microstructure of the gas diffusion layer and microporous layer of PEM fuel cell, and used LBM to simulate the coupling process of gas-liquid two-phase transport and oxygen reaction transport in the surface cracks. The simulation results show that the number of cracks has a great impact on the two-phase transport in porous media structure. When the cracks increase, liquid water covers the porous structure interface, which is not conducive to the transmission of liquid water.

This chapter briefly introduces the method of constructing microporous media structure by random algorithm and simulates gas-liquid two-phase transport in porous structure by SC pseudo-potential multiphase LBM. Then, it enumerates the influence of structural characteristics of PEM fuel cell GDL microporous media on two-phase transport and discusses the intrinsic relationship between its structural parameters and two-phase transport.

2. Construction of microporous media structure by random algorithm

The structure of porous media constructed by random algorithm can fully characterize its microstructure characteristics. At the same time, the model reconstructed by this method is completely based on program language or modeling software reconstruction. All structural parameters are flexible and controllable, and the cost is low. Therefore, the random algorithm is widely used in the study of simulating the heat and mass transfer process inside the microporous media.

GDL is a porous media material formed by a certain number of carbon fibers interlaced into a single layer and compressed to form a hydrophilic microstructure (50° ∼ 80°). To achieve more hydrophobic surface physical properties (100° ∼ 150°), GDL is usually placed in a certain mass fraction of PTFE solution for hydrophobic treatment. From the specific parameters of porous media structure, GDL is about 200 μm thick and carbon fibers are 7 μm in diameter, and its internal pore diameter varies from tens of micrometers to hundreds of micrometers, with a porosity of 60 ∼ 80%. This chapter takes carbon paper GDL as an example for reconstruction, and the idea of reconstructing GDL through random algorithms is as follows:

1. Generate randomly distributed straight lines in a plane with a specified area size to represent the centerline of carbon fibers;

2. Expand the fiber centerline in the first step into a cylinder in three dimensions, forming a single carbon fiber. The processing method is to determine whether the distance from the pore lattice around the centerline to the centerline is less than the radius of the carbon fiber. If so, the lattice will be transformed into a fiber lattice, otherwise, it will still remain a pore lattice;

3. Repeat the work of the first and second steps to generate a number of carbon fibers in the plane, and the carbon fibers are randomly interlaced until the porosity of the layer meets the requirements (only carbon fibers);

4. PTFE is randomly generated from the carbon fiber layer obtained in the third step of work, and the generation position is mainly at the position where the intersecting pores of carbon fibers and carbon fibers are small, until the PTFE content of the carbon fiber layer meets the requirements (the porosity of each layer of carbon fibers is ultimately determined by the PTFE content);

5. Repeat the first to fourth steps to generate the target number of carbon fiber layers and then overlay them to obtain the fiber structure of the target thickness.

Figure 1 shows the steps for randomly reconstructing the microstructure of GDL.

Figure 1.

(a) Distribution diagram and (b) flow diagram of GDL reconstruction.

3. Theory and application analysis of lattice Boltzmann method

LBM has better numerical stability and structural diversity in the study of multiphase transport processes in porous media structures. Generally, LBM can be divided into single relaxation time and multiple relaxation time according to the collision operator. Because the single relaxation time collision operator scheme has the advantages of simplicity and low computational cost, it is widely used in porous media flow simulation. This chapter mainly introduces the principle of Shan-Chen (SC) pseudo potential multiphase LBM with SRT collision operator and provides a theoretical basis for the subsequent numerical simulation of gas-liquid two-phase transport in GDL microporous media structure.

3.1 LBM content and introduction

For component k, the distribution function evolution equation for SRT and D2Q9 formats is:

$${f_{\mathrm{i}}}^{\mathrm{k}}\left(x+c_{\mathrm{i}}\mathit{\Delta t},t+\mathit{\Delta t}\right)-{f_{\mathrm{i}}}^{\mathrm{k}}\left(x,t\right)=-\frac{1}{\tau }\left[{f_{\mathrm{i}}}^{\mathrm{k}}\left(x,t\right)-{f_{\mathrm{i}}}^{\mathrm{eq},\mathrm{k}}\left(x,t\right)\right]$$E1

Where ${f_{\mathrm{i}}}^{\mathrm{k}}\left(x,t\right)$ is the density distribution function of component $k$ (oxygen and liquid water) in lattice direction $i$ at lattice position $x$ (unit lu) and lattice time $t$ (unit ts). In this chapter, the D2Q9 model, which is defined as a velocity model on a square, is used for two-dimensional simulation. The discrete velocity $c_i$ in 9 directions is:

$$c_i=\left\{\begin{array}{cc}0& i=0\\ \left(cos\left[\frac{\left(i-1\right)\pi }{2}\right],sin\left[\frac{\left(i-1\right)\pi }{2}\right]\right)& i=1,2,3,4\\ \sqrt{2}\left(cos\left[\frac{\left(i-5\right)\pi }{2}+\frac{\pi }{4}\right],sin\left[\frac{\left(i-5\right)\pi }{2}+\frac{\pi }{4}\right]\right)& i=5,6,7,8\end{array}\right.$$E2

For the DnQb model with single relaxation time (STR), the equilibrium distribution function ${f_{\mathrm{i}}}^{\mathrm{eq}}$ can be expressed as:

$${f_{\mathrm{i}}}^{\mathrm{eq}}={\omega }_{\mathrm{i}}\rho \left[1+\frac{1}{c^2}\left(c_{\mathrm{i}}\cdot u\right)+\frac{1}{2c^4}{\left(c_{\mathrm{i}}\cdot u\right)}^2-\frac{1}{2c^2}{u}^2\right]$$E3

Where $c_{\mathrm{s}}$ is the lattice sound velocity ($c_{\mathrm{s}}=\sqrt{\mathit{RT}}$), and for the D2Q9 model, the weight factor ${\omega }_{\mathrm{i}}$ is

$${\omega }_{\mathrm{i}}=\left\{\begin{array}{lll}\frac{4}{9}i& =& 0\\ \frac{1}{9}i& =& 1,2,3,4\\ \frac{1}{36}i& =& 5,6,7,8\end{array}\right.$$E4

In lattice units, $\tau$ (lattice unit ts) is defined as the relaxation time, which is related to the kinematic viscosity $\nu$ (lu²ts⁻¹):

$$\nu ={c_s}^2\left({\tau }_{\nu }-0.5\right)\mathit{\Delta t}$$E5

The macroscopic density $\rho$ and macroscopic velocity $u$ of fluid component $k$ can be obtained by the following equation:

$${\rho }^k=\sum {f_i}^k$$E6

$${\rho }^k{u}^k=\sum {f_i}^k{e}_i+{\tau }^k{F}^k$$E7

$$u^{\mathrm{k}}=u^{\text{'}}+\frac{{\tau }^{\mathrm{k}}{F}^{\mathrm{k}}}{{\rho }_{\mathrm{k}}}$$E8

For multi-component and multiphase flow problems, the forces acting on various fluid components $k$ are classified as surface tension $F_{1\mathrm{k}}$ between fluid and fluid, adhesion $F_{2k}$ between fluid and solid, and external force $F_{3k}$. $F^{\mathrm{k}}$ represents the total force acting on component $k$:

$$F^{\mathrm{k}}=F_{1\mathrm{k}}+F_{2\mathrm{k}}+F_{3\mathrm{k}}$$E9

Therefore, for the component $k$ particles at position $x$ subjected to the force of fluid particles from position $x^{\text{'}}$:

$$F_{1\mathrm{k}}\left(x\right)=-{\psi }_{\mathrm{k}}\left({\rho }_{\mathrm{k}}\left(x\right)\right)\sum _{\mathrm{x}\text{'}}\sum _{\overline{\mathrm{k}}}{G}_{\mathrm{k}\overline{\mathrm{k}}}\left(x,x^{\text{'}}\right){\psi }_{\mathrm{k}}\left({\rho }_{\mathrm{k}}\left(x\right)\right)\left(x^{\text{'}}-x\right)$$E10

Where $\psi \left({\rho }_{\mathrm{k}}\left(x\right)\right)$ is the pseudo-potential function (also representing the effective density), defined as:

$$\psi \left({\rho }_{\mathrm{k}}\left(x\right)\right)={\rho }_0\left[1-exp\left(-{\rho }_{\mathrm{k}}\right){\rho }_0\right]$$E11

Where $G_{\mathrm{k}\overline{\mathrm{k}}}$ represents the strength of the interaction between the fluid and the fluid, considering only the forces of the nearest neighboring particles. By controlling the value of $G_{\mathrm{k}\overline{\mathrm{k}}}$, the insolubility and surface tension of the two-phase fluid are controlled.

When the fluid $k$ particle at position $x$ is transported to position $x^{\text{'}}$, it is subjected to a force from adjacent solid particles:

$$F_{2\mathrm{k}}\left(x\right)=-{\rho }_{\mathrm{k}}\left(x\right)\sum _{\mathrm{x}\text{'}}{G}_{\mathrm{ks}}\left(x,x^{\text{'}}\right)n_{\mathrm{s}}\left(x^{\text{'}}\right)\left(x^{\text{'}}-x\right)$$E12

Where $G_{\mathrm{ks}}$ represents the intensity of fluid-solid interaction between fluid particles and solid surfaces, and the different wetting characteristics between fluid and solid surfaces can be achieved by adjusting the value of $G_{\mathrm{ks}}$. The $n_{\mathrm{s}}$ indicator function is used to distinguish each phase. When $n_{\mathrm{s}}=1$, it represents the solid phase, and when $n_{\mathrm{s}}=0$, it represents the liquid phase.

In this section, the SC pseudo-potential multiphase LB model is proposed to study the gas-liquid two-phase transport process in porous media and the influence of various structural design parameters on the dynamic transport process of liquid water.

3.2 LB method verification

Two numerical experiments were conducted in this section to validate the SC pseudo-potential multiphase LB method mentioned in Section 3.1 to simulate the physical phenomena of two-phase fluid flow in porous structures: one was a bubble experiment to obtain the gas-liquid interaction force $G_{\mathrm{k}\overline{\mathrm{k}}}$, and the other was a static contact angle experiment to determine the fluid-solid interaction force $G_{\mathrm{ks}}$, as shown in Eqs. (10) and (12).

3.2.1 Bubble test

Surface tension is caused by the adhesion of liquid molecules at the liquid gas interface. It plays a key role in determining the flow of two-phase fluid in porous media, and is closely related to capillary action and wetting phenomenon. For liquid droplets suspended in gas, according to the well-known Laplace law, the pressure difference at the liquid-gas interface is related to the radius of the bubble:

$$\mathit{\Delta P}=\frac{\sigma }{R}$$E13

Where $\sigma$ is the surface tension, and $\mathit{\Delta P}$ is the pressure difference between the inside and outside of the bubble, expressed in lm lu⁻¹ ts⁻².

A circular bubble with a radius of 20 lattices is initially located at the center of the 100 × 100 lattice domain, as shown in Figure 2. All outer sides are considered as periodic boundary conditions. Inside the droplet, the initial densities of air and water are set to 2.00 and 1.00 × 10⁻⁵, respectively. Outside the droplet, the initial densities of air and water are set to 2.00 and 0.00, respectively. The difference between the maximum and minimum density of a gas is a function of $G_{\mathrm{k}\overline{\mathrm{k}}}$. Figure 3 shows that when the value of $G_{\mathrm{k}\overline{\mathrm{k}}}$ reaches above 0.07, the density difference appears rapidly, reflecting the beginning of phase separation. Once $G_{\mathrm{k}\overline{\mathrm{k}}}$ exceeds 0.12, the density contrast will increase to the maximum value of 2.00. Therefore, $G_{\mathrm{k}\overline{\mathrm{k}}}=0.12$ is used in the simulation.

Figure 2.

Relationship between the internal and external pressure difference of the droplet and the reciprocal of the droplet radius.

Figure 3.

Variation in the difference between the maximum and minimum densities of the gas with different values of $G_{\mathrm{k}\overline{\mathrm{k}}}$.

The pressure difference ($\mathit{\Delta P}$) is at the liquid-gas interface (expressed in lattice units) and the bubble radius ($R$) at the final equilibrium state (expressed in lattice units). As shown in Figure 2, the results well comply with Laplace’s law.

3.2.2 Static contact angle test

The wettability of solid surfaces can be reflected by the contact angle, as described in Section 3.1. By adjusting the value of $G_{\mathrm{ks}}$, different wettability between fluid and solid surfaces can be achieved. In this simulation, semicircular droplets are placed on a horizontal solid plate to simulate the evolution of liquid shape on the solid wall at different $G_{\mathrm{ks}}$ values. The left and right sides of the calculation domain are considered as Periodic boundary conditions, while the upper and lower sides are set as no-slip boundary conditions. Similar to bubble testing, without any external force, the initial density of the fluid inside and outside the droplet is set to 2.00 and 1.0 × 10⁻⁵, respectively. After reaching a stable state, the shape of the droplet and the angle of the liquid-solid interface were recorded using the method proposed in the reference literature to obtain the contact angle. In Figure 4, an approximate linear relationship between the predicted contact angle and $G_{\mathrm{ks}}$ can be observed.

Figure 4.

Static contact angle as a function of $G_{\mathrm{ks}}$.

4. Results and analysis

4.1 Influence of surface characteristics of porous structures

The wettability of solid phase in porous structure is an important characteristic that affects its internal capillary pressure. In order to explore the influence of changes in the internal surface characteristics of porous media structure on gas-liquid two-phase transport, the wettability of porous structure is changed by changing the content of PTFE in GDL.

To investigate the effect of solid phase wettability on the gas-liquid two-phase transport inside GDL, this section first investigated the dynamic transport behavior of liquid water in traditional GDL with five different PTFE contents (0 wt., 5 wt., 10 wt., 15 wt., and 20 wt.%). The porosity of GDL without PTFE treatment (0 wt.%) was set to 60%, while the overall porosity of GDL decreased after PTFE solution treatment. The PTFE content and structural design parameters in GDL are shown in Table 1. Figure 5 shows three representative reconstructed two-dimensional microstructures of GDL with PTFE content of 0 wt., 10 wt., and 20 wt.%.

Case Number	PTFE content (wt.%)	Porosity (%)
Case 1	0 wt.%	60.00%
Case 2	5 wt.%	56.80%
Case 3	10 wt.%	53.60%
Case 4	15 wt.%	50.40%
Case 5	20 wt.%	47.20%

Table 1.

Conventional GDL with different PTFE contents.

Figure 5.

Three representative reconstructed of the two-dimensional microstructure of GDLs with PTFE contents of (a) 0 wt.%, (b) 10 wt.%, and (c) 20 wt.%.

This section discusses the dynamic behavior of two-phase fluid in five GDL porous structures with uniformly distributed PTFE content, and the structural design is shown in Table 1. Figure 6 shows a cross-sectional view of liquid water invading GDL with different PTFE contents in steady-state. In all cases shown in Figure 6, capillary fingering flow can be observed in liquid water, and capillary pressure is the main driving force for liquid transport in porous structures. Liquid water tends to invade pores with minimal capillary resistance, such as hydrophilic channels with larger pores. After breaking through GDL, liquid water will form a fixed liquid water flow channel, and subsequent liquid water will first break through this channel. From Figure 6(a) and (b), it can be seen that in untreated PTFE and GDL with low PTFE content, many large pores, even hydrophilic small pores, are occupied by liquid water. This can be attributed to the low capillary resistance caused by the hydrophilicity or low hydrophobicity of the solid phase, allowing liquid water to invade these pores. However, increasing the content of PTFE can lead to the need for liquid water to overcome greater capillary resistance, leading to the need for liquid water to break through larger pores rather than relatively smaller ones, which can effectively improve the transport of liquid water, as shown in Figure 6(b)–(d)).

Figure 6.

Profiles of liquid water invading conventional GDLs with different PTFE contents in the steady state.

From Figure 7(a), it can also be seen that the higher the PTFE content, the lower the liquid water saturation. As shown in Figure 7(b), as the PTFE content in GDL increases, the gas-liquid steady-state time decreases, indicating that the higher the hydrophobicity, the more favorable the drainage of GDL. However, the excessively hydrophobic porous structure significantly increases the entry pressure of liquid water, which may make the inlet side of GDL more susceptible to water flooding. It can be seen that evaluating the overall drainage performance of fuel cells solely based on liquid saturation is not enough. The main purpose of drainage is to enhance the gas mass diffusion of fuel cells. Therefore, effective porosity closely related to gas mass diffusion is used as another parameter for performance evaluation. Therefore, this study uses effective porosity as an additional criterion for evaluating GDL under different PTFE contents. As the PTFE content increases, the porosity of GDL decreases, but the effective porosity of GDL shows a trend of first increasing and then decreasing, as shown in Figure 7(a). When the PTFE content is 10 wt.%, the effective porosity of GDL is the highest, reaching 43.42%. Compared to GDL without PTFE and with a PTFE content of 20 wt.%, it increases by 14.80 and 6.90%, respectively.

Figure 7.

(a) Curves of effective porosity and liquid saturation and (b) steady-state time curves of gas and liquid for conventional GDLs with different PTFE contents.

Therefore, for the porous media structure with capillary pressure as the main driving force, the solid phase is too hydrophilic, which will lead to the retention of liquid water in the pores, affecting the transport and discharge of liquid, while the excessively hydrophobic will inhibit the invasion of liquid phase, and the pores are mostly occupied by gas.

4.2 Influence of pore distribution in porous structures

The distribution of pore and fixed phase wettability has a significant impact on its multiphase transport. This section will construct a porous microstructure with pore gradient distribution and discuss the dynamic changes in its internal two-phase transport.

In this chapter, a two-dimensional profile is obtained from the three-dimensional GDL porous media structure constructed using the random algorithm in the second section. The above random algorithm was used to design and reconstruct GDLs with PTFE gradient distribution, including PTFE bi-gradient and tri-gradient GDLs. The inlet and outlet regions of PTFE bi-gradient GDL have different PTFE contents, while the inlet, middle, and outlet regions of PTFE tri-gradient GDL have different PTFE contents. Figure 8 shows the representative reconstructed two-dimensional microstructure of bi-gradient and tri-gradient GDL with a total PTFE content of 10 wt.%. For bi-gradient GDL, as shown in Figure 8(a), the corresponding PTFE content in the inlet and outlet regions of GDL is 14 wt. and 6 wt.%. Similarly, Figure 8(b) shows the three-dimensional microstructure of GDL with different PTFE contents in the inlet, middle, and outlet regions.

Figure 8.

Representative reconstructed microstructure of bi-gradient and tri-gradient PTFE GDLs with a total PTFE content of 10 wt.%: (a) inlet regions: 14 wt.%, outlet regions: 6 wt.%, and (b) inlet regions: 16 wt.%, middle regions: 10 wt.%, outlet regions: 4 wt.%.

This section proposes six GDL structural designs for bi-gradient PTFE with a total PTFE content of 10 wt.% and constructs a GDL porous structure with gradient distribution of wetting characteristics. The design scheme and parameters of the bi-gradient GDL structure with PTFE content are shown in Table 2.

Case number	PTFE content (wt.%)	Region width (μm)	Porosity
Case 6; Case 7; Case 8; Case 9; Case 10; Case 11	in: 4, out:16	in: 50, out: 50	53.60%
	in: 6, out: 14	in: 50, out: 50	53.60%
	in: 8, out: 12	in: 50, out: 50	53.60%
	in: 12, out: 8	in: 50, out: 50	53.60%
	in: 14, out: 6	in: 50, out: 50	53.60%
	in: 16, out: 4	in: 50, out: 50	53.60%

Table 2.

Bi-gradient PTFE GDL with the total PTFE content of 10 wt.%.

Notes: “in” and “out” represent the inlet region and out region of the gradient PTFE GDL.

On this basis, the performance comparison simulation of PTFE bi-gradient GDL and traditional GDL (Case 3) was conducted using SC multiphase pseudo-potential LBM. Figure 9 shows the profile of PTFE content bi-gradient GDL (Case 6–11) and traditional GDL (Case 3) with PTFE content of 10 wt.% under stable state when liquid water invades two equally sized segmentation regions. It can be observed that liquid water invades and exits the GDL structure under capillary pressure. However, in Figure 9(a)–(c), the PTFE content in the treated bi-gradient GDL inlet area is relatively low, and the capillary resistance that liquid water needs to overcome when invading the pores of the GDL inlet area is low. The low hydrophobicity leads to relatively high adhesion, making liquid water occupy these pores and unable to be discharged. Therefore, compared to traditional GDL (Case 3), more liquid water occupies the pores of the inlet area. The liquid saturation in the inlet area is higher than that of traditional GDL structures. When liquid water begins to enter the pores of the outlet area with high PTFE content, the capillary resistance suddenly increases. In this case, the outlet area with high PTFE content acts as a capillary barrier, hindering the breakthrough of liquid water, which in turn increases the residual water content in the pores of the inlet area. Therefore, cases 6 and 7 have higher water saturation and lower effective porosity, and the breakthrough time of liquid water from GDL is longer. As shown in Figure 10(a) and (b), the liquid saturation in cases 6–8 is relatively high, and more liquid water occupies and invades the pores, thus hindering the gas mass diffusion ability, reducing the effective porosity. The structure is unfavorable to the drainage performance and gas diffusion performance of porous media.

Figure 9.

Profiles of liquid water invading the bi-gradient PTFE GDLs with two regions of the same size and the conventional GDL with a PTFE content of 10 wt.% at steady state.

Figure 10.

(a) Curves of effective porosity and liquid water saturation and (b) steady-state time curves of gas and liquid for the bi-gradient GDL with two regions of the same size and the conventional GDL with a PTFE content of 10 wt.%.

In Figure 9(e)–(g), compared to the traditional GDL structure design, the treated GDL has a relatively high PTFE content in the inlet area and a relatively high hydrophobic property. The capillary barrier formed in the inlet area has hindered the breakthrough and occupation of liquid water, resulting in less residual liquid water content in the pores of the inlet area. At the same time, the pores with lower hydrophobicity in the outlet area reduce the capillary resistance overcome by liquid water breakthrough, which is more conducive to the breakthrough and discharge of liquid water. As shown in Figure 9(a) and (b), compared to the performance parameters of traditional GDL (Case 3), PTFE with a higher content in the inlet area has a lower design liquid saturation and higher effective porosity. At the same time, the lower liquid equilibrium time in Figure 10(b) indicates that the time for liquid breakthrough and discharge of GDL is lower, and the drainage performance is improved. Therefore, a higher distribution of PTFE content in the inlet section is beneficial for the liquid transport of GDL. However, as shown in Case 11 in Figure 10(a), the high PTFE content in the inlet area means that the PTFE content in the outlet area is too low. In this case, its extremely low hydrophobicity leads to more liquid water occupying the pores of the outlet area, which is not conducive to the discharge of liquid water and the transportation efficiency of the gas. In addition, a high PTFE content in the inlet area can also lead to a high inlet pressure when liquid water enters GDL, hindering the liquid water from entering GDL from the inlet side and discharging into the fuel cell, thereby making the inlet side more susceptible to water flooding.

Similarly, this section discusses three GDL structures with tri-gradients of PTFE content, with higher PTFE content distributed near the inlet area. The PTFE content in the middle region remains at 10 wt.%. Therefore, when the average PTFE content is 10 wt.%, the PTFE content in the export area is determined by the PTFE content in the import area. The tri-gradient structure design is shown in Table 3.

Case number	PTFE content (wt.%)	Region width (μm)	Porosity
Case 12	in: 12, mid: 10, out: 8	in: 33.3, mid: 33.3, out 33.3	53.60%
Case 13	in: 14, mid: 10, out: 6	in: 33.3, mid: 33.3, out 33.3	53.60%
Case 14	in: 16, mid: 10, out: 4	in: 33.3, mid: 33.3, out 33.3	53.60%

Table 3.

Tri-gradient PTFE GDL with a total PTFE content of 10 wt.%.

Notes: “in,” “mid,” and “out” represent the inlet region, middle region, and outer region of the gradient PTFE GDL, respectively.

Figure 11 shows the GDL structure with three gradients of PTFE content when liquid water invades three regions of the same size during steady state, as well as a cross-sectional view of the traditional GDL of Case 3 and the GDL with bi-gradients of PTFE content of Case 10. According to the discussion and explanation in Section 4.2, a high PTFE content (14 wt.%) in the inlet area is beneficial for the drainage of GDL, while a low or high PTFE content in the inlet area is not conducive to the liquid transport of GDL. Similar conclusions can be applied to the GDL of PTFE with three gradients, as shown in cases 12 and 13 in Figure 11. However, for the PTFE content three gradient GDL, a reasonable PTFE content in the inlet area (14 wt.%), compared to the traditional GDL in Case 3 and the optimal PTFE bi-gradient GDL in Case 10, resulted in an increase in effective porosity of 5.80 and 2.00%, respectively, as shown in Figure 12. Due to the presence of an intermediate region with a PTFE content of 10 wt.%, the size of the outlet region with relatively low PTFE content is reduced, resulting in less residual water in the pores of the outlet region. In addition, the middle region of the PTFE tri-gradient GDL reduces the PTFE content gradient difference between the inlet region (14 wt.%) and the outlet region (6 wt.%). In this case, the PTFE tri-gradient GDL generates relatively smooth capillary resistance, which may be beneficial for the drainage of liquid water. As shown in Figure 12(b), the steady-state time curves of gases and liquids can also support this point.

Figure 11.

Profile of liquid water intrusion into GDL with PTFE tri-gradient under stable-state time (conventional GDL in case 3 and PTFE bi-gradient GDL in Case10).

Figure 12.

Steady-state time of (a) curves of effective porosity and liquid saturation degree and (b) gas and liquid for the different tri-gradient GDLs with three regions of the same sizes, the conventional GDL of case 3 and the bi-gradient PTFE GDL of case 10 at steady state.

Therefore, the gradient distribution of wettability will make the capillary pressure distribution in local regions uneven, or even form a capillary barrier, thus changing the two-phase transport state. It is of great significance to reasonably distribute the wettability of solid phase for improving the internal transport of porous media.

4.3 Effects of porous structure aging

The porous structure will cause structural damage and changes in the surface characteristics of the solid phase due to corrosion, shedding, and other conditions, thus affecting its internal multiphase flow state. In order to obtain the influence mechanism of the aging of porous media structure on two-phase transport, this section will conduct aging treatment and study on GDL structure.

To study the influence mechanism of porous media aging on liquid water transport behavior, the GDL of carbon fiber and PTFE degraded at the same time with different content of PTFE was reconstructed. Assuming that the final degradation amounts of carbon fiber and PTFE after 20 aging cycles in GDL are 30 and 15%, respectively, four different aging levels of GDL were reconstructed using exponential rate random degradation for 5, 10, 15, and 20 aging cycles, simulating liquid water transport in GDL under different aging levels. Figure 13(a)–(c) shows the reconstruction results of fresh GDL and aged GDL subjected to 10 cycles of aging. Figure 13(f) shows a large image of the local GDL component degradation. Compared to fresh GDL, aged GDL exhibits corrosion of carbon fibers and detachment of PTFE. The structural design scheme is shown in Table 4.

Figure 13.

(a)–(c) fresh GDL and GDL with different aging processes, and (f) enlarged view of local aging.

PTFE Content (wt.%)	Case Number	Carbon fiber degradation amount (%)	PTFE degradation amount (%)
5 wt.%	Case 15	0%	0.00%
	Case 16 (5-times)	2%	1.61%
	Case 17 (10-times)	6%	4.20%
	Case 18 (15-times)	14%	8.37%
	Case 19 (20-times)	30%	15.00%
10 wt.%	Case 20	0%	0.00%
	Case 21 (5-times)	2%	1.61%
	Case 22 (10-times)	6%	4.20%
	Case 23 (15-times)	14%	8.37%
	Case 24 (20-times)	30%	15.00%
15 wt.%	Case 25	0%	0.00%
	Case 26 (5-times)	2%	1.61%
	Case 27 (10-times)	6%	4.20%
	Case 28 (15-times)	14%	8.37%
	Case 29 (20-times)	30%	15.00%

Table 4.

Fresh GDL and aged GDL with different initial PTFE contents under exponential rate degradation.

Figure 14 shows the changes in the transport behavior of liquid water in fresh GDL and aged GDL as the aging degree increases when the initial PTFE content is 5 wt., 10 wt., and 15 wt.%. The increase in initial PTFE content can occupy some pores, resulting in a decrease in total porosity and an increase in capillary pressure. In addition, more carbon fibers are covered with PTFE, leading to an increase in the hydrophobic pathway of GDL and a further increase in capillary pressure. Therefore, GDL with a high initial PTFE content produces lower liquid water saturation than GDL with a low initial PTFE content. As GDL ages, higher PTFE content also leads to lower liquid water saturation in GDL compared to GDL structures with lower PTFE content, as shown in Cases 19, 24, and 29 in Figure 14.

Figure 14.

Profile of liquid water invading a fresh GDL and GDLs aged with initial PTFE contents of (a) 5 wt.%, (b) 10 wt.%, and (c) 15 wt.%.

From Figure 15(a), it can be seen that increasing the initial PTFE content can indeed reduce the liquid water saturation in GDL, but it will also occupy the overall porosity of GDL and reduce the effective porosity, thereby reducing the gas mass diffusion ability. In traditional GDL designs with different initial PTFE contents, the saturation of liquid water increases with increasing aging degree, while the corresponding effective porosity continues to decrease. However, during the entire aging process of GDL, GDL with an initial PTFE content of 10 wt.% showed the highest effective porosity, corresponding to the best overall performance. At the end of aging, the liquid saturation of GDL with a PTFE content of 10 wt.% increased by 10.20% compared to the initial GDL, but the effective porosity decreased by 4.20%. When the PTFE content is 15 wt.%, the liquid saturation increases by 9.1% and the effective porosity decreases by 4.30%. In addition, excessively hydrophobic GDL can lead to high inlet pressure of liquid water, resulting in liquid water generated by electrochemical reactions in CL remaining in CL and not being discharged from the fuel cell interior through GDL. From Figure 15(b), it can be seen that GDL with an initial PTFE content of 10 wt.% has a higher gas-liquid steady-state time and better drainage performance.

Figure 15.

Curves for (a) liquid saturation and effective porosity and (b) gas and liquid steady-state time for a fresh GDL and GDLs aged with initial PTFE contents of 5 wt.%, 10 wt.%, and 15 wt.%.

The above results show that aging will destroy the structure of porous media, reduce the overall capillary pressure, and even change the hydrophilic and hydrophobic properties of local solid phases, making liquid water accumulate in local pores and unable to be discharged, seriously affecting the multiphase transport efficiency inside porous media. Therefore, improving durability is an important way to ensure the internal transport of porous media structures.

5. Conclusion

This chapter introduces the construction of microporous media structure through random algorithm and uses the SC pseudo potential LB method to study multiphase transport in porous media. Taking the GDL structure of PEM fuel cell as an example, the structural characteristics of porous media, such as solid phase wettability, pore distribution, and structural aging, are deeply discussed, and structural aging, were discussed in depth, and the following conclusions were obtained:

1. The higher the PTFE content in GDL porous media structure, the shorter the liquid water penetration time, but the effective porosity of GDL decreases. When the PTFE content reaches 10 wt.%, GDL has the maximum effective porosity.

2. For gradient GDL, a reasonably high PTFE content near the inlet will reduce the liquid water penetration time, and the effective porosity of bi-gradient (14 wt.%/6 wt.%) and tri-gradient (14 wt.%, 10 wt.%, and 6 wt.%) GDL will increase by 4.20 and 5.60%.

3. The initial PTFE content in GDL will affect the anti-aging performance of GDL. Higher PTFE content can ensure that GDL has better anti-aging performance, but too high initial PTFE content will also affect the gas mass diffusion ability.

This chapter obtained the intrinsic relationship between structural characteristic parameters and multiphase fluid flow is obtained, this lays a theoretical foundation for studying gas-liquid transport in microporous structures.

Acknowledgments

We thank the supports for this project provided by National Natural Science Foundation of China (No. 52176084) and National Key Research and Development Program of China (No. 2022YFE0207600).

Conflict of interest

The authors declare no conflict of interest.

Appendices and nomenclature

$c_s$	speed of sound in lattice, lu ts^−1
$d$	pore diameter, lu
$e_i$	velocity of particles in lattice direction $i$, lu ts^−1
$F^k$	total force acting on $k$ component, lm lu ts^−2
$F_{1k}$	fluid/fluid interaction force acting on $k$ component, lm lu ts^−2
$F_{2k}$	fluid/solid interaction force acting on $k$ component, lm lu ts^−2
$F_{3k}$	other external forces acting on the $k$ component, lm lu ts^−2
$f_i^k$	density distribution function of $k$ component in lattice direction $i$
$f_i^{\mathit{eq},k}$	density equilibrium distribution function of $k$ component in lattice direction $i$
$G_{k\overline{k}}$	cohesion factor between components $k$ and $\overline{k}$
$G_{\mathit{ks}}$	adhesion factor between component $k$ and solid
$J\left(s\right)$	Leverett J function
$P$	pressure, lm lu^−1 ts^−2
$\mathit{\Delta P}$	pressure difference inside and outside the bubble, lm lu^−1 ts^−2
$R$	droplet radius, lu
$t$	lattice time, ts
$\mathit{\Delta t}$	lattice time step, ts
$u^{\prime }$	macroscopic velocity, lu ts^−1
$u^{\mathit{eq}}$	lattice velocity of component $k$, lu ts^−1
$u^k$	equilibrium velocity of component $k$, lu ts^−1
${\omega }_i$	weight factor in lattice direction $i$
$x$	lattice position, lu

Greek symbols

$\tau$	relaxation time, ts
$\nu$	kinematic viscosity, lu^2 ts^−1
$\rho$	density, lm lu^−3
$\psi$	effective density
$\sigma$	surface tension
$\theta$	static contact angle
$\epsilon$	porosity

Lattice unit

$\mathit{lm}$	lattice unit of mass
$\mathit{lu}$	lattice unit of length
$\mathit{ts}$	lattice time step