Nowadays, the study of Polymer Electrolyte Membrane Fuel Cell (PEMFC) as a power source in vehicles has attracted much attention from researchers. Reasons for this tendency include high energy conservation efficiency, low levels of harmful environmental pollutants, low operating temperatures, and relatively rapid onset. Due to the formation of liquid water in the catalyst layer on the cathode side of PEMFC, it is important to remove it from the gas channel. The accumulation of liquid water causes flooding in the gas channel and blocks the oxygen transport pathway, thereby reducing fuel cell efficiency.

Numerous experimental and numerical investigations have been carried out on the removal of liquid water and flooding in the gas channel of the fuel cell. Using neutron radiography measurements, Kramer et al. [1] investigated the two-phase flow phenomena in hydrogen-fueled Polymer Electrolyte Fuel Cells (PEFC) by successfully detecting liquid accumulation in the flow field and Gas Diffusion Layer (GDL) under various operating conditions. Hickner et al. [2] studied the effect of current density and temperature on the water content of an operating Proton Exchange Membrane Fuel Cell (PEMFC). Using neutron radiography, they found that the peak in water content of the fuel cell was observed at 60°C cell temperature.

Kumbur et al. [3] investigated theoretically and experimentally the effect of operational air flow rate on liquid droplet deformation at the interface of the Diffusion Media (DM) and the gas flow channel. At high flow rates, removal of liquid water from the DM surface occurs faster, which prevents local channel flooding. Yu et al. [4] experimentally simulated a cathode flow channel system with the high-sensitivity double parallel conductance probes inspecting system. Hussaini et al. [5] studied the effect of cathode flooding in the operation of a fuel cell. They introduced a new parameter, the wetting area ratio, to describe the channel flooding and liquid water coverage on a GDL. Zhang et al. [6] presented a Liquid Water Removal from a Polymer Electrolyte Fuel Cell. Their experimental results show that there are two major modes of liquid water removal from the GDL surface. One is the droplet detachment by shear forces of gas flow that leads to a mist flow in the gas channel, and the other is the capillary wicking onto the more hydrophilic channel, which results in the growth of the annular film flow and liquid slug flow in the channel.

Nilsson and Rothstein [7] present the effect of contact angle hysteresis on the dynamics of the coalescence of sessile drops. They show that droplet deformation and oscillation during the coalescence process are also more intense than in the single droplet flow. Considerable research has also been done on the numerical simulation of liquid water transport behaviors in PEM fuel cell gas channels that provide more details on the dynamics of water droplet movement in the fuel cell. The behavior of the liquid water inside a serpentine channel for a Proton Exchange Membrane (PEM) fuel cell was investigated by Quan et al. [8]. The numerical method used was the Volume of Fluid (VOF). Their results indicated that the bending region of the serpentine flow channel has a significant influence on the flow field and water transport, and that flooding may occur after the bend section and might block the transport of reactants inside the flow channel.

Similarly, liquid water transport in the micro-parallel-channels with PEM fuel cell stack inlet and outlet c was studied by Jiao et al. [9]. They found that excessive and unevenly distributed water in the PEMFC gas channel caused airflow to be blocked or unevenly distributed. Zho et al. [10] investigated numerically, through a VOF, the dynamical behavior of a water droplet that enters the GDL gas channel. They found that channel height plays an important role in the deformation and detachment of a water droplet. The use of the Lattice Boltzmann method as a powerful technique, especially for multiphase flow, has attracted many researchers in the last decade [11–14].

Ben Salah et al. [15] simulated with a multiphase free-energy Lattice Boltzmann model the behavior of a water droplet in a gas channel of a PEMFC. Their results show that, with increasing droplet velocity, the drainage efficiency in shallow channels is better than in deep channels. However, their efficiency decreases dramatically when the droplet hits the corner or above the channel's walls. Ben Salah et al. [16] employed the free energy LB model with a high density ratio to find the optimum GC geometry for water drainage. The results indicated that in a rectangular channel at a reasonable pressure drop, it has the best possible characteristics for water removal from the gas channel.

Hao and Cheng [17] used the Lattice Boltzmann method to study the emergence of a water droplet through a hydrophobic GDL in PEMFC. They realized that the droplet motion on the GDL surface was influenced by the drag shear force of the gas flow. Their results indicate that when the velocity of gas flow or contact angle increases, water droplet removal in facilities also increases.

Han and Meng [18] employed the LB method to simulate the removal of a water droplet from two different types of serpentine gas channels. The first one is a smooth U-shaped, and the second one has a sharp right-angle turning region. Their results indicate that the first type of serpentine gas channel is more suitable for water removal. Increasing the gas flow velocity or the surface contact angle of the channel can help the process.

Ben Amara and Ben Nasrallah [19] simulated with a single-component LB method the behavior of droplet motion in the gas channel of PEMFC. They observed that as the capillary number increased, the droplet deformation also increased. They also found that the average contact angle of the hysteresis increases along the droplet motion. Wu and Huang [20] investigated the dynamic behavior of liquid droplets on the GDL surface. They have realized that decreasing the viscosity ratio or increasing the Capillary number, the droplet detaches from the GDL surface easily.

One of the multiphase flow models in Lattice Boltzmann methods is the Shan-Chen pseudopotential method. The main advantages of the pseudopotential method are that the phases separate naturally via the effect of the lattice interaction force, realized by the density-based pseudopotential. Due to the fact that in the flow channel the capillary number is higher, the effect of high density and viscosity ratio is significant [21], and therefore it should be considered in modeling. Based on a single-component model with a large density ratio, Li et al. [14], [22] applied the pseudopotential Shan-Chen model and improved forcing scheme, suggesting the Multi Relaxation Time (MRT) LB model in order to achieve thermodynamic consistency and large density ratio in the model, Wu et al. [23]. optimized the method used by Li et al. [14], [22] to achieve a higher density ratio and higher Reynolds number.

In addition, a few studies have been conducted on multicomponent systems with high density ratios [24], [25]. Bao and Schaefer [26] introduced a new multicomponent multiphase LBE model that could be simulated with a high-density ratio. Zhu et al. [27] incorporated the Equation of State (EOS) into the pseudopotential model. They improved it by using the MRT model for approximating the collision operator and the correction coefficient of K.

Hou et al. [28] simulated the dynamic behavior of droplet coalescence in the gas channel of PEMFC using a multicomponent multiphase LMB method with a high-density ratio. Their results showed that the shear force strengthened during the droplet coalescence and helped in the droplet movement from the gas channel. In this paper, a multicomponent multiphase pseudopotential Shan-Chen LB model with MRT collision operator is developed to study the dynamic behavior of the coalescence process of liquid droplets in the gas channel of PEMFC. Also, the parameters such as surface contact angle, pressure differences, radius of droplet, and distance between two droplets are investigated.

In this paper, Shan and Chen's multicomponent multiphase lattice Boltzmann model with Multi-Relaxation Schemes (MRT) is applied because it has higher stability and accuracy than the single relaxation method. In the lattice Boltzmann, a pseudopotential function is used to consider inter-particle interactions. For component k, the collision operator for multicomponent multiphase flows in a 2D lattice can be expressed as

fak(x + eaΔt, t + Δt) − f ka(x, t) = −Ω[f ka(x, t) − fak,eq(x, t)] + ΔtSka(x, t), (1)

where is the Density Distribution Function (DDF), is its equilibrium distribution for Component k, and is the forcing term is the time step, is the collision matrix, t is the time, x is the spatial position, in the velocity space. is the discrete velocity along the a th direction, for the D2Q9, is defined as follows: In this paper, an MRT collision operator is used. The MRT schemes offer a higher stability, accuracy, low spurious velocity, and achievable density ratio [11], [14], [29–31]. Eq. (1) by the MRT method can be transformed to the following form:

where and are vectors of moments and the equilibrium vectors of moments for Component k, respectively. M is a transformation matrix. A linear transformation can perform the mapping between velocity and moment spaces.

For the D2Q9 lattice, the equilibria mk and mk, eq can be written as

is the relaxation matrix and is diagonal in the momentum space:

where the s1, s4, and s6 should be equal to each other and are considered to be one. s2=s3 and s8 = s9 are related to the bulk and kinematic viscosity, respectively:

where is the lattice sound velocity. s5=s7 and equal1.1. I is the unit tensor, which is the force term in momentum space. Since in the MRT method the relaxation time is related to the kinematic

viscosity, and at a temperature of the gas-to-liquid dynamic viscosity ratio of about 12, we set the

relaxation time to 1.2 for the gas phase, and from the ratio of kinematic viscosity, calculate the relaxation time for the liquid phase.

In this paper, to implement the external force into the LB framework, the Gue's method proposed by Gue et al. [32], [33] and McCracken et al. [31] was used, which is claimed to be derived directly from the Boltzmann equation. In this scheme, the force term in the MRT model should be

In this model, the density of the component and velocity can be expressed as follows [33]:

where is the total interaction force, acting on the component and expressed as:

where is fluid- fluid interaction force and defined:

where is the intra-molecular interaction force, and can be written as:

where the interaction potential, Gkk, is the interaction strength, and the weights are only

dependent on the length of the vector . It is given by w(1) = 1/3 and w(2) =1/12 for the nearest-neighbor interactions on the D2Q9 lattice [34]. The interaction force between two components is shown as [35]:

F12 = −G12�1(x) ∑ w(|eα|2)�2(x + eα)eα. (16)

where and are different from and expressed as:

�1(ρ2) = 1 − exp ( ). (18)

�2(ρ1) = a0 − exp ( ). (19)

The value s of, , and are critical for multicomponent multiphase systems [25], [35], [36]. The higher values of them lead to controlling the magnitude of the mutual diffusivity in the gas phase.

In this paper, we set the fluid-solid adhesion

force and calculated it by:

where the adhesion factor is determined, and the strength of the fluid–solid interaction for component k. S is equal to 1 or 0 for a solid and fluid domain node, respectively.

In this paper, component 1 is water and component 2 is air. The air is considered an ideal fluid and thus G22=0. The water is considered a non-ideal fluid following Carnahan-Starling (C-S) EOS, which is given by [37].

R Tc RTc

where a = 0.4963 P , b = 0.8727 P . The parameter a affects the interface thickness; reducing its value causes a thicker interface, which leads to a decrease in the spurious velocity and increases stability in high-density ratios. In the present study, we set b=4, R=1, a=0.5, and c=1.

For the water component, the effective mass is given by [38]:

In this paper, G11 is set to -1 [25]. Also, for the air component, the interaction potential is equal to its density. In order to adjust the coexistence densities and surface tension, Li et al. [14] improved the forcing term. Similarly, for the water component, the following relation was proposed for the force term [28].

In the MRT LB equation, surface tension cannot be tuned independently and is related to the density ratio. In order to solve this problem, Li and Laue [11] added a source term for a single component. In the present model, we adopted it for the water component as follows:

where the source term is expressed as:

The variable DNA is given by:

where the parameter K can be tuned to a surface tension with its value ranging from 0 to 1. In Eq. (27) , which may be less than zero, Wu et al. [23] employed a limiter function A1 with a smaller relaxation matrix in order to return an MRT collision to this point again, so for component 1, the new value is defined.

Subsequently, after A1for higher stability, limiter function A2 is also used as a reserve.

In this paper, for the relaxation matrix

In this section, for evaluating Laplace's law of a model, a circular droplet with a radiu s initially placed at the center of the domain is considered. Periodic boundary conditions are imposed in the x and y directions. The initial density fields are given as follows:

W here is the center of the domain, and W is the interface thickness, set to 2 lattice unit s. Moreover, the densities of the inside and outside of the droplet are, respectively.

The grid size is , and the temperature is set to , where for water and while for air component are set to 0.00055 and 0.00003 for inside and outside of the droplet. The pressure of the droplet is calculated with [25].

Laplace's law shows that the pressure difference between inside and outside of a droplet is linearly related to the inverse of the droplet radius, and is written as follows:

where Pin and Pout are respectively the fluid pressure inside and outside of the droplet , is the surface tension, and r is the radius of the droplet. Fig. 1 shows the distribution of density for Tc = 0.5 initial radius r0=40 in the lattice domain and k= 0.2. The smoothness in the phase change region and the absence of a minimum or maximum value at the interface boundaries indicate that the method used is stable and close to the physical state.

Fig. 2 illustrates the result of the simulation at T = 0.5 with different values of K for the validation of Laplace's law. There is a linear relationship between pressure difference and the inverse of droplet radius. The results demonstrate that by changing the value of k, the surface tension can be adjusted independently of the density

Fig. 3 shows the pressure distribution along the x direction of the computational domain for different values of k. It can be found that large fluctuations of the pressure near the vapor- liquid interface are also due to a sharp change of density and are ignored when calculating ΔP.

The contact angle determines the wettability of the surface. In fact, if the contact angle is less than 900, the surface is hydrophilic, and the liquid phase tends to expand as a film on the solid surface. In contrast, if the contact angle is greater than 900, the surface is non-wetting or hydrophobic, and the liquid phase forms as a droplet on the solid surface. In these simulations, the computational domain is 200 × 100. Periodic boundary conditions are applied on the left and right edges. A bounce-back boundary condition is applied on the bottom and top solid surfaces. After the simulation is converged, the contact angle is calculated using the expression:

where r is calculated by

where h and b are the height and base length of the droplet on the surface, respectively. Fig. 4 demonstrates a liquid droplet in different contact angles with a solid wall.

θ = 150 θ = 120 θ = 90 θ = 70

The computational domain is a PEMFC gas channel with 2000μm (length) ×200μm (Width). The bounce-back boundary condition is used for the top and bottom walls. For gas stream flow at the inlet and outlet of a rectangular channel, constant pressure introduced by Zou and He [39] is employed. For the water component, the Zou and He [39] constant pressure is used at the inlet, and a Convective Boundary Condition (CBC) [40] at the outlet of the channel. For this simulation, the computational domain is discretized with a 1000×200 lattice unit. The time step is set to , and we set k=0.125. In order to

calculate the inlet and outlet density for the air component, first, the pressure difference between the inlet and outlet of the gas channel must be calculated in the Lattice Boltzmann unit, so that the dimensionless parameter of the skin friction coefficient is equal to the real unit and the LBM unit.

ΔP ΔP

( ) = ( ). (34)

1 2 1 2 2 ρairu 2 ρairu

real

Considering the relationship between ideal gas and sound velocity, the relation between pressure and density can be written as follows:

2 ρLBM

PLBM = Cs ρLBM = . (35) Thus, by calculating the inlet and outlet of air density and considering the ideal gas relation for the air component, their relation is as follows:

ρin,LBM = ρout,LBM + 3ΔPLBM. (36)

In this way, by increasing the inlet pressure and density of the channel relative to the output, it will move the fluid forward. In this simulation, the behavior of motion and the coalescence process of two water droplets placed at the bottom of the solid wall of the PEM gas channel is discussed. Effect parameters such as pressure difference, surface contact angle, radius of droplets, and distance between two droplets are investigated. In this model, the effect of roughness is neglected, and a smooth surface is assumed [10]. Also, the gravity is neglected. Liquid water generated in the catalyst layer is passed through the porous GDL and forms as a droplet in the gas channel, which must be removed. Draining of droplets depends on multiple structural and operational parameters such as droplet size, wall contact angle, etc. Various pressure differences influence the dynamic motion of the droplets. In Fig. 5, the effect of the pressure difference on the behavior of the dynamic droplets' motion has been studied. It is known that increasing the pressure difference between the inlet and the outlet of the gas channel increases the gas flow velocity (Fig. 5), which increases the shear force, and eventually the droplets' movement is faster.

It is also evident that the effect of shear force on the leading droplet is greater than that of the trailing droplet, so its movement is faster. Eventually, the coalescence process between the two droplets happens. As shown in Fig. 5, the droplets merging is faster with increasing pressure difference but occurs farther from the upstream of the gas channel. It is also seen in Fig. 6 that increasing the pressure difference not only boosts the droplet motion but also causes its deformation.

Fig. 8 shows the effect of wall contact angle on the droplet movement in the gas channel. The larger contact angle increases the repulsive force between the water and the wall, and the separation of the droplet is facilitate.d The hydrophilic wall plays a role in resisting the force between the wall and the droplet, reducing the amount of shear force by decreasing its height; as a result, the droplet motion retard. In the hydrophobic wall, due to the repulsive force between the wall and the droplet and more shear force, the droplet motion is faster, and the coalescence process between two droplets occurs farther from the upstream of the gas channel. The size of the droplet formed in the gas channel is different, depending on the structural and operational parameters of the fuel cell. In Fig. 8, the effect of the droplet size on its movement along the channel is shown. Note that the radius of the leading droplet is assumed to be greater than or equal to the radius of the trailing droplet, but the sum of the volumes of the two droplets is equal. It is also evident that the shape of the first droplet before the coalescence process is similar in all three cases. When the size of the first droplet before the coalescence process is bigger, the smaller the distance between the droplet and the top wall, resulting in higher acceleration and maximum velocity, which makes the droplet move faster.

Due to the complex structure of the GDL, the droplets detaching are formed randomly at different locations, so the distance between them is different. In this section, the effect of the distance between two droplets is studied in Fig. 9, so that the distance between the two droplets is considered 300,400 and 500µm, while in all three cases the pressure difference is the same.

The results show that when the distance between the two droplets is 300 µm, the drag shear force on the trailing droplet is reduced, so the leading droplet hits the trailing droplet faster, and the coalescence process is performed. After that, the radius of the droplet becomes bigger as the gap between the droplet and the top wall of the gas channel decreases, so the fluid accelerates, and the velocity of the air component increases; thus, the drag shear force enhances, and as a result, the droplet movement is faster. In general, it can be found that the coalescence process is beneficial for the existing droplet in the channel. Fig. 9 illustrates the distribution of fluid velocity before and after the coalescence process. When the distance between the two droplets is 500 µm, the trailing droplet is more affected by the drag shear force of the air flow than in the other two cases, so its movement is faster, while the drag shear force applied to the trailing droplet is the same in all three cases. Eventually, the coalescence process does not occur, and both droplets drain out of the channel separately and at a higher time than the other two cases.

A numerical investigation on the interaction and coalescence process of two water droplets in the gas channel of PEMFC was presented using a multicomponent multiphase pseudopotential Lattice Boltzmann method with MRT collision operator.

In this way, two droplets with different distances between them were placed inside the gas channel, and their movement and coalescence process were studied under different parameters such as contact angle, pressure difference, radius of droplets, and distance between them. The important findings and conclusions are summarized in the following.

I. Increasing the pressure difference causes the shear force to increase; as a result, the droplet movement and merging from the gas channel are faster.

II. Our simulations demonstrate the effects of the GDL contact angle on the removal and coalescence process of the droplet, where a more hydrophobic GDL results in a higher repulsive force between the wall and the droplet, and the motion of the droplet is faster.

III. The effect of the distance between two droplets has been investigated, as decreasing it reduces the shear drag force applied to the trailing droplet. The leading droplet hits the trailing droplet faster and merges with it, after the coalescence process, because the shear stress is applied to the droplet, and the droplet velocity becomes faster, and the droplet movement is faster.

The authors declare that there are no conflicts of interest regarding the publication of this paper.

The datasets generated and/or analyzed during the current study are included in this article.

This study was carried out without any specific funding from public, commercial, or non-profit organizations.