Discrete-Particle Model to Optimize Operational Conditions of Proton-Exchange Membrane Fuel-Cell Gas Channels

Operation of proton-exchange membrane fuel cells is highly deteriorated by mass transfer loss, which is a result of spatial and temporal interaction between airflow, water flow, channel geometry, and its wettability. Prediction of two-phase flow dynamics in gas channels is essential for the optimization of the design and operating of fuel cells. We propose a mechanistic discrete particle model (DPM) to delineate dynamic water distribution in fuel cell gas channels and optimize the operating conditions. Similar to the experimental observations, the model predicts seven types of flow regimes from isolated, side wall, corner, slug, film, and plug flow droplets for industrial temporal and spatial scales. Consequently, two-phase flow regime maps are proposed. The results suggest that an increase in water accumulation in the channel is related to the increase in the water cluster density emerging from the gas diffusion layer rather than the increased water flow rate through constant water pathways. From a modeling perspective, the DPM replicated well volume-of-fluid channel simulation results in terms of saturation, water coverage ratio, and interface locations with an estimated 5 orders of magnitude increase in calculation speed.


Regime map translations
The presented Figure S2 (a) in this section contain the translation between water and air velocity to equivalent current density and Stoichiometry for the channel conditions specified for Figure 7. The region of current density and S2 Stoichiometry investigated in Figure 8 is also highlighed in red, although water velocity will be different due to different water cluster density used. Figure S2

Extraction of forces from CFD
Droplet detachment in a channel was simulated using the volume of fluid (VoF) method on a three-dimensional domain (625,000 cells, x = 20 µm) shown in Figure 2 (b). Forces acting on the droplet over time was extracted. This information was used as semi-validation, where experimental data for forces acting on discrete droplets is difficult to obtain. In VoF, both phases share the S3 same pressure field, separated by a interface with a specific volume thickness.
To remove the effect of the capillary pressure on the air pressure drop in the channel, the interface was constructed based on a contour 0.01 alpha (volume fraction water). The forces were extracted using Paraview using the following methods 60 .
The inertial force F p,x was found by integrating the pressure on the interface over the area normal to the direction of flow (x direction): The viscous shear force F µ,x was found by calculating the gradient of velocity in the entire domain, where at the interface between water and air the gradient of velocity was multiplied by the viscosity of air to find the shear stress acting on the surface in the two directions (⌧ xz , ⌧ xy ) normal to the direction of flow: The shear stress was integrated over the interface surface area normal to the gradient of velocity in the component direction to find the shear force: The adhesion force F ,x was found by extracting the solid-liquid-air contact line and integrating the normals with respect to the direction of flow, multiplied by S4 the contact angle ✓ with respect to the z direction (channel height direction).

Force comparison between DPM and VoF for capillary bridge
For three different volume water droplets in the capillary bridge regime, Figure   S3 shows the comparison of forces and interface shapes predicted DPM and VoF models. The shear force, inertial force, adhesion force and air pressure drop was extracted using the same method as above at one time step for each case. The magnitude of the forces predicted by each method are of similar values for each force. The three-dimensional representation used for the coalescence algorithm matches the dimensions of capillary bridges produced by the interface resolving VoF method at different wall contact angles (except form the spherical coordinates that is used for the capillary bridge in the DPM model). This analysis, combined with the comparison against transient VoF data in Figure 4 provided confidence in the predictions of momentum exchange between air and water droplets of different geometry.

Derivation of Force Balance
The equations used to estimate the forces presented in Section which act on different droplet geometry were approximated from first principles using the following approaches.

Inertial Force
The obstruction of air flow by droplets in the channel causes acceleration of flow, converting pressure energy into kinetic energy resulting in pressure loss.
Using the conservative form of the Bernoulli equation, the pressure drop was estimated between point 1 (before the droplet) and point 2 (at the droplet cross section) 61 : where p is the pressure, u is the velocity, z is the height coordinate and g is the gravitational acceleration. Equation (S5) is reformulated, ignoring gravity to obtain the pressure drop between points 1 and 2: The velocity directly above the droplet u 2 was determined using the continuity equation for air flow at a constant flow rate where u 1 is equal to the average velocity of air supplied at the inlet to the channel: The inertial pressure in eq. (S6) acts on the interface area of the droplet normal to the direction of flow A c :

Viscous Shear Force
The shear stress ⌧ zy acting on the acting on the air-water interface area A s is found by the viscosity of air multiplied by the gradient in velocity in the y-direction, evaluated at the droplet surface (y = 0): The flow of air in the channel section above the droplet is estimated as planar Poiseuille flow, with the average velocityū: where H is the height of the channel section between the droplet and the top channel surface. The shear stress on the droplet surface is therefore estimated as: The fluid-fluid interfacial area A s determined through analytical approximation is used to calculate the shear force where the average velocity above the droplet isū = u 2 : where b is the approximate thickness of the viscous boundary layer in laminar flow, estimated based based on the droplet configurations.

Derivation of Droplet Regimes
In this section, additional explanation and equations are presented for the regimes schematics presented in Figure 3.

Regime 1 -Emerging
Equation (13) is solved using an Newton iterative method to find the height of the droplet based on a known contact radius R p : The air medial axis position b shown in Figure 3 The area of water that covers the GDL surface, used to determine the water coverage ratio and the liquid-solid shear force was calculated using the radius of the pore: A surf = ⇡R 2 p (S15)

Regime 2 -Isolated
In contrast to regime 1, the contact radius on the GDL surface increases with droplet volume and was calculated using the radius of curvature from eq. (19): The height of the droplet was found using: The area of water covering the GDL A surf and air medial axis location b are calculated using eq. (S15) and (S14) respectively.

Regime 3 -Side Wall
Using trigonometric identities and equations for volume of spherical segments, the dimensions in Figure 3 (c) was estimated. The length from the geometric center to the channel wall x 1 is: The distance from the channel wall to the droplet interface x 2 is: The distance of the geometric center to the GDL surface x 3 is: The distance from the geometric center to the top of the GDL interface on the side wall of the channel x 4 is: The distance of the interface on the GDL surface to the channel wall x 5 is: The angle between the geometric center and the top interface point on the channel side wall ↵ w is: The angle at the geometric center, between the points at the top of the side wall and on the GDL surface ✓ c is: The volume of this shape is determined by truncating the volume of a spherical cap: where the height of the droplet h in this scenario is x 2 shown in Figure 3 (c) and is related to the radius of curvature using: The spherical cap is truncated by the GDL surface and can be approximated by the ratio of areas between a circular segment with an internal angle of ✓ c and a circle with radius R c , this ratio 1 is: Combining eq. (S25) and (S26) produces the droplet volume, expressed as a function of R c : The radius of curvature is found by multiplying eq. (S28) by the area ratio 1 in eq. (S27): The cross section area of the droplet was determined using equations for the circular and triangular cross-sections shown in Figure 3 (c): The surface area of the droplet was determined as the spherical cap surface area multiplied by the area ratio 1 . The area of the GDL covered by water is determined using the dimension x 5 as the height of a circular segment with central angle of 2✓ w . The radius corresponding to this circular segment is found S11 by: Therefore, the area of the GDL that the droplet covers was determined as: The air medial axis was approximated by the difference between x 2 and the

Regime 4 -Corner Droplet
The corner drop is approximated as a spherical cap with contact angle ✓ w , truncated by the channel corner. The dimensions in Figure 3 (d) are as follows.
The angle from the interface at the channel wall to the interface geometric center is: The distance from the geometric center to the channel wall is: The distance from the channel wall to the interface maximum across the channel width is: The volume of the corner droplet was determined by combining eq. (S27) and (S25): The radius of curvature R c was determined as a function of droplet volume and the contact angle of the wall shown. The air medial axis is found using eq. (S33).

Regime 5 -Truncated Capillary Bridge
The dimensions in Figure 3 (e) were used to derive equations for the two radius of curvature. The volume of this shape was found by the addition of two volumes, firstly the volume of a cylindrical shell between R 2 and R 3 , multiplied by the area fraction occupied by water at the cross section W z 1 . Secondly, the volume of a cylinder multiplied by the area ratio at the channel walls. The radius of curvature of the internal curvature R 1 is: The distance between R 3 and R 2 is: Due to the curvature of the interface, the area comprised of liquid within W z 1 is estimated using the area ratio 1 , where ✓ 1 = ⇡ 2✓ w : The next clipping stage involved the truncation of the volume of the cylinder V R 2 by the channel location, which is determined by the wall contact angles.
The secondary area ratio 2 is defined with a central segment angle of circle 2, The equation for the volume of a truncated capillary bridge was determined a function of the three cylindrical volumes: Following each cylinder has a volume V = ⇡R 2 W , eq. (S42) is: Using the schematic in Figure 3 (e) where R 3 = (R 2 z 1 ), the full equation becomes: This quadratic equation was solved using the quadratic formula to determine R 2 : The height of the droplet in the direction of the channel height h at the side walls (R 2 ) and in the middle of the channel (R 3 ) is determined:

Regime 6 -Film
This shape is comprised of a rectangular section shown in Figure 3 (f), with length L i spanning the width of the channel W and two truncated cylindrical sements. The rectangular volume, has an internal curvature determined by the distance between the two side walls as in Figure 3 (e): The internal radius of curvature is found using eq. (S38). The volume of the two end segments V R 2 is found by the circular segment area formed by the contact angle and ⇡ 4 (parallel to the channel walls): V 2 = (hx 1 + R 2 (✓ 2 sin ✓ 2 )) (S48) S15 where x 1 is the distance on the channel wall to the normal of the radius and ✓ 2 = 2✓ w ⇡ 2 . The radius R 2 is found by: Combining eq. (S47) and (S48), the length of the film is estimated as:  Figure S3: Comparison of shear force, inertial force, adhesion force and pressure drop for three different capillary bridge cases predicted by VoF and the DPM model. Model outputs for the interfaces of the VoF (blue) and DPM model (red) at different conditions provided.