Effect of Pore Shape and Spacing on Water Droplet Dynamics in Flow Channels of Proton Exchange Membrane Fuel Cells

Abstract

Effective water management increases the performance of proton exchange membrane fuel cells (PEMFCs). The liquid droplet movement mechanism in the cathode channel, the gas-liquid two-phase flow pattern, and the resulting pressure drop are important to water management in PEMFCs. This work employed computational fluid dynamics (CFD) with a volume of fluid (VOF) to simulate the effects of two operating parameters on the liquid water flow in the cathode flow channel: Gas diffusion layer (GDL) pore shape for water emergence, and distance between GDL pores. From seven pore shapes considered in this work, the longer the windward side of the micropore is, the larger the droplet can grow, and the duration of droplet growth movement will be longer. In the cases of two micropores for water introduction, a critical pore distance is noted for whether two droplets coalesce. When the micropore distance was shorter than this critical value, different droplets coalesce after the droplets grew to a certain extent. These results indicate that the pore shape and the distance between pores should be accounted for in future simulations of PEMFC droplet dynamics and that these parameters need to be optimized when designing novel GDL structures.

1. Introduction

Proton exchange membrane fuel cells (PEMFCs), as a clean, efficient, and environmentally friendly energy generation device, have been used in the fields of transportation, aerospace, and communication [1]. However, high cost and short life still limit the commercialization of fuel cells [2,3]. The transport and removal of water is critically important to the performance, stability, and durability of the fuel cell [4]. An important water management issue is the gas–liquid flow as water enters the reactant flow field channels, typically mini-channels [5].

The movement of liquid water in the cathode channels of fuel cells and the formation of gas–liquid two-phase flow are important to water management in PEMFCs. Researchers have studied droplet motion mechanisms of liquid water entering the gas channel through gas diffusion layers (GDLs) by experimental efforts and numerical simulation. Ye et al. [6] experimentally studied the effect of cell structure on the two-phase flow. A parallel channel flow system was machined from two acrylic blocks, and the GDL was sandwiched between the two blocks. Liquid water was injected into the channel through two 100 μm inlets. The airflow bypassed the channel blocked by the water. The resulting slug flow caused severe maldistribution and significant fluctuations in the pressure drop. Hussaini et al. [7] used an optically transparent fuel cell test platform to visually study the liquid water on the cathode side of the fuel cell. They determined the gas–liquid two-phase flow pattern in the fuel cell channel and developed the first flow regime map with slug, film, droplet, single-phase flows. Anderson et al. [8] proposed an additional flow pattern, the accumulation regime, from a large collection of experimental data from the literature. Bazylak et al. [9] conducted an experiment to study the dynamics of liquid water entering the gas channel through the GDL in a fuel cell. They found that the droplets tended to fix near the breakthrough position during the process of forming slug flow when the droplets detached from the GDL surface. This increased the liquid water content in the channel.

Although the flow of liquid water can be observed directly in experiments, the high cost and complexity of the experimental setup have made it important to optimize the water management of fuel cells through simulation [10,11]. Mathematical modeling and numerical simulation can theoretically explore the dynamic parameters of gas–liquid two-phase flow in the channel and clearly observe the flow area distribution in the flow channel. Quan et al. [12] first applied computational fluid dynamics simulation with the volume of fluid (VOF) model for interface tracking to study gas–liquid two-phase flow in a PEMFC and investigate the blockage of liquid water in curved areas. The VOF method possesses the capability of tracking a distinct interface between immiscible fluids by capturing the effect of surface tension, which is an important force under flow conditions relevant to fuel cell operations. The gas–liquid two-phase flow pattern formed in the flow channel has been studied. Cai et al. [13] studied the influence of wettability of the channel wall on the removal effect of liquid water. It is found that the wettability of the channel wall has a great impact on water transport in the channel, the distribution of water along the channel, and the removal time. Zhan et al. [14] found that a hydrophobic GDL and hydrophilic channel walls is the best combination, which is not only good for water drainage, but also for diffusion of the reaction gas from the channel to the catalytic layer. Moosa et al. [15] simulated vertical and horizontal fuel cells to study the influence of gravity on the gas–liquid two phase flow in a single-serpentine flow field. Ding et al. [16,17] used a VOF model combined with a 1D model to form a membrane electrode to study the effect of two-phase flow on the electrochemical performance. Comparing the single-phase flow model against two-phase flow models, it was found that the existence of slug flow reduces the output voltage due to the resultant limitation in mass transfer. Jiao et al. [18,19] designed gas channel models of parallel serpentine channels and straight channels to study the influence of different channel geometries on water drainage performance. The results showed that even a little water in both channels can cause an uneven distribution of the gas. Zhang et al. [20] used the Eulerian-Eulerian model and showed that the serpentine flow field is more favorable to the removal of liquid water than the parallel flow field, which improved the performance of fuel cells. Magesh et al. [21] designed a sinuous flow field by increasing the size of the PEMFC from 25 cm² to 100 cm². The results showed that the newly developed flow field was more beneficial to water removal than a serpentine flow field.

As discussed, many researchers have studied the influence of wall wettability, gas–liquid velocity and channel geometry on liquid water flow in the fuel cell’s channel. The microstructure of the GDLs also has a significant effect on the flow of liquid water in the channel. Real GDL micropores usually have an irregular geometry and are not simply circular or cylindrical; however, most simulation studies have introduced liquid water through a simplified single circular GDL micropore. In addition, most numerical simulations in literature only consider a simplified single droplet model. In the actual operation of fuel cells, the interaction between multiple droplets has an important impact on the liquid water flow. Understanding the emergence and interactions of multiple droplets from the GDL surface is still lacking in the literature. Therefore, in this work, the VOF model combined with CFD simulation was used to study the influence of the pore shape on the liquid water behavior in the fuel cells and the effect of the pore spacing (distance between pores). Liquid water introduced from two pores was studied with two water inlets set on the same straight line along the direction of airflow. The mechanism of the movement of the two droplets in the mini-channel and its influential factors were discussed.

2. Models and Methods

2.1. The VOF Model

The VOF multiphase flow model simulates two or more non-miscible fluids by solving a set of momentum equations and tracking the volume fraction of each fluid passing through the computational domain. It is assumed in the VOF model that two or multiple phases are not interspersed with each other, and each additional phase in the calculation region can be represented by phase volume fractions. In each governing equation, the total volume fraction of all phases is 1. This work uses air and water as the two fluids, obtains the instantaneous phase interface of arbitrary gas–liquid two phases by using the VOF model. The model assumes that air is a compressible ideal gas and water is an incompressible fluid. The VOF model equations are written as follows:

Continuity equation:

$$\frac{\partial {\alpha }_q}{\partial t}+\overrightarrow{v}\cdot \nabla {\alpha }_q=\frac{S_{{\alpha }_q}}{q},{\displaystyle \sum }_{q=1}^2{\alpha }_q=1$$

(1)

Mixture properties:

$$\rho ={\alpha }_2{\rho }_2+\left(1-{\alpha }_2\right){\rho }_1$$

(2)

$$\mu ={\alpha }_2{\mu }_2+\left(1-{\alpha }_2\right){\mu }_1$$

(3)

Momentum equation:

$$\frac{\partial }{\partial t}\left(\rho \overrightarrow{v}\right)+\nabla \cdot \left(\rho \overrightarrow{v}\overrightarrow{v}\right)=-\nabla p+\nabla \cdot \left[\mu \left(\nabla \overrightarrow{v}+\nabla {\overrightarrow{v}}^T\right)\right]+\rho \overrightarrow{g}+\overrightarrow{F}$$

(4)

Energy equation:

$$\frac{\partial }{\partial t}\left(\rho E\right)+\nabla \cdot \left(\overrightarrow{\nu }\left(\rho E+p\right)\right)=\nabla \cdot \left(k_{eff}\nabla T\right)+S_h$$

(5)

where α is volume fraction, ρ is the density (kg/m³), P is the pressure (Pa), μ is the viscosity (Pa·s), t is time (s), v is the velocity (m/s), k_eff is effective thermal conductivity (W·m⁻²·K⁻¹), and E is the energy (J).

A key advantage of the VOF model is that it can examine the effect of surface tension on fluid flow in microchannels. The VOF model includes two kinds of surface tension models: The continuum surface force model (CSF) and surface stress model (CSS) [22]. Brackbill et al. [23] proposed the CSF model to explain the effect of surface tension on the source term $\overrightarrow{F}$ of the momentum equation in the VOF model and is determined by the coefficient of surface tension (σ) and surface curvature (κ). It is a volume force expressed by the divergence theorem. The influence of the wall adhesion force was also considered in the VOF model. The equations of surface tension in the model are shown as follows. The curvature (κ) is defined in terms of the divergence of the unit normal:

$$\kappa =\nabla \cdot \left(\frac{n}{\left|n\right|}\right)$$

(6)

where n is surface normal, the gradient of α.

$$n=\nabla {\alpha }_q, p_2-p_1=\sigma \left(\frac{1}{R_1}+\frac{1}{R_2}\right)$$

(7)

Terms p₁ and p₂ are the pressures in the two fluids on either side of the interface and R (m) is the radii in the orthogonal direction. The forces from the CSF and CSS model are:

$$F\left(CSF\right)={\sigma }_{ij}\frac{\rho {\kappa }_i\nabla {\alpha }_i}{\frac{1}{2}\left({\rho }_i+{\rho }_j\right)}, F\left(CSS\right)=\nabla \cdot \left[\sigma \left(\left|\nabla \alpha \right|I-\frac{\nabla \alpha \otimes \nabla \alpha }{\left|\nabla \alpha \right|}\right)\right]$$

(8)

where σ is the surface tension (N/m) and κ is the curvature.

The equations of the wall adhesion force model are:

$$\widehat{n}={\widehat{n}}_{w}cos{\theta }_w+{\widehat{t}}_{w}sin{\theta }_w$$

(9)

where ${\widehat{n}}_w$ and ${\widehat{t}}_w$ are the unit vectors normal and tangential to the wall, and ${\theta }_w$ is the contact angle at the wall.

2.2. Model Settings and Boundary Conditions

In this work, the simplified geometric model was based on simplifying the cathode channel in the experimental device of Cheah et al. [24,25]. The simplified channel model is shown in Figure 1. The geometric size of the channel was 60.0 × 1.0 × 1.0 mm³, and the cross section was a regular quadrilateral. The water inlet micropore(s) were located at the center of the GDL, 20 mm from the gas inlet. In the case of two water inlets, the downstream micropore was located at a certain distance along the x-direction of the upstream micropore. GAMBIT software was used to divide the computing area. The grid diagram is shown in Figure 2.

The boundary conditions were set as follows: The gas velocity inlet was set at x = 0, the pressure-outlet was set at x = 60, the liquid velocity-inlet was marked at X₁, Y₁, and Z₁, the channel bottom was set with a hydrophobic GDL (135°), other surfaces were set as hydrophilic walls (45°).

2.3. Solution Procedure

Under the same simulation parameters, four grid sizes were tested: 1 × 5 × 5 (10⁻⁵ m), 1 × 5 × 10 (10⁻⁵ m), 10 × 10 × 10 (10⁻⁵ m), and 20 × 20 × 20 (10⁻⁵ m). The corresponding number of grids was 2.53 million, 1.55 million, 550,000, and 94,000, respectively. For the four mesh sizes, liquid phase fraction cloud maps were compared with different mesh precision. To obtain a clear phase interface and a suitable flow pattern and save computing resources as much as possible, the mesh size of 1 × 5 × 5 (10⁻⁵ m) was selected for this work. Under the unsteady state condition, the convergence of simulation can be obtained by selecting the time step that satisfies the calculation. According to the time step mentioned by Wu et al. [26], the maximum time step is 1.5 × 10⁻⁶ s under the condition that the maximum Courant number is 0.25 and the minimum grid size is 6.189 × 10⁻⁶ m. In the simulation process, we chose a time step of 10⁻⁶ s to meet the accuracy and time requirements of the simulation.

The commercial CFD software FLUENT 14.5 (Fluent, Inc., New York, NY, USA) was used to conduct a single-precision numerical simulation of gas-liquid two-phase flow on the 3D channel geometry model. The solution flow process is the laminar flow state, the solver is set as unsteady state mode, and the governing equation is the VOF model. The pressure-velocity coupling was obtained by using the PISO scheme. The pressure interpolation was obtained by using the PRESTO! Scheme, and the second-order upwind differencing scheme was used for the momentum equation. Geo-Reconstruct was used for the volume fraction calculation.

3. Results and Discussion

3.1. Effect of the Pore Shape

After the liquid water enters the channel through the pores, small spherical droplets form on the surface of the GDL and continue to grow. The droplet growth is affected by shear force, viscous force, inertial force, surface tension, pressure, and gravity [27,28]. Related dimensionless numbers that represent the ratio of surface tension to other forces are the Capillary number ($Ca=\left({\mu }_G{U}_G/\sigma \right)$), Weber number ($We={\rho }_G{U}_G^2{D}_h/\sigma$), and Bond number ($Bo=\left({\rho }_{L}g{D}_h^2/\sigma \right)$). Since the pore size was similar in this work, the water injection pore diameter 2R was used as the characteristic length D_h. The following dimensionless numbers were obtained: Ca = 2.46 × 10⁻⁴, We = 1.68 × 10⁻³, Bo = 1.34 × 10⁻³. From these results, the effect of surface tension is much greater than the other forces, which explains why liquid water forms near-spherical droplets on the surface of hydrophobic GDLs.

3.1.1. Effect of the Pore Shape on Droplet Volume and Cycle Time

The characterization of GDL micropores is generally based on the porosity measurement to obtain the pore size distribution, but the pore size distribution may not be the only factor affecting the performance of GDLs. Parikh et al. [29] characterized and compared the pore diameter, pore shape, direction, and distribution of three types of GDLs, and defined the roundness of the pore $S=4\pi A/P^2$, to express the pore shape more clearly, where A is the pore area and P is the pore circumference. Most existing research used simplified the pore into a circle or a square to study the flow of droplets in the cathode channel. In this section, the influence of the pore shape on the flow of a single droplet is discussed. Seven pore shapes for the micropores were investigated. The circular pore is 0.1 mm in diameter, and the area of the remaining pores was the same as that of the circular pore. Detailed simulation parameters for the seven pore shapes are provided in

Table 1. The pore shape parameter with area kept constant in all cases.

Case	Shape of Pore	A (mm2)	P (mm)	$\mathit{S} = 4\mathit{\pi }\mathit{A}/{\mathit{P}}^2$	${\mathit{L}}_{\mathit{w}}\left(\mathbf{mm}\right)$	UG (m/s)	UL (m/s)	$\mathit{\theta }(°)$
1		0.00785	0.3142	1.000	0.1571	1.0	0.1237	>90
2		0.00785	0.3300	0.906	0.1100	1.0	0.1273	120
3		0.00785	0.3380	0.863	0.1352	1.0	0.1273	126
4		0.00785	0.3544	0.785	0.0886	1.0	0.1273	90
5		0.00785	0.4038	0.604	0.1346	1.0	0.1273	120
6		0.00785	0.4430	0.502	0.1772	1.0	0.1273	90
7		0.00785	0.4430	0.502	0.0443	1.0	0.1273	90

(U_G and U_L represent the inlet rate of gas and liquid water; L_w represents the length of windward side of the pore; θ is the angle between the airflow and L_w).

Figure 3 shows critical droplet sizes for the seven cases. The droplet is considered the critical droplet when it is about to collide with the wall or detach from the pores and move along the GDL surface, and the droplet size is called critical crushing size. It can be seen in this figure that the critical volume of liquid droplets decreases in an order of case 4 > case 6 > case 7 > case 1 > case 5 > case 2 > case 3. The critical droplet size is closely related to the droplet movement cycle. Therefore, the droplet movement cycle can also reflect the critical volume of the droplet.

Figure 4 shows that the droplet movement cycle has the same trend as the critical droplet volume, and it decreases in the same order as discussed above. Therefore, the larger the critical size of the droplet, the longer the movement period. A droplet movement cycle refers to an entire droplet growth stage including emergence, growth, deformation, and fragmentation. At the beginning of droplet formation, the surface tension dominated and kept the droplet nearly spherical due to its small size. With the growth of the droplet, the blocking effect of the droplet on the airflow increased, causing the droplet to deform and detach under the action of the shear force. The droplet growth process is mainly affected by shear force, viscous force, inertial force, surface tension, pressure, and gravity [27,28,30,31]. As shown in

Table 2. Dimensionless number.

Name	Meaning	Formula	Value
Capillary number	Viscous force/surface tension	$Ca={\mu }_G{U}_G/\sigma$	$2.46\times {10}^{-4}$
Weber number	Inertial force/surface tension	$W{e}_G={\rho }_G{U}_G^2{D}_h/\sigma$	$1.68\times {10}^{-3}$
Bond number	Gravity/surface tension	$Bo={\rho }_{L}g{D}_h^2/\sigma$	$1.34\times {10}^{-3}$

, μ_G, U_G, σ, ρ_G, ρ_L, g, D_h are gas viscosity, gas inlet velocity (equal to gas apparent velocity and gas flow rate), surface tension coefficient, gas density, liquid density, acceleration of gravity, and the hydraulic diameter, using the dimensionless number to correlate these forces. Regarding the pore diameter 2R as the hydraulic diameter Dh, the calculation results in the table confirm that the surface tension in the smaller channels and GDL micropores plays a major role in the droplet size change. Since the critical droplet size is different in the above cases, it is concluded that the pore shape has a great effect on the droplet growth.

As shown in Table 1, the windward side length L_w of the water inlet pore decreases in the order of case 6 > case 1 > case 3 > case 5 > case 2 > case 4 > case 7. The angle $\theta$ between the windward side of the pore and the direction of the airflow, which is defined as the windward side slope, has a different increasing trend in the order of case 6 = case 4 = case 7 > case 1 > case 5 = case 2 > case 3. The results show that although the length of the windward side is the smallest in case 7, the critical volume of the droplet is larger than the non-quadrangular pore condition due to the larger slope of the windward side. This result indicates that the increasing trend of the critical volume of the droplet follows the law of windward slope and the important influence of the windward slope on droplet movement.

It can be seen that when the slope of the windward side is all 90 degrees, the difference of windward side length affects the critical volume of droplets. By comparing case 6 with case 7, the longer the windward side, the larger the critical droplet volume. When the angle of the windward side is all 120 degrees, the windward side length of case 5 is larger than that of case 2, which also supports the conclusion that the longer the windward side is, the larger the critical crushing volume of liquid droplets. Although an increase in the length of the windward side increases the adhesive force, it also makes the transverse diameter of the droplet larger. The obstruction effect on airflow increases and the drag force of airflow increases accordingly. When the drag force increase is greater than the adhesion force, the large length of the windward side will shorten the droplet movement cycle. This can explain why the critical droplet volume of case 4 is larger than that for cases 6 and 7.

It can also be observed from the schematic diagram of the water inlet pore that the slope of the windward side of the water inlet pore and the length of the windward side of the water inlet pore are determined by the direction of pore placement and the roundness of the pore. Therefore, it is concluded that the droplet motion period is closely related to the pore direction and roundness. From these simulations, when the pore shape is symmetrical and there is a corner rushing toward the flow direction, roundness can be used as a parameter to evaluate the critical volume of liquid droplets.

3.1.2. Effect of the Pore Shape on Downstream Liquid Flow Pattern

Figure 5 shows the liquid flow pattern in the channel at t = 1.0 s. With the change of the shape of the pore, the liquid flow pattern in the channel will change correspondingly. In the cases of square and rectangular pores, the longer droplet movement period and the larger critical droplet volume results in additional channel obstruction. The removed droplet in the last period moves a long distance downstream of the channel under the greater airflow drag force, so it will not merge with the broken droplet in this period. Due to this, a new droplet is formed after each droplet in the channel is removed from the pore. As noted previously, the liquid droplet movement period of the circular pore is shorter than that of a quadrilateral, and the critical volume and channel obstruction of the liquid droplet are also smaller. Therefore, the broken droplets in the last cycle move downstream for a relatively short distance under the relatively weak drag force of airflow, which coalesce with the detached droplets in this cycle and eventually accumulates to form a corner flow. In other cases, with the shortened droplet movement cycle and the decline of the critical droplet volume, the broken droplet will merge and accumulate into a larger droplet. This droplet moves downstream under the constant drag force of the air stream, becoming an independent droplet when the droplet moves beyond the range covered by the smaller droplet generated in the new cycle.

3.1.3. Effect of the Pore Shape on Channel Pressure Drop

Changes in the droplet morphology can also be identified by changes in the liquid water coverage on each wall and the pressure drop at the inlet and outlet of the channel. The coverage of liquid water on the GDL surface is an important parameter, as excess water coverage negatively impacts the performance of the fuel cell. The inlet and outlet pressure drop can reflect the blockage degree of liquid water to airflow. Low coverage of liquid water on the GDL surface and low-pressure drop are favorable for improving the cell performance [32].

Liquid water coverage on the GDL surface is closely correlated with the droplet motion cycle. According to Figure 4 and Figure 6, although case 4 and case 6 have the longest cycles, the maximum liquid coverage and average time on the GDL surface are smaller compared to other cases. It is known that in the case of large channel obstruction, to ensure sufficient passage of airflow, the contact surface between the droplet and the GDL surface will decrease to a minimum because of the large drag force of airflow, and the droplet height will increase. This phenomenon is mainly due to the channel blockage caused by the critical droplet volume under the condition of a quadrilateral pore being stronger than that of other shapes of pores.

The simulation calculates the final moment of liquid water coverage of SUM-wall (SUM-wall is the sum of the average liquid water coverage on the other three walls except for the GDL surface). The liquid water coverage on the wall in the case of a circular pore is much higher than that for other cases, while the rectangular and quadrilateral pores lead to walls with the least water coverage. This phenomenon indicates that as the number and size of the corner flow increases, liquid water coverage decreases.

The pressure drop between the inlet and outlet of the channel is shown in Figure 7. It can be seen in this figure that the pressure drop has the same trend as the droplet movement cycle. Long droplet movement cycles increase the accumulation of liquid water in the channel and the blockage of the channel, resulting in an increase in the overall channel pressure drop. For example, the droplet movement cycle under case 4 is the longest, the pressure drop is also the largest. The periodic appearance, growth, and breakage of the liquid droplet leads to fluctuations of the pressure drop as can be seen in Figure 7. Moreover, it can be observed that while the overall pressure drop behavior in the channel is similar, the periodic increment will increase as the droplet movement cycle increases.

3.2. Effect of the Pore Distance

The effect of the pore distance on liquid water flow in the channel was studied by varying the spacing of the two circular pores. The specific simulation parameters are as follows: R_a = R_b = 0.05 mm, U_G = 1.0 m/s, U_L = 0.1273 m/s (R_a represents the upstream pore radius; R_b represents the downstream pore radius; U_G and U_L represent the inlet rate of gas and liquid water). The contact angles for the GDL surface and the channel wall was set at 135° and 45°, respectively. Cases 8, 9, 10, 11, and 12 correspond to the pore spacing values of 0.3, 0.4, 0.6, 0.8, 1.2 mm, respectively.

3.2.1. Effect of the Pore Distance on Droplet Interactions

The changes in droplet morphology over time are shown in Figure 8. In the initial stage, the two droplets emerge and grow from the two micropores. When the volume of the two droplets grow to be in contact with each other, they will show different morphological changes due to the different pore distance. When the pore distance is 0.3 mm, as shown in case 8 in Figure 8, at 0.004 s, the two droplets merge into a larger droplet (e.g., case 8 (a)), which will remain nearly spherical and continue to grow (e.g., case 8 (b), case 8 (c)) and break from the surface (e.g., case 8 (d)) due to the surface tension.

When the pore distance is 0.4 mm, as shown in case 9 in Figure 8, the two small droplets merge into one large droplet (e.g., case 9 (a)) and rapidly move to the downstream micropore. This newly formed droplet continues to grow and coalesce with the new formed small droplet until it collides with the wall surface and breaks up. In addition, it is found that the combined, larger droplet in different cycles will transfer to the upstream and downstream micropores alternately and become the coalescent droplets that do not stop coalescing with small droplets on the other side. This is because, after the breakup of the coalescent droplet in the previous cycle, the small droplet on the micropore on the side of the coalescent droplet has grown to a certain volume, while the small droplet that has just appeared on the micropore on the side of the coalescent droplet has a small volume. Therefore, after coalescence, the center of gravity of the merged droplet is inclined to the micropore on the side of the coalescent droplet in the last cycle and transferred to the pores. This results in the phenomenon that the coalescent droplet appears alternately on both sides of the micropores.

When the pore distance is 0.6 mm, as shown in case 10 in Figure 8, the upstream droplets deform to the downstream under the shear force of convection, while the downstream droplets deform to the upstream under the action of reverse airflow. At some point, the two droplets will come into contact and merge into one droplet (e.g., case 10 (b)). Under the action of airflow and surface tension, the fusion droplets will leave the micropores, shrink to be nearly spherical, and move to the middle of the two pores. The droplet can then be broken due to collision with the wall of the channel or transfer to a side of the micropores for further growth or fusion (e.g., case 10 (c)). Fusion droplets can block the channel and cause a pressure difference between the upstream and downstream of the droplet. In this case, a countercurrent eddy is formed in the lower reaches to push the fusion droplets towards the upstream micropores (e.g., case 10 (d)). After the fusion droplet is transferred to the upstream, it will continue to fuse with the other side of the micropore droplet where it will separate and grow until broken. The main reason is that the gravity and the projected area at the bottom of the fused droplet are much larger than that of the small droplet in the other micropore. Therefore, the fusion droplets will continue to fuse with the small droplets on the micropores until they are broken.

When the pore distance is 0.8 mm, as shown in case 11 in Figure 8, the pores are widely spaced, and the droplets grow until they break up. The two droplets do not contact each other or coalesce during the entire process. However, since the airflow first contacts the upstream droplet and flows around the droplet, the drag force and deformation caused by the drag force on the upstream droplet is greater than that on the downstream droplet. Thus, the upstream droplet first breaks up and is removed from the GDL surface, while the downstream droplet will be broken when the droplet volume is large enough. As shown in case 11(d) of Figure 8, as the liquid water in the channel increases, slug flow forms in the channel and increases the channel blockage. Therefore, the gas flow has to pass through the gap between the upstream droplet and the slug flow. This exerts a stronger drag force on the downstream droplet, leading to the downstream droplet detaching from the micropore and moving along the GDL surface

When the pore distance continues to increase to 1.2 mm, as shown in case 12 in Figure 8, the two droplets have nearly no interaction between each other, so they grow and detach from the pores separately. This result is due to the large pore distance. Due to the large pore distance, the two droplets will not coalesce before detaching from the pores. However, the upstream droplet will indirectly affect the downstream droplet through the airflow during the growth process. This part will be discussed later.

3.2.2. Effect of the Pore Distance on Gas Velocity Profile, Flow Regime, and Pressure Drop

As can be seen from the above, the degree of interaction between the two droplets changes with an increase in the pore space. It is found that when the pore spacing is less than or equal to 0.6 mm, the droplets will merge and detach from the pores. When the pore distance exceeds 0.6 mm, the droplets will no longer merge, but indirectly interact with each other by affecting the airflow, and the droplets will move downstream from the pores along the GDL surface. Based on this observation, 0.6 mm was considered as the critical pore distance in these conditions—the maximum pore spacing for coalescing between two droplets. The influence of liquid droplets on airflow flow can be seen in Figure 9. Figure 9a is the velocity flow diagram of the x-z section at Y = 0.1 and Figure 9b is the velocity flow diagram of the x-y section at Z = 0.

In Figure 9, case 8 and case 9 are gas velocity flow fields when the two droplets start to merge. As can be seen from this figure, the droplet volume is small at this time, but it can still obstruct the airflow and cause the phenomenon of bypass flow, creating a void space downstream of the merged droplet. Due to the small pore distance at this time, the downstream liquid droplet is just in the void space of the upstream liquid droplet. The reverse eddy current formed at the downstream position of the upstream liquid drop will move to the downstream position of the downstream liquid drop because of the obstruction of the downstream liquid droplet, and there will be an area of no airflow between the two droplets. When the pore distance is 0.6 mm, the gas velocity flow field is shown in case 10 of Figure 9. Different from case 8 and case 9, with an increase in the pore distance, the volume of the two droplets will increase when they are about to merge, but the increase is less than that of the pore distance, so some airflow will pass through the area between the two droplets. As the pore distance continues to increase, coalescence no longer occurs between the two droplets. In this case, a large amount of airflow will pass through the space between the two droplets, as shown in the velocity flow diagram of case 11 and case 12 in Figure 9. Due to the large volume of liquid droplets at this time, the airflow will collide with the wall surface when it flows around the upstream liquid droplets and flows in the opposite direction to the area between the two liquid droplets. After that, it will change its direction again due to the influence of the downstream liquid droplet, resulting in irregular airflow through the space of the two liquid droplets. This change in the direction of the airflow causes the upstream and downstream droplets to collide with different sidewalls and break up, then spread along the wall surface. It can be seen that the droplets do affect each other by affecting the airflow flow even when the pore distance is large.

Figure 10 shows the flow pattern of liquid water in the channel when t = 1.0 s. When the spacing of pores is less than or equal to the critical spacing of 0.6 mm, most liquid water exists at the top of one side of the channel in the form of corner flow. When the spacing between the two pores is greater than 0.6 mm, liquid water tends to be evenly distributed on both sides of the top wall. Part of the upstream liquid water still exists in the form of corner flow on both sides, but the downstream liquid water will form slug flow. This result is because when the pore space is less than 0.6 mm, the droplet will break up in the form of a large droplet.

The pressure drop between the inlet and outlet of the channel is shown in Figure 11. As can be seen in this figure, when the pore distance is less than 0.6 mm, the pressure drop increases periodically with the droplet movement period, but the overall pressure drop is small. When the pore distance is greater than 0.6 mm, the initial pressure drop increases periodically and the overall pressure drop is small, but in the later stage, the pressure drop suddenly increases and fluctuates significantly. The change in the pressure drop is closely related to the change of liquid water in the channel. It is known that the pore distance will affect the droplet motion period, and 0.6 mm is the critical pore distance in this work. Taking the critical hole spacing of 0.6 mm as the demarcation point, the droplet movement period in the left and right sections decreases with the increase of the hole spacing, which makes the pressure drop changes in the two sections the same as the hole spacing. However, when the pore distance is greater than 0.6 mm, the pressure drop increases sharply because of the formation of slug flow downstream of the channel in the later period. This also indicates that slug flow has a much higher degree of obstruction to the channel than corner flow, and increases the pressure drop fluctuations of the channel.

3.3. Implications in PEMFC Applications

The gas diffusion layer is a key component of PEMFCs, providing a needed flow domain for reactant gas and product water. The real GDL is typically carbon paper made of interwoven carbon fibers, which is a porous medium structure with complex and random microstructures. Despite the geometric complexities of the GDL, simulated liquid water entering the gas channel through the GDL micropores in literature is mostly based on the shape of a single hole or a round hole. Therefore, changing the hole shape and hole spacing as done in this work can better approach the disordered fiber stacking situation encountered in actual fuel cell operation. The shape and spacing of the micropores simulated in this work show significant effects on the droplet motion, flow pattern distribution, and coverage of liquid water in the gas channel. The different critical droplet volumes caused by different pore shapes affect the accumulation of liquid water in the channel and the overall pressure drop of the channel, which in turn can affect the performance of PEMFCs. When considering multiple pores, although the existence of slug flow makes it easier for liquid water to move out of the channel, it increases the pressure drop and pressure drop fluctuation of the channel. This can adversely result in increased parasitic energy losses and decreased fuel cell performance. These effects should be considered in fuel cell simulations. The current results also provide a guide for effective water management of the proton exchange membrane fuel cell when considering new GDL materials.

4. Conclusions

CFD simulation coupled with a volume of fluid (VOF) method was employed to investigate liquid water droplet dynamics in flow channels for PEMFCs. The main parameters investigated include the pore shape (under the single water inlet condition) and the pore distance between two pores. Seven pore shapes and five distances between droplets were considered. These parameters strongly affect the droplet cycle, flow regime, and pressure drop, and should be considered in numerical simulations. The main findings from this work are as follows:

The pore shape showed a significant impact on the critical droplet volume and the growth period of the droplet. It was found that the larger the windward side slope and side length was, the larger the critical volume and the longer the growth cycle. It should be noted that the critical droplet volume would be reduced if the windward side of a pore is too long. The results showed that under the same pore area, the pentagon and the hexagon pore structure is more conducive to the discharge of liquid water in the channel and can improve the performance of the cell.

Under the range of conditions investigated in this work, the pore distance influences the air flow pattern, whether droplets coalesce, and the resulting flow regimes. In these conditions, a distance of 0.6 mm between droplets seemed to be the critical distance, beyond which droplet coalescence did not occur. When the micropore distance was smaller than this value, droplets coalesced when the droplets grew. In addition, the large pore distance caused the formation of slug flow in the downstream of the channel, leading to the increase of the overall pressure drop.