Influences of the perforation on effective transport properties of gas diffusion layers

International Journal of Heat and Mass Transfer 126 (2018) 243–255

Contents lists available at ScienceDirect

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

Influences of the perforation on effective transport properties of gas diffusion layers Wen-Zhen Fang a, Yu-Qing Tang a, Li Chen a, Qin-Jun Kang b, Wen-Quan Tao a,⇑ a b

Key Laboratory of Thermo-Fluid Science and Engineering, MOE, Xi’an Jiaotong University, Xi’an, China Computational Earth Science Group (EES-16), Los Alamos National Laboratory, Los Alamos, NM, USA

a r t i c l e

i n f o

Article history: Received 29 September 2017 Received in revised form 1 May 2018 Accepted 2 May 2018

Keywords: Gas diffusion layer Effective transport property Perforation Lattice Boltzmann method Saturation Multiple-relaxation-time

a b s t r a c t In this paper, the through-plane and in-plane effective transport properties, including permeability, diffusivity and thermal conductivity, of the perforated gas diffusion layer (GDL) are predicted using multiple-relaxation-time (MRT) lattice Boltzmann method (LBM) based on stochastic reconstructed microstructures. When predicting effective thermal conductivities of GDL, the effect of anisotropic conductive property of fibers is considered. The effective transport properties of dry perforated GDL are fitted as a function of perforation diameter and porosity. It is found that the permeability and effective diffusivity of GDL increase with perforation diameter and porosity while the effective thermal conductivity decreases. The two-phase LBM is adopted to simulate water distributions in perforated GDLs, and dependences of effective transport properties on saturation are then obtained. The results show that: the existence of the perforation significantly affects the water transport in hydrophobic perforated GDLs if its diameter is larger than the average pore size of GDL. The effective permeability and diffusivity of GDL decrease while effective thermal conductivity increases with saturation. The effective transport properties of perforated GDLs change less significantly with saturation than those of non-perforated GDL if the water droplet intruding into the perforation is displaced, while change more rapidly with saturation if the water droplet remains inside the perforation. Ó 2018 Published by Elsevier Ltd.

1. Introduction The proton exchange membrane fuel cell (PEMFC) is a promising power source for automotive and stationary applications due to its high power densities apart from its other advantages. The gas diffusion layer (GDL) plays an important role on the reactant gas transport, water and thermal management in PEMFCs [1–3]. The most widely adopted GDL material is the carbon paper which is composed of fibers randomly distributed in the plane, resulting in the anisotropy of GDLs with different transport properties in the in-plane and through-plane directions [4,5]. In recent years, the effective transport properties of GDL, including gas permeability [6–9], effective diffusivity [10–13] and thermal conductivity [14–17], are widely investigated by experimental measurements or numerical simulations. The presence of liquid water in GDLs impedes reactant gas transport to the active sites, resulting in a significant loss in the performance of PEMFCs, especially at the high current density. Therefore, the dependence of effective

⇑ Corresponding author. E-mail address: [email protected] (W.-Q. Tao). https://doi.org/10.1016/j.ijheatmasstransfer.2018.05.016 0017-9310/Ó 2018 Published by Elsevier Ltd.

transport properties on the saturation of GDL have also been a focus and widely studied by many researchers [10,13,18–20]. To avoid the water flooding problem in the PEMFC, Gerteisen et al. [21,22] first proposed a perforated GDL in which the laserperforation (80 lm in diameter) is made in the through-plane direction. The test cell with the perforated GDL shows a higher current density than the cell with original GDL. After that, many scholars started to investigate influences of the perforation on water transport in GDL. Manahan et al. [23,24] conducted neutron radiography experiments to observe the changes in liquid water behaviors due to the perforation. Markötter et al. [25] performed the in-situ X-ray radiography experiments to observe the water distribution in perforated GDL. Subsequently, his cooperators Haubmann et al. [26] verified the improvement of water transport in perforated GDL and achieved the highest performance with the perforation diameter of 60 lm among three different diameters (30 lm, 60 lm and 150 lm). Okuhata et al. [27] also investigated the effect of perforation structure on liquid water transport by optically visualized technique. However, due to the high costs of experiments, the above studies only focused on some kinds of perforated GDLs, and could not obtain the correlations of effective transport properties with the perforation diameter and porosity.

244

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

water distribution in perforated GDL. Finally, the effective transport properties as a function of porosity, perforation diameter and water saturation are obtained. To reveal the effects of perforation, comparisons of effective transport properties of perforated GDL and non-perforated GDL are also made in this paper. This paper is organized as follows. In Section 2, the LBM including single-phase and two-phase flows and transport models with specified boundary conditions are introduced. In Section 3, the accuracy of the models is verified. In Section 4.1, the effective transport properties of dry perforated GDL are obtained. The water droplet distributions are simulated by two-phase LB model in Section 4.2, and the effective transport properties of partially saturated perforated GDL are obtained in Section 4.3. Finally, some conclusions are drawn in Section 5.

In this aspect, numerical simulations have a great advantage over the experimental studies. During the past few decades, the lattice Boltzmann method (LBM) has developed as an alternative and powerful technique to simulate fluid flows and transport processes, especially in applications involving complex geometries and interfacial dynamics [28–30]. The LBM allows to simulate pore-scale single-phase flow [7,8], mass or heat transports [17,31], and two-phase (water-gas) flow [32–35] in GDL considering the real structures of GDL. It is reported that with the singlerelaxation-time LBM to deal with fluid-solid boundaries, the predicted effective transport coefficients are viscosity dependent. Then, the multiple-relaxation-time (MRT) LBM is developed to overcome this drawback and improve the numerical stability [8,36]. In most previous studies, when predicting the effective thermal conductivity of GDL, the solid fibers are assumed to be isotropic [14–17]. However, the fiber is anisotropic with different thermal conductivities in axial and radius directions. A D3Q7 MRT LB model with off-diagonal components in the relaxation time matrix was developed to predict the effective thermal conductivity of heterogeneous materials with anisotropic components [37,38]. However, as indicated in Ref. [30], a reduced LB model (D3Q7) is insufficient for the complex porous medium, and therefore a D3Q19 model is extended in this paper. Although many empirical correlations of permeability, effective diffusivity and thermal conductivities of GDLs are proposed, they usually ignore the effects of real structure and water droplet distributions on the estimations of effective transport properties, which further reduces the accuracy of PEMFC modeling. As for the perforated GDL, there exist few correlations for predicting effective transport properties. In the present paper, the effective transport properties of perforated GDL, including permeability, diffusivity, thermal conductivity, are predicted by MRT D3Q19 LBM accounting for the microstructure of GDL. Based on a large number of predicted data, the correlations of effective transport properties of dry perforated GDL are fitted as a function of porosity and perforation diameter. Then, the two-phase LBM is adopted to simulate the

2

1;

1;

1;

1;

1;

1;

1;

1;

6 30; 11; 11; 11; 11; 11; 11; 8; 6 6 6 12; 4; 4; 4; 4; 4; 4; 1; 6 6 0; 1; 1; 0; 0; 0; 0; 1; 6 6 6 0; 4; 4; 0; 0; 0; 0; 1; 6 6 0; 0; 0; 1; 1; 0; 0; 1; 6 6 6 0; 0; 0; 4; 4; 0; 0; 1; 6 6 0; 0; 0; 0; 0; 1; 1; 0; 6 6 6 0; 0; 0; 0; 0; 4; 4; 0; 6 6 2; 2; 1; 1; 1; 1; 1; M ¼ 6 0; 6 6 0; 4; 4; 2; 2; 2; 2; 1; 6 6 0; 0; 1; 1; 1; 1; 1; 6 0; 6 6 0; 0; 0; 2; 2; 2; 2; 1; 6 6 6 0; 0; 0; 0; 0; 0; 0; 1; 6 6 0; 0; 0; 0; 0; 0; 0; 0; 6 6 6 0; 0; 0; 0; 0; 0; 0; 0; 6 6 0; 0; 0; 0; 0; 0; 0; 1; 6 6 4 0; 0; 0; 0; 0; 0; 0; 1; 0;

0;

0;

0;

0;

0;

0;

0;

2. Lattice Boltzmann method 2.1. MRT lattice Boltzmann method In the present paper, a D3Q19 MRT LB model is adopted to simulate the gas flow, diffusion and heat transport in the GDL. The evolution equation of the particle distribution function is [31]:

f i ðx þ ei dt; t þ dtÞ  f i ðx; tÞ ¼ ðM1 SMÞij ½f j  f j  eq

where fi is the density distribution function; x denotes the position; t is the time; dt is the time step, ei is the discrete velocity along the ith direction:

2

3 0; 1; 1; 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 1; 0; 0; 0; 0 6 7 ½ei  ¼ 4 0; 0; 0; 1; 1; 0; 0; 1; 1; 1; 1; 0; 0; 0; 0; 1; 1; 1; 1 5 0; 0; 0; 0; 0; 1; 1; 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 1

ð2Þ

In Eq. (1), M is an orthogonal transformation matrix, defined as:

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

8;

8;

8;

8;

8;

8;

8;

8;

8;

8;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

0;

0;

0;

1; 1; 1; 1; 1; 1;

1; 0;

1; 0;

1; 0;

1; 0;

0; 1;

0; 1;

0; 1;

0;

0;

1;

ð1Þ

1; 1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

0;

0;

1;

1;

1; 1;

1;

1;

1; 1;

1;

1;

1;

1;

2; 2; 2;

1;

1;

1;

2; 2; 2;

0;

0;

0;

0;

0;

0;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1;

1; 1; 1; 1;

0;

0;

0;

1; 1;

1; 1; 1; 1; 0; 0; 0; 0;

0; 0;

0; 0;

0; 0;

1; 1; 1; 1; 0;

0;

0;

0;

0;

0;

0;

1;

1;

1;

1;

1;

0;

0;

0;

1;

1;

1;

0;

0;

0;

1;

1;

1;

0;

0;

0;

0;

1;

1;

1;

0;

0;

0;

1;

1;

1; 1; 1; 1;

1;

1; 1;

0;

0;

1; 1;

0;

1; 1;

1

3

8 7 7 7 1 7 7 0 7 7 7 0 7 7 1 7 7 7 1 7 7 1 7 7 7 1 7 7 7 2 7 7 2 7 7 7 0 7 7 0 7 7 7 0 7 7 1 7 7 7 0 7 7 0 7 7 7 1 5 1

ð3Þ

245

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255 eq

and S is the relaxation time coefficient matrix; f i is the equilibrium distribution function. Their definitions are described below in detail for different LB models.

where dij is the Kronecker symbol and Dij is the diffusivity matrix. The local macroscopic quantity of the scalar variable and its partial derivative can be obtained by:

2.1.1. Fluid flow LB model The simulation of gas flow in GDL is performed by the LBM. For fluid flow LB model, the equilibrium distribution function is defined as:

/¼

"

eq fi

ei  u ðei  uÞ2 u2 ¼ xi q 1 þ 2 þ  2 cs 2c4s 2cs

#

ð4Þ

pffiffiffi where x0 ¼ 1=3, x16 ¼ 1=18, x718 ¼ 1=36, and cs ¼ 1= 3. The relaxation time coefficient matrix, S, is given by [36]:

S ¼ diagðs0 ; s1 ; . . . ; s17 ; s18 Þ

ð5Þ

X @/ 1 X f i ; sij ¼ 2 f ei @xj cs dx i i i

ð12Þ

By means of asymptotic analysis [39], the recovered macroscopic diffusion equation is:

  @/ @ @/ Dij ¼ @t @xi @xj

ð13Þ

If the mass or heat diffusion is isotropic,

sxx ¼ syy ¼ szz and

sij ¼ 0 ði–jÞ. 2.2. Multicomponent multiphase LB model

with

s0 ¼ s3 ¼ s5 ¼ s7 ¼ 0; s1 ¼ s2 ¼ s915 ¼ ¼ s1618 ¼ 8

1

s

; s4 ¼ s6 ¼ s8

2s  1 8s  1

ð6Þ

where s is related to the kinematic viscosity:

s¼

m c2s dt

þ 0:5

ð7Þ

The multicomponent multiphase model proposed by Shan and Chen (SC model) is adopted to simulate the two-phase flow in GDL. The SC model introduces k distribution functions for a fluid mixture composed of k components. Each distribution function represents a component and satisfies the evolution equation. The evolution equation can be expressed as [40,41]: k

k

f i ðx þ ei dt; t þ dtÞ  f i ðx; tÞ ¼

The density and velocity are determined by:

q¼

X

f i ; qu ¼

i

X f i ei

2.1.2. Mass and heat transport LB model The mass diffusion and heat transfer simulations in GDL are performed to determine the effective diffusivity and thermal conductivity of GDL. For the mass or heat transport LB model, the equilibrium distribution function is: eq

sk

ð14Þ

k;eq

ð8Þ

i

f i ¼ xi /

i 1 h k;eq k f i ðx; tÞ  f i ðx; tÞ

ð9Þ

where / represents the concentration or temperature. The relaxation time coefficient matrix has 18  18 elements and the nonzero elements are shown as follows:

where f i is the equilibrium distribution function of kth component, and has the same expression as Eq. (4), in which the equilibrium velocity is shifted by considering the interaction between multiple components:

u ¼ u0 þ

sk Fk qk

ð15Þ

where u0 is the common velocity of all fluid components, defined as:

u0 ¼

X q uk k k

sk

,

Xq

k

k

sk

ð16Þ

In Eq. (15), Fk is the total force acting on the kth component including fluid-fluid interaction, fluid-solid adhesion and body force. The fluid-fluid interaction can be obtained by:

XX Fkf ¼ wk ðqk ðxÞÞ Gkk ðx; x0 Þwk ðqk ðxÞÞðx0  xÞ x0

ð17Þ

k

where wk ðqk ðxÞÞ is the potential function (or effective density) which is used to account for the forces between neighboring lattices and often take the value of local density in the multiphase simulation; G is defined as:

8 jx  x0 j ¼ 1; > < g kk pffiffiffi 0 Gkk ðx; x Þ ¼ g kk =2; jx  x0 j ¼ 2 > : 0; otherwise

ð18Þ

where g kk controls the interaction strength between fluids, and will be illustrated in the discussion later.The fluid-solid adhesion can be described as:

X Fks ¼ wk ðqk ðxÞÞ g kw sðx0 Þðx0  xÞ

ð19Þ

x0

ð10Þ

sij is related to the macroscopic mass or heat diffusivity:

sij ¼

Dij 1 þ dij c2s dt 2

ð11Þ

where sðx0 Þ is the indicator function which equals 0 for fluids and 1 for solids; g kw controls the interaction strength between fluid and the solid. The interaction strength g kw is positive for the nonwetting surface while negative for the wetting surface. The body force can be simply introduced by:

Fkb ¼ qk g while g is the body force per unit mass.

ð20Þ

246

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

Fig. 1. Microstructure of the carbon paper GDL.

Considering the inter-particle force between neighboring nodes, the relation between the pressure and density is:

p ¼ c2s

X k

qk þ 3

X g kk wk wk

ð21Þ

k;k

the continuities of temperature and heat flux can be satisfied by assuming that [37,38]:

ðqcp Þf ¼ ðqcp Þs

where subscripts f and s represent the fluid and solid, respectively. The permeability of GDL can be obtained by Darcy’s law:

2.3. Computational domain and reconstructed microstructure In the present paper, a 120  120  120 grid system is adopted and each grid step represents 2 lm. The carbon paper GDL is composed of randomly distributed carbon fibers in the in-plane direction, resulting in the different transport properties between in-plane and through-plane directions. In the present paper, the reconstructed non-perforated GDL is assumed to be the stack of several carbon fiber layers, and each layer is composed of straight fibers with a fixed diameter of 8 lm randomly distributed in the plane. Therefore, there exist contacting points between two fiber layers which are favor of heat removal. The reconstructed GDL with a porosity of 0.72 is shown in Fig. 1(a). The gaps between fibers are filled with the air. To avoid water flooding, a perforated GDL is adopted to enhance liquid water transport performance. The diameter of perforation is an important parameter to control the liquid invasion process. Several values of perforation diameters are chosen to investigate their influences on water transport performance. Fig. 1(b) shows the perforated GDL with a perforation diameter of 80 lm. 2.4. Boundary conditions and calculations of effective transport properties In the present paper, the in-plane and through-plane permeability, effective diffusivity and thermal conductivity of GDL are predicted by the single-phase LB model (introduced in Sections 2.1.1 and 2.1.2). To predict the effective transport properties of GDLs in one certain direction, gradients of pressure, concentration, temperature are applied at the two end surfaces of GDL in that direction while other surfaces are assumed to be periodic. The solid fibers in the computational domain are treated to be impermeable with no slip and no flux boundary conditions on their surfaces when predicting permeability and effective diffusivity of GDL, while the continuities of heat flux and temperature should be satisfied at fluid-solid interfaces when predicting the effective thermal conductivity of GDL. With the ‘half lattice division scheme’,

ð22Þ

K¼

lu rp

ð23Þ

where l is the dynamic viscosity; u is the average superficial velocity and rp is the pressure gradient between the top and bottom surfaces. The effective diffusivity of GDL can be obtained by Fick’s law:

Deff ¼ 

j

rc

ð24Þ

where j is the average diffusion flux and rc is the concentration gradient between the top and bottom surfaces. The effective thermal conductivity of GDL can be obtained by Fourier’s law:

keff ¼ 

q

rT

ð25Þ

where q is the average heat flux and rT is the concentration gradient between the top and bottom surfaces. The two-phase LB model (introduced in Section 2.2) is adopted to obtain the water distribution in GDL. In the present paper, two kinds of boundary conditions are applied to investigate water transport behaviors in GDL. The pressure gradient boundary condition [32] is adopted to investigate the water intrusion process in GDL, while a body force is used as the driving force [35] to obtain the water distribution in GDL. It is because the capillary pressure gradient results in a saturation gradient along the flow direction, especially the water flooding problem in the inlet. 3. Validation The accuracy of the fluid flow LB model has been validated in our previous work [30], and is not presented in this paper again for brevity. In this Section, accuracies of the D3Q19 anisotropic heat transfer LB model (introduced in Section 2.1.2) and two-phase LB model (introduced in Section 2.2) are validated.

247

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

(T-Tmin ) /(Tmax-Tmin )

1.0

LBM Fluent

0.8 0.6 0.4 0.2 0.0 0.0

0.2

0.4

z

0.6

0.8

0.12

0.14

1.0

Fig. 2. Comparison of temperature distributions.

0.12 0.10

LBM simulation linear fit

Δp

0.08 0.06 0.04 0.02 0.00 0.00

0.02

0.04

0.06 0.08 1/R

0.10

Fig. 3. Bubble test.

3.1. Anisotropic heat transfer LB model validation The effective thermal conductivity of a layered dual-component composite material is evaluated with thermal conductivity matrices:

2

4 1

6 k1ij ¼ 4 1 8

2

3

7 1 5; 2 1 10

2

10 1

6 k2ij ¼ 4 1 4

4

3

7 1 5 1 12 8

ð26Þ

The temperature contours and temperature distribution on the vertical middle line obtained by the present LB model are compared with that obtained by Ansys Fluent 14.0, as shown in Fig. 2. It can be seen that the temperature contours from two computations are in good agreement. Furthermore, predicted effective thermal conductivities of the composite material along z direction with the grid 30  30  30 are 10.41 W/(mK) and 10.34 W/(mK), respectively. The deviation is only 0.68%, which validates the accuracy of the present model. 3.2. Two-phase LB model validation The bubble test and static droplet test are conducted to determine the fluid-fluid interaction parameter and fluid-solid

parameter in the two-phase SC model. The bubble test consists of a spherical bubble located at the center of an 80  80  80 grid system. The value of g kk controls the fluid-fluid interfacial tension. In SC model, two immiscible fluids separate with each other only if the value of g kk is larger than a critical value [42]. We set q1 ¼ 2, q2 ¼ 0:01 inside the bubble with an initial radius of 20 while q1 ¼ 0:01, q2 ¼ 2 outside the bubble. After a series of tests, g kk is determined to be 0.08, and the steady state bubble is show in Fig. 3(a). The pressure jump across the interface is given by the Laplace’s law:

pi  po ¼ r=R

ð27Þ

where pi and po represent the pressure inside and outside the bubble, respectively. r is the surface tension. Fig. 3(b) shows the pressure drop as a function of the inverse radius. It can be seen that the pressure drop is proportional to the inverse radius, which agrees well with Laplace’s law. The three-phase contact angle of the static droplet on the wall is determined by the value of g kw (k equals 1 or 2 to represent liquid and gas, respectively, and w stands for the wall). The value of g 1w is varied to obtain the static droplet with different contact angle while maintains g 1w ¼ g 2w . The steady state of a droplet with

248

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

140

Contact angle/

120 100 80 60 40 -0.06

-0.04

-0.02

0.00

g1w

0.02

0.04

0.06

g1w Fig. 4. The static droplet test.

g 1w to be 0.03 is shown in Fig. 4(a), while the static contact angle varying with g 1w is shown in Fig. 4(b).

4.1.1. Diffusivity The effective diffusivity of non-perforated GDL is a function of porosity:

Deff ¼ D0 f ðeÞ

4. Results and discussion In the following, the effective transport properties of dry GDL will be presented first, followed by the discussion on wetted GDL within which water droplet drainage happens, with focus being put on the determination of the effective transport properties of partially saturated GDL.

ð28Þ

A famous empirical correlation was proposed to calculate the effective diffusivity of the non-perforated GDL [10]:

 Deff ¼ D0 e

e  ep 1  ep

a ð29Þ

where ep is a percolation threshold, and a is an empirical parameter. ep was found to be 0.11; while a was found to be 0.521 and 0.785 for 4.1. Effective transport properties of dry GDL The in-plane and through plane effective transport properties of the dry perforated GDL are determined as a function of porosity and perforation diameters. The porosity varies from 0.58 to 0.83; while the diameter of perforations varies from 0 to 96 lm.

in-plane and through-plane direction, respectively.In the following, the effective diffusivity of dry perforated GDL is determined by LBM introduced in Section 2.1.2. The predicted through-plane and in-plane effective diffusivities varying with the porosity and perforation diameter are shown in Fig. 5(a) and (b), respectively. In the figure, the scatter data are the results of LBM, while the curved

Fig. 5. Effective diffusivity of the perforated GDL.

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

249

surfaces are the fitting results. The formula of the fitting function is given by:

  e  ep a c Deff ¼ D0 e ðbðd=dav e Þ þ 1Þ 1  ep

ð30Þ

where a, b, c, ep are the fitting parameters; d is the diameter of the perforation in GDL and dave is the average diameter of the pore size of GDL, calculated as:

dav e ¼

pffiffiffi e pffiffiffi df 1 e

ð31Þ

where df is the diameter of fibers. Compared with Eq. (29), the modified term in Eq. (30) represents the contribution of the perforation on the effective diffusivity of GDL. The fitting process is conducted by First Optimization (1stOpt) Software. First, the value of ep is determined by fitting the predicted results of non-perforated GDL with Eq. (29) (the value of d in Eq. (30) is set to be zero). With the given value of ep, we finally determine the values of a, b, c. The fitting results show that ep = 0.25, while a = 0.855, b =0.049, c =1.524 and a = 0.443, b =0.046, c =0.943 for the through-plane and in-plane directions, respectively. The fitting values of ep and a different from Eq. (29) may result from the significant microstructure difference among GDLs. It can be seen that the perforation significantly enhances the through-plane and in-plane effective diffusivities of the GDL. The increase tendency in the throughplane direction is larger than that in the in-plane direction. 4.1.2. Permeability The widely adopted empirical correlation to predict permeability is the Kozeny-Carman equation [10]: 2

K=df ¼

e3 16kK ð1  eÞ2

ð32Þ

where kK is the KC constant (usually, 0.5 < e < 0.7), and it is often taken as 0.6. The gas permeability of dry perforated GDL is determined by LBM introduced in Section 2.1.1. The predicted throughplane and in-plane gas permeability varying with the porosity and perforation diameter is shown in Fig. 6(a) and (b), respectively. It can be seen that both through-plane and in-plane permeability increase with the porosity and perforation diameter. However, the through-plane permeability increases more intensely with perforation diameter than that of in-plane permeability. It is because the perforation is made in the through-plane direction. Due to

Fig. 7. Effective thermal conductivity of the non-perforated GDL.

different variation tendencies of the through-plane and in-plane directions, their corresponding fitting functions are 2

K=df ¼

e3 c þ beðd=df Þ 16ð1 þ eÞa ð1  eÞ2

ð33Þ

e3 c ðbðd=dL Þ þ 1Þ 16ð1 þ eÞa ð1  eÞ2

ð34Þ

and 2

K=df ¼

respectively. The fitting parameters are a = 3.027, b = 3.767  104 and c = 4 for through-plane direction while a = 2.513, b = 2.698 and c = 2.102 for in-plane direction. Compared with the KC equation, the modified terms in Eqs. (33) and (34) represent the effect of the perforation. Note that kK is no longer a constant but a function of porosity, i.e. ð1 þ eÞa , which may result from the wider range of porosity studied and significant microstructure difference among GDLs. 4.1.3. Thermal conductivity The gas thermal conductivity kg is set to be 0.03 W/(mK) in this paper. The fibers in GDL are anisotropic with different thermal conductivities in axial and radial directions. The reported axial thermal conductivity of fiber is 120 W/(mK). However, the radial thermal conductivity of fiber is seldom reported. In the present paper, the variation range of the radial thermal conductivity is taken from

Fig. 6. Permeability of the perforated GDL.

250

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

Fig. 8. Effective thermal conductivity of the perforated GDL.

12 to 120 W/(mK), and the effective thermal conductivity of nonperforated GDL is obtained as a function of the radial thermal conductivity of fiber. Fig. 7 shows the ratio of kTe =kg of non-perforated GDL varying with porosity and radial thermal conductivity of fibers.In Fig. 7, the curved surface is fitted by the following function:

kTe g þ ðn  1Þ þ ðn  1Þðg  1Þð1  eÞ ¼ þ agb ð1  eÞc kg g þ ðn  1Þ þ ð1  gÞð1  eÞ

ð35Þ

where g ¼ krf =kg is the ratio of radial thermal conductivity of fibers and gas thermal conductivity; a, b, c are the fitting parameters, and a = 7.768, b = 0.347 and c = 1.435. The first term at the right side of Eq. (35) is the famous empirical Hamilton model [43], and n is the shape factor of dispersed phase (n = 6 for fiber). The second term at the right side of Eq. (35) represents the contribution part of connection point of fibers between neighboring layers. Then, the influences of perforation diameter and porosity on the effective thermal conductivity of perforated GDL are studied with radius thermal conductivity of fiber being 12 W/(mK). The predicted through-plane and in-plane effective thermal conductivities varying with the porosity and perforation diameter are shown in Fig. 8(a) and (b), respectively. It can be seen that both throughplane and in-plane effective thermal conductivities decrease with the increase in porosity and perforation diameter, and the effect of the porosity is more significant. The effective thermal conductivity in the in-plane direction is significantly larger than that in the through-plane direction. In Fig. 8(a), the function format of fitting curved surface is expressed as follows:

  kTe g þ ðn  1Þ þ ðn  1Þðg  1Þð1  eÞ ¼ þ agb ð1  eÞc kg g þ ðn  1Þ þ ð1  gÞð1  eÞ b

ð1  a2 ðd=dL Þ 2 Þ

ð36Þ

where the part in the square brace is obtained from the results of non-perforated GDL, and the remaining part represents the effect of perforation. The values of fitting parameter equals 0.745 and 1.912 for a2 and b2 , respectively. Similar to the fitting process of through-plane effective thermal conductivities, the fitting function of in-plane effective thermal conductivities can be obtained as follows:

kIe b ¼ ½1 þ agð1  eÞb ð1  a2 ðd=dL Þ 2 Þ kg

ð37Þ

with a = 0.101, b = 1.455 and a

2

= 1.166, b

2

= 1.789.

4.2. Water transport in GDL As indicated above in the present study, the SC model is adopted to investigate water transport behaviors in GDL. The pressure difference between inlet and outlet is adopted to study the water intrusion process in GDL (primary drainage simulations), while a body force is adopted to obtain the water distribution in GDL (steady state flow simulations). As indicated in previous studies [32,34], the water flow in GDL is dominated by capillary pressure and the effect of density and viscosity ratio of liquid versus gas can be neglected. In the present paper, all the density ratio is around unity and s1 = s2 = 1. 4.2.1. Primary drainage The predicted water intrusion process in the perforated hydrophobic GDL with a contact angle of 113.5° is presented in this section. A constant water pressure and zero gas pressure zone (with five lattice width) are applied at the inlet while a constant gas pressure and zero water pressure (with five lattice width) are applied at the outlet to maintain constant pressure drop. The influences of pressure difference and perforation diameter on the water transport behavior are shown in Fig. 9. The inlet water density qw;in is maintained to be 2.0 while the outlet gas density qg;out is set to be 1.9 and 1.8 for Fig. 9(a) and (b), respectively. The diameters of perforation of the left and right figures are 48 lm and 80 lm, respectively. It can be seen that water invades into GDL more easily when the pressure drop increases. At low pressure difference, the water only invades into GDL at preferential locations due to the resistance of capillary pressure. The capillary resistance is proportional to the inverse diameter of pores, and therefore water can invade a larger size pore more easily. The perforation diameter will have an effect on water transport behavior of GDL if it is larger than the average pore size. It is the reason that the perforation diameter of 80 lm can significantly enhance the water transport process than that of 48 lm. Hereafter, the studied perforation diameter is 80 lm if not mentioned. 4.2.2. Steady state flow In this section, the body force g is set to be 4  106 to replace the fluid pressure gradient, and the periodic boundary conditions are applied at all six directions. Such treatment is also adopted in

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

(1) d = 48

251

(2) d = 80 (a) ρ w,in = 2 , ρg,out = 1.9

(1) d = 48

(2) d = 80 (b) ρ w,in = 2 , ρg,out = 1.8 Fig. 9. Water intrusion in perforated hydrophobic GDL.

Refs. [32,34,35] to obtain the water distribution in GDL. Considering the effect of consideration on the water distribution in GDL remains a tough problem, which is neglected in this paper. An initial random distribution of liquid droplet is placed within the microstructure of GDL until reaches the desired water saturation. Here, the saturation is defined as the ratio of the water volume and void volume. Two immiscible fluids can flow simultaneously and the corresponding saturation and fluid flow rate can be directly calculated if the steady state is reached. The dynamic liquid water distribution in perforated hydrophobic GDL is shown in Fig. 10 with an initial saturation of 6.5%. The static contact angle of hydrophobic GDL is 113.5°. From the yellow dashed circle in Fig. 10, it can be seen that the liquid water will intrude into the 80 lm-diameter perforation due to its smaller capillary resistance. At time t = 8000dt (in Fig. 10(e)), the intruding liquid droplet merges with another liquid droplet to form a larger size droplet. With the existence of perforation, the liquid water droplet close to the perforation will intrude into the perforation and soon goes away, and therefore reduces the water flooding situation inside the GDL.

4.3. Effective transport properties of partially saturated GDL The permeability, effective diffusivity and thermal conductivity of partially saturated GDL are related to the saturation. The water distribution in GDL is obtained by the two-phase SC model. The liquid droplets occupying in GDL pores are treated to be rigid. For predicting effective permeability and diffusivity of GDL, the liquid is assumed to be impermeable solid; and for predicting effective thermal conductivity, it is assumed to be a conductive component with the thermal conductivity of 0.6 W/(mK). In this section, we modify the previous correlations of effective transport properties of dry GDL (introduced in Section 4.1) by multiplying a saturation term to gain the effective transport properties of partially saturated GDL:

Us ðe; d; sÞ ¼ Udry ðe; dÞf ðsÞ

ð38Þ

where U represents effective diffusivity, permeability or thermal conductivity; s is the saturation;f ðsÞ ¼ Us =Udry is a function of saturation, which represents the effect of saturation on effective

252

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

Fig. 10. Dynamic liquid water distributions in perforated hydrophobic GDL.

Fig. 11. Influences of saturation on effective diffusivity of GDL.

transport properties of GDL. The quantitative relations of effective transport properties as the function of saturations will be presented in the following discussions with the porosity of non-perforated

GDL being 0.72. It should be noted that as shown in Fig. 10, the water droplets intruding into the perforation remain in the perforation by reentering into the bottom side surface because of the peri-

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

253

Fig. 12. Influences of saturation on effective permeability of GDL.

Fig. 13. Influences of saturation on effective thermal conductivity of GDL.

odic boundary condition and then merged to be a larger water droplet, resulting in the block of perforation. Actually, water droplets will not reenter into the perforation but be discharged into GDL into gas channel. The block of perforation will significantly reduce the effective diffusivity and permeability. Therefore, we artificially remove the water droplets out of the perforation after the steady state condition is reached and then obtain the effective transport properties of GDL as a function of the initial saturation. 4.3.1. Diffusivity The influences of saturation on the through-plane and in-plane effective diffusivities of GDL are shown in Fig. 11(a) and (b), respectively. In the figure, the discrete points are the simulation results, while the straight lines are the fitting results. The existence of perforation has an impact on the water distribution in perforated GDL, resulting in different variation tendencies with saturation of the non-perforated and perforated GDLs. The format of fitting function is:

f ðsÞ ¼ ð1  sÞa

ð39Þ

where a is the fitting parameter. For through-plane direction, the fitting parameter a are 1.40, 1.58 and 0.51 for non-perforated GDL, perforated GDL and perforated GDL with water discharged from the perforation (hereafter, it is denoted as ‘draining perforated GDL’), respectively; while for the in-plane direction, the fitting parameter a are 1.06, 1.20 and 0.30, respectively. It can be seen that the decrease tendency of effective diffusivities of GDL in the inplane direction is larger than that in through-plane direction when saturation increases. Besides, when the saturation increases, effective diffusivities of draining perforated GDL decrease slowest while those of perforated GDL decrease most rapidly if the water droplets are still inside the perforation (red and circular points in the figures). The water droplets accumulated inside the perforation (before the moment that water droplets are discharged) occupy the space of perforation, resulting in the rapid decrease of the effective diffusivity. On the contrary, if the water droplets in the perforation are discharged, the decrease tendency of effective diffusivity of GDL will reduce compared with that of nonperforated GDL. It is because the discharged water droplets are out of GDL, leaving more free paths for gas to diffuse.

254

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255

4.3.2. Permeability The influences of saturation on the through-plane and in-plane effective permeability of GDL are shown in Fig. 12(a) and (b), respectively. The effective permeability of GDL is a function of saturation:

 f ðsÞ ¼

1s 1 þ bs

a ð40Þ

where a and b are the fitting parameters. In the through-plane direction, a = 0.48, b = 9.74 for non-perforated GDL, a = 0.53, b = 156.47 for perforated GDL, and a = 0.45, b = 1 for draining perforated GDL; while in the in-plane direction, a = 0.48, b = 9.74 for non-perforated GDL, a = 0.53, b = 156.47 for perforated GDL, and a = 0.45, b = 1 for draining perforated GDL. The effective permeability of draining perforated GDL decreases slowest with the increase in saturation while that of perforated GDL decreases most rapidly. The reason is similar to the diffusivity occasions. 4.3.3. Thermal conductivity The influences of saturation on the through-plane and in-plane effective thermal conductivity of GDL are shown in Fig. 13(a) and (b), respectively. The effective thermal conductivity of GDL is a function of saturation:

f ðsÞ ¼ ð1 þ sÞa

ð41Þ

where a is the fitting parameter. For the through-plane direction, the fitting parameter a are 2.84, 3.15 and 1.20 for non-perforated GDL, perforated GDL and draining perforated GDL, respectively; while for in-plane direction, the fitting parameter a are 0.35, 0.52 and 0.21, respectively. On the contrary, the effective thermal conductivity of GDL increases with saturation. Besides, the effective thermal conductivity of draining perforated GDL increases slowest while that of perforated GDL increases most rapidly. The reason is opposite to the diffusivity occasions. 5. Conclusion In the present paper, the through-plane and in-plane effective transport properties of perforated GDLs are estimated by the pore-scale LBM based on the reconstructed microstructure, and the dependences of effective transport properties on saturation are obtained. The major results are as follows: (1) The effective thermal properties of dry perforated GDL can be fitted as a function of porosity and perforation diameter; the permeability and effective diffusivity of GDL increase while the effective thermal conductivity decreases with the increase in perforation diameter and porosity; (2) The effective permeability and diffusivity of partially saturated GDL decrease while thermal conductivity increases with saturation; the perforation can significantly influence the water transport in a hydrophobic perforated GDL if its diameter is larger than the average pore size of the GDL; the effective transport properties of the non-perforated GDL change more significantly with saturation than thoes of perforated GDL with a hole diameter of 80 lm if the water droplets intruding into the perforation are discharged out of GDL while change less rapidly if the water droplets remain inside the perforation. Conflict of interest The authors declare that there is no conflict of interest.

Acknowledgement This research was supported by the 13-5 National Key Research Project of China (2017YFB0102702) and the Key Project of International Joint Research of National Natural Science Foundation of China (51320105004). Q.K. acknowledges the support from LANL’s LDRD Program. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.ijheatmasstransfer. 2018.05.016. References [1] M.H. Shojaeefard, G.R. Molaeimanesh, M. Nazemian, et al., A review on microstructure reconstruction of PEM fuel cells porous electrodes for pore scale simulation, Int. J. Hydrogen Energy 41 (44) (2016) 20276–20293. [2] G.R. Molaeimanesh, H.S. Googarchin, A.Q. Moqaddam, Lattice Boltzmann simulation of proton exchange membrane fuel cells–a review on opportunities and challenges, Int. J. Hydrogen Energy 41 (47) (2016) 22221– 22245. [3] N. Zamel, X. Li, Effective transport properties for polymer electrolyte membrane fuel cells–with a focus on the gas diffusion layer, Prog. Energy Combust. Sci. 39 (2013) 111–146. [4] J. Pharoah, K. Karan, W. Sun, On effective transport coefficients in PEM fuel cell electrodes: anisotropy of the porous transport layers, J. Power Sources 161 (2006) 214–224. [5] A.G. Yiotis, M.E. Kainourgiakis, G.C. Charalambopoulou, A.K. Stubos, Microscale characterisation of stochastically reconstructed carbon fiber-based Gas Diffusion Layers; effects of anisotropy and resin content, J. Power Sources 320 (2016) 153–167. [6] J.T. Gostick, M.W. Fowler, M.D. Pritzker, M.A. Ioannidis, L.M. Behra, In-plane and through-plane gas permeability of carbon fiber electrode backing layers, J. Power Sources 162 (2006) 228–238. [7] L. Hao, P. Cheng, Lattice Boltzmann simulations of anisotropic permeabilities in carbon paper gas diffusion layers, J. Power Sources 186 (2009) 104–114. [8] Y. Gao, X. Zhang, P. Rama, Y. Liu, R. Chen, H. Ostadi, K. Jiang, Modeling Fluid Flow in the Gas Diffusion Layers in PEMFC Using the Multiple Relaxation-time Lattice Boltzmann Method, Fuel Cells 12 (2012) 365–381. [9] A. Tamayol, F. McGregor, M. Bahrami, Single phase through-plane permeability of carbon paper gas diffusion layers, J. Power Sources 204 (2012) 94–99. [10] J.H. Nam, M. Kaviany, Effective diffusivity and water-saturation distribution in single-and two-layer PEMFC diffusion medium, Int. J. Heat Mass Transfer 46 (2003) 4595–4611. [11] L.M. Pant, S.K. Mitra, M. Secanell, Absolute permeability and Knudsen diffusivity measurements in PEMFC gas diffusion layers and micro porous layers, J. Power Sources 206 (2012) 153–160. [12] P.K. Das, X. Li, Z.S. Liu, Effective transport coefficients in PEM fuel cell catalyst and gas diffusion layers: beyond Bruggeman approximation, Appl. Energy 87 (2010) 2785–2796. [13] 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. [14] J. Ramousse, S. Didierjean, O. Lottin, D. Maillet, Estimation of the effective thermal conductivity of carbon felts used as PEMFC Gas Diffusion Layers, Int. J. Therm. Sci. 47 (2008) 1–6. [15] D. Veyret, G. Tsotridis, Numerical determination of the effective thermal conductivity of fibrous materials. Application to proton exchange membrane fuel cell gas diffusion layers, J. Power Sources 195 (2010) 1302–1307. [16] A. Radhakrishnan, Z. Lu, S.G. Kandlikar, Effective thermal conductivity of gas diffusion layers used in PEMFC: measured with guarded-hot-plate method and predicted by a fractal model, ECS Trans. 33 (2010) 1163–1176. [17] J. Yablecki, A. Nabovati, A. Bazylak, Modeling the effective thermal conductivity of an anisotropic gas diffusion layer in a polymer electrolyte membrane fuel cell, J. Electrochem. Soc. 159 (2012) B647–B653. [18] J.T. Gostick, M.A. Ioannidis, M.W. Fowler, M.D. Pritzker, Pore network modeling of fibrous gas diffusion layers for polymer electrolyte membrane fuel cells, J. Power Sources 173 (2007) 277–290. [19] G. Xu, J. LaManna, J. Clement, M. Mench, Direct measurement of through-plane thermal conductivity of partially saturated fuel cell diffusion media, J. Power Sources 256 (2014) 212–219. [20] T. Bednarek, G. Tsotridis, Calculation of effective transport properties of partially saturated gas diffusion layers, J. Power Sources 340 (2017) 111–120. [21] D. Gerteisen, T. Heilmann, C. Ziegler, Enhancing liquid water transport by laser perforation of a GDL in a PEM fuel cell, J. Power Sources 177 (2008) 348–354. [22] D. Gerteisen, C. Sadeler, Stability and performance improvement of a polymer electrolyte membrane fuel cell stack by laser perforation of gas diffusion layers, J. Power Sources 195 (2010) 5252–5257.

W.-Z. Fang et al. / International Journal of Heat and Mass Transfer 126 (2018) 243–255 [23] M. Manahan, M. Hatzell, E. Kumbur, M. Mench, Laser perforated fuel cell diffusion media. Part I: Related changes in performance and water content, J. Power Sources 196 (2011) 5573–5582. [24] M. Manahan, J. Clement, A. Srouji, S. Brown, T. Reutzel, M. Mench, Laser modified fuel cell diffusion media: engineering enhanced performance via localized water redistribution, J. Electrochem. Soc. 161 (2014) F1061–F1069. [25] H. Markötter, R. Alink, J. Haußmann, K. Dittmann, T. Arlt, F. Wieder, C. Tötzke, M. Klages, C. Reiter, H. Riesemeier, Visualization of the water distribution in perforated gas diffusion layers by means of synchrotron X-ray radiography, Int. J. Hydrogen Energy 37 (2012) 7757–7761. [26] J. Haußmann, H. Markötter, R. Alink, A. Bauder, K. Dittmann, I. Manke, J. Scholta, Synchrotron radiography and tomography of water transport in perforated gas diffusion media, J. Power Sources 239 (2013) 611–622. [27] G. Okuhata, T. Tonoike, K. Nishida, S. Tsushima, S. Hirai, Effect of perforation structure of cathode GDL on liquid water removal in PEFC, ECS Trans. 58 (2013) 1047–1057. [28] W.Z. Fang, L. Chen, Q.J. Kang, et al., Lattice Boltzmann modeling of pool boiling with large liquid-gas density ratio, Int. J. Therm. Sci. 114 (2017) 172–183. [29] T. Koido, T. Furusawa, K. Moriyama, An approach to modeling two-phase transport in the gas diffusion layer of a proton exchange membrane fuel cell, J. Power Sources 175 (2008) 127–136. [30] A. Xu, W. Shyy, T.S. Zhao, Lattice Boltzmann modeling of transport phenomena in fuel cells and flow batteries, Acta Mech. Sin. (2017) 1–20. [31] L. Li, C. Chen, R. Mei, J.F. Klausner, Conjugate heat and mass transfer in the lattice Boltzmann equation method, Phys. Rev. E 89 (2014) 043308. [32] P.P. Mukherjee, Q. Kang, C.Y. Wang, Pore-scale modeling of two-phase transport in polymer electrolyte fuel cells—progress and perspective, Energy Environ. Sci. 4 (2011) 346–369.

255

[33] G.R. Molaeimanesh, M.H. Akbari, A three-dimensional pore-scale model of the cathode electrode in polymer-electrolyte membrane fuel cell by lattice Boltzmann method, J. Power Sources 258 (2014) 89–97. [34] P.K. Sinha, P.P. Mukherjee, C.Y. Wang, Impact of GDL structure and wettability on water management in polymer electrolyte fuel cells, J. Mater. Chem. 17 (2007) 3089–3103. [35] H. Li, C. Pan, C.T. Miller, Pore-scale investigation of viscous coupling effects for two-phase flow in porous media, Phys. Rev. E 72 (2005) 026705. [36] C. Pan, L.S. Luo, C.T. Miller, An evaluation of lattice Boltzmann schemes for porous medium flow simulation, Comput. Fluids 35 (2006) 898–909. [37] W.Z. Fang, L. Chen, J.J. Gou, W.Q. Tao, Predictions of effective thermal conductivities for three-dimensional four-directional braided composites using the lattice Boltzmann method, Int. J. Heat Mass Transfer 92 (2016) 120–130. [38] W.Z. Fang, J.J. Gou, H. Zhang, Q. Kang, W.Q. Tao, Numerical predictions of the effective thermal conductivity for needled C/C-SiC composite materials, Numer. Heat Tr. A-Appl. 70 (2016) 1101–1117. [39] H. Yoshida, M. Nagaoka, Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation, J. Comput. Phys. 229 (2010) 7774–7795. [40] Q. Kang, D. Zhang, S. Chen, Displacement of a three-dimensional immiscible droplet in a duct, J. Fluid Mech. 545 (2005) 41–66. [41] Q. Kang, D. Zhang, S. Chen, Displacement of a two-dimensional immiscible droplet in a channel, Phys. Fluids 14 (2002) 3203–3214. [42] H. Huang, D.T. Thorne Jr, M.G. Schaap, M.C. Sukop, Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models, Phys. Rev. E 76 (2007) 066701. [43] M. Wang, N. Pan, Predictions of effective physical properties of complex multiphase materials, Mat. Sci. Eng. R-Rep. 63 (2008) 1–30.