Proton Exchange Membrane Fuel Cells (PEMFCs) have gained popularity in recent decades as an alternative energy source for several applications, including providing traction power in FCEVs (fuel cell electric vehicles) [1], [2]. Some of the major benefits of the PEMFC include low operating temperature, high power density, high efficiency, less noise, and lower pollution [3]. Gas Diffusion Layer (GDL) of the PEMFC, which is made up of porous materials such as carbon fibers, facilitates mass transfer and conductance in the CL (catalyst layer) [4]. The compression force exerted on the GDL does not destroy carbon fibers and does not reduce their volume but only reduces the pore space. Moreover, porosity beneath the rib is reduced due to the compression force. Hence, the compression force alters its diffusivity and conductivity [1]. For the cell to achieve optimal performance, it must be able to remove water droplets formed in the CL from the GDL. Compression and deformation of the GDL play a significant role in the ejection of droplets from the GDL and control of water in the GDL [5].

However, it should be taken into consideration that, when the produced water is not removed at a proper rate from the electrode, it will accumulate in the form of a water film that fills up the pores in GDL and leads to flooding, making the transportation of reactants difficult [6], [7]. Several measures have been taken to mitigate the flooding problem by increasing the hydrophobicity of GDL and CL through the introduction of highly hydrophobic material such as Polytetrafluoroethylene (PTFE), insertion of Microporous Layer (MPL) between GDL and CL, and advanced designs of the flow field [3].

One of the famous and efficient flow fields that has been introduced is the interdigitated flow field introduced by Nguyen [8]. In this flow field, the creation of the in-plane pressure gradient in the GDL causes the gas to flow in-plane under the land. This flow helps in the transportation of the reactants' supplies towards the reactive sites, and also the transport of the water through the alternate mechanism from the diffusion to the forced convection in the GDL; consequently, this flow field has solved the flooding problem at the cathode site. As a result, the interdigitated flow field increases the performance of the PEMFC in comparison to the parallel flow field [9–13].

The performance of the interdigitated PEMFCs has been investigated by several studies using different physical parameters. The first parameter under consideration is the influence of the GDL compression. Different studies have been performed in this respect. For instance, Su et al. [14] performed a numerical simulation to investigate the influence of the compression forces on the PEMFC performance. They tested the performance of three different types of GDL and concluded that in the case of non-homogeneous GDL, there would be a marked improvement in the predicted cell performance. By means of numerical simulation, Shin and Wang [15] demonstrated that the influence of the compression forces would become significant in the high current density region. It was found that in the case of ignoring the assembly compression, there would be an overprediction of the fuel cell performance.

Lin et al. [16] investigated the Effect of compaction on two types of GDL fuel cells. The results showed that in the high compression region, the performance and the limiting current of the fuel cell decreased and accordingly, an optimal value for the fuel cell performance was determined. Chipper et al. [17] indicated that the non-uniform compression of the GDL and intrusion of the GDL layer would significantly influence the increase in the degree of non-uniformity in the liquid saturation, the oxygen concentration, and current density distribution. By means of the Lattice Boltzmann Method (LBM), Molaeimanesh and Nazemian [18] proved that the application of compression leads to poorer oxygen distribution and more water vapor accommodation near the catalyst layer surface. Moreover, they obtained the optimal value for the compression ratio for the maximum average current density. Ramiar et al. [19] investigated the effects of the inhomogeneous GDL compression on the performance of PEMFC with an interdigitated flow field. All of the above-mentioned efforts indicate that increasing the compression from 10% to 35% leads to an increase in the flooding of GDL, limited transfer of reactants, and a decrease in the efficiency of the PEMFC.

Shi et al. [20] investigated the Effect of the GDL deformation. It was found that the GDL deformation causes lower oxygen concentration at the reaction sites and an increase in the saturated liquid water level along the flow direction.

Ihonan et al. [21] examined the effects of four commercial GDL materials. They concluded that increasing the clamping pressure leads to flooding of the fuel cell because of the reduction in porosity and the temperature gradient between the GDL and the current collector. In this study, the presence of water is assumed to be a continuous liquid. In contrast, the continuous movement of liquid water through the pores of the GDL structure is not reported as an in situ observation. On the contrary, the environmental scanning electron microscope observations show that there are small water droplets that exist within the fibrous structure of GDL [22], [23].

LBM can be considered a strong numerical method that enables one to simulate pore-scale multiphase flows such as flow in porous media as well as mass transport in the real structures of GDL [24], [25]. There are many studies that have been done concerning the application of the LBM method in the investigation of liquid water and the transfer process of reactants and the structure of GDL [26–31].

Hao and Cheng [31] analyzed the Effect of GDL wettability on liquid water transport in GDL. The results obtained in this study showed that the employment of hydrophilic passages leads to the improvement in the removal of liquid water from the GDL. As a result, it reduces the possibility of flooding in the GDL and CL. Li Chen et al. [32] studied the Effect of land width and GDL contact angle on the liquid water removal time and found that the narrower land reduces the liquid water removal time. Higher contact angle enhances the removal and decreases the residual saturation.

Kim et al. [30] investigated the liquid water transfer taking into account the MPL and GDL. The authors concluded that there was a reduction in the liquid content of PTL when applying a thicker MPL and/or increasing the hydrophobicity of the solid surface. Molaeimanesh et al. [33] analyzed the Effect of GDL compression on liquid water removal. It was revealed that in the case of hydrophobic GDLs, the Effect of increased compression on the removal process does not depend on the compression amount. On the other hand, compression may not always be effective for hydrophilic GDLs.

Deng et al. [34] applied the LBM method to study the Effect of GDL compression on liquid water removal in PEMFC with an interdigitated flow field. It was found that there is an optimal level of compression for the enhanced liquid water removal. Deng et al. [34] examined the Effect of the way of contact angle definition, PTFE distribution, and pressure difference at the interface of gas diffusion and MPL s. It was found that for the sample of GDL that dried under vacuum conditions, the PTFE distribution inside the whole GDL is almost uniform, thus requiring less time to overcome the water barrier.

The primary objective of the current research work is to examine the effect of liquid-water injection location on pore-scale water transport, drainage process, and water saturation in compressed Porous Transport Layers (PTLs) of an interdigitated PEMFCs. In order to accomplish the above-mentioned objective, a pore-scale multiphase MRT-LBM based on the Shan-Chen pseudopotential approach is established for the simulation of gas and liquid transport in the PTLs.

Initially, the validity of the proposed numerical model is proved by comparing the results of the model with Laplace's law, droplet contact angle behavior, and existing correlations of the single-phase flow through porous media, showing that the model can precisely describe the physical mechanism of the process. Then, after validation of the model, the effects of three representative liquid-water injection locations, such as under the gas inlet, under the rib region, and under the gas outlet, on the liquid-water transport, drainage properties, and water saturation time evolution are examined using the validated model.

In this paper, Shan and Chen's multicomponent multiphase LBM with Multi-Relaxation Scheme (MRT) is adopted because of its high stability and accuracy compared to that of the single relaxation scheme. Furthermore, a pseudopotential function is applied to simulate the interparticle interaction in lattice Boltzmann. Considering the component k, the collision operator for multicomponent multiphase flows in a 2D lattice can be written as follows:

fak(x + eaΔt, t + Δt) − f ka(x, t) = −Ω[f ka(x, t) − fak,eq(x, t)] + ΔtSka(x, t), (1) where fak is the Density Distribution Function (DDF), fak, eq is its equilibrium distribution for component k, Ω is the collision matrix, t is the time, x is the spatial position, Δt is the time step, Ska is the forcing term in the velocity space. ea is the discrete velocity along the ath direction; for D2Q9, ea is defined as follows:

0 a = 0, (a − 1)π (a − 1)π (cos [ ] , sin [ ]) a = 1,2,3,4.

ea = 4 4 (2)

(a − 1)π (a − 1)π √2 (cos [ ] , sin[ ]) a = 5,6,7,8.

{ 4 4 The MRT shames offer higher stability, accuracy, low spurious velocity, and achievable density ratio [35–37]. Eq. (1) by the MRT method can be transformed to the following form:

m∗,k = mk − Λ[mk(x, t) − mk,eq(x, t)]| + Δt(I − 0.5Λ)Fkm, (3)

where mk(x, t) and mk,eq(x, t) are vectors of moments and the equilibrium vectors of moments for component k, respectively. M is the transformation matrix. A linear transformation can perform the mapping between velocity and moment spaces.

mk = Mf k, f k = M−1mk. (4)

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

mk = (ρ, e, ε, jx, qx, jy, qy, pxx, pyy)T. (5)

mk,eq = ρ(1, −2 + 3(uk2 + vk2), 1 − 3(uk2 + vk2), uk, −uk, v, −v, u2 − v2, uv)T. (6)

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

Λ = diag(S1, S2, S3, S4, S5, S6, S7, S8, S9)

−1 −1 −1 −1 −1 −1 −1 −1 −1 (7)

= diag(τρ τe τς τj τq τj τq τ τ ), where S1, S4, and S6 should be equal to each other and are considered to be one. S5 and S7 equal 1.1, and S2=S3 and S8 = S9 are related to the bulk and kinematic viscosity, respectively:

λ = ( − 0.5) Cs . (8)

ν = (τ − 0.5)Cs 2, (9)

where Cs = 1⁄√3 is the lattice sound velocity, and τ = S ; also, I is the unit tensor, and (I − 0.5Λ)Fm 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 ratio of Tc = 0.5, the gas-to-liquid dynamic viscosity ratio is about 25, we set relaxation time τg to 1.0 for the gas phase, The relaxation time for the liquid phase is calculated using the ratio of kinematic viscosity.

To implement the external force into the LB framework, Gue's method, proposed by Gue et al. [38] is incorporated. McCracken et al. [37], which is claimed to be derived directly from the Boltzmann equation. In this scheme, the force term in the MRT model is derived as follows:

Fkm = [0,6(uFkx + vFyk), −6(uFkx + vFky), Fkx, −Fkx, Fky, −Fky, 2(uFxk − vFky),2(uFkx

T (10)

− vFky)] .

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

ρk = ∑ fa. (11)

∑k (∑a faea + Δt )

Vk = 2 , (12)

∑k ρk where F = (Fx, Fy) is the total interaction force acting on the component and is expressed as:

Fk = Fkcoh + Fkadh, (13)

where Fkcoh is the fluid-fluid interaction force and is defined in the following:

Fkcoh = Fkk + Fkj, (14)

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

Fkk = −Gkkψk(x)Cs2 ∑ w(|eα|2)ψk(x + eα)eα. (15)

In which ψk(x) is the interaction potential, Gkk is the interaction strength, and w(|eα|2) are the weights only depending on the length of the vector eα. The weights are set as w(1) = 1/3 and w(2) =1/12 for the nearest-neighbor interactions on the D2Q9 lattice [39].

The Fkj is the interaction force between two components and is written as [40]:

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

F21 = −G21ϕ2(x) ∑ w(|eα|2)ϕ1(x + eα)eα, (17)

where ϕ1(x) and ϕ2(x) are different from ψk(x) and are expressed as:

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

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

The values of ρ10,ρ20,a0, G12 and G21 are critical for a multicomponent multiphase system [40–43]; the higher value of which leads to limiting the magnitude of the mutual diffusivity in the gas phase.

−0.0008

In this paper we set a0 = 0.005, ρ10 = log(a ) ,ρ20 = 0.0003 [44], and G12 = G21 = 0.0001.

Fkadh is the fluid-solid adhesion force and is calculated by the following equation:

Fkadh(x) = −Gkadhψk(x) ∑ w(|eα|2)S(x + eαΔt)eα. (20)

In which Gkadh is the adhesion factor that determines the strength of the fluid-solid interaction for component k, S is equal to 1 and 0 for solid and fluid domain nodes, respectively.

Component 1 is water and component 2 is air. Air is considered an ideal fluid; thus, G22 = 0. Water is considered a non-ideal fluid that follows the Carnahan-Starling (C-S) Equation of State (EOS) that is given by [24]:

1 + b ρ⁄4 + (b ρ⁄4)2 − (b ρ⁄4)3 2

PEos = ρRT 3 − aρ . (21)

(1 − b ρ⁄4)

R Tc RTc a

where a = 0.4963 P , b = 0.8727 P , and the critical temperature is Tc = 0.3773 bR. In the present study, we set b=4, R=1, and c=1. Also, we set a=0.5 and Tc ≈ 0.047, the same as Li et al. [45]. The parameter a affects the interface thickness; reducing its value results in a thicker interface, leading to decreased spurious velocity and increased stability in high density ratios. For the water component, the effective mass is given by [46]:

1 2(PEOS − ρcs 2)

ψ (x) = √ 2 . (22)

G11cs G11 is set to -1, as well as for the air component, the interaction potential is equal to its density [47].

To adjust the coexistence densities and surface tension, Li et al. [45] improved the forcing term. Similarly, for the water component, the following relation is proposed for the force term [48]:

12σ|F|2 6(uFx + vFy) + 2 ψ δt(τe − 0.5) 12σ|F|2 −6(uFx + vFy) − 2 ψ δt(τς − 0.5)

F1m = Fx , (23)

−Fx Fy −Fy 2(uFx − vFy) [ (uFy + vFx) ] where |F|2 = Fx 2 + Fy 2 and σ are used to achieve good thermodynamic consistency and to change the mechanical stability condition [48].

In the MRT LB equation, surface tension cannot be tuned independently as it is related to the density ratio. Hou et al. [49] added a source term Eq. (25) to the right-hand side of the MRT collision equation for a single component. In the present model, we added it to the collision MRT equation for the water component as written in the following:

∗,1 1 1 1,eq Λ 1

m = m − Λ[m (x, t) − m (x, t)]| + Δt (I − ) Fm + ΔtC. (24)

The source term in the above equation is expressed as:

1.5τe−1(Qxx + Qyy) −1.5τe−1(Qxx + Qyy)

C = 0 . (25)

−τv−1(Qxx − Qyy) [ −τv−1Qxy ] The variables Qxx ,Qyy, and Qxy are given by:

KG11

Q = 1 2 1 1 . (26)

2ψ (x) ∑ w(|eα| ) (ψ (x + eα) − ψ (x)) eα eα

In which the parameter K can be tuned to surface tension, with its value ranging from 0 to 1 [50]. It is important to note that in Eq. (24), m∗,1 can be less than zero. Ahmed and Sung [51] employed a limiter function A1 with a smaller τe−1 and τς−1 in the relaxation matrix Λ1 in order to return an MRT collision to this point again, so for component 1 the new f1,∗a(x, t) is defined as:

A1: m∗,1(x, t 1 1 1,eq Λ 1,∗ m − Λ[m (x, t) − m (x, t)]| + Δt (I − ) S + ΔtCiff (x, t + Δt) ≥ 0,

2 (27)

+ Δt) { 1 1 1,eq Λ1 m − Λ1[m (x, t) − m (x, t)]| + Δt (I − )S + ΔtC otherwise.

Subsequently, after A1, for higher stability, limiter function A2 is also used as a reserve, as follows:

f1,∗(x, t + Δt)iff ∗(x, t + Δt) ≥ 0,

A2: f1,∗(x, t + Δt) { 1 1,∗ ∗ (28)

max(0, ⁄9 ∑ f (x + eaΔt, t + Δt)) iff (x, t + Δt) < 0.

For relaxation matrix Λ1 we have set τe−1 = τς−1 = 0.02.

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

ρin + ρout ρin − ρout 2(r − r0)

ρ(x, y) = − tanh ( ). (29)

2 2 W

In which r = √(x − xc)2 + (y − yc)2.

The center of the domain is (xc, ) yc, and W is the interface thickness and is set to 2 lattice units. ρin and ρout are the densities inside and outside the droplet, respectively.

The grid size is 200 × 200, and the temperature is set to 0.5Tc(T = 500C), where for water ρin(liquid) = 0.456 and ρout(vapor) = 0.000055, while for the air component, it is set to 5 × 10−7 and 0.000567 for inside and outside of the droplet.

After achieving a steady state, the pressure of the system is calculated as follows [43]:

2 1 2 σ 2 1 2 σ σ̄

Psys = cs ∑ ρσ + 2 c ∑ Gσσ [ψ (x)] + 2 cs ∑ Gσσ̄ ϕ (x) ϕ (x). (30)

σ σ σ≠σ̄

Laplace's law gives the system pressure difference between inside and outside of the droplet, which is linearly related to the inverse of the droplet radius. For two-dimensional simulation, it is written as follows:

(Pin − Pout ) = , (31)

where Pin and Pout are respectively the pressure inside and outside of the droplet, γ is the surface tension, and R is the radius of the droplet.

As is seen in Fig. 1, for T/Tc=0.5 and the initial radius of r0=40, the distribution of density in a 200×200 lattice domain with k= 0.68. The lack of fluctuations in the phase change region and the absence of fluctuations at the boundaries of the interface indicate the stability of the method used and its proximity to the physical one.

Fig. 2 is devoted to the visualization of the results obtained by the simulation at T/T_c =0.5 with different values of k for the confirmation of Laplace's law. It is clear from the figure that there is a linear relationship between the pressure difference and the inverse of the droplet radius. From the results, it is observed that by changing the value of k, the surface tension could be varied independently of the density ratio. For this study, we take k=0.68.

⁄ = . .

Fig. 3 shows the pressure distribution in the x direction of the computational domain for different values of k. It is clear from the figure that the pressure fluctuations near the vapor-liquid interface happen due to a sharp variation of density, which is ignored in the calculation of ΔP.

The contact angle characterizes the wettability of the surface. In fact, if the contact angle is less than 900, then the surface is hydrophilic, and the liquid phase expands like a film on the solid surface. On the other hand, if the contact angle is more than 900, then the surface is non-wetting or hydrophobic, and the liquid phase is formed like a droplet on the solid surface.

In the simulations, the computational domain is 200 × 100. Periodic boundary conditions have been applied on the left and right edges. Bounce-back boundary conditions have been employed on the bottom and top solid surfaces. After convergence of the simulation, the contact angle has been evaluated according to the following expression:

arcsin( ) θ ≤ 90 ,

θ = { (32)

π − arcsin( ) θ > 90 ,

where r is defined as:

4h2 + b2

r = , (33)

where h and b are the height and the base length of the droplet on the surface, respectively. Fig. 4 is about liquid droplets in different contact angles on a solid wall, indicating an almost linear relationship between the adhesion factor and the contact angle.

A standard PEMFC electrode utilizing the interdigitated gas distributor as designed by Nguyen [8] is depicted in Fig. 5. According to Fig. 5a, air enters the multilayer hydrophobic PTLs (GDL and MPL) and is directed through the lands into the outlet channel fingers. In Fig. 5b, the computational domain used within the present study is illustrated. This domain encompasses half the inlet channel, half the outlet channel, one land, and PTFL. The CL is taken to function like an extremely thin interface attached to the bottom of the PTLs. In Fig. 5b, the GDL (region "A") layer is formed by randomly positioning circular solids of average diameter 16 μm in the layer with thickness 150 μm. The porosity of the GDL layer is approximately equal to 0.8 prior to analyzing the PLTs, and its contact angle is taken to be 115-135 degrees. After that, an MPL layer with a thickness of 25-70 μm (region "B") is created through the random positioning of square particles (4 μm edge length). MPL averages 0.76 with respect to porosity, whereas its contact angle equals 150 degrees.

The standard PEMFC electrode having the structure of an interdigitated gas distributor as invented by Nguyen [8] is shown in Fig. 5. From Fig. 5a, it can be seen that air is fed into the multilayer hydrophobic PTLs (GDL and MPL) and is carried via the lands into the outlet channel fingers. In Fig. 5b, the computational domain utilized in the current study is shown. The computational domain includes half the inlet channel, half the outlet channel, one land, and PTFL. It is considered that the CL acts like a very thin interface connected to the bottom of the PTLs.

In Fig. 5b, the GDL (region "A") layer is generated by random placement of circles having an average diameter of 16 μm in a layer having a thickness of 150 μm. The porosity of the GDL layer is about 0.8 before analysis of the PLTs, and its contact angle is 115-135 degrees. Then the MPL layer, having a thickness of 25-70 μm (region "B"), is generated by random placement of square solids having an edge length of 4 μm. MPL has a porosity of about 0.76, and its contact angle is 150 degrees. The values of the parameters used in the computational domain are given in Table 2.

The location of liquid-water injection into the PTLs is a critical parameter governing water management in PEMFCs. The injection location directly affects the liquid-water transport pathway, its accumulation within the porous structure, the spatial distribution of water saturation, and the time required for the system to reach a quasi-steady state.

In general, when liquid water is generated or injected at a location farther from the outlet, it must travel a longer distance through the porous media before being removed. This extended transport pathway promotes liquid-water accumulation within the PTLs, resulting in higher water saturation and an increased risk of flooding.

In the present study, the Effect of the liquid-water injection location on water transport through the PTLs was investigated under identical operating and structural conditions. The compression ratio was set to 85%, the capillary number was Ca = 0.0004, and the MPL thickness was 70 μm. The initial porosities of the GDL and MPL were set to 82% and 79%, respectively, while the inlet air velocity was maintained at 3 m/s. Five water injection points with an average diameter of approximately 5 μm were considered at different locations within the computational domain.

Based on their relative positions with respect to the gas inlet, rib, and outlet, three representative injection configurations were examined: 1) injection beneath the gas inlet region, 2) injection in the middle region beneath the rib, and 3) injection beneath the gas outlet region. The corresponding injection locations are summarized in Table 3.

In this study, a pore-scale multiphase MRT-LBM model was developed to investigate liquid-water transport and accumulation within the PTLs of an interdigitated PEMFC. The model was validated against Laplace's law, droplet contact-angle behavior, and permeability correlations, showing good agreement with analytical and empirical results. The simulations demonstrated that the liquid-water injection location strongly affects water accumulation and drainage within the PTLs. Injection beneath the gas inlet resulted in the highest average water saturation (approximately 0.34–0.35), whereas injection beneath the outlet produced the lowest values (approximately 0.12–0.15). The outlet-side injection configuration reduced the average water saturation by approximately 57–65% compared with the inlet-side configuration and by 42–54% compared with the middle-region configuration.

Overall, positioning the liquid-water injection location closer to the outlet promotes more efficient water removal, reduces water accumulation, and helps maintain open gas transport pathways within the PTLs. The developed pore-scale MRT-LBM framework provides a useful tool for analyzing multiphase transport and optimizing water management in PEMFCs.

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.