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

Proton exchange membrane fuel cell (PEMFC) reaction byproduct water is eventually discharged from the cathode channel through the surface of the diffusion layer (GDL). Although there are many studies on the kinetic behavior of liquid water in the cathode flow channel, there is still a lack of litera-ture on water transport behavior in the flow channel considering the pore structure distribution on the GDL surface. In this paper, the effect of different size liquid water inlet arrangements on the GDL surface on the water management performance of the cathode runner is investigated. The volume of fluid (VOF) model is used to simulate the gas-liquid two-phase flow phenomenon in the channel. The results show that the droplet size discharged from the cathode flow channel is basically the uniformity when the liquid water inlet is distributed in gas flow from the largest to the smallest in size. At the same time, the maximum total liquid water content achievable in the flow channel, the maximum GDL surface coverage is the smallest, and the cathode flow channel has the best water management performance. The effect of pore structure distribution on the surface of GDL should be considered when conducting the study of water transport behavior in the cathode flow channel to reveal further the transport process of liquid water in the cathode channel and the removal mechanism.

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 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).

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

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.
is the velocity of liquid water, once α 2 is known, α 1 can be found by the following equation.
The continuity equation is formulated as follows.
where S m is the mass source phase of the two phases.
The momentum equation is given by: where SV F  is the equivalent volume force form, P is the static pressure, g  and is the acceleration of gravity.
The air-water density ρ and the kinetic viscosity μ are each obtained by the following equation.
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 SV F  to the momentum equation, which was calculated as: 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.
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: where n w is the unit surface normal vector, t w is the unit surface tangential vector, and θ is the contact angle.

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 Journal of Power and Energy Engineering 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.

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 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.

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: 2 WV dV α = ∫∫∫ (10) where α 2 indicates the volume fraction of liquid water. Journal of Power and Energy Engineering  (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   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.

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 direc-tion 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.