Effects of surface wettability on two-phase flow in the compressed gas diffusion layer microstructures

International Journal of Heat and Mass Transfer 151 (2020) 119370

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/hmt

Effects of surface wettability on two-phase flow in the compressed gas diffusion layer microstructures Xia Zhou, Lizhen Wu, Zhiqiang Niu, Zhiming Bao, Xiaoyan Sun, Zhanrui Liu, Yanan Li, Kui Jiao, Zhi Liu∗, Qing Du∗ State Key Laboratory of Engines, Tianjin University, 135 Yaguan Rd, Tianjin 300350, China

a r t i c l e

i n f o

Article history: Received 10 October 2019 Revised 28 December 2019 Accepted 11 January 2020

Keywords: Surface wettability Gas diffusion layer Compression Two-phase flow

a b s t r a c t The surface wettability plays a significant role in two-phase flow in the gas diffusion layer (GDL). However, few reported researches have considered directional single gradient polytetrafluoroethylene (PTFE) distribution, PTFE immersion depth, and the GDL compression. In this study, effects of various surface wettability distribution schemes on two-phase flow in the GDL microstructures are investigated. In addition, the GDL deformation is also considered via clamping pressure simulation based on the finite element method (FEM). Two-phase flow in the compressed GDL microstructures is modeled using the volume of fluid (VOF) model. The results show that the direction of single gradient of the surface wettability is influential in two-phase flow, especially for the compressed GDL. Moreover, thicker PTFE immersion depth contributes to water removal from the GDL. Finally, novel surface wettability distribution schemes are first proposed, and they are demonstrated as important references for controllable water transport in the GDL. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction Proton exchange membrane fuel cell (PEMFC) is a promising candidate as the power source in vehicles [1-2]. As an integral part of the PEMFC, the gas diffusion layer (GDL) plays a key role in the distribution of reactant gases and the removal of product water [3,4]. The accumulation of excess liquid water in the porous GDL may hinder the transport of reactant gases, which in turn worsens the performance of PEMFC [5,6]. Water management is one of the major issues that restricts the commercial development of PEMFC. In order to ensure effective sealing, current collection and reactant delivery, proper clamping pressure should be guaranteed in the process of assembling PEMFC [7]. Clamping pressure can deform the porous GDL [8], which may influence the material transport properties of the GDL such as diffusivity and permeability [9]. The surface wettability of the GDL, which relates to the content and distribution of polytetrafluoroethylene (PTFE), is another significant issue contributing to water transport in the GDL. The process of the PTFE coating on the GDL fibers is always non-uniform, which has a great influence on water drainage in the GDL [10]. Thus, it is meaningful to further explore the effects of various PTFE

∗

Corresponding authors. E-mail addresses: [email protected] (Z. Liu), [email protected] (Q. Du).

https://doi.org/10.1016/j.ijheatmasstransfer.2020.119370 0017-9310/© 2020 Elsevier Ltd. All rights reserved.

distributions on water transport behaviors in the GDL, taking compression of the GDL into account. Numerous experiments [11–13,26,33,34] and numerical methods [12–25,27–32,35–40] have been adopt to investigate the twophase flow in the GDL. For instance, Afra et al. [11] experimentally studied the capillary fingering mechanisms in the GDL by using the images of high-speed camera. They observed that the lateral flow of the catalyst layer (CL) can be reduced by adding micro porous layer (MPL). Satjaritanun et al. [12] adopted direct-modeling-based Lattice Boltzmann Method (LBM) combined with in-situ flow visualization to conduct a fundamental study of water transport in the GDL. Niu et al. [14] utilized volume of fluid (VOF) method to study the effects of perforation on water transport and found that the water content in the GDL can be greatly reduced by perforation. Zhan et al. [15] adopted pore network model (PNM) to investigate the effect of graded porosity on liquid water distribution and oxygen diffusion inside the GDL. With these methods, the effects of surface wettability on water transport in the GDL were investigated under various conditions. Molaeimanesh et al. [23] explored the difference of droplet removal in the cross-flow field of four different PTFE-distributed GDLs. The results showed that the interfacial force between untreated carbon fiber and water droplets would dominate the droplet removal process. In view of the actual hydrophobic treatment of the GDL, surface wettability of the GDL is always inhomo-

2

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370 Table 1 Previous literatures on the effects of PTFE (surface wettability) on water transport in GDLs [17, 23-33, 39]. (EXP: experiment, N: no, Y: yes). Literatures

Items

Dimension

Uneven PTFE

Compression

Method

Ahmad et al. [27] Molaeimanesh et al. [23] Yin et al. [24] Chen et al. [25] Molaeimanesh and Akbari [28] Tötzke et al. [33] Sakaida et al. [31] Niu et al. [39] Shakerinejad et al. [29] Niu et al. [32] Jinuntuya et al. [26] Deng et al. [17] Yu et al. [30] This paper

100°, 120°, 140° 20 (wt% PTFE) 90°, 120°, 150° 10, 18.1,25.6 (wt% PTFE) 80°, 105°, 150° 10, 20 (wt% PTFE) 50°, 130° 90°, 120° 115°, 145° 0, 10, 20 (wt% PTFE) 60°, 80°, 90°,100°, 120°,140° 90°, 115°, 130° 0, 11, 16, 20, 27 (wt% PTFE) 60°–150°

2D 2D 3D 3D 3D 3D 3D 3D 2D 3D 3D 3D 3D 3D

N N Y Y N Y Y N N N N Y Y Y

N N N N N Y N N N N N N N Y

VOF LBM VOF LBM LBM EXP. LBM VOF LBM VOF LBM LBM LBM VOF

geneous. Yin et al. [24] studied the effect of variable contact angles in one GDL on liquid water transport characteristics and found that the water intrusion process in the GDL with variable contact angles was similar to that of the corresponding average fixed contact angle. However, Chen et al. [25] set the GDL sub-region neighboring to gas channel (GC) as the region with high PTFE content and found that the content and distribution of PTFE had an important influence on the liquid-gas flow and the relationship between relative permeability and phase saturation. By changing the contact angle and pressure difference, Jinuntuya et al. [26] found that the transition between stable displacement and capillary fingering was not a gradual process. The hydrophilic contact angle showed the characteristics of stable displacement, while the hydrophobic contact angle showed the characteristics of capillary fingering. Typical researches about effects of surface wettability on water transport inside GDLs are summarized in Table 1. Since the above studies obtained different results of liquid water migration behaviors in the GDL with inhomogeneous surface wettability, it is urgent to further understand their underlying mechanism. Reviewing the mentioned researches, there are still some obvious limitations, such as ignoring the effects of the GDL compression. Considering the GDL compression is an inevitable research trend due its great impact on water transport. Refs [8,33– 36,38,44,45] have investigated the effects of compression on water transport in the GDL. To study the two-phase flow in GDLs with different PTFE contents, Tötzke et al. [33] used the synchrotron Xray computed tomography (CT) to observe the GDL compression, and a dedicated compression device to reproduce the inhomogeneous compression in PEMFC. Zenyuk et al. [34] adopted X-ray CT to study the effects of geometrical land and channel on the distribution of spatial liquid water in the GDL at different compression ratios (CRs). Fazeli et al. [35] used the pore networks extracted from compressed GDLs micro-computed tomography images based on synchronous dispersion. They simulated the liquid water transport of the GDL materials within a certain clamping pressure range. Molaeimanesh et al. [36] adopted LBM to study the effect of the GDL compression on the removal of a water droplet from an electrode of a cell with interdigitated flow field. To the best of the author’s knowledge, few reported studies take both the GDL compression and wettability into account simultaneously. Moreover, the effects of non-uniform PTFE distribution on water transport have not been fully understood. Therefore, it is meaningful to further investigate the effects of various surface wettability distribution schemes on two-phase flow in the compressed GDL microstructures. In this study, a stochastic algorithm [37] is adopted to reconstruct the uncompressed GDL microstructures. The compression process of the GDL is simulated based on the finite element

method (FEM). A three-dimensional VOF model is employed to simulate the two-phase flow in the GDL. The simulated results are verified with experimental study [38]. This study aims to explore the influence of various PTFE distribution schemes on water transport in compressed GDLs. The directional single gradient PTFE distribution, PTFE immersion depth, and the GDL compression are investigated. Moreover, the novel PTFE distribution schemes are first proposed and discussed. This study provides new ideas for the future design of hydrophobic treatment of the GDL for further improvement of water management.

2. Numerical model 2.1. Computational domains and boundary conditions As shown in Fig. 1, the computational domain includes a threedimensional (3D) fibrous GDL, two ribs, and a channel. It can be seen that the non-uniform deformation of the GDL is obvious. A stochastic algorithm [37] is employed to reconstruct the uncompressed GDL microstructures, the whole process is shown in the supplementary materials. More details of the GDL reconstruction have been elaborated in our previous study [46]. In the FEM model, the normal displacement of bipolar plate (BP) and the GDL is limited along the z-axis. The direction of clamping pressure is set along the negative z-axis. The VOF model is adopted to simulate the two-phase flow in the GDL. The fed velocity of liquid water is assumed to be 0.01 m s−1 [27]. All walls are set with a contact angle 90° and the pressure outlets are set to 7 kPa [39]. Details of other basic geometric parameters and boundary conditions can be found in the Fig. 1(a). The relevant material characteristic parameters of the FEM model and the physical parameters of the VOF model are shown in Table 2. The considered surface wettability distribution schemes in this study are as follows: (1) single gradient treatment of surface wettability along the through-plane direction, as shown in Figs. 1(b) (the more hydrophobic GDL on the channel side and the less hydrophobic or even hydrophilic GDL on the CL side) and 1(c) (the less hydrophobic or even hydrophilic GDL on the channel side and the more hydrophobic GDL on the CL side); (2) sandwich treatment of surface wettability along the through-plane direction, as shown in Figs. 1(d) (the 15% thickness on each side of the GDL is more hydrophobic) and 1(e) (the 30% thickness on each side of the GDL is more hydrophobic); (3) New designs of surface wettability distribution of the GDL, as shown in Fig. 1(f) (under-land GDL: less hydrophobic or hydrophilic) and Fig. 1(g) (under-land GDL: more hydrophobic, other part: less hydrophobic or hydrophilic).

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

3

Fig. 1. Schematics of numerical model based on GDL microstructures in the proton exchange membrane fuel cell (PEMFC): (a) The computational domains and boundary conditions; (b) Single gradient scheme 1; (c) Single gradient scheme 2; (d) Through-plane sandwich scheme 1; (e) Through-plane sandwich treatment 2; (f) In-plane sandwich scheme 1; (g) In-plane sandwich scheme 2.

2.2. Governing equations 2.2.1. FEM model In this study, the elastic deformation of the GDL is considered. The related balanced differential equations, geometric equations and physical equations can be used to describe the relationship between strain and stress components, shown as follows [44]:

(1) Balanced differential equations:

⎧ ∂σ ⎪ ⎨ ∂ xx +

∂ τxy ∂ τxz ∂ y + ∂ z + Fx = 0 ∂ σy ∂ τyx ∂ τyz ∂ y + ∂ x + ∂ z + Fy = 0 ⎪ ⎩ ∂ σz + ∂ τzx + ∂ τzy + F = 0 z ∂z ∂x ∂y

(1)

4

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

from 0 to 1 indicates that the cell is at air-water interface. Once α is solved successfully, the average volume density ρ and dynamic viscosity μ for the two-phase mixture in the each cell can be updated according to the α [14]:

Table 2 List of the parameters and their values in the study. Physical parameters Water viscosity (Pa•s) Surface tension (N/m) Young’s modulus (MPa) Poisson’s ratio Mesh parameters Type Mesh size Part 1 VOF FEM

3.0 μm 2.0 μm

0.65 0.07 BP: 13,000, fiber: 300 BP: 0.26, fiber: 0.256

[40] [41] [42,43], Assumed [42,43]

Thickness

Porosity

Mesh number (Million)

200 μm

0.7

1.55 4.25

Compression Ratio (CR)

μ=μl α +μg (1−α )

(7)

(1) Continuity equation:

− →

∇·U =0

Values

Mesh number (Million)

0% 10% 20% 30% 40%

1.98 1.85 1.70 1.55 1.43

− → − → − → U = α U l + (1 − α ) U g

(10)

− → where the “compression velocity” U r indicates the relative velocity of liquid and gas at the interface, the subscript r stands for relative velocity.

(11)

(1) Momentum equation:

−

→ ∂ ρU

(2)

∂t

+

−

−

−

→− → → → ∇ · ρ U U − ∇ · μ∇ U − ∇ U · ∇μ − → − → = −∇ pd − g · x ∇ ρ + σ k∇ α

(12)

where σ and k are the surface tension coefficient and mean inter− → − → face curvature respectively, g denotes gravity vector and x indicates position vector. The modified pressure pd is defined to simplify the definition of boundary conditions:

(3)

→ − → − pd = p − ρ g · x

(13)

Given the effects of surface tension at the gas-liquid interface, the continuum surface force (CSF) model is implemented in this VOF model by adding a force to Eq. (12), which is defined as follow:

fσ = σ k∇α

(1) Physical equations:

(14)

where the mean interface curvature k can be calculated by:

Normal stresses are calculated by:

 σx = λe + 2Gεx σy = λe + 2Gεy σz = λe + 2Gεz

− −

→ → ∂α + ∇ · U α + ∇ · U r α (1 − α ) = 0 ∂t

− → − → − → Ur= Ul− Ug

whereɛx , ɛy andɛz are the normal strains, γ xy ,γ xz andγ yz are the shear strains.

 γxy = τxy /G γxz = τxz /G γyz = τyz /G

(9)

(1) Phase conservation equation:

(1) Geometric equations:

⎧ 2 ∂ 2ε ∂ 2 γxy ∂ ε ⎪ =0 ⎪ 2x + 2y − ⎪ ∂ x∂ y ∂y ∂x ⎪ ⎪ ⎪ 2 2 2 ⎪ ∂ εy ∂ εz ∂ γyz ⎪ ⎪ + − =0 ⎪ ⎪ ∂ y∂ z ∂ z2 ∂ y2 ⎪ 2 2 ⎪ 2 ∂ γ ∂ εz ∂ εx zx ⎪ ⎪ ⎨ ∂ x2+ ∂ z2 − ∂ z∂ x = 0  ∂ ∂ γxy ∂ γxz ∂ γyz ∂ 2 εx + − −2 =0 ⎪ ⎪ ∂x  ∂z ∂y ∂x  ∂ y∂ z ⎪ ⎪ ⎪ ⎪ ∂ 2 εy ∂ ∂ γyx ∂ γyz ∂ γxz ⎪ ⎪ + − −2 =0 ⎪ ⎪ ∂ y ∂ z ∂x ∂y  ∂ x∂ z ⎪ ⎪ ⎪ ⎪ ∂ ∂ γzx ∂ γzy ∂ γxy ∂ 2 εz ⎪ ⎩ + − −2 =0 ∂z ∂y ∂x ∂z ∂ x∂ y

(8)

− → where U represents the effective velocity shared by the two phases throughout the flow domain, which can be calculated by:

whereσ x ,σ y andσ z are the normal stresses, τ xy ,τ xz andτ yz are the shear stresses, Fx , Fy and Fz are the body forces.

− → k = −∇ · n = −∇ · (4)



∇α |∇ α|

(5)

2(1+μ )

where λ is the Lame constant, G is the shear elastic modulus, E is the Young’s modulus, μ is the Poisson ratio. 2.2.2. VOF model In order to capture the interface between two immiscible fluids, the VOF model is adopted to resolve the advection of phase volume fraction α . The value of α is 1 means that the cell is completely filled with liquid phase, and conversely, when it’s filled with gas phase, α is equal to 0. The volume fraction α ranging



(15)

− → where the n is the surface unit normal, which is adjusted in the cells close to the wall according to the following equation:

− → − → − → n = n w cos θ + t w sin θ

e, λ and G can be calculated by:

⎧ ⎨e = εx + εy + εz μE λ= ⎩G = (1+μE )(1−2μ)

(6)

where subscripts l and g represent the liquid and gas phases respectively. The governing equations required for the VOF model are presented as follows [16]:

Part 2: Formal cases Type Objective VOF

ρ =ρl α +ρg (1−α )

(16)

− → − → where n w and t w represent the unit vector normal to the wall and unit vector tangential to the wall respectively. The subscript w denotes the wall. In this paper, θ is the contact angle. 2.3. Model assumptions In this study, assumptions are outlined as follows [24,39,45]: (1) The difference in mechanical behavior between cathode and anode of the PEMFC is ignored; (2) Two-phase flow inside the GDL is laminar and incompressible; (3) No phase change between liquid water and water vapor occurs in the process; (4) The effect of PTFE treatment on volume of pores is ignored; (5) Effect of an MPL can be neglected; (6) Large deformation theory is ignored.

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

5

Fig. 2. Comparisons between the simulated result and experimental results [38,45,46]: (a) quantitative results [45]; (b) qualitative results [46]; (c) mesh independency test [46]; (d) added validation case in this paper (water breakthrough behaviors in the GDL with compression ratio 30%, contact angle: 120 °) [38]; (e) equivalent view.

2.4. Numerical procedures The hexahedral mesh is used to discretize both the computational domains of FEM and VOF models, as shown in Fig. 1(a). The FEM model is used to compress the GDL (solid computational domain: fibers and ribs). Hexahedral meshes are realized with the help of the software called "Hypermesh". Then the compressed GDL microstructures are captured and the fluid computational domain (the channel and pores between fibers) is extracted in reverse from the solid domain. Hexahedral meshes of the fluid computational domain are realized with the help of mesh tools in "Open FOAM". The adopted mesh size and mesh number are shown in Table 2. In order to avoid sliding, the contact surfaces between the components are set bonded due to the existence of bolts. Each case for the FEM model is processed in parallel by using 32 Intel Xeon @3.50 GHz processors. It takes approximately 12 h for each case. The discrete scheme of governing equation of VOF model is second order. The pressure and velocity coupling solutions are solved by a pressure implicit split operator (PISO) scheme coupling the semi-implicit method. Owing to the excellent parallelism of open-source software Open FOAM, the computing domains are decomposed with 140 Intel Xeon @2.93 GHz processors. Each case takes about 36 h. 3. Results and discussions In this paper, water breakthrough time and local water accumulation are the targets which are used to compare the GDL performances between various PTFE distribution schemes. Breakthrough time indicates the time needed for water removal from the GDL, while local water accumulation means the water distribution in the GDL. Faster water removal shows that the GDL is less possible

to suffer from water flooding, whereas that does not mean the GDL is also suitable for reactant gas transport from the channel to the catalyst layer. It is expected that the paths for water (from the catalyst layer to the channel) and gas (from the channel to the catalyst layer) should be preferably separated. Hence, local water accumulation is important to evaluate whether the GDL is suitable for both the water and gas transport. From the authors’ perspective, these two factors are equally impactful on fuel cell performance.

3.1. Model validation Both quantitative (compression-stress curve, porosity curve, and water saturation curve) and qualitative (the morphology of the compressed GDL and dynamic water transport behaviors) comparisons have been conducted in our previous works [45,46], as shown in Fig. 2(a) and (b). The simulated results showed a good agreement with the reported experimental results. Mesh independency tests have also been carried out in Ref. [46] and the results showed that the mesh sizes in Ref. [46] (same with the mesh sizes used in this paper) are accurate enough, as shown in Fig. 2(c). This paper further considers one more case to validate the model. The simulated results (position of water breakthrough, water droplet morphology) show a good agreement with the experimental results reported in Ref. [38], as shown in Fig. 2(d). That means this model is reliable and shows enough accuracy for exploring the water transport in the compressed the GDL. The blue surface shown in Fig. 2(e) is the gas-liquid phase interface, and Fig. 2(e) shows how the blue surface corresponds to water volume in the GDL. The whole manuscript adopts the form of the blue surface to clearly show the water distribution in the GDL.

6

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

Fig. 3. Water dynamic behaviors emerging from the bottom of the GDL (initial porosity: 0.7, fiber diameter: 8 μm, and initial thickness: 200 μm) with different compression ratios (CR): (a) CR = 0%; (b) CR = 10%; (c) CR = 20%; (d) CR = 30%; (e) CR = 40%; and (f) the corresponding slices at y = 10 μm (CPN: compression neglected, CP10: compression ratio 10%).

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

7

Fig. 4. Water dynamic behaviors emerging from the bottom of the GDL (CR = 30%) treated with the single PTFE gradient design 1 (1-3): (a) 150°, 130°, 110°; (b) 130°, 110°, 90°; (c) 110°, 90°, 70°(d) 90°, 70°, 50°

3.2. Effect of compression on water transport in the GDL In this section, the effects of compression on water transport in the GDL are studied. Fig. 3(a) shows water transport behaviors in the uncompressed GDL. There is no significant accumulation of water when it emerges from the GDL bottom to the channel (see t = 1.8 ms). As time goes by, the water chooses a dominant path, and then water breakthrough takes place at the position near the land-channel interface (see t = 9.6 ms). This phenomenon is consistent with the results reported by Jeon et al. [38]. The reason is that there are some preferential paths with lower resistance for water to flow in the porous the GDL microstructures. Fig. 3(b) shows the water transport behaviors in the GDL with the compression ratio 10%. Different from the phenomenon of water breakthrough in Fig. 3(a), the time of water breakthrough is earlier and the position of breakthrough moves away from the landchannel interface (see t = 7.4 ms) in Fig. 3(b). In addition, there is a special phenomenon that less water distribution is found in the GDL under the land-channel interface. This is because com-

pression deforms the GDL microstructure and changes the dominant path. In general, resistance under the land is greater than that under the channel, and compression exacerbates this trend. Water tends to accumulate in the under-channel GDL, so as to select the dominant path for breakthrough. In Fig. 3(b) and (c), it can be observed that there is minor water transport difference in GDLs with the 10% and 20% compression ratios. As shown in Fig. 3(d), there is much less water in the GDL under the land-channel interface when compression ratio is 30%. The water breakthrough occurs in the middle of the flow channel and breaks through in advance (see t = 5.6 ms). Due to the higher compression ratio of 40% in Fig. 3(e), it can be seen that the number of water breakthrough position increases and the time is earlier (see t = 5.0 ms). Fig. 3(f) shows the slice of water transport behaviors after the breakthrough event at y = 10 μm. It is found that once the preferential path is formed when the breakthrough occurs, water will follow this path continuously which keeps a stable water distribution in the GDL. Furthermore, the under-land GDL is drier when compression ratio increases.

8

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

Fig. 5. Water dynamic behaviors emerging from the bottom of the GDL (CR = 30%) treated with the single PTFE gradient design 2 (1-3): (a) 110°, 130°, 150°; (b) 90°, 110°, 130°; (c) 70°, 90°, 110°; (d) 50°, 70°, 90°.

In summary, it can be summarized that compression reduces the time required for water breakthrough. The higher the compression rate, the more obvious this phenomenon. It is worth noting that the microstructures of the GDL is changed by compression, which has a great impact on the mechanism of water transport in the GDL. Therefore, it is necessary to consider the effects compression on studies relating to the two-phase flow in the GDL. In this study, considering the rationality of the breakthrough location of water, the compression ratio 30% is finally selected for subsequent researches. 3.3. Effects of surface wettability distribution on water transport in the compressed GDL In this section, effects of various surface wettability distribution schemes on water transport inside the compressed GDL microstructures are investigated. Directional single gradient surface wettability distribution schemes are discussed in Section 3.3.1. The effects of various PTFE immersion depths on two-phase flow in the

GDL are discussed in Section 3.3.2. Novel surface wettability distribution schemes are proposed and analyzed in Section 3.3.3. All simulations mentioned above have taken the GDL compression into account. 3.3.1. Effects of directional single gradient surface wettability distribution For the PTFE treated GDL, its surface wettability may change along the through-plane direction, resulting in different two-phase flow transport characteristics [24]. As mentioned in the introduction, although recent reported works have investigated the effects of single surface wettability gradient on water transport in the GDL by VOF and LBM methods, they ignored the effects of the gradient direction of PTFE on water transport. Especially for the compressed GDL, the gradient direction is meaningful for water transport owing to non-uniform deformation of the GDL. 3.3.1.1. Hydrophobicity decreases from the GDL top to bottom. In this section, the relative positions among channel, CL and the GDL

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

9

Fig. 6. Water dynamic behaviors emerging from the bottom of the GDL (CR = 30%) treated with the PTFE immersion design 1 (1-3): (a) 90°, 90°, 90°; (b) 120°, 90°, 120°; (c) 140°, 80°, 140°; (d) 150°, 60°, 150°.

are considered. The upper layer is close to the channel side (top side) and the lower layer is close to the CL side (bottom side). Fig. 4 shows the water transport behaviors when the contact angle decreases with a constant gradient of 20° from top to bottom, which is called the single gradient treatment of the GDL. As shown

in Fig. 4(a), it can be observed that early water transport at the bottom of the GDL is relatively flat (see t = 4.6 ms). Two water breakthrough positions (see t = 6.6 ms) are found close to the left interface between the channel and land. When water in the underchannel GDL breaks through, it only reaches the middle layer of

10

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

Fig. 7. Water dynamic behaviors emerging from the bottom of the GDL (CR = 30%) treated with the PTFE immersion design 2 (1-3): (a) 120°, 90°, 120°; (b) 140°, 80°, 140°; (c) 150°, 60°, 150°

the GDL. Compared with Fig. 4(a) and (b) shows that water breakthrough is later and the position number is reduced to one (see t = 8.6 ms). The whole process of water flow moving upward is smoother in Fig. 4(b) and water under the land goes up a little bit during water breakthrough. In Figs. 4(c) and (d), the GDL is also subjected to gradient treatment from top to bottom, but the difference is that its surface wettability tends to be hydrophilic. An interesting phenomenon is found that water transport behaviors in Fig. 4(c) and (b) are

similar. However, there is more water under the land in Fig. 4(c). The hydrophobicity of the GDL decreases and even changes to hydrophilicity, which delays the water breakthrough (see Fig. 4(c) t = 9.0 ms). It can be seen that early water transport in Fig. 4(d) is similar with that in Fig. 4(c), while the phenomena before (see t = 6.0 ms) and after (see t = 9.0 ms) water breakthrough show obvious differences in Fig. 4(c). Water under the land reaches the interface between land and the GDL when water breaks through. Due to the nearly-hydrophilic layer (upper layer 90°) of the GDL,

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

relatively smooth layered water film (see t = 6.0 ms) appears in the GDL and the state of the water after breakthrough looks like two steamed buns. Breakthrough occurs after water accumulates, and the flatter transport process is not conducive to this behavior. The smooth water transport process is probably caused by a decrease in hydrophobicity. The hydrophobicity of the GDL decreases to be less hydrophobic or even hydrophilic, which is not conducive to the water flowing behavior and will trap water. Therefore, under-land water increases because of the weak lateral flow. 3.3.1.2. Hydrophobicity increases from the GDL top to bottom. In Fig. 5, the GDL is treated with the same gradient surface wettability, but the more hydrophilic layer is close to the channel. In Fig. 5(a), water flows from both sides to the middle inside the GDL and begins to converge (see t = 2.0 ms). Eventually, water breakthrough occurs (see t = 3.0 ms). It can be seen that there are some differences between Fig. 5(a) and (b), but the overall trend is similar. There is more water in the GDL under the land. The reason causing the delay of water breakthrough is that 90° is a nearhydrophobic contact angle for more difficult water breakthrough. In Fig. 5(c), some water flows and converges horizontally at the bottom, and other water flows upward in the GDL under the land on both sides at the same time (see t = 3.0 ms). The two flows also converge at the interface of the GDL and channel, but it is difficult to see water breakthrough (see t = 7.0 ms). This is because the hydrophilic GDL easily traps water and impedes its flow. As shown in Fig. 5(d), the water inside the GDL under the channel no longer accumulates from the lower layer (see t = 4.0 ms). Water from left side and right side converges at the middle of interface of the GDL and the channel (see t = 6.4 ms). It’s interesting that the state of water breakthrough is lamellar water film, which hindering the supply of reactant gas to the CL. As expected, two-phase flow in Figs. 4 and 5 shows obvious differences due to opposite gradient PTFE distribution directions. Take Figs. 4(b) and 5(b) for example, uniform water distribution is found in Fig. 4(b), but water is found mainly accumulating in the middle of under-channel GDL and at the margin of the under-land GDL in Fig. 5(b). Moreover, water tends to move horizontally from both the under-land sides to the under-channel region. The more the hydrophobic GDL bottom, the more non-uniform the water distribution. Moreover, water breakthrough occurs at 4.0 ms in Fig. 5(b), whereas the time is 8.6 ms in Fig. 4(b). In summary, the above results mean that directional single gradient of surface wettability highly influences the water distribution and the time of water removal. More hydrophobic GDL near the CL is preferred for better water management. Thus, investigation of the PTFE gradient on water transport in the GDL must consider its relative position with the channel and the CL, especially for the condition that the GDL compression is considered. 3.3.2. Effects of various PTFE immersion depths In the hydrophobic treatment process of the GDL, the GDL is always put into PTFE solution for full immersion, and then dried and sintered. It is easy to understand that the internal hydrophobicity of the GDL is always smaller than the external hydrophobicity [24]. However, the immersion depths are seldom considered in the previous studies. In this section, two-phase flow in the GDL with different immersion degree (30% and 60% total thickness) of PTFE solution is studied. 3.3.2.1. More hydrophobic treatment of 30% total thickness. The hydrophobicity of the upper, middle and lower layers of the GDL in Fig. 6(a) remains consistent. Water at the bottom of the GDL moves upward synchronously in a relatively smooth manner. As shown in Fig. 6(b), there is a large dry area, but some water still exists in the middle layer of the GDL. Eventually, water breaks through

11

at both the left and right sides of the channel (see t = 9.6 ms). Fig. 6(c) shows a more exaggerated water stratification than that shown in Fig. 6(b). It can be observed that water tends to move upward inside the GDL under the land firstly, and then moves horizontally due to changes in surface wettability when it reaches the interface between the middle layer and the upper layer (see t = 3.0 ms). Eventually, the two water flows meet in the middle (see t = 6.0 ms), and move up for breakthrough (see t = 9.2 ms). The first two small breakthrough droplets form a large droplet (see t = 10.2 ms). In theory, water flows more easily through the hydrophobic GDL and is trapped in the hydrophilic GDL. However, the lower layer is so thin that water does not have time to move horizontally, which prevents the water from accumulating in the GDL under the channel. The middle layer of the GDL is hydrophilic, and the compression under land may cause the water to tend to gather here. The difference between Figs. 6(c) and 10(d) is only reflected in the time of water breakthrough (see Fig. 6(d) t = 8.6 ms). The reason is that the larger contact angle difference may amplify the effects which is similar to those shown in Fig. 6(c), and faster convergence leads to earlier water breakthrough. In this section, two-phase flow in the sandwich treatment the GDL is further understood. Obvious water stratification occurs due to the different contact angles between the middle layer and the adjacent two layers. Therefore, water tends to move upward from the left and right side, and then moves horizontally for converging. The larger contact angle difference will amplify the effects. 3.3.2.2. More hydrophobic treatment of 60% total thickness. The hydrophobic treatment of the upper and lower layers is further increased to 30% of the total thickness separately depending on the degree of immersion. The water stratification of the middle layer in Fig. 7(a) is still obvious, but the intersection position of water is shifted to the left (see t = 4.0 ms) when compared with that in Fig. 7(b). The number of water breakthrough position is also reduced to one. Although water in the hydrophilic region does not tend to move towards the center, the hydrophobic region reflects a stronger role due to the decreased thickness of the middle layer. In addition, the distribution of porosity is non-uniform under the compression, which changes the preferential path of water transport and characteristics of water breakthrough. In Fig. 7(b), water converges earlier compared with that in Fig. 6(c) and the intersection position of water is also shifted to the left (see t = 4.0 ms). Compared with Fig. 6(d), Fig. 7(c) shows that water converges at the interface of the middle layer and the upper layer earlier. It is obvious that the water breakthrough in Fig. 7 has been advanced in all the cases and the difference between direct upward transport of water and upward transport from left and right side becomes insignificant, when compared with that in Fig. 6. The reason is that a larger proportion of the hydrophilic layer makes the location of water convergence move upward (the interface between the middle layer and the upper layer). In Fig. 6, the thickness of the middle layer and the lower layer of the GDL together accounts for 85%, while it is 70% in Fig. 7. In Section 3.3.2, it is easy to conclude that water stratification effect occurs when various PTFE immersion depths are considered. Water breakthrough occurs earlier and more proportion of the under-land GDL keeps dry when deeper PTFE immersion is adopted. Increasing the contact angle difference between the middle layer and the two adjacent layers will not magnify this effect. Thicker hydrophobic layer contributes to water removal from the GDL. 3.3.3. Novel PTFE distribution schemes In this section, new surface wettability distribution schemes are designed. The surface wettability of the left, middle and right parts of the GDL is treated by sandwich distribution, as shown in Figs. 8

12

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

Fig. 8. Water dynamic behaviors emerging from the bottom of the GDL (CR = 30%) treated with the novel PTFE design 1 (1-3): (a) 90°, 120°, 90°; (b) 80°, 140°, 80°; (c) 60°, 150°, 60°.

and 9. The difference of contact angle between the under-land GDL and the under-channel GDL increases gradually, which can help us fully understand two-phase flow in the GDL with the new designed surface wettability distribution schemes. 3.3.3.1. The more hydrophobic under-channel and less hydrophobic under-land design. In this section, more hydrophobic under-

channel and less hydrophobic under-land GDL is designed. In Fig. 8(a), water in the left and right parts of the GDL firstly moves upward, and then breaks through at the intersection of left and right channel-land the GDL (see t = 6.2 ms). Water under the channel builds up later and begins to flow upward (see t = 7.0 ms). It is an interesting question why water does not firstly flow into

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

13

Fig. 9. Water dynamic behaviors emerging from the bottom of the GDL (CR = 30%) treated with the novel PTFE design 2 (1-3): (a) 120°, 90°, 120°; (b) 140°, 80°, 140°; (c)150°, 60°, 150°.

14

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

the more hydrophobic part under the channel. Water is more likely to be trapped in the hydrophilic GDL compared with the hydrophobic GDL. When the GDL is compressed, the porosity of the GDL under land is reduced compared with that under the channel, and water is more inclined to accumulate in the middle part. However, when the left and right side are treated as hydrophilic, water does not move into the middle part. That means the effect of hydrophilic treatment on water movement here is more obvious than that of porosity reduction. Little horizontal water transport is observed at t = 6.2 ms. Water in the middle part does not accumulate, which makes it difficult for water to break through in the middle part. As shown in Figs. 8(b) and (c), the increased contact angle difference (60°, 90°) between the middle part and the side parts makes the water breakthrough earlier (see Fig. 8(b) t = 5.8 ms and Fig. 8(c) t = 3.4 ms), and the phenomenon of water accumulation in the middle part is less obvious (see Fig. 8(b) t = 7.0 ms). The reason is that larger contact angle difference contributes to the trend shown in Fig. 8(a). Compared with Fig. 3(d), Fig. 8(a) shows that although water breaks through relatively late, the GDL under the channel is relatively dry, which is beneficial for reactant gas arriving the catalytic layer (CL) for electrochemical reaction. Larger adjacent contact angle difference in Fig. 8(c) makes the time of water breakthrough shorter than that in Fig. 3(d). In this section, the novel schemes that more hydrophobic under-channel GDL and less hydrophobic (or even hydrophilic) under-land GDL is designed. Water is found mainly existing in the under-land GDL. Dry under-channel GDL is first realized. This novel design is demonstrated beneficial for better water distribution, which can supply reactant gas to the CL effectively. 3.3.3.2. The less hydrophobic under-channel and more hydrophobic under-land design. Revert to the new design in Section 3.3.3.1, the middle part of the GDL is designed to be less hydrophobic or hydrophilic and the left and right parts are designed to be more hydrophobic in this section. As shown in Fig. 9(a), water in the middle part firstly flows upward and breaks through at the interface of the channel and the GDL (see t = 7.0 ms). It can be inspired by the cases in Fig. 8 that water itself prefers to stay in the less hydrophobic or hydrophilic region. In addition, the porosity under the channel is not significantly affected by compression and porosity of the GDL under the land decreases due to compression. Thus, the porosity and contact angle in Fig. 9(a) are more favorable for water transport in the middle part. Water looks like steamed buns scattered on both sides (see t = 8.4 ms) after breakthrough. This is because the nearly-hydrophilic contact angle and the hydrophilic contact angle are not conducive for water breakthrough. In Fig. 9(b), water breakthrough occurs in the right half of the interface between the channel and the GDL (see t = 9.0 ms). However, the state of the water after the breakthrough is water film rather than droplets. In addition, there are no other obvious differences between Figs. 9(a) and (b). In Fig. 9(c), the phenomenon of water transport is not significantly different from the previous cases. However, water breaks through in the form of a smooth water film, which covers the entire channel (see t = 6.6 ms). The last two cases in Fig. 9 can better illustrate that water tends to flow into hydrophilic GDL, but it is difficult for it to flow out of this region. What’s more, the greater the contact angle difference between the middle part and the adjacent part, the stronger the effect. Compared with Fig. 8, Fig. 9 shows that the hydrophilicity of the middle part in the GDL is not only unfavorable for water breakthrough, but it also may cover the channel and the GDL interface, causing difficulties in the supply of reactant gas. In a word, in order to avoid poor water management, this kind of the GDL design is not expected to appear in the future design of the GDL.

Through the cases in Figs. 8 and 9, it can be known that the new design proposed in this study can realize the macro-control of water transport inside the GDL. By special PTFE treatment on the GDL, the control of water flow path in the GDL can be realized. This is of great guiding significance to optimize the water drainage performance of the GDL in the future.

3.4. In-plane water saturation In this section, in-plane water saturation distribution in the GDLs are considered to further investigate the effects of various PTFE distribution schemes on two-phase flow. The results are shown in Fig. 10. Fig. 10(a) shows the in-plane (along x-axis) water saturation distribution in GDLs at various compression states. It is easy to find that water saturation in the left (negative x-axis side) under-land GDL is much higher than that in the under-channel and right (positive x-axis side) when compression ratio is lower than 20%. That may indicate the pore size of the right under-land GDL is much smaller than that in the left under-land GDL, which makes water prefer to flow into the left side. As the compression ratio increases, water saturation shows an obvious increase and is significantly higher than that in both under-land GDL. Furthermore, in-plane water saturation keeps a stable state when water breakthrough occurs, which is same with the phenomenon found in Fig. 3(f). In Fig. 10(b), four typical schemes (GS1-GS4) are selected for comparison to discuss the effects of directional single PTFE gradient on in-plane water saturation distribution. The results indicate that water saturation in the under-channel GDL is much lower than that in the under-land GDL when gradient scheme 2 (hydrophobicity increases from the GDL top to bottom) is used. As far as the gradient scheme 1 (hydrophobicity decreases from the GDL top to bottom) is concerned, water saturation in the under-channel GDL is much higher than that in the under-land GDL when the whole GDL is hydrophobic, whereas they are almost similar when the whole GDL is hydrophilic. Fig. 10(c) shows the corresponding in-plane water saturation. It is easy to find that under-channel water saturations are truly lower than that under the land. In addition, water breakthrough time is earlier when the contact angle difference increases. It is interesting to find that water saturation in the right under-land GDL is higher than that in both under-channel and left underland GDL, which is completely different from the phenomenon observed in Fig. 10(a). The PTFE immersion design (S2-S4) shows lower under-channel water saturation than that in the uniform design (S1). Fig. 10(d) shows the in-plane water saturation in the GDL when 30% PTFE immersion thickness is considered. The results show that under-channel water saturation is similar when different contact angles are used. For the S2 design, water saturation in both under-land GDLs is almost same. As the angle difference between layers increases, right under-land water saturation increases whereas left under-land water saturation decreases. Furthermore, time of water breakthrough is earlier when larger angle difference is used. Fig. 10(e) shows the in-plane water saturation distribution in the GDL treated with the above two novel PTFE schemes. It is found that the scheme S3 (60°,150°,60°) has the earliest water breakthrough. Different from that found in previous results, underchannel water saturation does not show a stable distribution when the scheme NS1 (more hydrophobic under-channel GDL design) is investigated. This may be caused by water accumulation near the land-channel interface. In a word, the above discussions show that the direction of single gradient of the surface wettability is influential in dynamic water transport behaviors and water distribution, especially for the compressed GDL. Furthermore, thicker PTFE immersion depth contributes to water removal from the GDL. Finally, novel PTFE distribution schemes are demonstrated as important references for controllable water transport in the GDL.

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

15

Fig. 10. In-plane water saturation distribution in the GDL with different compression ratios and PTFE distribution schemes.

In this study, the heat and electron transfer processes are not considered and relevant assumptions have been clarified in Section 2.3. If heat and electron transfer is considered, water accumulation in the under-land GDL, from the authors’ prospective, may contribute to heat and electron transfer between the GDL and the CL, as thermal conductivity and conductivity of water are much higher than that of gas in the pores between fibers and gaps between the GDL and ribs. Thus, the first optimal design (the more hydrophobic under-channel the GDL and less hydrophobic under-

land GDL) may be expected to have more potential to improve the cell performance when heat and electron transfer is considered. In addition, as far as some PTFE distribution schemes (Figs. 5–8) are concerned, most water flows in the under-land GDL (or at early stage), which means water condensation may not influence the results here. As for the other PTFE distribution scheme (Fig. 9), water mainly flows in the under-channel GDL. Thanks to strong in-plane water transport ability caused by wettability difference (from the authors’ aspect), the under-land water may be

16

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370

expected to flow into the under-channel GDL. That means, from the authors’ prospective, the results in this study may be directly applicable to fuel cell conditions. However, because the effect of evaporation and condensation is influential in two-phase flow in some GDLs, it should be considered in the future study. 4. Conclusion In this study, the GDL microstructures are reconstructed using a stochastic algorithm. The deformation of the GDL microstructures is also considered based on the finite element method (FEM). Twophase flow in the compressed the GDL is investigated by a threedimensional (3D) volume of fluid (VOF) model. According to the actual non-uniformity of PTFE treatment (along the through-plane direction) in the GDL, the directional gradient treatment schemes and PTFE immersion depth are considered. In addition, new PTFE treatment schemes (along the in-plane direction) are proposed to further improve water management inside the GDL. The main conclusions are summarized as follows: (1) Directional single gradient of surface wettability highly influences the water distribution and the time of water removal. The more hydrophobic GDL near the CL is preferred for better water management. (2) Water breakthrough occurs earlier and more proportion of the under-land GDL keeps dry when deeper PTFE immersion is adopted. However, excessive PTFE immersion depth may destroy the water stratification effects. (3) The new surface wettability scheme that more hydrophobic under-channel GDL and less hydrophilic under-land GDL is preferred for effective reactant gas transport. The novel PTFE treatment design provides a vital reference for the macro-control of water transport inside the GDL. Declaration of Competing Interest None. CRediT authorship contribution statement Xia Zhou: Conceptualization, Methodology, Software, Data curation, Writing - review & editing. Lizhen Wu: Writing - original draft, Investigation, Data curation. Zhiqiang Niu: Methodology, Software. Zhiming Bao: Methodology, Software. Xiaoyan Sun: Data curation, Investigation, Validation. Zhanrui Liu: Data curation, Investigation, Software. Yanan Li: Data curation, Investigation. Kui Jiao: Supervision, Funding acquisition, Resources. Zhi Liu: Supervision, Project administration. Qing Du: Supervision, Funding acquisition, Resources. Acknowledgments This work is supported by the National Natural Science Foundation of China (Grant No. 51976138), and the National Natural Science Foundation of Tianjin(China) for Distinguished Young Scholars (Grant No. 18JCJQJC46700) Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ijheatmasstransfer. 2020.119370. References [1] G. Zhang, L. Fan, J. Sun, K. Jiao, A 3D model of PEMFC considering detailed multiphase flow and anisotropic transport properties, Int. J. Heat Mass Transfer 115 (2017) 714–724.

[2] S. Wang, Y. Wang, Investigation of the through-plane effective oxygen diffusivity in the porous media of PEM fuel cells: Effects of the pore size distribution and water saturation distribution, Int. J. Heat Mass Transfer 98 (2016) 541–549. [3] L. Hao, P. Cheng, Capillary pressures in carbon paper gas diffusion layers having hydrophilic and hydrophobic pores, Int. J. Heat Mass Transfer 55 (1-3) (2012) 133–139. [4] Y. Wang, S. Cho, R. Thiedmann, V. Schmidt, W. Lehnert, X.H. Feng, Stochastic modeling and direct simulation of the diffusion media for polymer electrolyte fuel cells, Int. J. Heat Mass Transfer 53 (5-6) (2010) 1128–1138. [5] W.Z. Fang, Y.Q. Tang, L. Chen, Q.J. Kang, W.Q. Tao, Influences of the perforation on effective transport properties of gas diffusion layers, Int. J. Heat Mass Transfer 126 (2018) 243–255. [6] N.W. Lee, S.I. Kim, Y.S. Kim, S.H. Kim, B.K. Ahn, M.S. Kim, An effective discharge method for condensed water inside the GDL using pressure gradient of a PEM fuel cell, Int. J. Heat Mass Transfer 85 (2015) 703–710. [7] J. Millichamp, T.J. Mason, T.P. Neville, N. Rajalakshmi, R. Jervis, P.R. Shearing, D.J.L. Brett, Mechanisms and effects of mechanical compression and dimensional change in polymer electrolyte fuel cells–A review, J. Power Sources 284 (2015) 305–320. [8] D. Kanda, H. Watanabe, K. Okazaki, Effect of local stress concentration near the rib edge on water and electron transport phenomena in polymer electrolyte fuel cell, Int. J. Heat Mass Transfer 67 (2013) 659–665. [9] G.R. Molaeimanesh, M. Nazemian, Investigation of GDL compression effects on the performance of a PEM fuel cell cathode by lattice Boltzmann method, J. Power Sources 359 (2017) 494–506. [10] T. Chen, S. Liu, J. Zhang, M. Tang, Study on the characteristics of GDL with different PTFE content and its effect on the performance of PEMFC, Int. J. Heat Mass Transfer 128 (2019) 1168–1174. [11] M. Afra, M. Nazari, M.H. Kayhani, M. Sharifpur, J.P. Meyer, 3D experimental visualization of water flooding in proton exchange membrane fuel cells, Energy 175 (2019) 967–977. [12] P. Satjaritanun, J.W. Weidner, S. Hirano, Z. Lu, Y. Khunatorn, S. Ogawa, S.E. Litster, A.D. Shum, I.V. Zenyuk, S. Shimpalee, Micro-scale analysis of liquid water breakthrough inside gas diffusion layer for PEMFC using X-ray computed tomography and Lattice Boltzmann Method, J. Electrochem. Soc. 164 (11) (2017) E3359–E3371. [13] M. Sabharwal, J.T. Gostick, M. Secanell, Virtual liquid water intrusion in fuel cell gas diffusion media, J. Electrochem. Soc. 165 (7) (2018) F553–F563. [14] Z. Niu, Y. Wang, K. Jiao, J. Wu, Two-Phase flow dynamics in the gas diffusion layer of proton exchange membrane fuel cells: volume of fluid modeling and comparison with experiment, J. Electrochem. Soc. 165 (9) (2018) F613–F620. [15] N. Zhan, W. Wu, S. Wang, Pore network modeling of liquid water and oxygen transport through the porosity-graded bilayer gas diffusion layer of polymer electrolyte membrane fuel cells, Electrochim. Acta 306 (2019) 264–276. [16] Z. Niu, J. Wu, Z. Bao, Y. Wang, Y. Yin, K. Jiao, Two-phase flow and oxygen transport in the perforated gas diffusion layer of proton exchange membrane fuel cell, Int. J. Heat Mass Transfer 139 (2019) 58–68. [17] H. Deng, Y. Hou, K. Jiao, Lattice Boltzmann simulation of liquid water transport inside and at interface of gas diffusion and micro-porous layers of PEM fuel cells, Int. J. Heat Mass Transfer 140 (2019) 1074–1090. [18] Y. Wang, S. Wang, S. Liu, H. Li, K. Zhu, Three-dimensional Simulation of a PEM Fuel Cell with Experimentally Measured Through-plane Gas Effective Diffusivity Considering Knudsen Diffusion and the Liquid Water Effect in Porous Electrodes, Electrochim. Acta 318 (20) (2019) 770–782. [19] G. Molaeimanesh, M.H. Akbari, Water droplet dynamic behavior during removal from a proton exchange membrane fuel cell gas diffusion layer by Lattice-Boltzmann method, Korean J. Chem. Eng. 31 (4) (2014) 598–610. [20] F.C. Cetinbas, R.K. Ahluwalia, A.D. Shum, I.V. Zenyuk, Direct Simulations of Pore-Scale Water Transport through Diffusion Media, J. Electrochem. Soc. 166 (7) (2019) F30 01–F30 08. [21] P. Carrere, M. Prat, Liquid water in cathode gas diffusion layers of PEM fuel cells: Identification of various pore filling regimes from pore network simulations, Int. J. Heat Mass Transfer 129 (2019) 1043–1056. [22] P.A. García-Salaberri, G. Hwang, M. Vera, A.Z. Weber, J.T. Gostick, Effective diffusivity in partially-saturated carbon-fiber gas diffusion layers: Effect of through-plane saturation distribution, Int. J. Heat Mass Transfer 86 (2015) 319–333. [23] G.R. Molaeimanesh, M.H. Akbari, Impact of PTFE distribution on the removal of liquid water from a PEMFC electrode by lattice Boltzmann method, Int. J. Hydrogen Energy 39 (16) (2014) 8401–8409. [24] Y. Yin, T. Wu, P. He, Q. Du, K. Jiao, Numerical simulation of two-phase cross flow in microstructure of gas diffusion layer with variable contact angle, Int. J. Hydrogen Energy 39 (28) (2014) 15772–15785. [25] W. Chen, F. Jiang, Impact of PTFE content and distribution on liquid–gas flow in PEMFC carbon paper gas distribution layer: 3D lattice Boltzmann simulations, Int. J. Hydrogen Energy 41 (20) (2016) 8550–8562. [26] F. Jinuntuya, M. Whiteley, R. Chen, A. Fly, The effects of gas diffusion layers structure on water transportation using X-ray computed tomography based Lattice Boltzmann method, J. Power Sources 378 (2018) 53–65. [27] Z.Y. Ahmad, S. Didari, J. Moon, T.A.L. Harris, Computational fluid dynamics of water droplet formation and detachment from gas diffusion layer, ECS Trans. 45 (23) (2013) 89–100. [28] G.R. Molaeimanesh, M.H. Akbari, Role of wettability and water droplet size during water removal from a PEMFC GDL by lattice Boltzmann method, Int. J. Hydrogen Energy 41 (33) (2016) 14872–14884.

X. Zhou, L. Wu and Z. Niu et al. / International Journal of Heat and Mass Transfer 151 (2020) 119370 [29] E. Shakerinejad, M.H. Kayhani, M. Nazari, M. Nazari, A. Tamayol, Increasing the performance of gas diffusion layer by insertion of small hydrophilic layer in proton-exchange membrane fuel cells, Int. J. Hydrogen Energy 43 (4) (2018) 2410–2428. [30] J. Yu, D. Froning, U. Reimer, W Lehnert, Polytetrafluorethylene effects on liquid water flowing through the gas diffusion layer of polymer electrolyte membrane fuel cells, J. Power Sources 438 (2019) 226975. [31] S. Sakaida, Y. Tabe, T. Chikahisa, K. Tanaka, M. Konno, Study on PEFC Gas Diffusion Layer with Designed Wettability Pattern Tolerant to Flooding, ECS Trans. 86 (13) (2018) 111–118. [32] Z. Niu, Z. Bao, J. Wu, Y. Wang, K. Jiao, Two-phase flow in the mixed-wettability gas diffusion layer of proton exchange membrane fuel cells, Appl. Energy 232 (2018) 443–450. [33] C. Tötzke, G. Gaiselmann, M. Osenberg, T. Arlt, H. Markötter, A. Hilger, A. Kupsch, B.R. Müller, V. Schmidt, W. Lehnert, I. Manke, Influence of hydrophobic treatment on the structure of compressed gas diffusion layers, J. Power Sources 324 (2016) 625–636. [34] I.V. Zenyuk, D.Y. Parkinson, G. Hwang, A.Z. Weber, Probing water distribution in compressed fuel-cell gas-diffusion layers using X-ray computed tomography, Electrochem. Commun. 53 (2015) 24–28. [35] M. Fazeli, J. Hinebaugh, Z. Fishman, C. Tötzke, W. Lehnert, I. Manke, A. Bazylak, Pore network modeling to explore the effects of compression on multiphase transport in polymer electrolyte membrane fuel cell gas diffusion layers, J. Power Sources 335 (2016) 162–171. [36] G.R. Molaeimanesh, M.H. Shojaeefard, M.R. Moqaddari, Effects of electrode compression on the water droplet removal from proton exchange membrane fuel cells, Korean J. Chem. Eng. 36 (1) (2019) 136–145.

17

[37] L. Chen, H.B. Luan, Y.L. He, W.Q. Tao, Pore-scale flow and mass transport in gas diffusion layer of proton exchange membrane fuel cell with interdigitated flow fields, Int. J. Therm. Sci. 51 (2012) 132–144. [38] D.H. Jeon, H. Kim, Effect of compression on water transport in gas diffusion layer of polymer electrolyte membrane fuel cell using lattice Boltzmann method, J. Power Sources 294 (2015) 393–405. [39] Z. Niu, K. Jiao, Y. Wang, Q. Du, Y. Yin, Numerical simulation of two-phase cross flow in the gas diffusion layer microstructure of proton exchange membrane fuel cells, Int. J. Energy Res. 42 (2) (2018) 802–816. [40] L. Korson, W. Drost-Hansen, F.J. Millero, Viscosity of water at various temperatures, J. Phys. Chem. 73 (1) (1969) 34–39. [41] N.B. Vargaftik, B.N. Volkov, L.D. Voljak, International tables of the surface tension of water, J. Phys. Chem. Ref. Data 12 (3) (1983) 817–820. [42] C. Carral, P. Mele, A numerical analysis of PEMFC stack assembly through a 3D finite element model, Int. J. Hydrogen Energy 39 (9) (2014) 4516–4530. [43] T.J. Mason, J. Millichamp, T.P. Neville, A. El-kharouf, G.P. Bruno, J.L.B. Daniel, Effect of clamping pressure on ohmic resistance and compression of gas diffusion layers for polymer electrolyte fuel cells, J. Power Sources 219 (2012) 52–59. [44] M. Movahedi, A. Ramiar, A.A. Ranjber, 3D numerical investigation of clamping pressure effect on the performance of proton exchange membrane fuel cell with interdigitated flow field, Energy 142 (2018) 617–632. [45] X. Zhou, Z. Niu, Z. Bao, J. Wang, Z. Liu, Y. Yin, Q. Du, K. Jiao, Two-phase flow in compressed gas diffusion layer: Finite element and volume of fluid modeling, J. Power Sources 437 (2019) 226933. [46] X. Zhou, Z. Niu, Y.N. Li., X.Y. Sun, Q. Du, J. Xuan, K. Jiao, Investigation of two-phase flow in the compressed gas diffusion layer microstructures, Int. J. Hydrogen Energy 44 (2019) 26498–26516.