Numerical Study of Gas-Liquid Two-Phase Flow in PEMFC Cathode Flow Channel with Different Diffusion Layer Surface Structure

Show more

1. Introduction

Proton exchange membrane fuel cells (PEMFC) have a wide range of promising applications in automotive, aerospace, and portable electronic devices due to their efficient energy conversion efficiency and environmental friendliness [1] [2]. PEMFC is an energy conversion device that converts chemical energy directly into electrical energy, mainly consisting of a bipolar plate, a polymer electrolyte membrane, two catalyst layers, and two gas diffusion layers. In particular, GDL acts as a platform for the transport of multiphase reactants and reaction products [3]. As shown in Figure 1(a), the reaction by-product water needs to be discharged by the cathode runner through the surface of the diffusion layer, resulting in a relatively complex flooding problem in the cathode runner in PEMFC. Therefore, a deeper understanding of liquid water’s transport and removal mechanism in the flow channel is needed.

Due to the small flow channel size, experiments are used to study the kinetic behavior of liquid water in the flow channel, and the experiments require a very high resolution of the imaging device [4]. Consequently, most of the research works use simulation to study the transport process of liquid water within the flow channel. Many fuel cell visualization techniques are available to show the quantitative and qualitative behavior of liquid water, among which the numerical simulation of PEMFC two-phase flow using the volume of fluid function method (VOF) is the most widely used [5]. Numerous factors are affecting the water transport characteristics of the cathode flow channel; the design of the flow channel geometry is crucial, of which the most representative is the three-dimensional fine mesh flow field. Bao [6] used optical microscope images to reconstruct a three-dimensional flow field model and explored single-phase

Figure 1. (a) Schematic diagram of PEMFC; (b) Simplified model of cathode flow channel; (c) Liquid water inlet arrangement on GDL surface.

and two-phase flow characteristics. Zhu [7] investigated the sensitivity of liquid water behavior to the cross-sectional shape of the cathode channel in a low-temperature fuel cell by VOF simulations. Six cross-sectional shapes such as triangle, semicircle, variable-length ratio rectangle, trapezoid, inverted trapezoid, and arch were investigated numerically. The results show that the cathode runner geometry significantly affects the droplet detachment time, GDL surface coverage, and water saturation. In addition to the geometry of the gas channel, influences of the wettability of the channel walls are essential for the two-phase flow pattern. Yue [8] observed the water distribution in fuel cell GDL by x-ray imaging to study the effect of cathode GDL hydrophobicity on fuel cell water distribution. Cao [9] used a VOF model to qualitatively analyze the performance of wall wettability on removing liquid water from different cross-sectional flow channels. Moosal [10] simulated vertical and horizontal fuel cells and investigated the effect of gravity on the gas-liquid two-phase flow in a single serpentine flow field. Although current research has focused on two-phase flow phenomena influenced by macroscopic flow channel structure, flow channel surface properties, and airflow size and direction, little attention has been paid to the effect of GDL surface pore structure distribution on liquid water transport behavior within the flow channel.

In this paper, an air-water two-phase flow model in the PEMFC cathode flow channel is developed using the VOF method. The effect of different sizes of liquid water inlet arrangement on the GDL surface on the water management performance of the PEMFC cathode runner is investigated, as shown in Figure 1(c). The results are qualitatively described by tracking the 3D gas-liquid interface corresponding to each case through the VOF model. The removal characteristics of the different configurations were quantified by post-treatment parameters such as two-phase pressure drop (∆P), the total liquid water content in the flow channel (WV), and liquid water coverage of the GDL surface (GWCR).

2. Mathematical Model

2.1. Computational Domain and Model Assumptions

As the PEMFC undergoes a reduction reaction at the cathode, the appearance of liquid water in the cathode flow channel is preemptive. Therefore, this study only considers the dynamic behavior of liquid water in the cathode flow channel. Wang [11] proposed that the cross-sectional dimension of the flow channel for PEMFC to have the best output performance is 0.535 × 0.535 mm^{2}; Fan [12] research shows that the minimum water inlet distance to ensure that liquid water does not coalesce in the flow channel is 0.6 mm; based on the above research, the specific size of the cathode flow channel in this study is determined to be 0.535 × 0.5353 mm^{3}, of which the length of the cathode flow channel is 3 mm (L = 3 mm). As shown in Figure 1(b), the left side of the flow channel is the air inlet, and the right side is the outlet. The lower surface of the flow channel is the GDL surface, the two sides are the sidewalls of the flow channel, and the top wall is the top wall of the flow channel. The liquid water inlet is set at 0.5 mm (Le = 0.5 mm). As shown in Figure 1(c), the liquid water inlet apertures are respectively φ50 um, φ30 um, and φ20 um, which are sequentially distributed along the centerline of the cathode flow channel along the gas flow direction, and the distance between the holes is 0.6 mm. Finally, the geometric area and boundary of the cathode flow channel together constitute the entire calculation domain of the VOF model. Based on the VOF model, several assumptions are made for the gas-liquid two-phase flow in the cathode flow channel.

1) Gas flow in the flow channel is laminar (Re = 338 < 2300), transient, incompressible gas.

2) Isothermal working conditions.

3) Absence of heat generation and heat transfer in the flow channel.

4) Condensation and evaporation in the cathode runner are ignored.

5) Surface tension and fluid properties remain constant.

2.2. Mathematical Methods

The kinetic behavior of liquid water in the cathode flow channel is studied based on different microstructure distributions on the surface of the diffusion layer, which belongs to the air-water mutually immiscible two-phase interface capture problem. Consequently, a surface tracking technique, the VOF model, is used to model the two-phase flow. In this methodology, air and water share a standard set of momentum equations, and the tracing of the gas-liquid interface in the computational domain is achieved by introducing the variable of phase volume fraction α, where α_{q} represents the ratio of the volume of a phase to the volume of the grid in which it is located. α_{q} = 1, when the grid is filled with the liquid phase; α_{q} = 0, when the grid is entirely in the gas phase; When the grid is a gas-liquid two-phase coexistence, 0 ≤ α_{q} ≤ 1; There are only two phases of air (α_{1})-water (α_{2}) in this research problem, and the interface of liquid water can be traced by solving the volume fraction equation for the liquid phase α_{2}, which is given by the following equation.

$\frac{\partial}{\partial t}\left({\alpha}_{2}\right)+\nabla \cdot \left({\alpha}_{2}{\stackrel{\to}{V}}_{2}\right)=0$ (1)

where
${\stackrel{\to}{V}}_{2}$ is the velocity of liquid water, once α_{2} is known, α_{1} can be found by the following equation.

${\alpha}_{1}+{\alpha}_{2}=1$ (2)

The continuity equation is formulated as follows.

$\frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho \stackrel{\to}{V}\right)={S}_{m}$ (3)

where S_{m} is the mass source phase of the two phases.

The momentum equation is given by:

$\frac{\partial (\rho \stackrel{\to}{V})}{\partial t}+\nabla \cdot \left(\rho \stackrel{\to}{V}\stackrel{\to}{V}\right)=-\nabla P+\nabla \left[\mu \left(\nabla \stackrel{\to}{V}+\nabla {\stackrel{\to}{V}}^{T}\right)\right]+\rho \stackrel{\to}{g}+\stackrel{\to}{{F}_{SV}}$ (4)

where $\stackrel{\to}{{F}_{SV}}$ is the equivalent volume force form, P is the static pressure, $\stackrel{\to}{g}$ and is the acceleration of gravity.

The air-water density ρ and the kinetic viscosity μ are each obtained by the following equation.

$\rho =\left(1-{\alpha}_{2}\right){\rho}_{1}+{\alpha}_{2}{\rho}_{2}$ (5)

$\mu =\left(1-{\alpha}_{2}\right){\mu}_{1}+{\alpha}_{2}{\mu}_{2}$ (6)

where the aerodynamic viscosity μ_{1} is 2.0425 × 10^{−5} Pa∙s, the dynamic viscosity of liquid water μ_{2} is 4.06 × 10^{−5} Pa∙s, the density of air ρ_{1} is 1.29 kg/m^{3}, the density of liquid water ρ_{2} is 1 × 10^{3} kg/m^{3}.

Due to the apparent effect of surface tension on the two-phase flow, the continuous medium surface force (CSF) model [13] was used by adding source phase $\stackrel{\to}{{F}_{SV}}$ to the momentum equation, which was calculated as:

${F}_{SV}=\frac{2\sigma \rho \kappa \nabla {\alpha}_{2}}{{\rho}_{1}+{\rho}_{2}}$ (7)

where σ is the surface tension coefficient, which is taken as σ = 0.0735 N/m in this study; κ is the curvature of the gas-liquid interface.

$\kappa =\nabla \cdot \stackrel{\to}{n}$ (8)

Considering the influence of liquid water by surface adhesion, the normal unit vector of the two-phase interface in the grid is determined according to the wall contact angle, which is calculated as:

$\stackrel{^}{n}={\stackrel{^}{n}}_{w}\mathrm{cos}\theta +{\stackrel{^}{t}}_{w}\mathrm{sin}\theta $ (9)

where n_{w} is the unit surface normal vector, t_{w} is the unit surface tangential vector, and θ is the contact angle.

2.3. Boundary Conditions and Mesh Independence Verification

The boundary conditions of the governing equation are specified at the air inlet, water inlet, outlet, and channel wall in Figure 1(b). According to the working conditions of the fuel cell: the activation area is 30 cm^{2}, the current density is 0.8 A∙cm^{2}, the cathode stoichiometric ratio is 2, and the airflow rate at the inlet is 10 m/s [2] [14]. The inlet boundary condition is specified as a uniform velocity distribution. The research of Jon [15] showed that the hydrophilic contact angle on the top wall and sidewall of the cathode flow channel combined with the hydrophobic angle of the GDL surface, the flow channel has efficient drainage performance; Ding [16] research pointed out that the liquid water injection rate changes Two orders of magnitude, the two-phase flow pattern remains unchanged; Quan and Lai [17] showed that since the air velocity is several orders of magnitude higher than the liquid water production rate on the GDL surface, the liquid water injection velocity will not significantly deviate from the two-phase flow pattern in the cathode flow channel. Based on the above research, this research stipulates that the air inlet velocity is 10 m/s; the water inlet velocity is 1 m/s, to approximate the emergence speed of liquid water on the GDL surface in the cathode flow channel; the boundary condition of the cathode flow channel outlet is outflow; The contact angle of the GDL surface is set to 140, which represents the wettability of the carbon paper surface treated with polytetrafluoroethylene (PTFE) [18]; the contact angle between the sidewall of the cathode flow channel and the top wall is 60; Lu [19] Experimental studies have shown that PEMFC operates normally under the conditions of temperature, operating pressure and gravitational acceleration of 298 k, 1 bar, and 9.81 m/s^{2} respectively. Therefore, these operating conditions are used in this simulation for solution calculation.

This research uses GAMBIT 2.4.6 to divide the computational domain, and the number of grids in the computational domain is 85,368, 96,325, 116,073, 124,122, 131,404. The pressure drop of the two-phase flow in the cathode channel is selected to verify the grid independence of the calculation domain. As shown in Figure 2, the pressure drop of the two-phase flow in the cathode channel is selected to verify the grid independence of the calculation domain. After the number of grids is 116,073, the pressure drop of the two-phase flow basically does not change when the number of grids is increased. Therefore, this study calculates the number of domain grids is 116,073.

2.4. Solving Method

In this research, the powerful CFD commercial software ANSYS 19.0 was employed to simulate the gas-liquid two-phase flow using the pressure solver for the explicit VOF model. A Green-Gauss cell-based algorithm for gradient evolution and a first-order implicit time discretization method is used to capture the liquid-phase interface dynamically. The momentum equation is discretized using the second-order windward format; the pressure-velocity coupling equation is discretized using the SIMPLE algorithm, and the pressure term is discretized using the Body Force Weight. The air-water actual shape is obtained by solving

Figure 2. Grid independence verification.

the volume fraction equation using the Geo-Reconstruct method. The convergence condition of the continuity equation and the momentum equation is specified as 10^{−4}. In order to achieve time accuracy and ensure numerical stability, the Courant number of each grid cell must be less than 1. In this investigation, the time step is set to 10^{−6} s.

3. Results and Discussion

3.1. Effect of Different Sizes of Liquid Water Inlet Arrangement on the Two-Phase Flow Pattern

This section discusses the effect of the distribution model of the liquid water inlet on the surface of GDL with pore diameters of φ50 um, φ30 um, and φ20 um on the flow of gas-liquid two phases, respectively, and to simplify the research problem, let all three be on the centerline of the flow channel. As shown in Figure 2, the liquid water flow pattern in the corresponding cathode runner of each case is gradually discharged in the form of liquid beads. In case 1 and case 2, φ50 um inlet is the closest to the air inlet, and the droplets emerging here are the first to break away under the action of airflow shear force and start to move toward the outlet of the flow channel, and the separated droplets merge with the growing droplets at φ30 um and φ20 um inlets successively. Moreover, the merged droplets were discharged from the flow channel at 5.7 ms and 5.64 ms, respectively. Compared with the two, the liquid water inlet is distributed sequentially along the direction of gas flow in order of size from large to small, which is more conducive to the rapid discharge of liquid water in the flow channel. In case 3, the droplet detached at the φ30 um inlet, and the growing droplet at the φ50 um inlet merge at the 6.9 ms moment. In case 4 the droplets detached at the φ30 um inlet merge with the growing droplets at the φ20 um and φ50 um inlets, and the merged droplets catch up with the detached droplets at the φ50 um inlet at 12 ms. It can be seen that both case 3 and case 4 tend to form fat droplets in the flow channel, which will increase the risk of the cathode flow channel being blocked by liquid water. From case 5 and case 6, it can be observed that the smaller the size of the inlet is, the closer it is to the inlet, and many tiny broken droplets will be trapped in the flow channel, which will cause the accumulation of liquid water on the surface of the GDL and lead to the blockage of the porous cathode electrode, thus causing the PEMFC to work less efficiently.

3.2. Quantification of Water Management Performance of Cathode Runners by Different Sizes of Liquid Water Inlet Arrangement

The total content of liquid water (WV) in the cathode runner at any given moment is shown in Figure 3(a), where WV characterizes the total volume of liquid water resting in the cathode runner. Specifically, it is defined as:

$\text{WV}={\displaystyle \iiint {\alpha}_{2}dV}$ (10)

where α_{2} indicates the volume fraction of liquid water.

Figure 3. Quantification curves of liquid water in the cathode runner for different sizes of water inlet arrangement. (a) Total liquid water content (b) GDL surface coverage.

From Figure 3(a), it can be concluded that each case corresponds to the maximum total liquid water content (Max. WV) in the flow channel in the order of smallest to largest, which are 1.41 × 10^{−2} mm³, 1.4 × 10^{−2} mm^{3}, 1.57 × 10^{−2} mm³, 1.89 × 10^{−2} mm^{3}, 1.35 × 10^{−2} mm^{3}, 1.48 × 10^{−2} mm^{3}. A sudden drop in the total liquid water content in the cathode runner occurs at a particular time interval, corresponding to the discharge of liquid water from the runner one by one in a bead-like manner. In which case 1 and case 2 correspond to the total liquid water content in the flow channel decreasing in a specific time interval with basically the same value, indicating that the droplet volume in the discharge flow channel is uniform, which will be conducive to the stable discharge of liquid water in the cathode flow channel.

The degree of blockage of the reactive gas channels in the GDL by liquid water is characterized by the magnitude of the diffusion layer surface coverage (GWCR), which is specifically defined as:

$\text{GWCR}=\frac{{\displaystyle \iint {\alpha}_{2}dA}}{{A}_{surface}}$ (11)

The smaller the Value of GWCR indicates that the less liquid water covers the surface of GDL, which keeps more reaction gas channels open and can avoid the oxygen starvation phenomenon in PEMFC. As shown in Figure 3(b), the cases are arranged in the order from most negligible to largest, corresponding to the maximum coverage of the GDL surface by liquid water in the cathode runner (Max. GWCR) in the following order: 4.58%, 4.53%, 5.1%, 6%, 5%, and 4.81%. The liquid water coverage of the GDL surface is minimal in case 2, thus revealing that case 2 is more favorable to keep the reaction gas channels open. In the meantime, the maximum pressure drop (Max.∆p) required to discharge the liquid water in the flow channel in each case is 57.6 Pa, 55 Pa, 56.4 Pa, 85 Pa, 55.7 Pa, 51.8 Pa, in descending sequence. It can be noticed that case 4 corresponds to the immense value of pressure drop in the cathode flow channel, which is attributed to the partial droplet discharge along the sidewalls of the flow channel during transport, as shown in Figure 4 (t = 19 ms), where the adsorption of droplets on the hydrophilic walls increases the droplet discharge resistance. Case 6 corresponds to the smallest value of pressure drop in the cathode runner, and case 2 is the second, indicating that the extra pumping power required to

Figure 4. Corresponding flow pattern of each case with different size of inlet arrangement in the cathode runner.

Figure 5. Comprehensive comparison of water management performance of cathode runners for different sizes of water inlet arrangements.

discharge the droplets in the cathode runner is small, which helps to improve the efficiency of the PEMFC.

In summary, the pore structure distribution on the GDL surface can significantly affect the kinetic behavior of liquid water in the cathode flow channel. The maximum pressure drop value, the maximum liquid water content in the runner, and the maximum coverage of the GDL surface by liquid water are used as evaluation indexes, and the smaller the values are, the less likely the cathode runner flooding phenomenon is to occur. As shown in Figure 5, a comparative analysis of each case shows that case 2 has the best water management performance with smaller values for each corresponding cathode runner.

4. Conclusions

The PEMFC operates as an electrical stack under assembly pressure, and the components are subject to different degrees of deformation, especially the GDL pore structure. This paper simulates the gas-liquid two-phase flow phenomenon in the cathode flow channel by the volume of fluid (VOF) method based on different GDL surface microstructures. The preliminary study investigates the effect of different sizes of liquid water inlet arrangement on the GDL surface on the water management performance of the cathode runner. Its main conclusions are as follows.

1) The PEMFC cathode runner has the best water management performance when the liquid water inlets are distributed sequentially along the gas flow direction from largest to smallest in size.

2) The smaller the size of the liquid water inlet closer to the inlet, the greater the number of small liquid beads of different sizes that stay in the cathode runner, resulting in an unstable mode of liquid water transport in the runner.

These findings will guide the construction of the surface structure of the diffusion layer.

References

[1] Zhang, X., Higier, A., Zhang, X. and Liu, H. (2019) Experimental Studies of Effect of Land Width in PEM Fuel Cells with Serpentine Flow Field and Carbon Cloth. Energies, 12. https://doi.org/10.3390/en12030471

[2] Samad, J. (2018) Numerical Investigation of Gas-Liquid Two-Phase Flow inside PEMFC Gas Channels with Rectangular and Trapezoidal Cross Sections. Energy Weekly News.

[3] Mehdi, M. and Kazuya, T. (2014) In-Plane Microstructure of Gas Diffusion Layers with Different Properties for PEFC. Journal of Fuel Cell Science and Technology, 11. https://doi.org/10.1115/1.4025930

[4] Bazylak, A. (2009) Liquid Water Visualization in PEM Fuel Cells: A Review. International Journal of Hydrogen Energy, 34.
https://doi.org/10.1016/j.ijhydene.2009.02.084

[5] Anderson, R., Zhang, L., Ding, Y., Blanco, M., Bi, X. and Wilkinson, D.P. (2009) A Critical Review of Two-Phase Flow in Gas Flow Channels of Proton Exchange Membrane Fuel Cells. Journal of Power Sources, 195.
https://doi.org/10.1016/j.jpowsour.2009.12.123

[6] Zhiming, B., Zhiqiang, N. and Kui, J. (2019) Analysis of Single- and Two-Phase Flow Characteristics of 3-D Fine Mesh Flow Field of Proton Exchange Membrane Fuel Cells. Journal of Power Sources, 438.
https://doi.org/10.1016/j.jpowsour.2019.226995

[7] Zhu, X., Liao, Q., Sui, P.C. and Djilali, N. (2009) Numerical Investigation of Water Droplet Dynamics in a Low-Temperature Fuel Cell Microchannel: Effect of Channel Geometry. Journal of Power Sources, 195.
https://doi.org/10.1016/j.jpowsour.2009.08.021

[8] Like, Y., Shixue, W., Takuto, A., Yoshio, U. and Yulin, W. (2021) Effect of Water Distribution in Gas Diffusion Layer on Proton Exchange Membrane Fuel Cell Performance. International Journal of Hydrogen Energy, 46.
https://doi.org/10.1016/j.ijhydene.2020.06.184

[9] Yan, C., et al. (2021) PEM Fuel Cell Cathode-Side Flow Field Design Optimization Based on Multi-Criteria Analysis of Liquid-Slug Dynamics. Journal of Industrial and Engineering Chemistry, 98. https://doi.org/10.1016/j.jiec.2021.03.024

[10] Ashrafi, M. and Shams, M. (2017) The Effects of Flow-Field Orientation on Water Management in PEM Fuel Cells with Serpentine Channels. Applied Energy, 208. https://doi.org/10.1016/j.apenergy.2017.09.044

[11] Wang, X.-D., Yan, W.-M., Duan, Y.-Y., Weng, F.-B., Jung, G.-B. and Lee, C.-Y. (2009) Numerical Study on Channel Size Effect for Proton Exchange Membrane fuel Cell with Serpentine Flow Field. Energy Conversion and Management, 51.
https://doi.org/10.1016/j.enconman.2009.11.037

[12] Mengying, F., et al. (2021) Effect of Pore Shape and Spacing on Water Droplet Dynamics in Flow Channels of Proton Exchange Membrane Fuel Cells. Energies, 14.
https://doi.org/10.3390/en14051250

[13] B. J.U, K. D.B, and Z. C (1992) A Continuum Method for Modeling Surface Tension. Academic Press, Vol. 100. https://doi.org/10.1016/0021-9991(92)90240-Y

[14] Hossain, M., Islam, S.Z., Colley-Davies, A. and Adom, E. (2013) Water Dynamics inside a Cathode Channel of a Polymer Electrolyte Membrane Fuel Cell. Renewable Energy, 50. https://doi.org/10.1016/j.renene.2012.08.041

[15] Owejan, J.P., Gagliardo, J.J., Sergi, J.M., Kandlikar, S.G. and Trabold, T.A. (2008) Water Management Studies in PEM Fuel Cells, Part I: Fuel Cell Design and in Situ Water Distributions. International Journal of Hydrogen Energy, 34.
https://doi.org/10.1016/j.ijhydene.2008.12.100

[16] Ding, Y., Bi, H.T. and Wilkinson, D.P. (2010) Three-Dimensional Numerical Simulation of Water Droplet Emerging from a Gas Diffusion Layer Surface in Micro-Channels. Journal of Power Sources, 195.
https://doi.org/10.1016/j.jpowsour.2010.05.059

[17] Quan, P. and Lai, M.-C. (2006) Numerical Study of Water Management in the Air Flow Channel of a PEM Fuel Cell Cathode. Journal of Power Sources, 164.
https://doi.org/10.1016/j.jpowsour.2006.09.110

[18] Gopalan, P. and Kandlikar, S.G. (2014) Effect of Channel Materials and Trapezoidal Corner Angles on Emerging Droplet Behavior in Proton Exchange Membrane Fuel Cell Gas Channels. Journal of Power Sources, 248.
https://doi.org/10.1016/j.jpowsour.2013.09.070

[19] Lu, Z., Rath, C., Zhang, G. and Kandlikar, S.G. (2011) Water Management Studies in PEM Fuel Cells, Part IV: Effects of Channel Surface Wettability, Geometry and Orientation on the Two-Phase Flow in Parallel Gas Channels. International Journal of Hydrogen Energy, 36. https://doi.org/10.1016/j.ijhydene.2011.04.226