Investigation of two-phase flow in the compressed gas diffusion layer microstructures

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/he

Investigation of two-phase flow in the compressed gas diffusion layer microstructures Xia Zhou a, Zhiqiang Niu a, Yanan Li a, Xiaoyan Sun a, Qing Du a, Jin Xuan b,**, Kui Jiao a,* a b

State Key Laboratory of Engines, Tianjin University, 135 Yaguan Rd, Tianjin 300350, China Department of Chemical Engineering, Loughborough University, Loughborough, United Kingdom

highlights
Three-dimensional microstructures of GDLs are digitally reconstructed.
Two-phase flow in uncompressed and compressed GDLs is considered.
Effects of multiple parameters on two-phase flow in GDLs are investigated.
Mechanisms of water transport in GDLs are discussed.

article info

abstract

Article history:

This study aims to investigate how multiple parameters affect the two-phase flow in

Received 27 April 2019

compressed gas diffusion layer (GDL). A stochastic model is adopted to reconstruct the GDL

Received in revised form

microstructures. Solid mechanics simulations on the reconstructed GDL microstructures

29 July 2019

are performed, based on the finite element method (FEM). Various pore morphologies and

Accepted 12 August 2019

distributions of compressed GDLs are observed. Two-phase flow in GDL is simulated using

Available online 7 September 2019

a volume of fluid (VOF) model. Corner droplet (on the GDL surface) and water flow (emerging from GDL bottom) are considered. It is found that two-phase flow in the GDL is

Keywords:

highly influenced by compression, fiber diameter, porosity, and GDL thickness. The results

Multiple parameters

indicate that a larger fiber diameter or higher porosity contributes to the water transport

Two-phase flow

due to larger average pore size. Furthermore, water removal from a thicker GDL is more

Gas diffusion layer

difficult, whereas water transport in the lower part of a compressed thick GDL is easy.

Compression

© 2019 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved.

Introduction Proton exchange membrane fuel cell (PEMFC) is a promising energy-conversion technology for vehicles owing to its high efficiency and power density [1]. PEMFC operation is highly influenced by gas diffusion layer (GDL), because GDL is a

crucial component in PEMFC and provides necessary pathways for reactants and products of the electrochemical reactions [2]. Excessive liquid water trapped in the GDL will block the transport of reactant gases and worsens cell performance [3,4]. PEMFC is always assembled under compressive loads to avoid gas leakage and to reduce contact resistance. In the clamping process, the GDL is compressed

* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (J. Xuan), [email protected] (K. Jiao). https://doi.org/10.1016/j.ijhydene.2019.08.108 0360-3199/© 2019 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

causing a variation of its microstructure with respective to pore size and pore distribution [5]. Due to the complex geometries (fibers with several parameters) and multiple working conditions (especially various clamping pressure conditions) involved, two-phase flow in GDL is not yet fully understood. Hence, further understanding of water transport in various GDLs is highly beneficial to the performance of PEMFCs.

Literature review During the last decade, numerous experimental studies have been conducted to investigate two-phase flow in gas diffusion layers (GDLs). Obeisun et al. [6] studied water droplet dynamics on the surface of a GDL by wettability and thermal characterization. Optical microscopy [7,8], neutron imaging [9e13], and X-ray computed tomography (CT) [14e21] are popular experimental techniques developed to investigate effects of various parameters on water transport in undeformed and deformed GDLs. Hinebaugh et al. [22] used a variety of imaging methods, including optical microscopy, microscale computed tomography, and scanning electron microscopy to characterize seven commercially available GDLs. Fishman et al. [7] observed side view images of liquidegasesolid interfaces during the evaporation of liquid water droplets on various commercial GDLs. Cooper et al. [13] adopted neutron radiography to examine the through-plane liquid water distribution of an operating PEMFC. Zenyuk et al. [17] used X-ray CT to investigate geometrical and channel effects on spatial liquid water distribution in compressed GDL and found water distribution was highly influenced by compression. Banerjee et al. [23] used X-ray Computed Tomography (X-CT) to study the impact of the s rib-channel compression on porosity profiles. Although several experimental researches have studied water transport in the GDL, they are inconvenient to investigate various parameters affecting two-phase flow in the GDL. In addition, the experimental techniques are always costly and have high requirement for facilities. As an alternative, numerical methods have been conducted using the volume-of-fluid (VOF) approach, lattice Boltzmann (LBM) method, and pore network model (PNM). Andersson et al. [24] employed a VOF model to investigate the interface resolved two-phase flow behavior inside the gas channel including the GDL surface. But the morphology of the GDL microstructures was neglected. In fact, achieving a realistic and accurate pore geometry is vital for a trustworthy pore scale simulation [25]. Effects of various material parameters such as PTFE distribution (wettability or contact angle), porosity, thickness, and surface morphology of GDLs on liquid water transport performed by these methods are mostly outlined in Table 1 [26e51]. PTFE content or PTFE distribution (wettability or contact angle) is a popular topic related to water transport in GDLs. Molaeimanesh et al. [26,27], Chen et al. [28], Sakaida et al. [30,33], Shakerinejad et al. [31], Jinuntuya et al. [34], and Sadeghi et al. [35] employed the LBM model to investigate twophase flow in GDLs when various wettability was considered. Molaeimanesh et al. [52] reviewed LBM simulation of PEMFCs, and discussed opportunities and challenges for this method.

26499

The effects of PTFE wettability are also investigated via the PNM [38] and VOF [42e45] modeling, taking mixed wettability and variable contact angle into account. The effects of GDL thickness [29,38], water droplet size [27], pressure difference [34], porosity [35,40], and rib width [33,37] on water transport are also included. Effects of several parameters on water transport have been studied in the above studies, but they did not take the clamping pressure into account. Clamping pressure affects the performance and structure of PEMFCs [53]. Jeon et al. [46], Molaeimanesh et al. [47], Bao et al. [50], and Han et al. [51] investigate the effects of compression ratio on water transport in the GDL microstructures. But simplified GDL structures were applied, such as two-dimensional fibers. X-ray CT is an accurate technology to reconstruct actual GDL microstructures, especially for compressed GDL microstructures. It was adopted in Refs. [48,49] to capture the compressed GDL microstructures for two-phase flow simulation. However, this technology is always costly and not convenient to reconstruct GDLs with various parameters (porosity, diameter, thickness, compression ratio etc.). Most of the published work in the literature is related to one or two of factors, which does not accurately reflect the realistic water transport process in the GDL. There are still some obvious limitations in the reported studies [26e51]: 1). Systematic investigation on the effects of multiple material parameters on two-phase flow in compressed GDLs is still needed; 2). Clear effects of various compressions on water transport and water distribution in the GDL need to be concluded; 3). Simplified stochastic method or high spatial resolution X-ray CT was adopted to reconstruct GDL microstructures, which is of high cost and lack of facility availability.

Present study In this study, a stochastic algorithm [54,55] is employed to reconstruct the uncompressed GDL microstructures. Solid mechanics simulations, which are based on the finite element method (FEM), are applied for capturing the compressed GDL microstructures. Two-phase flow in the compressed GDL is simulated using a volume of fluid (VOF) model. Both models are compared with experimental results [46]. Because functional requirements of the GDL are highly dependent on morphological and mechanical characteristics [56], systematic key parameters including carbon fiber diameter, porosity, GDL thickness, compression ratio and pore location at the GDL bottom where water emerges are investigated. To the best of the authors’ knowledge, this work represents the first study on the effects of multiple parameters on two-phase flow in the compressed GDL microstructures.

Model formulation Domain and boundary conditions Each computational domain consists of a three-dimensional (3D) fibrous GDL, a rib from bipolar plate (two ribs for validation sample), and two channels in the bipolar plate (one

26500

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Table 1 e The effects of various parameters on water transport in GDLs using LBM, PNM, and VOF simulations [26e51]. (SR: stochastic reconstruction; XR: X-ray reconstruction; SFS: simplified fiber structure). Literatures and Published year

Research Parameters

1. Compression neglected (LBM, PNM, and VOF) (1) LBM Molaeimanesh et al. [26] (2014) PTFE distribution Molaeimanesh and Akbari [27] (2016) Wettability and water droplet size Chen and Jiang [28] (2016) PTFE content and distribution Gao et al. [29] (2017) Porosity and thickness Sakaida et al. [30] (2017) Hydrophobic-hydrophilic pattern Shakerinejad et al. [31] (2018) Wettability of GDL Yu et al. [32] (2018) Water breakthrough Sakaida et al. [33] (2018) Wettability and the rib width Jinuntuya et al. [34] (2018) Contact angles and pressure differences Sadeghi et al. [35] (2018) Porosity, and contact angle Liu et al. [36] (2019) Heterogeneous multi-transport (2) PNM Lee et al. [37] (2014) Width of the rib, thickness of the GDL Shahraeeni and Hoorfar [38] (2014) GDL hydrophobicity and thickness Fazeli et al. [39] (2015) Boundary condition of water inlet Huang et al. [40] (2019) Network structure (uniform and gradient porosity) (3) VOF Chen et al. [41] (2013) GDL surface Ahmad et al. [42] (2013) PTFE loading Yin et al. [43] (2014) Variable contact angle Niu et al. [44] (2017) Liquid water at four locations, contact angle Niu et al. [45] (2018) Mixed wettability 2. Compression considered (LBM, PNM, and VOF) (4) LBM Jeon and Kim [46] (2015) Compression ratio Molaeimanesh et al. [47] (2019) Removal of a droplet from GDL Satjaritanun et al. [48] (2018) Water evolution, water saturation, and breakthrough pressure (5) PNM Fazeli et al. [49] (2016) Compression value (6) VOF Bao et al. [50] (2014) Droplet position and contact angles Han et al. [51] (2019) Air velocity and wettability Compression and multiple parameters None

channel for validation sample), as shown in Fig. 1. Other details of the FEM and VOF model will be discussed in the following subsections FEM modeling and VOF modeling.

FEM modeling The FEM modeling [57] is used to simulate the elastic deformation of fibrous GDLs. Compared with the large deformation of GDLs, the deformations of the catalyst layer (CL) and the membrane layer (ML), are always neglected due to their relatively larger elastic modulus [58]. Normal displacements (along negative z-axis direction) of bipolar plate (BP), GDL, and their interface are restricted. The value of clamping pressure exerted on the GDL depends on the selected parameters. Taking the default GDL (fiber diameter 8 mm, porosity 0.78, thickness 110 mm) as an example, clamping pressures 1.4 MPa, 2.8 MPa, 4.2 MPa and 5.6 MPa are adopted, corresponding to compression ratios 10%, 20%, 30% and 40%. The clamping pressures are applied downwards (along negative z-axis direction), modeling the fabrication stacking process. The mechanical parameters employed in this study are based on the measurements of the stress-strain relationships given by Mason et al. [59], as shown in Table 2.

Microstructure reconstruction

SR SR SR SR SR SR SR SR XR SR SR SFS SFS XR XR SFS SFS SR SR SR

SFS (2D) SFS (2D) XR XR SFS (plane) SFS (2D) SR

VOF model Interdigitated flow field [63] is one of the popular and efficient designs for fuel cells and has no continuous channels. Air is forced into the GDL and flows through the under-land GDL into lower-pressure channel, as shown in Fig. 1e. In addition, for serpentine flow field, two neighboring channels may also have such pressure difference causing the cross flow. In this study, the 3D computational domain is a small part (one rib, one higher-pressure channel and one lower-pressure channel) of an interdigitated flow filed or a serpentine flow field, as shown in Fig. 1. In fact, the liquid water is initially generated from two main processes: generation at the cathode CL owing to electrochemical reaction and vapor condensation on cold ribs. Correspondingly, water flow (emerging from GDL bottom) and initial corner droplet (on the GDL surface) are considered in the VOF model. The scenarios of droplet location and pore location are same with Ref. [44]. To simulate water emergence and growth of liquid water in the GDLs, the water is injected from a rectangular puncture with a cross section of 50 mm  50 mm. Three positions of the water injection hole at GDL bottom are considered in this study, as shown in Fig. 1f. Parameters used in the VOF model are gathered in Table 2.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

26501

Conservation equations e ¼ εx þ εy þ εz

(5)

FEM modeling Elastic deformation of the GDL is considered in this investigation. The relationship between the strain and stress components consists of balanced differential equations, geometric equations and physical equations. Relative equations are as follows [64]:

(1)

(2) Geometric equations:

(2)

(7)

The VOF method with interface reconstruction is strictly mass conservative. It is able to model immiscible fluids by solving a single set of momentum equations and then captures the volume fraction of each fluid throughout the computational domains. In the VOF model, water phase fraction a is introduced as the main variable for the model to solve. The value of a is zero at any region occupied by gas phase, and otherwise unity. Phase fraction between 0 and 1 means the cell contains gas, as well as liquid water. The phase indicator a is vital to calculate averaged volume density (r) and dynamic viscosity (m) for the two-phase mixture [44]: r ¼ rl a þ rg ð1  aÞ

(8)

m ¼ ml a þ mg ð1  aÞ

(9)

! V, U ¼0

(10)

Phase conservation equation:

(3)

va ! ! þ V , ð U aÞ þ V , ½ U r að1  aÞ ¼ 0 (11) vt ! where the symbol U represents the effective velocity vector which is shared by both two phases throughout the process. ! U r named ‘compression velocity’ indicates the relative velocity of liquid and gas at the interface and subscript rdenotes ! ! “relative velocity”. Both U and U r are calculated respectively by: ! ! ! U ¼ a U l þ ð1  aÞ U g

(12)

! ! ! Ur ¼ Ul  Ug

(13)

Momentum equation:

(3) Physical equations:

! vðr U Þ !! ! ! þ V , ðr U U Þ  V , ðmV U Þ  ðV U Þ , Vm ¼  Vpd  ! g ,! x Vr vt þ skVa

Normal stresses are calculated by: 8 > > < sx ¼ le þ 2Gεx sy ¼ le þ 2Gεy > > : sz ¼ le þ 2Gεz

E 2ð1 þ mÞ

G¼

where subscripts l and g indicate the liquid water and gas respectively. All the fundamental governing equations required for the VOF model here are shown as follows [45,65]: Continuity equation:

whereεx , εy andεz are the normal strains, gxy , gxz andgyz are the shear strains. 8  > > < gxy ¼ txy G gxz ¼ txz  =G > > : gyz ¼ tyz G

(6)

VOF model

wheresx , sy andsz are the normal stresses, txy , txz andtyz are the shear stresses, Fx , Fy and Fz are the body forces.

8 > 2 > > v2 εx v2 εy v gxy > > ¼0 þ 2  > 2 > vxvy > vy vx > > > > 2 > v2 ε > v2 εz v gyz > y > > > 2 þ 2  vyvz ¼ 0 > vz vy > > > > 2 > > > v ε v2 ε v2 gzx > > 2z þ 2x  ¼0 < vzvx vx vz   > > v vgxy vgxz vgyz v2 εx > > þ  2 ¼0 > > vx vz > vy vx vyvz > > > >   > > v2 εy v vgyx vgyz vgxz > > þ  2 ¼0 > > > vy vz vx vy vxvz > > > >   > > v vgzx vgzy vgxy v2 εz > > > þ  2 ¼0 > : vz vy vx vz vxvy

mE ð1 þ mÞð1  2mÞ

where l is the Lame constant, G is the shear elastic modulus, E is the Young’s modulus, mis the Poisson ratio.

(1) Balanced differential equations:

8 > vt > vs vt > > x þ xy þ xz þ Fx ¼ 0 > > vx vy vz > > > > < vs vtyx vtyz y þ þ þ Fy ¼ 0 > vy vx vz > > > > > vsz vtzx vtzy > > > þ þ þ Fz ¼ 0 > : vz vx vy

l¼

(4)

(14) wherea, s and k denote the phase fraction, surface tension coefficient and mean curvature of the phase interface

26502

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Fig. 1 e Schematics of numerical model based on GDL microstructures in the proton exchange membrane fuel cell (PEMFC): (aed) The computational domains and boundary conditions; (e) Typical interdigitated and serpentine flow filed of a PEMFC; (f) water droplet location (on the GDL surface) and pore location (at the GDL bottom).

respectively. pd is the modified pressure for simplifying boundary conditions, and it is determined by: g ,! x pd ¼ p  r!

(15)

The subscript d means “dynamic”. Here ! g indicates the gravity vector and ! x represents the position vector.

The continuum surface force (CSF) model is implemented in this model to take the effects of the surface tension at the gas-liquid interface into account by adding a force fs to Eq. (14), which is defined as follows [44]: fs ¼ skVa wherek is the mean curvature of the phase interface:

(16)

26503

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Table 2 e List of the parameters and their values in the study. Parameter name Operating parameters Water inlet velocity (m/s) Gas velocity (m/s) Higher pressure of gas channel (kPa) Lower pressure of gas channel (kPa) Compression ratio Geometric parameters Fiber diameter (mm) Porosity Thickness (mm) Toray H-060 Physical parameters Water viscosity (Pa$s) Surface tension (N/m) Contact angle Young’s modulus (MPa) Poisson’s ratio

  Va k¼ V,! n ¼  V, jVaj

Value

References

1 10 water droplet: 10, water emerging flow: 7 0 0%, 10%, 20%, 30%, 40%

[44] [44] [44] [44] [46]

7, 8, 9 0.7(Assumed), 0.78, 0.86 110, 190

[60] [28], Sigracet 25BA Toray H-030,

0.65 0.07 Fiber: 90 , 120 ; wall: 90 BP: 13000, fiber: 300 BP: 0.26, fiber: 0.256

[61] [62] [44] [57,59], Assumed [57,59]

(17)

the surface unit normal ! n is adjusted in the cells close the wall according to the equation as follow ! ! (18) n ¼! n w cos q þ t w sin q ! where ! n w and t w are the unit vector normal to the wall and unit vector tangential to the wall, the subscript w indicates the wall. Symbol q is the contact angle, and it is assumed constant here.

Model assumptions In this study, assumptions for the FEM and VOF models are outlined as follows: 1) The components in both cathode and anode of the PEMFC have the similar mechanical behavior; 2) Elastic deformation of GDLs is assumed throughout this investigation; 3) Two-phase flow inside GDLs are laminar and incompressible; 4) Transport properties of the fluids keep constant and no phase change occurs in the process.

split operator) PISO scheme coupling the semi-implicit method is adopted for the pressure and velocity coupling solution. Owing to the outstanding parallelism of the popular open-source software Open FOAM, the computational

Table 3 e Features of the meshes employed to study the effect of multiple parameters. Type

Mesh Size

Mesh Number (Million)

Part 1: Model Validation FEM 2.0 mm 200 mm 0.7 4.25 VOF 3.0 mm 1.56 Part 2: Mesh Independency Test FEM 1.3 mm 200 mm 0.7 10.17 3.0 mm 1.81 VOF 2.0 mm 0.64 4.0 mm 4.02 Part 3: Formal Cases (Default parameters: carbon fiber diameter 8 mm, porosity 0.78, thickness 110 mm, adopted mesh size FEM-2.0 mm/VOF-3.0 mm) (CR: compression ratio, All: water droplet and water flow)

Type FEM

Objective

Details

Compression Fiber diameter

6.80 million 7 mm (7.10 million) 9 mm (6.66 million) 0.86 (5.50 million), 0.7 (7.11 million) 190 mm (9.60 million)

Porosity Thickness

Two-phase Values Flow

Numerical procedures Both computational domains of the FEM the VOF models are discretized with hexahedral meshes, as shown in Fig. 1. Details of the mesh size and mesh number adopted are gathered in Table 3. Due to the existence of the bolts, the contact surfaces between all components are regarded as bonded to avoid slide. The commercial software ANSYS is used for the FEM modeling. Preconditioned conjugate gradient (PCG) method is applied to solve the conservation equations. To validate the FEM model, the GDL with thickness 200 mm and porosity 0.7 is reconstructed. The governing equations for the VOF model are discretized by a second order scheme. A (pressure implicit

Thickness Porosity

VOF Compression

Water droplet 0% 10% 20% 30% 40% Water flow 0% 20% 40% Fiber diameter All 7 mm 9 mm Porosity 0.7 0.84 Thickness 190 mm

Mesh Number (Million) 8.24 7.93 7.67 7.34 7.06 8.24 7.67 7.06 7.01 6.99 6.77 7.16 8.26

26504

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Appropriate compression of the GDL is a trade-off between having a low electrical contact resistance and effective transport of reactants and products [66]. In this section, the effects of four commonly adopted compression ratios 10%, 20%, 30%, and 40% on two-phase flow in GDLs are investigated. Fig. 4 shows the compressed GDL microstructures and the detailed pore morphology of the various GDLs at the slice of y ¼ 300 mm.

channel. No separated small droplet occurs throughout the process due to the sufficient porosity. The whole droplet reaches the lower-pressure channel at 0.9 ms. Fig. 5b shows the dynamic behavior of a droplet flowing through the under-land GDL with compression ratio 10% (CR ¼ 10%). Different from Fig. 5a, smaller droplets are separated from the initial droplet by the fronted fibers. It is caused by the relatively large size of the initial droplet compared with small pores in the compressed GDL. In addition, the initial droplet is forced to move in the direction of the airflow (positive y-axis) at the early stage (see t ¼ 0.9 ms). This is attributed to the difficulty of the large droplet entering the GDL in time due to a sudden decrease of the pore size close to the land. It is worth noting that the main part of the droplet does not reach the lowerpressure channel until 5.7 ms. The water dynamic behavior of a droplet in the GDL with compression ratio 20% (CR ¼ 20%) is depicted in Fig. 5c. It can be observed that apparent two-phase flow characteristics in both cases with compression 10% and 20% are similar, whereas more portions of the droplet are trapped in the pores between fibers under the land. Therefore, after 5.7 ms, some water still remains in the under-land GDL. This is because these small droplets spread on the surface of carbon fibers. The air from the higher-pressure channel can hardly move them when the pore size decreases to a certain level owing to compression. Fig. 5d shows the droplet behavior in the GDL with compression ratio 30% (CR ¼ 30%). It is observed that the water droplet no longer moves with airflow (positive y-axis) in the higher-pressure channel (see t ¼ 1.8 ms). After flowing through the under-land GDL, the smaller droplets accumulate in the location near center of lower-pressure channel (see t ¼ 8.1 ms). Comparing Fig. 5d with Fig. 5e, it is easy to find that less time is needed for the main part of the droplet passing through the under-land GDL for 40% compression (see t ¼ 4.5 ms in Fig. 5d and e). Moreover, the droplet stops in the center of the lower-pressure channel. Because of the reduction of the cross-sectional area under land, both capillary pressure and air velocity will increase. Higher compression of GDL leads to higher average velocity of air flow traveling through the under-land GDL, and then faster removal process of droplet may be expected. When the growth rate of positive force from the air flow is larger than that of negative force from the capillary pressure, the droplet removal process from the higher-pressure channel will be faster as compression increases. The quicker droplet transport behavior for CR ¼ 40% is caused by the higher air velocity. The results above indicate that compression does hinder the water droplet transport between two adjacent channels, whereas this hindering effect will be weakened when compression reaches a critical value (30% in this study).

Corner droplet on the GDL surface

Water emerging from the bottom of GDL

The dynamic behaviors of liquid water droplets with 100 mm initial radius during removal from higher-pressure channel are investigated. The initial position of the droplet can be found in Fig. 1f. As Fig. 5a shows, the water droplet flows through the uncompressed under-land GDL (CR ¼ 0%). At the initial stage, liquid droplet is forced into the pores between carbon fibers owing to the pressure difference between two adjacent channels, and then flows to the lower-pressure

This section considers three different pore locations of liquid water emerging (x ¼ 150 mm, 300 mm, and 450 mm) from the bottom of GDLs with compression ratios 0%, 20%, and 40% respectively. The schematic of three locations is shown in Fig. 1f.

domains are decomposed with 140 Intel Xeon @ 2.93 GHz processors. Adjustable time step is used in VOF simulation. The final time step maintains around 6  108 s. Two hours are required for reconstructing each kind of uncompressed GDL microstructures. Each case for the FEM modeling takes about 24 h by using 32 Intel Xeon @3.50 GHz processors in parallel. As for VOF model, each water droplet (on the GDL surface) case takes about 72 h, and each water flow (emerging from GDL bottom) case takes 48 h, by using 140 Intel Xeon @ 2.93 GHz processors.

Results and discussions The following sections are presented by discussing compression (section Effect of compression) and various parameters of GDLs (i.e. carbon fiber diameter (section Effect of carbon fiber diameter), porosity (section Effect of porosity), and thickness (section Effect of GDL thickness) affecting two-phase flow in uncompressed and compressed GDLs.

Model validation In this study, emphasis is placed on the effects of multiple parameters on the behavior of liquid water transport in the GDL. The validations of the FEM and VOF models are carried out in a qualitative way here. The simulated deformation of GDL microstructures and water emerging behavior in the GDL with 40% compression are compared with the experimental results reported in Ref. [46]. Both employed samples for validation are set the same geometrical and operating conditions with the experimental samples. The detailed geometric parameters of the samples are outlined in Fig. 1. Clearly, both simulated results are in good agreement with the experimental results, as shown in Fig. 2. To guarantee sufficient grids adopted, grid dependency tests are also performed by testing three mesh sizes for FEM and VOF models respectively. As shown in Fig. 3, mesh sizes adopted in this paper indicate grid independence.

Effect of compression

Middle location of under-channel GDL bottom. Fig. 6a shows the water dynamic behavior emerging from middle location of

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

26505

Fig. 2 e Comparisons between the simulated results and experimental results [40]. (a) compressed GDL microstructures (CR: compression ratio, CR ¼ 10%, 20%, 30%, and 40%); (b) water breakthrough behavior in GDL with compression ratio 40% (CR ¼ 40%).

26506

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Fig. 3 e Mesh independency tests for FEM modeling and VOF model. (a) Compressed GDL microstructures (CR ¼ 10%, 20%, 30%, and 40%); (b) water behaviors before and after breakthrough in GDL with compression ratio 40% (CR ¼ 40%); (3) compression.

under-channel GDL bottom (x ¼ 150 mm) in the uncompressed GDL. It can be observed that most water is discharged into the higher-pressure channel and only a little water flows through under-land GDL and reaches the lower-pressure channel. This phenomenon is same as the reported results concluded by Ref. [44]. However, this phenomenon is no longer exists when compression effects are considered, as shown in Fig. 6b. All of the water accumulates in the GDL under higher-pressure channel. This is caused by the increase of the resistance from higher capillary pressure due to the decrease of porosity. At the later stage, the liquid water cluster grows up and attaches to the left wall of the higher-pressure channel. When the GDL is compressed by 40%, the liquid water behavior is similar with that of CR ¼ 20%, as shown in Fig. 6c. This is because the compression only affects the under-land GDL

while has minor effects on the GDL under the channel. Once the carbon fibers under the land are compressed to a critical level which is enough to prevent the liquid water pass through the pores between under-land fibers, the increase of compression has minor effects on the liquid water behavior in the higher-pressure channel.

Location of GDL bottom at land-channel interface. Fig. 7a shows the water dynamic behavior emerging from location of uncompressed GDL bottom at land-channel interface (x ¼ 300 mm). It can be observed that water is totally discharged into the lower-pressure channel and no water flows through the GDL under the higher-pressure channel. However, some water flows through the pores between carbon fibers under the higher-pressure channel and forms water film adhering to

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

26507

Fig. 4 e Schematics of compressed GDL microstructures: (a) Compressed GDL with compression 10%, 20%, 30%, and 40% for formal cases; (b)e(j) detailed pore morphology of the compressed GDLs with multiple parameters at the slice of y ¼ 300 mm.

the wall near land-channel interface when compression effects are considered, as shown in Fig. 7b. Moreover, water breakthrough is found in both higher-pressure channel and lower-pressure channel instead of only the lower-pressure channel. This is caused by the higher transport resistance in the under-land GDL. Most water emerging from the channel side is unable to pass through the under-land GDL and accumulates in the higher-pressure GDL. When the under-land GDL is compressed by 40%, an interesting phenomenon is found that water emerging from the pore location at the landchannel interface prefers flowing into the higher-pressure channel, as shown in Fig. 7c. More water is found in the higher-pressure channel than that in the lower-pressure channel (see t ¼ 2.7 ms). This is because the dominant role which controls the water flow behavior changes between capillary pressure and pressure difference between two adjacent channels. For the uncompressed GDL, the pressure difference is the dominant role forcing water to flow into the lower-pressure channel. As the GDL is compressed, the decrease of porosity increases the capillary pressure hindering the water flow to the lower-pressure channel. Then some water accumulates in the GDL under higher-pressure channel and flows out to the higher-pressure channel.

Middle location of under-land GDL bottom. Fig. 8 shows the dynamic behaviors of liquid water emerging from the middle location of the under-land GDL bottom (x ¼ 450 mm). It can be observed that in the uncompressed GDL, the water flows to

the lower-pressure channel totally due to sufficient pressure difference between two adjacent channels, as shown in Fig. 8a. Moreover, the water tends to accumulate in the middle of the lower-pressure channel at the early stage (0e1.5 ms), and attaches to the side of the land forming water film in the later period. Water breakthrough is found in the position near the center of the lower-pressure channel. However, water tends to attach to the side of land once it reaches the lowerpressure channel in the GDL with compression ratio 20% (CR ¼ 20%), as shown in Fig. 8b. Water breakthrough is observed at the land/channel interface. This is because the GDL retreats the porous area as compression ratio increases. Therefore, the dominant path of water breakthrough is changed. As expected, more water is trapped in the pores between fibers under the land than that of CR ¼ 0%. Fig. 8c shows the similar water dynamic behavior with Fig. 8b, whereas more water is found to accumulates in the pores between fibers under the land. Moreover, the water is found spreading toward the place close to the higher-pressure channel. The results of section Water emerging from the bottom of DGL indicate that the effects of compression on the behavior of water emerging from the GDL bottom are vital to investigate the water distribution in the GDL. This is because the GDL retreats the porous region as the compression ratio increases, and hence the pore size and pore distribution are changed in the GDL microstructures. Water emerging from the middle location of the under-channel GDL bottom (x ¼ 150 mm) is

26508

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Fig. 5 e Dynamic behaviors of a liquid water droplet transport from the higher-pressure channel to the lower-pressure channel (initial porosity: 0.78, fiber diameter: 8 mm, and initial thickness: 110 mm) with different compression ratios: (a) CR ¼ 0%; (b) CR ¼ 10%; (c) CR ¼ 20%; (d) CR ¼ 30%; (e) CR ¼ 40%.

trapped in the higher-pressure channel totally when compression effects are considered. Water emerging from location of GDL bottom at land-channel (higher-pressure channel) interface (x ¼ 300 mm) tends to flow to higherpressure channel when compression effects are considered. Partial water emerging from the middle location of under-land GDL bottom (x ¼ 450 mm) tends to flow towards the place close to the higher-pressure channel when compression effects are considered. In addition, the wet condition of the under-land GDL is highly influenced by compression. The results of this subsection will be meaningful to further understand the role of compression in the removal of water generated from the reaction surface on the catalyst layer under actual working conditions.

The results of section Effect of compression indicate that compression should be considered when two-phase flow is investigated in the GDL. In general, compression will hinder the water transport between two adjacent channels resulting in a different wet condition of the under-land GDL. The results are useful to provide an important reference for optimizing the design of the flow field. Especially for the interdigitated flow filed, assembly pressure should be carefully controlled to avoid excessive compression of GDL. As mentioned in the introduction, GDL is a crucial component of PEMFC and material design of GDL is important for water management. Further understanding the physical properties of the GDL is also important for developing a GDL microstructure model, and these properties constitute the variable parameters in

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

26509

Fig. 6 e Dynamic behaviors of liquid water emerging from middle location of under-channel GDL bottom (x ¼ 150 mm) (initial porosity: 0.78, fiber diameter: 8 mm, and initial thickness: 110 mm) with different compression ratios: (a) CR ¼ 0%; (b) CR ¼ 20%; (c) CR ¼ 40%.

Fig. 7 e Dynamic behaviors of liquid water emerging from location of GDL bottom at land-channel interface (x ¼ 300 mm) (initial porosity: 0.78, fiber diameter: 8 mm, and initial thickness: 110 mm) with different compression ratios: (a) CR ¼ 0%; (b) CR ¼ 20%; (c) CR ¼ 40%.

26510

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Fig. 8 e Dynamic behaviors of liquid water emerging from middle location of under-land GDL bottom (x ¼ 450 mm) (initial porosity: 0.78, fiber diameter: 8 mm, and initial thickness: 110 mm) with different compression ratios: (a) CR ¼ 0%; (b) CR ¼ 20%; (c) CR ¼ 40%.

constructing a 3D GDL model. Several material parameters of the GDL are also included in this study to provide a further understanding of water transport in the compressed GDL. The results can provide an important reference when undesirable consequences of GDL design are considered. In the following sections, the effects of carbon fiber diameter, porosity, and thickness on two-phase flow in compressed GDLs with compression ratio 40% are investigated.

Effect of material parameters Effect of carbon fiber diameter Carbon fiber structures plays an important role in electrical and thermal resistances [67]. Increment of fiber diameter will decrease the total resistance [68]. However, studies on effect of carbon fiber diameter to water transport behavior in GDL are scarce. The effects of three different carbon fiber diameters (7 mm, 8 mm, and 9 mm) on two-phase flow are investigated. Fig. 4 (e)e(g) show the slice of GDL microstructures with three carbon fiber diameters at y ¼ 300 mm. It is observed that smaller pores are presented in the GDL with smaller carbon fiber diameter. Since the porosity is set constant, more pores exist in the GDL with smaller fiber diameter. Two-phase flow in the GDL with carbon diameter 8 mm has been discussed in section Effect of compression, as shown in Figs. 5e and 7c respectively. In terms of corner droplet on the GDL surface, larger carbon fiber diameter is found beneficial to water transport and less time is required in the case with fiber diameter 9 mm, as shown in Figs. 9a, b, and 5e. Less water is found trapped in the pores between the carbon fibers under

the land as the carbon fiber diameter increases. Moreover, once the droplet first reaches the lower-pressure channel, its appearance is also different for various fiber diameters. This is caused by the smaller pore size in the GDL with smaller fiber diameter. Initial droplet tends to be separated into small droplets when a smaller fiber diameter is adopted. However, it should be noted that the apparent droplet transport behaviors for three fiber diameters are similar. As for water emerging from the location of GDL bottom at land-channel interface (x ¼ 300 mm), it is observed that more water is forced into the higher-pressure channel as the carbon fiber diameter decreases, as shown in Figs. 9c, d, and 7c. In addition, it is found that water breakthrough first occurs in the higher-pressure channel for diameters of 7 mm and 8 mm, whereas the first water breakthrough is found in the lowerpressure channel for the diameter of 9 mm. Wet condition of the under-land GDL differs in these three carbon fiber diameters, and the water tends to diffuses around in the underland GDL as the diameter increases. This is because the larger averaged pore size between fibers with larger diameter is formed. As reported in Ref. [69], the capillary pressure pc is defined as pc ¼ pa  pw ¼

4s cosqf dpore

(19)

where pa and pw are the air and water pressures; dpore and s are the diameter of GDL pore and the surface tension coefficient of water. As the pore diameter increases, the capillary pressure decreases. Hence, water transport resistance in the under-

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

26511

Fig. 9 e Two-phase flow in the GDL with different carbon fiber diameter (initial porosity: 0.78, initial thickness: 110 mm, CR ¼ 40%, x ¼ 300 mm): (a) Water droplet in the GDL with fiber diameter7 mm; (b) water droplet in the GDL with fiber diameter 9 mm; (c) water emerging flow in the GDL with fiber diameter7 mm; (d) water emerging flow in the GDL with fiber diameter 9 mm.

land GDL is smaller for the larger carbon fiber diameter. When the porosity is set constant, more pores are found in the GDL with the fiber diameter 7 mm, because the GDL with the fiber diameter 7 mm has smaller averaged pores than that with the fiber diameter 9 mm, as shown in Fig. 4. Two-phase flow behaviors mentioned above also indicate that pore size has a more dominant effect than pore number on water transport in the GDL with same porosity. It is worth noting that carbon fiber diameter is also an important material parameter for water transport in the GDL besides porosity.

Effect of porosity The porosity of the GDL affects its other physical properties [70] such as conductivity, thermal conductivity and permeability of the GDL. In this section, the effects of three porosities of 0.7, 0.78, and 0.84 on water transport behavior and water distribution in the compressed GDL with compression ratio 40% are investigated. Figs. 4e, h, and 4i show the slices of GDL microstructures with different porosity at y ¼ 300 mm. It can be observed directly that the average pore size is larger when a higher porosity is adopted. In addition, the fiber distribution under the land is denser as the porosity decreases. Two-phase

flow in the GDL with porosity 0.78 is previously completed in section Effect of compression, as shown in Figs. 5e and 7c respectively. In terms of corner droplet on the GDL surface, only a relatively small portion of the water droplet passes through the pores under the land and reaches the lower-pressure channel when a porosity of 0.7 is adopted, as shown in Fig. 10a. This is because a high capillary pressure is achieved due to porosity reduction. The pressure difference between two adjacent channels is not able to overcome the resistance from the capillary pressure. When a higher porosity is adopted, the water droplet easily reaches the lower-pressure channel and less time is required for the transport process. In addition, it is observed that higher porosity contributes to droplet accumulating near the center of lower-pressure channel. This is caused by the lower resistance in the underland GDL. The droplet has a faster speed to pass through the under-land GDL and reaches the location away from the land/ channel interface. Since water transport is important for flow field, enough porosity should be guaranteed for sufficient water removal requirement. In addition, the compression should be fully considered when the commercial GDL with

26512

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Fig. 10 e Two-phase flow in the GDL with different porosity (fiber diameter: 8 mm, initial thickness: 110 mm, CR ¼ 40%, x ¼ 300 mm): (a) Water droplet in the GDL with porosity 0.7; (b) water droplet in the GDL with porosity 0.86; (c) water emerging flow in the GDL with porosity 0.7; (d) water emerging flow in the GDL with porosity 0.86.

appropriate porosity is chosen. Otherwise, water accumulated in the higher-pressure channel will cause water flooding which will worsen the cell performance. In terms of water emerging from the location of GDL bottom at land-channel interface (x ¼ 300 mm), the water tends to spread around in the under-land GDL and less water is able to reach the lower-pressure channel when a lower porosity is adopted as shown in Fig. 10c, d, and 7c. More water is observed trapped in the pores between fibers under the land as the porosity decrease. For the porosity of 0.86, water is discharged into lower-pressure channel totally, whereas water breakthrough occurs in both higher-pressure and lower-pressure channel for porosities of 0.7 and 0.78. This is caused by the higher resistance in the under-land GDL due to the lower porosity. These facts indicate that for the GDL with same carbon fiber diameter, higher porosity will contribute to improve the water transport between the adjacent channel and reduces the water accumulation in the higher-pressure channel. Moreover, porosity has a vital effect on water distribution in the GDL.

Effect of GDL thickness One of the important factors in designing and selection of the GDL is the thickness which can affect water

management and reactant gases delivery to the active area [71]. The microstructures of the compressed GDLs with original thickness 110 mm and 190 mm have been displayed in Fig. 4e and j. Comparing Fig. 11a with Fig. 11b, water droplet is found hindered in the cross-flow process. The droplet reaches the lower-pressure channel after 5.4 ms for the 190 mm thick GDL, while only 4.5 ms is required for the 110 mm thick GDL. This may be caused by the fact that droplet will first move downward and then pass through the pores in the transition region due to uneven GDL deformation and pressure difference between two adjacent channels. Herein, for the thicker GDL sample, more time is needed for the whole transport process. Moreover, it is observed that almost no droplets are trapped in the pores between the under-land fibers in the thicker GDL, as shown in Fig. 11b. This may be caused by the fact that the lower part of the thicker GDL is sufficient for the whole droplet to pass through the under-land GDL. Water droplets are also found to accumulate near the center of the lower-pressure channel for the thinner GDL, while accumulate at the interface between the lower-pressure channel and the land for the thicker GDL sample. This may be caused by the lower pressure gradient in the thicker GDL microstructures. Lower pressure gradient means lower velocity of the droplet, so that the

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

26513

Fig. 11 e Two-phase flow in the GDL with different thickness (porosity 0.78, carbon diameter: 8 mm, CR ¼ 40%, x ¼ 300 mm): (a) Water droplet in the GDL with thickness 110 mm; (b) water droplet in the GDL with thickness 190 mm; (c) water emerging flow in the GDL with thickness 110 mm; (d) water emerging flow in the GDL with thickness 190 mm.

droplets have insufficient time to reach the middle of the lower-pressure channel. Comparing Fig. 11c with Fig. 11d, the water emerging from the bottom of the thicker GDL tends to diffuse in pores between fibers and the water breakthrough occurs on the side wall away from the land. A larger portion of the GDL is filled with water while no water film is found at the interface between the lower-pressure channel and the land. It is observed that water breakthrough in the thicker GDL is more difficult than that in the thinner GDL, and more time is required for the breakthrough process in the thicker GDL. Moreover, water reaches the right wall (away from the land) of the lowerpressure channel and water breakthrough occurs in this region. This is different from the phenomenon found in the thinner GDL. The result is attributed to that water emerging from the thicker GDL bottom encounters smaller resistance due to thicker lower part of the GDL at the early stage. It is easier for the airflow to drive the water flow. Both the underland and under-channel areas of the GDL are found filled with large amount of water for the thicker GDL, whereas most water is trapped in the under-land region for the thinner GDL. The results in this section indicate that the thickness has vital effects on the water transport in the GDL. Water

redistribution occurs with the change of the thickness. A thicker GDL is not beneficial to water removal from the GDL. In terms of water droplet in the higher-pressure channel, dryer under-land GDL is found in the thicker GDL. In addition, water droplet tends to accumulate in the channel-land interface (lower-pressure channel). As for water emerging from the bottom of the GDL, thicker GDL reduces the formation of water film at the channel-land interface, whereas makes the whole GDL more humid. This water distribution difference between water droplet and water emerging flow may be caused by the different contact angles set for GDLs. In addition, water flow in the emerging process has an original inlet velocity which may make up the pressure gradient drop in the thicker GDL. In general, a thicker GDL is not beneficial to water removal from the GDL, and the wet condition of the GDL is highly influenced by GDL thickness.

Conclusion In this study, a stochastic algorithm is adopted to reconstruct the GDL microstructures. The deformation of the microstructures is simulated via a finite element method (FEM).

26514

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

Two-phase flow in the uncompressed and compressed GDLs is investigated using a volume of fluid (VOF) method. The effects of compression and multiple parameters (carbon fiber diameter, porosity, and GDL thickness) on two-phase flow in GDLs are investigated. Both initial corner stagnant water droplet (on the GDL surface) and water flow (emerging from the GDL bottom) are considered. The main conclusions are summarized as follows: (1) Two-phase flow in the GDL is highly influenced by compression. Compression hinders water transport between two adjacent channels and more water is found trapped in the under-land GDL as compression ratio increases. Moreover, water redistribution occurs when compression effects are considered. (2) As fiber diameter increases, larger average pore size is observed in the GDL. Larger fiber diameter contributes to water transport in the GDL. (3) As porosity increases, larger average pore size and more pores are observed in the GDL. Higher porosity contributes to water transport between two adjacent channels. Sufficient porosity should be guaranteed to avoid electrode water flooding. (4) As for a thick GDL, water transport in the lower part is easier than that in the lower part of a thin GDL when the same compression ratio is considered. But a thick GDL is undesirable for entire process of water removal from the GDL.

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

Acknowledgements This work is supported by the National Key Research and Development Program of China (Grant No. 2017YFB0102703), the China-UK International Cooperation and Exchange Project (Newton Advanced Fellowship) jointly supported by the National Natural Science Foundation of China (Grant No. 51861130359) and the UK Royal Society (Grant No. NAFyR1y180146), and the National Natural Science Foundation of Tianjin (China) for Distinguished Young Scholars (Grant No. 18JCJQJC46700).

[12]

[13]

[14]

Appendix A. Supplementary data [15]

Supplementary data to this article can be found online at https://doi.org/10.1016/j.ijhydene.2019.08.108. [16]

references [17] [1] Wang Y, Chen KS, Mishler J, et al. A review of polymer electrolyte membrane fuel cells: technology, applications, and needs on fundamental research. Appl Energy 2011;88(4):981e1007. [2] Hinebaugh J, Gostick J, Bazylak A. Stochastic modeling of polymer electrolyte membrane fuel cell gas diffusion layersePart 2: a comprehensive substrate model with pore

[18]

size distribution and heterogeneity effects. Int J Hydrogen Energy 2017;42(24):15872e86. Liu X, Peng F, Lou G, Wen Z. Liquid water transport characteristics of porous diffusion media in polymer electrolyte membrane fuel cells: a review. J Power Sources 2015;299:85e96. Omrani R, Shabani B. Gas diffusion layer modifications and treatments for improving the performance of proton exchange membrane fuel cells and electrolysers: a review. Int J Hydrogen Energy 2017;42(47):28515e36. Xu Y, Qiu D, Yi P, Lan S, Peng L. An integrated model of the water transport in nonuniform compressed gas diffusion layers for PEMFC. Int J Hydrogen Energy 2019;44(26):13777e85. Obeisun OA, Finegan DP, Engebretsen E, et al. Ex-situ characterisation of water droplet dynamics on the surface of a fuel cell gas diffusion layer through wettability analysis and thermal characterisation. Int J Hydrogen Energy 2017;42(7):4404e14. Fishman JZ, Leung H, Bazylak A. Droplet pinning by PEM fuel cell GDL surfaces. Int J Hydrogen Energy 2010;35(17):9144e50. Jiao K, Park J, Li X. Experimental investigations on liquid water removal from the gas diffusion layer by reactant flow in a PEM fuel cell. Appl Energy 2010;87(9):2770e7. Boillat P, Oberholzer P, Kaestner A, Siegrist R, Lehmann EH, Scherer GG, Wokaun A. Impact of water on PEFC performance evaluated by neutron imaging combined with pulsed helox operation. J Electrochem Soc 2012;159(7):F210e8. Park J, Li X, Tran D, Abdel-Baset T, Hussey DS, Jacobson DL, Arifc M. Neutron imaging investigation of liquid water distribution in and the performance of a PEM fuel cell. Int J Hydrogen Energy 2008;33(13):3373e84. Hussey DS, Spernjak D, Weber AZ, Mukundan R, Fairweather J, Brosha EL, Davey J, Spendelow JS, Jacobson DL, Borup RL. Accurate measurement of the through-plane water content of proton exchange membranes using neutron radiography. J Appl Phys 2012;112(10):104906. Park H. Effect of the hydrophilic and hydrophobic characteristics of the gas diffusion medium on polymer electrolyte fuel cell performance under non-humidification condition. Energy Convers Manag 2014;81:220e30. Cooper NJ, Santamaria AD, Becton MK, Park JW. Neutron radiography measurements of in-situ PEMFC liquid water saturation in 2D & 3D morphology gas diffusion layers. Int J Hydrogen Energy 2017;42(25):16269e78. € tter H, George MG, Liu H, Ge N, Lee J, Ince UU, Marko Bazylak A. Effects of compression on water distribution in gas diffusion layer materials of PEMFC in a point injection device by means of synchrotron X-ray imaging. Int J Hydrogen Energy 2018;43(1):391e406. Sassin MB, Garsany Y, Gould BD, Swider-Lyons K. Impact of compressive stress on MEA pore structure and its consequence on PEMFC performance. J Electrochem Soc 2016;163(8):F808e15. Zenyuk IV, Parkinson DY, Connolly LG, Weber AZ. Gasdiffusion-layer structural properties under compression via X-ray tomography. J Power Sources 2016;328:364e76. Zenyuk IV, Parkinson DY, Hwang G, Weber AZ. Probing water distribution in compressed fuel-cell gas-diffusion layers using X-ray computed tomography. Electrochem Commun 2015;53:24e8. Satjaritanun P, Weidner JW, Hirano S, Lu Z, Khunatorn Y, Ogawa S, Shimpalee S. Micro-scale analysis of liquid water breakthrough inside gas diffusion layer for PEMFC using Xray computed tomography and Lattice Boltzmann Method. J Electrochem Soc 2017;164(11):E3359e71.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

€ tzke C, Gaiselmann G, Osenberg M, Arlt T, Marko € tter H, [19] To Hilger A, Manke I. Influence of hydrophobic treatment on the structure of compressed gas diffusion layers. J Power Sources 2016;324:625e36. [20] Lee J, Yip R, Antonacci P, Ge N, Kotaka T, Tabuchi Y, Bazylak A. Synchrotron investigation of microporous layer thickness on liquid water distribution in a PEM fuel cell. J Electrochem Soc 2015;162(7):F669e76. [21] Sabharwal M, Gostick JT, Secanell M. Virtual liquid water intrusion in fuel cell gas diffusion media. J Electrochem Soc 2018;165(7):F553e63. [22] Hinebaugh J, Bazylak A. Stochastic modeling of polymer electrolyte membrane fuel cell gas diffusion layersePart 1: physical characterization. Int J Hydrogen Energy 2017;42(24):15861e71. [23] Banerjee R, Hinebaugh J, Liu H, Yip R, Ge N, Bazylak A. Heterogeneous porosity distributions of polymer electrolyte membrane fuel cell gas diffusion layer materials with ribchannel compression. Int J Hydrogen Energy 2016;41(33):14885e96. [24] Andersson M, Beale SB, Reimer U, Lehnert W, Stolten D. Interface resolving two-phase flow simulations in gas channels relevant for polymer electrolyte fuel cells using the volume of fluid approach. Int J Hydrogen Energy 2018;43(5):2961e76. [25] Shojaeefard MH, Molaeimanesh GR, Nazemian M, Moqaddari MR. A review on microstructure reconstruction of PEM fuel cells porous electrodes for pore scale simulation. Int J Hydrogen Energy 2016;41(44):20276e93. [26] Molaeimanesh GR, Akbari MH. Impact of PTFE distribution on the removal of liquid water from a PEMFC electrode by lattice Boltzmann method. Int J Hydrogen Energy 2014;39(16):8401e9. [27] Molaeimanesh GR, Akbari MH. Role of wettability and water droplet size during water removal from a PEMFC GDL by lattice Boltzmann method. Int J Hydrogen Energy 2016;41(33):14872e84. [28] Chen W, Jiang F. Impact of PTFE content and distribution on liquidegas flow in PEMFC carbon paper gas distribution layer: 3D lattice Boltzmann simulations. Int J Hydrogen Energy 2016;41(20):8550e62. [29] Gao Y, Montana A, Chen F. Evaluation of porosity and thickness on effective diffusivity in gas diffusion layer. J Power Sources 2017;342:252e65. [30] Sakaida S, Tabe Y, Chikahisa T. Study on gas diffusion layer structure tolerant to flooding in PEFC by scale model experiment and LBM simulation. ECS Trans 2017;80(8):123e31. [31] Shakerinejad E, Kayhani MH, Nazari M, Nazari M, Tamayol A. Increasing the performance of gas diffusion layer by insertion of small hydrophilic layer in proton-exchange membrane fuel cells. Int J Hydrogen Energy 2018;43(4):2410e28. [32] Yu J, Froning D, Reimer U, Lehnert W. Liquid water breakthrough location distances on a gas diffusion layer of polymer electrolyte membrane fuel cells. J Power Sources 2018;389:56e60. [33] Sakaida S, Tabe Y, Chikahisa T, Tanaka K, Konno M. Study on PEFC gas diffusion layer with designed wettability pattern tolerant to flooding. ECS Trans 2018;86(13):111e8. [34] Jinuntuya F, Whiteley M, Chen R, Fly A. The effects of gas diffusion layers structure on water transportation using Xray computed tomography based Lattice Boltzmann method. J Power Sources 2018;378:53e65. [35] Sadeghi R, Shadloo MS, Hopp-Hirschler M, Hadjadj A, Nieken U. Three-dimensional lattice Boltzmann simulations of high density ratio two-phase flows in porous media. Comput Math Appl 2018;75(7):2445e65.

26515

[36] Liu J, Shin S, Um S. Comprehensive statistical analysis of heterogeneous transport characteristics in multifunctional porous gas diffusion layers using lattice Boltzmann method for fuel cell applications. Renew Energy 2019;139:279e91. [37] Lee KJ, Kang JH, Nam JH. Liquid water distribution in hydrophobic gas-diffusion layers with interconnect rib geometry: an invasion-percolation pore network analysis. Int J Hydrogen Energy 2014;39(12):6646e56. [38] Straubhaar B, Pauchet J, Prat M. Pore network modelling of condensation in gas diffusion layers of proton exchange membrane fuel cells. Int J Heat Mass Transf 2016;102:891e901. [39] Fazeli M, Hinebaugh J, Bazylak A. Investigating inlet condition effects on PEMFC GDL liquid water transport through pore network modeling. J Electrochem Soc 2015;162(7):F661e8. [40] Huang X, He Y, Zhou W, Deng D, Zhao Y. Pore network modeling of fibrous porous media of uniform and gradient porosity. Powder Technol 2019;343:350e61. [41] Chen L, He YL, Tao WQ. Effects of surface microstructures of gas diffusion layer on water droplet dynamic behaviors in a micro gas channel of proton exchange membrane fuel cells. Int J Heat Mass Transf 2013;60:252e62. [42] Ahmad ZY, Didari S, Moon J, Harris TAL. Computational fluid dynamics of water droplet formation and detachment from gas diffusion layer. ECS Trans 2013;45(23):89e100. [43] Yin Y, Wu T, He P, Du Q, Jiao K. Numerical simulation of twophase cross flow in microstructure of gas diffusion layer with variable contact angle. Int J Hydrogen Energy 2014;39(28):15772e85. [44] Niu Z, Jiao K, Wang Y, Du Q, Yan Y. Numerical simulation of two-phase cross flow in the gas diffusion layer microstructure of proton exchange membrane fuel cells. Int J Energy Res 2018;42(2):802e16. [45] Niu Z, Bao Z, Wu J, Wang Y, Jiao K. Two-phase flow in the mixed-wettability gas diffusion layer of proton exchange membrane fuel cells. Appl Energy 2018;232:443e50. [46] Jeon DH, Kim H. Effect of compression on water transport in gas diffusion layer of polymer electrolyte membrane fuel cell using lattice Boltzmann method. J Power Sources 2015;294:393e405. [47] Molaeimanesh GR, Shojaeefard MH, Moqaddari MR. Effects of electrode compression on the water droplet removal from proton exchange membrane fuel cells. Korean J Chem Eng 2019;36(1):136e45. [48] Satjaritanun P, Hirano S, Shum AD, Zenyuk IV, Weber AZ, Weidner JW, Shimpalee S. Fundamental understanding of water movement in gas diffusion layer under different arrangements using combination of direct modeling and experimental visualization. J Electrochem Soc 2018;165(13):F1115e26. € tzke C, Lehnert W, [49] Fazeli M, Hinebaugh J, Fishman Z, To Manke I, Bazylak A. Pore network modeling to explore the effects of compression on multiphase transport in polymer electrolyte membrane fuel cell gas diffusion layers. J Power Sources 2016;335:162e71. [50] Bao N, Zhou Y, Jiao K, Yan Y, Du Q, Chen J. Effect of gas diffusion layer deformation on liquid water transport in proton exchange membrane fuel cell. Engineering applications of computational fluid mechanics 2014;8(1):26e43. [51] Han C, Chen Z. Numerical simulations of two-phase flow in a proton-exchange membrane fuel cell based on the generalized design method. Energy Sources, Part A Recovery, Util Environ Eff 2019;41(10):1253e71. [52] Molaeimanesh GR, Googarchin HS, Moqaddam AQ. Lattice Boltzmann simulation of proton exchange membrane fuel

26516

[53]

[54]

[55]

[56]

[57]

[58]

[59]

[60]

[61] [62]

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 4 ( 2 0 1 9 ) 2 6 4 9 8 e2 6 5 1 6

cellseA review on opportunities and challenges. Int J Hydrogen Energy 2016;41(47):22221e45. Dafalla AM, Jiang F. Stresses and their impacts on proton exchange membrane fuel cells: a review. Int J Hydrogen Energy 2018;43(4):2327e48. Chen L, Luan HB, He YL, Tao WQ. Pore-scale flow and mass transport in gas diffusion layer of proton exchange membrane fuel cell with interdigitated flow fields. Int J Therm Sci 2012;51:132e44. Zamel N, Li X, Shen J. Numerical estimation of the effective electrical conductivity in carbon paper diffusion media. Appl Energy 2012;93:39e44. Turkmen AC, Celik C. The effect of different gas diffusion layer porosity on proton exchange membrane fuel cells. Fuel 2018;222:465e74. Carral C, Mele P. A numerical analysis of PEMFC stack assembly through a 3D finite element model. Int J Hydrogen Energy 2014;39(9):4516e30. Tang Y, Karlsson AM, Santare MH, Gilbert M, Cleghorn S, Johnson WB. An experimental investigation of humidity and temperature effects on the mechanical properties of perfluorosulfonic acid membrane. Mater Sci Eng A 2006;425(1e2):297e304. Mason TJ, Millichamp J, Neville TP, El-kharouf A, Pollet BG, Brett DJ. Effect of clamping pressure on ohmic resistance and compression of gas diffusion layers for polymer electrolyte fuel cells. J Power Sources 2012;219:52e9. Molaeimanesh GR, Nazemian M. Investigation of GDL compression effects on the performance of a PEM fuel cell cathode by lattice Boltzmann method. J Power Sources 2017;359:494e506. Korson L, Drost-Hansen W, Millero FJ. Viscosity of water at various temperatures. J Phys Chem 1969;73(1):34e9. Vargaftik NB, Volkov BN, Voljak LD. International tables of the surface tension of water. J Phys Chem Ref Data 1983;12(3):817e20.

[63] Nguyen TV. A gas distributor design for proton-exchangemembrane fuel cells. J Electrochem Soc 1996;143(5):L103e5. [64] Movahedi M, Ramiar A, Ranjber AA. 3D numerical investigation of clamping pressure effect on the performance of proton exchange membrane fuel cell with interdigitated flow field. Energy 2018;142:617e32. [65] Niu Z, Wu J, Bao Z, Wang Y, Yan Y, Jiao K. Two-phase flow and oxygen transport in the perforated gas diffusion layer of proton exchange membrane fuel cell. Int J Heat Mass Transf 2019;139:58e68. [66] Li S, Sunden B. Effects of gas diffusion layer deformation on the transport phenomena and performance of pem fuel cells with interdigitated flow fields. Int J Hydrogen Energy 2018;43(33):16279e92. [67] Fadzillah DM, Rosli MI, Talib MZM, Kamarudin SK, Daud WRW. Review on microstructure modelling of a gas diffusion layer for proton exchange membrane fuel cells. Renew Sustain Energy Rev 2017;77:1001e9. [68] Hwang CM, Ishida M, Ito H, Maeda T, Nakano A, Hasegawa Y, Yoshida T. Influence of properties of gas diffusion layers on the performance of polymer electrolyte-based unitized reversible fuel cells. Int J Hydrogen Energy 2011;36(2):1740e53. [69] Jiao K, Li XG. Water transport in polymer electrolyte membrane fuel cells. Prog Energy Combust Sci 2011;37(3):221e91. [70] El-Kharouf A, Mason TJ, Brett DJL, Pollet BG. Ex-situ characterisation of gas diffusion layers for proton exchange membrane fuel cells. J Power Sources 2012;218:393e404. [71] Omrani R, Shabani B. Review of gas diffusion layer for proton exchange membrane-based technologies with a focus on unitised regenerative fuel cells. Int J Hydrogen Energy 2019;44(7):3834e60.