- Email: [email protected]
Impact of PTFE content and distribution on liquid–gas flow in PEMFC carbon paper gas distribution layer: 3D lattice Boltzmann simulations
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Impact of PTFE content and distribution on liquidegas flow in PEMFC carbon paper gas distribution layer: 3D lattice Boltzmann simulations Wang Chen, Fangming Jiang* Laboratory of Advanced Energy Systems, Guangdong Key Laboratory of New and Renewable Energy Research and Development, CAS Key Laboratory of Renewable Energy, Guangzhou Institute of Energy Conversion, Chinese Academy of Sciences (CAS), Guangzhou, 510640, China
article info
abstract
Article history:
A 3D multiphase lattice Boltzmann model (LBM) was established and employed to study
Received 20 November 2015
the impact of polytetrafluoroethylene (PTFE) content and distribution on the liquidegas
Received in revised form
transport in proton exchange membrane fuel cell (PEMFC) carbon-paper gas distribution
2 February 2016
layer (GDL). In the computer-generated carbon paper GDL, part of the fibers’ surface was
Accepted 9 February 2016
randomly specified as PTFE to achieve a desired PTFE content and distribution. For non-
Available online 19 April 2016
uniform PTFE distribution cases, a GDL sub-region neighboring to gas channel (GC) was set to have higher PTFE content. PTFE content and distribution show important influence
Keywords:
on the liquidegas flow and on the relationship of relative permeability versus phase
Proton exchange membrane fuel cell
saturation. The GDL sub-region of more PTFE acts like a capillary barrier and raises the
Gas distribution layer
liquid saturation level in the rest part of GDL. The simulated liquid saturation level in GDL
Lattice Boltzmann method
diminishes with the increase of PTFE content in GDL. The liquid phase relative perme-
Polytetrafluoroethylene
ability is larger in the GDL of lower PTFE content if the (overall) liquid saturation is lower
Liquidegas flow
than a threshold value, whereas it becomes smaller if the liquid saturation is higher than this threshold value. The liquid phase relative permeability in the GDL with non-uniform PTFE distribution is reduced when the liquid saturation in GDL is lower than a certain value, but when the liquid saturation is higher than this value, it is larger compared with that in the GDL with uniform PTFE distribution. The gas phase relative permeability in the GDL of lower PTFE content is larger than tat in the GDL of higher PTFE content and the nonuniform distribution of PTFE in GDL slightly reduces the gas phase relative permeability value. © 2016 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved.
* Corresponding author. Tel.: þ86 20 87057656. E-mail address: [email protected] (F. Jiang). http://dx.doi.org/10.1016/j.ijhydene.2016.02.159 0360-3199/© 2016 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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
Introduction Fuel cells directly convert chemical energy of fuels into electricity and thus have high energy conversion efficiency. Among many types of fuel cells, the H2/air proton exchange membrane fuel cell (PEMFC) may be the most promising to be used as alternative power engine of future vehicles due to advantages like low operating temperature (usually lower than 100 C), zero emission and high power density. Though tremendous progresses on the development of PEMFC engines have been achieved in the past two decades, challenges still remain for large-scale commercialization of passenger PEMFC cars. Further improving the performance (e.g. durability) and reducing the cost of PEMFC may be the most critical challenges. A thorough understanding to water transport in PEMFC is helpful to tackle these challenges. A PEMFC is usually layer-structured, consisting of two gas distribution layers (GDLs), two catalyst layers (CLs), a membrane, two gas flow channels, two bi-polar plates, and some coolant channels. The GDL, as one key component of PEMFC, is located between the gas flow channel and the catalyst layer (CL), and is normally highly porous. The pores in GDL provide pathways for the transport of gaseous reactants and the removal of product water. During practical operations of PEMFC, liquid water may emerge inevitably. Excess liquid water accumulating in GDL may block the pores and hinder the transport paths of reactant gases, a major mechanism limiting the fuel cell high current density operation. A deep knowledge of the multiphase flow in GDL is essential to achieving an increased PEMFC performance. Numerical modeling may be the most efficient method to study the complex multidisciplinary, multiphase transport in PEMFC. To date, a few theoretical models have been established for the simulation of liquidegas flow in PEMFC GDL, e.g., multiphase mixture model (M2) [1e3], two-fluid model (F2) [4e6], and volume of fluid (VOF) [7] model. The M2 model treats the multiphase fluid as a mixture when solving the momentum and continuity equations; it defines the saturation of liquid and considers the interactions between different phases (including phase change) when solving the species transport equations; the relative velocity between phases can be obtained via post processing of the simulated data. The F2 model solves two sets of conservation equations for the gas and liquid phase, respectively, while the interactions between phases are approximated empirically or semi-empirically. Different from the M2 and F2 models, the VOF model solves an additional scalar transport equation for the scalar of liquid volume fraction; it reconstructs the interface of gas and liquid phases every time step and thus explicitly differentiates the interface between phases, which enables more accurate implementation of interactions between phases. However, the VOF method is very computationally intensive, as it generally requires very small time step and very fine mesh to resolve the complicated geometry of liquidegas interface and the complex interactions between phases, which restricts the VOF model being only successfully used for modeling multiphase flows in PEMFC gas channels (GC) [8,9] or GDLs with artificially simplified pore structure [10,11].
8551
Besides the above-reviewed macroscopic models, the lattice Boltzmann (LB) method, based on the mesoscopic kinetic theory, has been employed to simulate the liquidegas flow in GC and porous GDL in recent years [12e15]. Unlike the conventional computational fluid dynamics (CFD) models, the LB model first requires resolving the mesoscopic pore structures in GDL and solves the Boltzmann equation being simplified for a lattice. Therefore, the LB model can give mesoscopic insight into the gaseliquid flow in GDL. The GDL is made of fibrous porous medium, either carbon paper or carbon cloth. To enhance the hydrophobicity of GDL, the polytetrafluoroethylene (PTFE) is usually used to treat the GDL [16,17]. After the treatment, about 5%e40% of PTFE by weight (i.e. 5 wt.% e 40 wt.% PTFE) may be coated on the surfaces of GDL fibers [18]. The PTFE coating, particularly in the cathode GDL, may lose due to hot-pressing and mechanical compression when assembling PEMFC [19], or could be washed off by flow of fluid during operations [18,19], or diminishes during sub-zero temperature cold start cycles [20]. How does the content of PTFE in GDL affect the PEMFC performance? Despite some works (mainly experimental) conducted, no common cognition has been reached. Paganin et al. [21] reported an optimal content of 15 wt.% whereas the result of Lim and Wang [22] was 10 wt.%. Moreover, when doing PTFE treatment to GDL, the PTFE may not totally penetrate the GDL and accumulates toward the GDL surface. The distribution of PTFE in GDL must have influence on the involved liquidegas flow, which has been rarely studied except a recent work by Maleimanesh and Akbari [23]. Their two-dimensional (2D) LB simulations revealed even a small uncoated (by PTFE) region in GDL would impede the removal of liquid water. The relative permeability is an important parameter characterizing the multiphase flow in porous media, defined as the ratio of the effective permeability to the absolute permeability of a fluid phase. In PEMFC, the relative permeability (Krl for liquid phase and Krg for gas phase) reflects the obstruction of one phase acts on the other. It is generally hard to directly measure the Krl and Krg in the GDL, instead, porescale numerical modeling is widely used as an effective tool to determine these parameters. Gostick [24] by pore network model, Niu et al. [25], Hao et al. [26] and Koido et al. [27] by LB models, simulated the two-phase flow in GDL and calculated the relative permeability for each phase. Niu et al. [25] found that the relationship between Krl and liquid saturation sl is approximately amenable to the van Genuchten function [28]. Koido et al. [27] predicted the Kr for each phase by a combination of the single phase and the multicomponent multiphase LB method and they validated their results by comparing the predicted gas relative permeability values against experimental data. The above numerical modeling works to some extent revealed the underlying mechanisms of the two-phase transport in fibrous GDL. However, the impact of PTFE's content and distribution were ignored or not directly simulated in these works. The present work will establish a three-dimensional (3D) LB model for modeling the gaseliquid flow in PEMFC carbon paper GDL. The major purpose is to reveal the effects of PTFE content and distribution on the involved multiphase transport. Specially, the relative permeability, Krl and Krg, will be
8552
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
derived with respect to carbon paper GDLs of various PTFE content and distribution scenarios.
3 9 eq u $! u þ ð! e i $! x i ; tÞ ¼ ui Di þ ui r 3! e i $! u ! u Þ2 : fi ð! 2 2
LB modeling A few LB models for simulating multicomponent multiphase flow have been developed, such as, color-gradient model [29], pseudo potential model [30, 31] and free energy model [32, 33]. The forces that may act on the liquidegas flow in PEMFC GDL include the inertial force, viscous force, gravity and surface tension. The flow regime is governed by a set of dimensionless parameters: i) the capillary number (Ca) which is defined as uh/s, representing the ratio of viscous force to surface tension, ii) the Reynolds number (Re) defined as uD/y, representing the ratio of inertial force to viscous force, and iii) the Bond number (Bo) defined as gjrl rgjD2/s, representing the ratio of gravity to surface tension. In normal PEMFC operations, the Re and Bo are less than the order of 104 and the Ca is less than 105, meaning that the flow is very dominated by the surface tension. Huang et al. [34] evaluated that the LB free energy model is relatively more suitable for the simulation of capillary force dominating multiphase flow in porous medium. In the present work, the LB free energy model improved by Zheng et al. [35] is employed to simulate the liquidegas flow in PEMFC carbon paper GDLs. In CFD schemes, the dynamic behaviors of two-phase flow can be described by the NaviereStokes equations together with an interface evolution equation (CahneHilliard equation), as vr! u ! þ V$r! u! u ¼ VP þ hV2 ! u þ F b; vt
(1)
v4 þ V$4! u ¼ MV2 m4 ; (2) vt where, r, h and ! u denote the fluid density, viscosity, and ! vector velocity, respectively; P is the pressure; t is the time; F b is a general body force; 4 is an order parameter, which identifies the phases and varies between 1 and 1; m4 is the chemical potential; the diffusivity M is chosen as M ¼ q(t4q 0.5)DtG, with q ¼ 1/(t4 þ 0.5), t4 being the dimensionless relaxation time, and G a parameter used to control the interfacial diffusion between phases. By simplification of Boltzmann equation in a lattice, Eq. (1) is reformulated as, x þ! e t Dt; t þ DtÞ ¼ fi ð! x ; tÞ þ Ui ; fi ð!
is the lattice speed of sound (in this work c2s ¼ 1/3); ui is the eq weight coefficient in direction i; fi is the equilibrium distribution function, which is defined as
(3)
where, fi denotes the density distribution function in direction i; ! x is the spatial position; ! e i is the particle velocity in direction i; Dt is the time step; Ui is the collision operator in direction i, which has fi ð! x ; tÞ ð! e i $! 1 ui ! ! uÞ! þ 1 u Þ þ e ð e Ui ¼ i i 2tr c2s c2s tr ! ! m4 V 4 þ F b Dt: eq fi ð! x ; tÞ
(4) In Eq. (4), tr represents the relaxation time which is required to be larger than 0.5 (tr is set to be 0.6 in this work); cs
(5)
The present work solves Eq. (3) in terms of the D3Q19 lattice scheme, in which Di and ui are given as follows. 8 1 > > 3r 6 4m r ; i¼0 þ > 4 < 3 ; Di ¼ > 1 > > : 3 4m4 þ r ; is0 3
8 1 > > > ; i¼0 > > 3 > > < 1 ui ¼ ; i¼16 : > 18 > > > > > > : 1 ; i ¼ 7 18 36
(6)
The macroscopic local density and momentum are obtained from fi by r¼
18 X i¼0
X ! 1 ! fi ; r! u ¼ fi e i þ m4 V4 þ F b : 2 i¼0 18
(7)
The LB equation with a set of D3Q7 distributed gi is solved to recover Eq. (2). x þ! e i Dt; t þ DtÞ ¼ gi ð! x ; tÞ þ ð1 qÞ gi ð! x þ! e i Dt; tÞ gi ð! ! g ð x ; tÞ þ U ; i
i
(8)
where, the collision operator Ui has a definition as Ui ¼
eq gi ð! x ; tÞ gi ð! xÞ t4
(9)
In this work t4 is set to be 1.5. The equilibrium distribution function of gi is formulated as eq gi ¼ Ai þ Bi 4 þ Ci 4! e i $! u:
(10)
The coefficients in Eq. (10) are given as follows.
A0 ¼ 3Gm4 ; Ai ¼ Gm4 2ðis0Þ; B0 ¼ 1; Bi ¼ 0ðis0Þ; Ci ¼ 1=2q: (11) Similarly, the order parameter 4 is obtained from the distribution functions gi by 4¼
6 X
gi :
(12)
i¼0
In this free energy model, the pressure term in Eq. (1) is related to the surface tension, and it can be written as [36]: VP ¼ 4Vm4 Vp0 , where p0 ¼ rc2s . A free energy function in a closed volume with mixture of two fluids is adopted [36], i.e. I k r lnr ; (13) J ¼ dV jð4Þ þ ðV4Þ2 þ 2 3 where, V is a control volume; the bulk free energy density per unit mass for the homogeneous system is jð4Þ ¼ að42 4*2 Þ2 , with a being an amplitude parameter to control the free energy of interactions between the two phases; the square of the gradient term is associated with variations of the density and contributes to the free energy excess in the interfacial region; k is a parameter related to surface tension. The chemical potential is calculated by [36]
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
m4 ¼ a 443 44*2 4 kV2 4;
(14)
and the corresponding pressure tensor can be written as [36] pxy
" # 2 v4 v4 v4 ¼ pdxy k dxy ; vx vx vy
(15)
where, dxy is the Kronecker operator; the subscripts x, y represent Cartesian coordinates; p is calculated by ðV4Þ2 r þ : p ¼ a 344 24*2 42 4*4 k4V2 4 þ 3 2
Static contact angle test
(17)
18 X ui ½4ð! x þ! e i DtÞ 4ð! x Þ V 4¼2 : 2 2 ðDtÞ c i¼1 s 2
(18)
In Eq. (17), m denotes the dimension, m ¼ x, y, z for Cartesian coordinates system. To specify a contact angle in the model, a surface energy jð4s Þ for solid nodes should be added into the free energy calculation and it is assumed to be a simple linear function of 4s [38] (i.e. jð4s Þ ¼ g4s ). A natural boundary constraint for 4s is assumed [38], as ! n $ðV4s Þ ¼ g=k;
(19)
where g is the wetting potential; ! n represents the unit vector with direction pointing to the nearest or a next-nearest fluidside node from the solid node. The 4s and its second order derivatives (V2 4s ) are approximated by [25] P 1
4s ¼
s 4f;1 d v4 vn N
V2 4s y
P h ¼
1
i 4f;1 þ d g=k ; N
P s 2 1 4f;1 4s d v4 vn N
i P h 2 1 4f;1 4s þ d g=k ; ¼ N
61 61 61 lu liquid domain with periodic boundary conditions applied to all the six boundary surfaces, was simulated. We recorded the pressure difference (DP, in lattice unit) across the liquid/gas interface and the radius (R, in lattice unit) of the bubble at the final equilibrium state. Fig. 1 displays the obtained relationship between DP and R. It is obvious that DP is directly proportional to the inverse of R, perfectly amenable to Laplace's Law.
(16)
The first and second derivatives of 4 are calculated by [37] 18 v4 X ui 4ð! x þ! e i DtÞ ¼ eim ; 2 Dt vm c s i¼1
8553
A small volume of liquid on a solid surface tends to spread as film if the surface is hydrophilic or forms a droplet if the surface is hydrophobic. Following the analysis in Refs. [38], the contact angle q of liquid sticking to the solid surface is related to the wetting potential g by pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi 3 3 ~ ~ ; 1g cos q ¼ ð1=2Þ 1þg
pffiffiffiffiffiffiffiffi ~ is a dimensionless parameter, g ~ ¼ g= 2ka. where, g Considering a liquid volume placed on a flat solid surface, we simulated the evolution of the shape of this liquid volume and recorded the shape at steady state. The contact angle q was then easily obtained. Different contact angles were obtained by adjusting the wetting potential related parameters, i.e. k and a. Fig. 2 displays the obtained contact angle q as a function of wetting potential g. The simulated results agree well with the analytical data. In Fig. 2, the simulated liquid droplet shape at final steady state for three representative solid surfaces: hydrophilic (q ¼ 50 ), neutral (q ¼ 90 ) and hydrophobic (q ¼ 130 ) solid surface, are also included.
Numerical reconstruction of carbon paper GDL (20)
Toray TGP-H-060 carbon paper is a typical PEMFC GDL, which is of 0.78 porosity and 7 mm fiber diameter. Its SEM images are displayed in Fig. 3. Morphological anisotropy between the
(21)
where, d is the unit grid space; 4f,l is the order parameter on the nearest or a next-nearest neighboring fluid-side node; N represents the total number of the nearest and next-nearest fluid-side nodes.
Model validation Bubble test For a gas bubble suspended in liquid, the pressure difference (DP) across the liquid/gas interface is related to the radius R of the bubble by Laplace's Law, i.e. DP ¼
2s : R
(23)
(22)
where, s is the surface tension. To verify the model's capability of dealing with the surface tension, the dynamic evolution of a bubble volume (not in spherical shape initially) positioned at the center of a
Fig. 1 e Pressure difference across liquid/gas interface versus bubble radius.
8554
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
direction (i.e. x), the GDL is divided into two regions, one of which has lower PTFE content and the other (neighboring to the gas channel) has higher PTFE content. The sectional images displayed in Fig. 4(c) and (d) give more detailed information about PTFE distribution inside the GDLs: one uniform, the other non-uniform. The reconstructed GDL contains 69 79 79 lattice elements; one lattice element represents a 1.5 1.5 1.5 mm cube, indicating the thickness of the reconstructed GDL is about 103.5 mm.
Results and discussion Influence of PTFE content and distribution on liquid water removal
Fig. 2 e Static contact angle as a function of wetting potential.
through-plane and in-plane direction is clearly seen. Wu and Jiang [39] numerically reconstructed the 3D micro-pore structure of TGP-H-060 carbon paper with the reconstruction method proposed by Schulz et al. [40]. This method was actually adapted from a non-woven structure stochastic generation technique [41]. The present work follows the same method used in Wu and Jiang [39] to reconstruct the Toray TGP-H-060 GDL. The PTFE content in GDL is within 5e40 wt.% and 1.0e13.2% by volume fraction. [43] It is assumed in this work that the PTFE coated on the GDL fibers’ surface has uniform thickness of 0.6 mm, very small compared with the fiber diameter. When distributing the PTFE on the fiber surfaces, we first identified all the numerical elements that represent the surfaces of fibers and randomly specified a certain number of these numerical elements as PTFE. The area integral of all the PTFE elements times the thickness (i.e. 0.6 mm, assumed) of PTFE film gives the total volume of PTFE added to the carbon paper GDL. Fig. 4a shows the micrograph of the reconstructed GDL of 10.0 wt.% PTFE uniformly distributed all through the GDL volume. Fig. 4 (b) shows the reconstructed GDL overall of the same amount of PTFE, 10.0 wt.%, but non-uniformly distributed. Along the GDL thickness (through-plane)
We simulated the liquidegas flow in the computer-generated GDLs (see in Fig. 4 for 3D microstructure) with the LB model described in section (LB Modeling). The PTFE is distinguished from the GDL fibers due to a different contact angle specified. The simulated domain contains 119 79 79 lattice elements, which consists of two sub-domains: one has 69 79 79 lattice elements, and is actually the computer-generated GDL, and the other is an open or empty space. The interface (x ¼ 69 lu) of the two sub-domains mimics the GDL/GC interface. Initially, all the pore space in the GDL and the open space are filled with gas. Liquid is flowing into the GDL from the x ¼ 0 boundary surface (mimics the GDL/CL interface) at a fixed velocity, 2.0 104 lu ts1. The x ¼ 119 lu boundary surface is a fixed pressure boundary with zero pressure value. At all the other boundary surfaces of the simulated domain, the symmetric boundary condition is assumed. Note that all the variables are in the lattice units. In terms of the liquid velocity (2.0 104 lu ts1) applied at the x ¼ 0 boundary surface, the capillary number (Ca) of the fluid flow is calculated to be around 1.08 105, indicating the fluid flow approximates to the capillary fingering flow in PEMFC GDLs [44]. The multiphase LB models generally have difficulty when simulating multiphase flow with a large density ratio between phases, e.g. the air/water system considered in the present work; the simulations are vulnerable to numerical instability [34]. As discussed in Section (LB Modeling), the multiphase flow in PEMFC GDL is normally capillary-fingering flow, dominated by the capillary force. That is to say, the density ratio between phases has relatively little effects on the fluid flow. Therefore, for better numerical stability, the density ratio of the two phases is artificially set to be 1.0 for all the LB
Fig. 3 e SEM images of Toray TGP-H-060 carbon paper GDL [42].
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
8555
Fig. 4 e The computer-generated Toray TGP-H-060 carbon paper GDLs, (a) uniformly distributed PTFE; (b) non-uniformly distributed PTFE, (c) 2D sectional images for the GDL with PTFE uniformly distributed, (d) 2D sectional images for the GDL with PTFE non-uniformly distributed. Red denotes fibers and blue PTFE; in Fig. 4 (d) the gray-colored plane indicates the interface that divides the GDL into two regions of different PTFE contents. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
simulations in the present work. The same treatment to the air/ water density ratio was taken by Hao and Cheng [45] when modeling multiphase transport in PEMFC GDL by LBM. Fig. 5 illustrates the invasion process of liquid water entering into the GDL. The simulated GDL contains 10.0 wt.% PTFE uniformly distributed. The assumed contact angles for the PTFE and carbon fiber are 150 and 92 , respectively. The liquid water is seen to gradually seep into the GDL, As the GDL is hydrophobic and the flow is dominated by the capillary force, the liquid selectively occupies the larger pores firstly. When it arrives to 240,000 time steps (ts), the liquid water is seen to break through the GDL at a few preferential flow paths. Further, Fig. 6 displays the liquid distribution images on three x-sectional planes (with x ¼ 10, 38, and 66 lu, respectively) at three differing time instants: 80,000, 160,000, and
240,000 ts. The liquid water phase is clearly seen to expand toward the y and z directions, occupying local pore space when it penetrates the GDL along the x direction. Fig. 7 presents the liquid saturation profile along the GDL through-plane direction at the final steady state. We justify the onset of steady-state flow according to two criteria: one is the relative deviation of the overall liquid flow rates at the GDL/CL interface and the GDL/GC interface is lower than 0.01, the other is the relative maximum local change of the order parameter 4 in the whole GDL region at two time instants with 200Dt intervals is smaller than 0.01. The displayed liquid saturation values in Fig. 7 are x-section-averaged values at designated x positions. Two cases are considered for comparison. The two cases differ from each other due to the specified PTFE distributions. One has a GDL of 10.0 wt.%
8556
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
Fig. 5 e Illustration of the invasion process of liquid water entering into GDL. Green: liquid water; blue: fibers. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
uniformly distributed PTFE while the other case considers a GDL with non-uniformly distributed PTFE but the overall PTFE content is 10.0 wt.%, the same as the former case. For the latter case, due to the non-uniform distribution of PTFE, the GDL is divided into two regions, indicated by “A” and “B” in Fig. 7. In region A, the PTFE content is 6.5 wt.%, while in region B, it is 21.8 wt.%. For both cases, the contact angle assumed at the GDL fiber surface is 92 while at the PTFE surface is 150 . Seen from Fig. 7, if the GDL contains non-uniformly distributed PTFE, more liquid water is accumulated in region A, harmful to the liquid water removal from CL to GC. For the case with GDL of non-uniformly distributed PTFE, the region B GDL that contains more PTFE acts like a “capillary barrier”, preventing liquid water from penetrating the GDL. This raises the liquid saturation level in region A. The experimental work by Song et al. [46] found that the existence of a less hydrophobic region in the central region of GDL would impede the water removal, qualitatively according with the observations got from Fig. 7. Checking the original data used to draw the two curves in Fig. 7, we find the largest relative difference of sectionaveraged liquid saturation between the two cases is more than 10%. Four insets of liquid water distribution images on typical x-sectional planes are specially depicted in Fig. 7, which show intuitively the differences of the liquid water distribution inside GDL between the two cases.
It is seen from Fig. 7 also that at the final steady state, partial region of the GDL close to the GDL/CL interface is full of liquid water for both cases. When the liquid water flows through the GDL, it will expand toward the in-plane (i.e. y and z) directions and occupy all the pores along its flow path (refer to Fig. 6) if the flow resistance in the through-plane (i.e. x) direction is larger than the capillary force the local smallest hydrophobic pore imposes on the liquid water. We deduce from Fig. 7 that the liquid flow resistance of ~40 mm thick GDL is comparable to the entry pressure of the smallest hydrophobic pore in the GDL. To investigate the influence of PTFE content on the liquid removal, cases with GDL of 10.0 wt.%, 18.1 wt.% and 25.6 wt.% uniformly distributed PTFE are considered for comparison. Fig. 8 gives the liquid saturation profile along the GDL throughplane direction (x) at the final steady state for all the three cases. The displayed liquid saturation values are x-sectionaveraged values. The liquid saturation level is generally seen to decrease if the PTFE content in GDL is higher. Lowering PTFE content means to reduce the hydrophobicity of GDL. A less hydrophobic GDL means the capillary resistance for liquid water enters into a pore is lowered. Hence, if the GDL contains less PTFE, not only larger pores, but also more pores of smaller pore size in the GDL will be occupied by the liquid water; liquid water dwells in more pores after it enters into the GDL, raising
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
8557
Fig. 6 e Images of liquid water distribution on x-sectional planes at three differing time instants. Green: liquid water; blue: fibers. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
the liquid saturation level in the GDL. Qualitatively similar results have been reported by Hao and Cheng [45]. Increasing PTFE content in GDL makes the GDL become more hydrophobic. The liquid water needs to overcome larger capillary force to expand in the in-plane directions and occupy less relatively smaller pores. We see from Fig. 8 that a smaller part of the GDL is full of liquid water if the PTFE content in GDL is increased from 10.0 wt.% to 18.1 wt.%; for the case with GDL of 25.6 wt.% PTFE content, even no GDL x-cross section is full of liquid water except the inlet GDL/CL interface. As described in the Introduction section, there is still lack of a general cognition about the optimum PTFE content in GDL. From the simulation results presented in Fig. 8, the liquid saturation level in GDL diminishes with the increase of PTFE content in GDL. It seems that there is no an optimum PTFE content. We speculate that determining the optimum PTFE content in GDL should comprehensively consider more involved factors, not solely rely on the resultant liquid water saturation level in GDL. As an example, an additional inset in Fig. 8 gives the calculated liquid pressure drops across the GDL for all the three cases. The more hydrophobic GDL significantly elevates the entry pressure of liquid water entering into
the GDL. This may lead to the CL more vulnerable to water flooding. The LB simulation results are dependent on the lattice resolution, which was examined by the Richardson extrapolation technique [47]. For the same computer-generated GDL (microstructure displayed in Fig. 4b), two additional LB simulations with lattice resolutions of 2.5 mm and 1.0 mm respectively were performed. Considering the obtained sectionaverage liquid saturation values on the GDL/GC interface, we do Richardson extrapolation and find the extrapolated accurate value that is from the ideal case with infinite small lattice resolution deviates from the calculated value based on the simulated case with 1.5 mm lattice resolution by only about 4%. This further corroborates the numerical accuracy of the LB simulations carried out in the present work.
Influence of PTFE content and distribution on relative permeability As described in Section (Influence of PTFE content and distribution on liquid water removal), PTFE content and distribution affect the liquidegas flow, which is bound to affect
8558
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
applied in the through-plane direction. The body force applied is well controlled to guarantee the capillary number of fluid flow is within the order of 105. When the fluid flow reaches steady-state, we calculate the average liquid (or gas) velocity in the through-plane direction (ux ) and then derive the relative permeability Kr according to Darcy’ law, as Kr ¼
Fig. 7 e Liquid saturation profile along the GDL thickness direction: effect of PTFE distribution. The insets display the liquid distribution images on x-sectional planes; liquid phase is green-colored. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 8 e Liquid saturation profile along the GDL thickness direction: effect of PTFE content.
the relative permeability of each phase. A series of cases are simulated to derive the relationship between relative permeability and phase saturation. The simulated domain is solely the computer-generated GDLs (illustrated in Fig. 4); periodic boundary condition is considered for the through-plane (x) direction while symmetric boundaries are assumed for the boundary surfaces in the other two directions. Initially, a certain amount of liquid is randomly positioned in the pore space of the GDL. The fluid flow is driven by a body force (Fb)
ux h KFb
(24)
where, Fb is the applied body force; K is the absolute permeability. K is determined to be 4.05 1012 m2 from single phase LB simulations performed in the present work, slightly deviates from the range of K values measured by Mathias et al. [43], 5 1012 m2 e 10 1012 m2. Changing the initially positioned liquid amount in the GDL, a case of different liquid water saturation level in GDL is set. Fig. 9 presents the simulated steady-state liquid water distribution in GDLs of different PTFE contents. The specified contact angles for GDL fiber surface and PTFE surface are 92 and 150 , respectively. The two cases considered have the same phase distribution initially, i.e. the same overall liquid water saturation (sl ¼ 0.42) in the GDL. As evident by Fig. 8, the two cases give significantly different results; the liquid water distribution for the case with GDL of higher PTFE content is more fragmented at the final steady-state, simply due to the more hydrophobic GDL. Fig. 10 presents the derived relationship of relative permeability versus phase saturation. Results from two sets of cases are included in Fig. 10 for comparison. The two sets of cases differ from each other due to the PTFE content. For all the cases, the specified contact angle on GDL fiber surface is 92 , while the contact angle specified on PTFE surface is 150 . It is seen that the relative permeability increases with the increasing phase saturation for the liquid and gas phases both. A natural explanation for this is a higher phase saturation level makes the corresponding fluid phase more easily form connective flow paths all-through the GDL. Inspection on the results of Krl versus sl shown in Fig. 10 finds when the sl is less than ~0.5, the Krl is larger in the GDL of lower PTFE content, whereas when sl is more than ~0.5, the Krl is larger in the GDL of higher PTFE content. There may exist two underlying mechanisms [48]: i) liquid flow in the GDL of higher PTFE content is more vulnerable to fragmentation due to the higher hydrophobicity of GDL; ii) the liquid in more hydrophobic GDL prefers dwelling in larger pores, which have smaller specific surface area interfacing with the solid GDL fibers, leading to smaller liquid flow resistance. It is speculated that the former mechanism dominates the liquid flow when the liquid saturation level in GDL is low. However, when the liquid saturation level in GDL is sufficiently high (say, >0.5), most of the pore space in the GDL is occupied by the liquid phase and there should always have some connective liquid flow paths all-through the GDL. Hence, the latter mechanism takes over the dominating role when the liquid saturation level in GDL is high enough, i.e. sl > 0.5 in the present work. Seen from Fig. 10 also, the PTFE content in GDL has influence on the relationship of Krg versus sg, the Krg in the GDL of lower PTFE content is larger than the one in the GDL of higher PTFE content. This can be explained by the improved
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
8559
Fig. 9 e Simulated liquid water distribution in GDL of (a) 10.0 wt.% and (b) 18.1 wt.% PTFE content. For both cases the overall liquid saturation in GDL is 0.42; on the 2D section cuts, green denotes liquid water. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
connectivity for the gas phase in the GDL of lower PTFE content according to the above-described mechanism i). Both the liquid and gas phase relative permeability values decrease to zero when the phase saturation reduces to a certain low level. For the cases with more hydrophobic GDL, the Krg is seen to decrease to zero when the gas phase saturation is lowered down to about 0.2. Qualitatively similar results were reported in Refs. [26 and 27]. For multi-phase flow in porous medium, if one phase is in very low saturation level, it can be totally fragmented and there is thus no possibility for this phase to form an effective all-through flow path and then gets stuck in the porous medium. This is to some extent similar to the concept of “percolation threshold” of porous media [47]. The inset in Fig. 10 summarizes the relative permeability values predicted by LB models available in open literature. Refs. [25 and 26] used the free energy LB model and the GDL contact angles were considered to be 105 and 145 , respectively; Ref. [27] used the SeC model and 135 GDL contact angle was assumed. The results from different literature vary in similar
trend though differences exist. All the literature data qualitatively agree with the results obtained in the present work. Fig. 11 compares the derived Kr-s relationships for GDL of uniformly distributed PTFE with those for GDL of nonuniformly distributed PTFE. The PTFE content and distribution for both the two sets of cases considered here have been detailed in Section (Influence of PTFE content and distribution on liquid water removal) in relation with Fig. 7. For the GDL of non-uniform PTFE distribution, the region with higher PTFE content acts like a “capillary barrier”, impeding the liquid flow. Seen from Fig. 11, the non-uniform PTFE distribution in GDL does reduce the Krl value when the sl is not more than ~0.55, but with the increase of sl (i.e., when sl is greater than ~0.55), the Krl in the GDL with non-uniform PTFE distribution exceeds the value of Krl in the GDL with uniform PTFE distribution. This can also be explained by the aforementioned two mechanisms. The PTFE distribution shows relatively smaller effect on the Krg; the Krg in the GDL of uniform PTFE distribution is slightly larger than that in the
8560
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
Genuchten function particularly when the phase saturation is low; the van Genuchten function underestimates the values of Krl especially when the non-wetting liquid phase saturation takes an intermediate value.
Summary and conclusions
Fig. 10 e Effects of PTFE content on the liquid and gas relative permeability of GDL.
Fig. 11 e Effects of PTFE distribution on the liquid and gas relative permeability of GDL.
GDL of non-uniform PTFE distribution, indicating the nonuniform distribution of PTFE in GDL may somehow prevents the gas phase transport. Suitable measures taken to avoid non-uniform distribution of PTFE in GDL may relieve the mass transport limitation of PEMFC at high current density operations. In Figs. 10 and 11, the LBM-calculated Kr-s data for both wetting and non-wetting phases are also compared with the van Genuchten function [28], which is formulated as 2ðj1Þ=j ðj1Þ=j ; Kr ¼ s1=3 1 ð1 sÞ
(25)
where the value of j ranges from 2 to 5. It is seen that the Kr-s data for the wetting gas phase generally better fit the van
A 3D multiphase LB model was developed and employed to study the impact of PTFE content and distribution on the liquidegas transport in PEMFC carbon-paper GDLs. The capability and reliability of the developed LB model were first tested by simulating the dynamic evolution of a gas bubble suspended in liquid and of a liquid volume put on a solid surface. The simulation results accord very well with the corresponding analytical results. The 3D micro-pore structure of TGP-H-060 carbon paper was numerically reconstructed and part of the fibers’ surface was randomly specified as PTFE. Two classes of cases were considered. One has PTFE uniformly distributed all through the GDL volume. The other has non-uniformly distributed PTFE: in the GDL sub-region neighboring to the GC, the PTFE content is higher, which mimics the PTFE accumulation toward the GDL/GC interface due to incomplete penetration of PTFE into GDL in practical PTFE treatment processes. The simulated liquid flow in GDL shows evident capillary force dominated characteristics, gradually seeping into the GDL and selectively occupying larger pores first. A GDL with higher PTFE content is more hydrophobic. The GDL sub-region containing more PTFE acts like a “capillary barrier”, which prevents liquid water from passing the GDL and raises the liquid saturation level in the remaining GDL region adjoining to the CL. The simulation results about liquid saturation level in GDL do not support that there exists an optimum PTFE content in GDL. Three cases with GDL of differing PTFE contents, 10.0 wt.%, 18.1 wt.% and 25.6 wt.%, respectively, were simulated. The simulated liquid saturation level in GDL diminishes with the increase of PTFE content in GDL. It is found from simulation results also that the liquid pressure drop across the GDL is higher if the case has more PTFE in GDL. We deduce determining the optimum PTFE content in GDL must consider more involved factors, and should not rely solely on the resultant liquid saturation level in GDL. LB simulations enable the calculation of relative permeability of liquid and gas phase in the carbon paper GDL. We derived the relationship of relative permeability versus phase saturation for both phases from a series of additional simulation cases. The relative permeability is found to increase with the increasing phase saturation for both phases. PTFE content and distribution are found to have important influence on the relationship of liquid phase relative permeability versus phase saturation. If the (overall) liquid saturation is low (say, lower than ~0.5 for the cases considered in the present work), the liquid phase relative permeability is larger in the GDL of lower PTFE content, whereas it is smaller if the liquid saturation is high (say, higher than ~0.5). The liquid phase relative permeability in the GDL with non-uniform PTFE distribution is found to be reduced when the liquid saturation in GDL is low (say, lower than ~0.55), but when the liquid saturation is high (say, higher than ~0.55), it is larger compared with that in the GDL
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
with uniform PTFE distribution. Two underlying mechanisms take effects here: i) liquid flow in the GDL of higher PTFE content is more vulnerable to fragmentation due to the higher hydrophobicity of GDL; ii) the liquid in more hydrophobic GDL prefers dwelling in larger pores, which have smaller specific surface area interfacing with the solid GDL fibers, leading to smaller liquid flow resistance. Mechanism i) dominates the liquid flow when the liquid saturation level in GDL is low, while mechanism ii) takes over the dominating role when the liquid saturation level in GDL is sufficiently high. The PTFE content and distribution in GDL also show influence on the relationship of gas phase relative permeability versus gas saturation. The gas phase relative permeability in the GDL of lower PTFE content is found to be larger than that in the GDL of higher PTFE content and the non-uniform distribution of PTFE in GDL is found to slightly reduce the gas phase relative permeability value. It is found as well that the LBM-predicted relationship of gas phase relative permeability versus phase saturation generally better accord with the van Genuchten function particularly when the phase saturation is low, whereas the van Genuchten function underestimates the liquid phase relative permeability especially when the liquid phase saturation takes an intermediate value.
Acknowledgments Financial support received from the Key Scientific Development Project of Guangdong Province (2015A030308019), the Guangzhou Scientific and Technological Development Plan (2014J4100217), and the CAS “100 talents” Program (FJ) is gratefully acknowledged.
Nomenclature a, k Ai, Bi, Ci, Bo Ca cs d D ! ei fi eq fi Fb g gi eq
gi
j K Kr lm lu M P
parameters related to surface tension Di coefficients in equilibrium distribution function Bond number Capillary number lattice speed of sound, lu ts1 unit grid space, lu characteristic length, m particle speed, lu ts1 density distribution function in the ith direction equilibrium density distribution function body force, lm lu2 ts2 gravity, m/s2 distribution function for the order parameter in the ith direction equilibrium distribution function of the order parameter exponent in the van Genuchten function absolute permeability, m2 relative permeability unit of mass in LBM unit of length in LBM mobility of fluid pressure of fluid, lm lu1 ts2
p0 p Q q R Re s T ts u ! u V ! x
8561
thermodynamic pressure, lm lu1 ts2 bulk pressure, lm lu1 ts2 flow rate, lm lu2 ts1 constant coefficient related to the relaxation time t4 droplet radius, lu Reynolds number saturation of a fluid phase the rate of the order parameter variation; time, ts unit of time in LBM velocity, lu ts1 velocity vector, lu ts1 controlling volume, lu3 spatial position
Greek symbols G parameter used to control the fluid mobility g wetting potential ~ g a dimensionless parameter in relation with Eq. (23) Kronecker operator dxy h viscosity, lm lu2 ts1 m4 chemical potential r density, lm lu3 s surface tension, lm ts2 tr, t4 dimensionless relaxation times, ts y kinematic viscosity, m2/s 4 order parameter 4* order parameter in equilibrium state j bulk free energy density J free energy function weight coefficient ui collision operator in the ith direction Ui Subscripts/superscripts x the direction x in Cartesian coordinates y the direction y in Cartesian coordinates l liquid g gas i the ith discrete distribution direction s solid in inlet out outlet
references
[1] Pasaogullari U, Wang CY. Liquid water transport in gas diffusion layer of polymer electrolyte fuel cells. J Electrochem Soc 2004;151:399e406. [2] Wang CY. Fundamental models for fuel cell engineering. Chem Soc Rev 2004;104:4727e66. [3] Jiang FM, Wang CY. Numerical modeling of liquid water motion in a polymer electrolyte fuel cell. Int J Hydrogen Energy 2014;39:942e50. [4] He W, Yi JS, Nguyen TV. Two-phase flow model of the cathode of PEM fuel cells using interdigitated flow field. Am Inst Chem Eng 2000;46:2053e64. [5] Natarajan D, Nguyen TV. A two-dimensional, two-phase, multicomponent, transient model for the cathode of a proton exchange membrane fuel cell using conventional gas distributors. J Electrochem Soc 2001;148:1324e35.
8562
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
[6] Wu H, Li XG, Berg P. On the modeling of water transport in polymer electrolyte membrane fuel cells. Electrochimica Acta 2009;54:6913e27. [7] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39:201e25. [8] Theodorakakos A, Ous T, Gavaises M. Dynamics of water droplets detached from porous surfaces of relevance to PEM fuel cells. J Colloid Interface Sci 2006;300:673e87. [9] Zhu X, Sui PC, Ned D. Three-dimensional numerical simulations of water droplet dynamics in a PEMFC gas channel. J Power Sources 2008;181:101e15. [10] Park JW, Jiao K, Li XG. Numerical investigations on liquid water removal from the porous gas diffusion layer by reactant flow. Appl Energy 2010;87:2180e6. [11] Aimy B, David S, Ned D. Dynamic water transport and droplet emergence in PEMFC gas diffusion layers. J Power Sources 2008;176:2180e6. [12] Mohamed EA, Ben A, Sass BN. Numerical simulation of droplet dynamics in a proton exchange membrane (PEMFC) fuel cell micro-channel. Int J Hydrogen Energy 2015;40:1333e42. [13] Park J, Li X. Multi-phase micro-scale flow simulation in the electrodes of a PEM fuel cell by lattice Boltzmann method. J Power Sources 2008;178:248e57. [14] Gao Y, Zhang XX, Rama P, Chen R, Ostadi H, Jiang K. Lattice Boltzmann simulation of water and gas flow in porous gas diffusion layers in fuel cells reconstructed from microtomography. Comput Math Appl 2013;65:891e900. [15] Chen L, Luan HB, He YL, Tao WQ. Numerical investigation of liquid water transport and distribution in porous gas diffusion layer of a proton exchange membrane fuel cell using lattice Boltzmann method. Russ J Electrochem 2012;48:712e26. [16] Bevers D, Rogers R, Bradke MV. Examination of the influence of PTFE coating on the properties of carbon paper in polymer electrolyte fuel cells. J Power Source 1998;63:193e201. [17] Lin GY, Nguyen TV. Effect of thickness and hydrophobic polymer content of the Gas diffusion layer on electrode flooding level in a PEMFC. J Electrochem Soc 2005;152:1942e8. [18] Kandlikar SG, Garofalo ML, Lu Z. Water management in a PEMFC: water transport mechanism and material degradation in gas diffusion layers. Fuel Cells 2011;11:814e23. [19] Wua JF, Martina JJ, Orfino FP, Wang HJ, Legzdins C, Yuan XZ, et al. In situ accelerated degradation of gas diffusion layer in proton exchange membrane fuel cell part I: effect of elevated temperature and flow rate. J Power Sources 2010;195:1888e94. [20] Lee C, Merida W. Gas diffusion layer durability under steadystate and freezing conditions. J Power Sources 2007;164:141e53. [21] Paganin VA, Ticianelli EA, Gonzalez ER. Development and electrochemical studies of gas diffusion electrodes for polymer electrolyte fuel cells. J Appl Electrochem 1996;26:297e304. [22] Lim C, Wang CY. Effects of hydrophobic polymer content in GDL on power performance of a PEM fuel cell. Electrochimica Acta 2004;49:4149e56. [23] 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:8401e9. [24] Gostick JT. Random pore network modeling of fibrous PEMFC gas diffusion media using voronoi and delaunay tessellations. J Electrochem Soc 2013;160:731e43. [25] Niu XD, Munekata T, Hyodo SA, Suga K. An investigation of water-gas transport processes in the gas-diffusion-layer. J Power Sources 2007;172:542e52. [26] Hao L, Cheng P. Pore-scale simulations on relative permeabilities of porous media by lattice Boltzmann method. Int J Heat Mass Transf 2010;53:1908e13.
[27] Koido T, Furusawa T, Moriyama K. An approach to modeling two-phase transport in the gas diffusion layer of a proton exchange membrane fuel cell. J Power Sources 2008;175:127e36. [28] van Genuchten MT. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J 1980;44:892e8. [29] Gunstensen AK, Rothman DH, Zaleski S, et al. Lattice Boltzmann models for nonideal gases. Phys Rev E 2000;62:4982. [30] Shan, X.; Chen, H, Lattice Boltzmann model for simulating flows with multiple phases and components. Phys Rev E 1993; 47; 1815. [31] Shan, X.; Chen, H, Simulation of nonideal gases and liquidgas phase transitions by lattice Boltzmann equation. Phys Rev E 1994; 49; 2941. [32] Swift MR, Osborn WR, Yeomans JM. Lattice Boltzmann simulations of nonideal fluids. Phys Rev Lett 1995;75:830. [33] Swift MR, Orlandini E, Osborn WR, Yeomans JM. Lattice Boltzmann simulations of liquid-gas and binary fluid systems. Phys Rev E 1996;54:5041. [34] Huang HB, Wang L, Lu XY. Evaluation of three lattice Boltzmann models for multiphase flows in porous media. Comput Math Appl 2011;61:3606e17. [35] Zheng HW, Shu C, Chew YT. A lattice Boltzmann model for multiphase flows with large density ratio. J Comput Phys 2006;218:353e71. [36] Kendon VM, Gates ME, Pagonabarraga I, Desplat JC, Bladon P. Inertial effects in three-dimensional synodal decomposition of a symmetric binary fluid mixture: a lattice Boltzmann study. J Fluid Mech 2001;440:147e203. [37] Lee T, Lin CL. A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio. J Comput Phys 2005;206:16e47. [38] Briant AJ, Papatzacos P, Yeomans JM. Lattice Boltzmann simulations of contact line motion in a liquidegas system. Philos Trans R Soc Lond A 2002;360:485e95. [39] Wu W, Jiang FM. Microstructure reconstruction and characterization of PEMFC electrodes. Int J Hydrogen Energy 2014;39:15894e906. [40] Schulz VP, Becker J, Wiegmann A, Mukherjee PP, Wang CY. Modeling of two-phase behavior in the gas diffusion medium of PEMFCs via full morphology approach. J Electrochem Soc 2007;154:419e26. [41] Katja S, Stefanie P, Doris RB, Andreas W, Joachim O. Design of acoustic trim based on geometric modeling and flow simulation for non-woven. Comput Mater Sci 2006;38:56e66. [42] Bazylak A, Sinton D, Liu ZS, Njilali N. Effect of compression on liquid water transport and microstructure of PEMFC gas diffusion layers. J Power Sources 2007;163:784e92. [43] Mathias MF, Roth J, Feleming J, Lehnert W. In: Lietsich W, Lamm A, Gasteiger HA, editors. In hand book of fuel cells: fundamentals, technology and applications, Vol. 3. Chichester: John Wiley & Sons; 2003. Part 1. [44] Sinha PK, Wang CY. Pore-network modeling of liquid water transport in gas diffusion layer of a polymer electrolyte fuel cell. Electrochimica Acta 2007;52:7936e45. [45] Hao L, Cheng P. Lattice Boltzmann simulations of water transport in gas diffusion layer of a polymer electrolyte membrane fuel cell. J Power Sources 2010;195:3870e81. [46] Song W, Yu HM, Shao ZG, Yi BL, Lin J, Liu N. Effect of polytetrafluoroethylene distribution in the gas diffusion layer on water flooding in proton exchange membrane fuel cells. Chin J Catal 2014;35:468e73. [47] Jiang FM, Oliveira MSA, Sousa ACM. Mesoscale SPH modeling of fluid flow in isotropic porous media. Comput Phys Commun 2007;176:471e80. [48] Li HN, Pan CX, Miller CT. Pore-scale investigation of viscous coupling effects for two-phase flow in porous media. Phys Rev 2005;72:026705.
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Impact of PTFE content and distribution on liquidegas flow in PEMFC carbon paper gas distribution layer: 3D lattice Boltzmann simulations Wang Chen, Fangming Jiang* Laboratory of Advanced Energy Systems, Guangdong Key Laboratory of New and Renewable Energy Research and Development, CAS Key Laboratory of Renewable Energy, Guangzhou Institute of Energy Conversion, Chinese Academy of Sciences (CAS), Guangzhou, 510640, China
article info
abstract
Article history:
A 3D multiphase lattice Boltzmann model (LBM) was established and employed to study
Received 20 November 2015
the impact of polytetrafluoroethylene (PTFE) content and distribution on the liquidegas
Received in revised form
transport in proton exchange membrane fuel cell (PEMFC) carbon-paper gas distribution
2 February 2016
layer (GDL). In the computer-generated carbon paper GDL, part of the fibers’ surface was
Accepted 9 February 2016
randomly specified as PTFE to achieve a desired PTFE content and distribution. For non-
Available online 19 April 2016
uniform PTFE distribution cases, a GDL sub-region neighboring to gas channel (GC) was set to have higher PTFE content. PTFE content and distribution show important influence
Keywords:
on the liquidegas flow and on the relationship of relative permeability versus phase
Proton exchange membrane fuel cell
saturation. The GDL sub-region of more PTFE acts like a capillary barrier and raises the
Gas distribution layer
liquid saturation level in the rest part of GDL. The simulated liquid saturation level in GDL
Lattice Boltzmann method
diminishes with the increase of PTFE content in GDL. The liquid phase relative perme-
Polytetrafluoroethylene
ability is larger in the GDL of lower PTFE content if the (overall) liquid saturation is lower
Liquidegas flow
than a threshold value, whereas it becomes smaller if the liquid saturation is higher than this threshold value. The liquid phase relative permeability in the GDL with non-uniform PTFE distribution is reduced when the liquid saturation in GDL is lower than a certain value, but when the liquid saturation is higher than this value, it is larger compared with that in the GDL with uniform PTFE distribution. The gas phase relative permeability in the GDL of lower PTFE content is larger than tat in the GDL of higher PTFE content and the nonuniform distribution of PTFE in GDL slightly reduces the gas phase relative permeability value. © 2016 Hydrogen Energy Publications LLC. Published by Elsevier Ltd. All rights reserved.
* Corresponding author. Tel.: þ86 20 87057656. E-mail address: [email protected] (F. Jiang). http://dx.doi.org/10.1016/j.ijhydene.2016.02.159 0360-3199/© 2016 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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
Introduction Fuel cells directly convert chemical energy of fuels into electricity and thus have high energy conversion efficiency. Among many types of fuel cells, the H2/air proton exchange membrane fuel cell (PEMFC) may be the most promising to be used as alternative power engine of future vehicles due to advantages like low operating temperature (usually lower than 100 C), zero emission and high power density. Though tremendous progresses on the development of PEMFC engines have been achieved in the past two decades, challenges still remain for large-scale commercialization of passenger PEMFC cars. Further improving the performance (e.g. durability) and reducing the cost of PEMFC may be the most critical challenges. A thorough understanding to water transport in PEMFC is helpful to tackle these challenges. A PEMFC is usually layer-structured, consisting of two gas distribution layers (GDLs), two catalyst layers (CLs), a membrane, two gas flow channels, two bi-polar plates, and some coolant channels. The GDL, as one key component of PEMFC, is located between the gas flow channel and the catalyst layer (CL), and is normally highly porous. The pores in GDL provide pathways for the transport of gaseous reactants and the removal of product water. During practical operations of PEMFC, liquid water may emerge inevitably. Excess liquid water accumulating in GDL may block the pores and hinder the transport paths of reactant gases, a major mechanism limiting the fuel cell high current density operation. A deep knowledge of the multiphase flow in GDL is essential to achieving an increased PEMFC performance. Numerical modeling may be the most efficient method to study the complex multidisciplinary, multiphase transport in PEMFC. To date, a few theoretical models have been established for the simulation of liquidegas flow in PEMFC GDL, e.g., multiphase mixture model (M2) [1e3], two-fluid model (F2) [4e6], and volume of fluid (VOF) [7] model. The M2 model treats the multiphase fluid as a mixture when solving the momentum and continuity equations; it defines the saturation of liquid and considers the interactions between different phases (including phase change) when solving the species transport equations; the relative velocity between phases can be obtained via post processing of the simulated data. The F2 model solves two sets of conservation equations for the gas and liquid phase, respectively, while the interactions between phases are approximated empirically or semi-empirically. Different from the M2 and F2 models, the VOF model solves an additional scalar transport equation for the scalar of liquid volume fraction; it reconstructs the interface of gas and liquid phases every time step and thus explicitly differentiates the interface between phases, which enables more accurate implementation of interactions between phases. However, the VOF method is very computationally intensive, as it generally requires very small time step and very fine mesh to resolve the complicated geometry of liquidegas interface and the complex interactions between phases, which restricts the VOF model being only successfully used for modeling multiphase flows in PEMFC gas channels (GC) [8,9] or GDLs with artificially simplified pore structure [10,11].
8551
Besides the above-reviewed macroscopic models, the lattice Boltzmann (LB) method, based on the mesoscopic kinetic theory, has been employed to simulate the liquidegas flow in GC and porous GDL in recent years [12e15]. Unlike the conventional computational fluid dynamics (CFD) models, the LB model first requires resolving the mesoscopic pore structures in GDL and solves the Boltzmann equation being simplified for a lattice. Therefore, the LB model can give mesoscopic insight into the gaseliquid flow in GDL. The GDL is made of fibrous porous medium, either carbon paper or carbon cloth. To enhance the hydrophobicity of GDL, the polytetrafluoroethylene (PTFE) is usually used to treat the GDL [16,17]. After the treatment, about 5%e40% of PTFE by weight (i.e. 5 wt.% e 40 wt.% PTFE) may be coated on the surfaces of GDL fibers [18]. The PTFE coating, particularly in the cathode GDL, may lose due to hot-pressing and mechanical compression when assembling PEMFC [19], or could be washed off by flow of fluid during operations [18,19], or diminishes during sub-zero temperature cold start cycles [20]. How does the content of PTFE in GDL affect the PEMFC performance? Despite some works (mainly experimental) conducted, no common cognition has been reached. Paganin et al. [21] reported an optimal content of 15 wt.% whereas the result of Lim and Wang [22] was 10 wt.%. Moreover, when doing PTFE treatment to GDL, the PTFE may not totally penetrate the GDL and accumulates toward the GDL surface. The distribution of PTFE in GDL must have influence on the involved liquidegas flow, which has been rarely studied except a recent work by Maleimanesh and Akbari [23]. Their two-dimensional (2D) LB simulations revealed even a small uncoated (by PTFE) region in GDL would impede the removal of liquid water. The relative permeability is an important parameter characterizing the multiphase flow in porous media, defined as the ratio of the effective permeability to the absolute permeability of a fluid phase. In PEMFC, the relative permeability (Krl for liquid phase and Krg for gas phase) reflects the obstruction of one phase acts on the other. It is generally hard to directly measure the Krl and Krg in the GDL, instead, porescale numerical modeling is widely used as an effective tool to determine these parameters. Gostick [24] by pore network model, Niu et al. [25], Hao et al. [26] and Koido et al. [27] by LB models, simulated the two-phase flow in GDL and calculated the relative permeability for each phase. Niu et al. [25] found that the relationship between Krl and liquid saturation sl is approximately amenable to the van Genuchten function [28]. Koido et al. [27] predicted the Kr for each phase by a combination of the single phase and the multicomponent multiphase LB method and they validated their results by comparing the predicted gas relative permeability values against experimental data. The above numerical modeling works to some extent revealed the underlying mechanisms of the two-phase transport in fibrous GDL. However, the impact of PTFE's content and distribution were ignored or not directly simulated in these works. The present work will establish a three-dimensional (3D) LB model for modeling the gaseliquid flow in PEMFC carbon paper GDL. The major purpose is to reveal the effects of PTFE content and distribution on the involved multiphase transport. Specially, the relative permeability, Krl and Krg, will be
8552
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
derived with respect to carbon paper GDLs of various PTFE content and distribution scenarios.
3 9 eq u $! u þ ð! e i $! x i ; tÞ ¼ ui Di þ ui r 3! e i $! u ! u Þ2 : fi ð! 2 2
LB modeling A few LB models for simulating multicomponent multiphase flow have been developed, such as, color-gradient model [29], pseudo potential model [30, 31] and free energy model [32, 33]. The forces that may act on the liquidegas flow in PEMFC GDL include the inertial force, viscous force, gravity and surface tension. The flow regime is governed by a set of dimensionless parameters: i) the capillary number (Ca) which is defined as uh/s, representing the ratio of viscous force to surface tension, ii) the Reynolds number (Re) defined as uD/y, representing the ratio of inertial force to viscous force, and iii) the Bond number (Bo) defined as gjrl rgjD2/s, representing the ratio of gravity to surface tension. In normal PEMFC operations, the Re and Bo are less than the order of 104 and the Ca is less than 105, meaning that the flow is very dominated by the surface tension. Huang et al. [34] evaluated that the LB free energy model is relatively more suitable for the simulation of capillary force dominating multiphase flow in porous medium. In the present work, the LB free energy model improved by Zheng et al. [35] is employed to simulate the liquidegas flow in PEMFC carbon paper GDLs. In CFD schemes, the dynamic behaviors of two-phase flow can be described by the NaviereStokes equations together with an interface evolution equation (CahneHilliard equation), as vr! u ! þ V$r! u! u ¼ VP þ hV2 ! u þ F b; vt
(1)
v4 þ V$4! u ¼ MV2 m4 ; (2) vt where, r, h and ! u denote the fluid density, viscosity, and ! vector velocity, respectively; P is the pressure; t is the time; F b is a general body force; 4 is an order parameter, which identifies the phases and varies between 1 and 1; m4 is the chemical potential; the diffusivity M is chosen as M ¼ q(t4q 0.5)DtG, with q ¼ 1/(t4 þ 0.5), t4 being the dimensionless relaxation time, and G a parameter used to control the interfacial diffusion between phases. By simplification of Boltzmann equation in a lattice, Eq. (1) is reformulated as, x þ! e t Dt; t þ DtÞ ¼ fi ð! x ; tÞ þ Ui ; fi ð!
is the lattice speed of sound (in this work c2s ¼ 1/3); ui is the eq weight coefficient in direction i; fi is the equilibrium distribution function, which is defined as
(3)
where, fi denotes the density distribution function in direction i; ! x is the spatial position; ! e i is the particle velocity in direction i; Dt is the time step; Ui is the collision operator in direction i, which has fi ð! x ; tÞ ð! e i $! 1 ui ! ! uÞ! þ 1 u Þ þ e ð e Ui ¼ i i 2tr c2s c2s tr ! ! m4 V 4 þ F b Dt: eq fi ð! x ; tÞ
(4) In Eq. (4), tr represents the relaxation time which is required to be larger than 0.5 (tr is set to be 0.6 in this work); cs
(5)
The present work solves Eq. (3) in terms of the D3Q19 lattice scheme, in which Di and ui are given as follows. 8 1 > > 3r 6 4m r ; i¼0 þ > 4 < 3 ; Di ¼ > 1 > > : 3 4m4 þ r ; is0 3
8 1 > > > ; i¼0 > > 3 > > < 1 ui ¼ ; i¼16 : > 18 > > > > > > : 1 ; i ¼ 7 18 36
(6)
The macroscopic local density and momentum are obtained from fi by r¼
18 X i¼0
X ! 1 ! fi ; r! u ¼ fi e i þ m4 V4 þ F b : 2 i¼0 18
(7)
The LB equation with a set of D3Q7 distributed gi is solved to recover Eq. (2). x þ! e i Dt; t þ DtÞ ¼ gi ð! x ; tÞ þ ð1 qÞ gi ð! x þ! e i Dt; tÞ gi ð! ! g ð x ; tÞ þ U ; i
i
(8)
where, the collision operator Ui has a definition as Ui ¼
eq gi ð! x ; tÞ gi ð! xÞ t4
(9)
In this work t4 is set to be 1.5. The equilibrium distribution function of gi is formulated as eq gi ¼ Ai þ Bi 4 þ Ci 4! e i $! u:
(10)
The coefficients in Eq. (10) are given as follows.
A0 ¼ 3Gm4 ; Ai ¼ Gm4 2ðis0Þ; B0 ¼ 1; Bi ¼ 0ðis0Þ; Ci ¼ 1=2q: (11) Similarly, the order parameter 4 is obtained from the distribution functions gi by 4¼
6 X
gi :
(12)
i¼0
In this free energy model, the pressure term in Eq. (1) is related to the surface tension, and it can be written as [36]: VP ¼ 4Vm4 Vp0 , where p0 ¼ rc2s . A free energy function in a closed volume with mixture of two fluids is adopted [36], i.e. I k r lnr ; (13) J ¼ dV jð4Þ þ ðV4Þ2 þ 2 3 where, V is a control volume; the bulk free energy density per unit mass for the homogeneous system is jð4Þ ¼ að42 4*2 Þ2 , with a being an amplitude parameter to control the free energy of interactions between the two phases; the square of the gradient term is associated with variations of the density and contributes to the free energy excess in the interfacial region; k is a parameter related to surface tension. The chemical potential is calculated by [36]
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
m4 ¼ a 443 44*2 4 kV2 4;
(14)
and the corresponding pressure tensor can be written as [36] pxy
" # 2 v4 v4 v4 ¼ pdxy k dxy ; vx vx vy
(15)
where, dxy is the Kronecker operator; the subscripts x, y represent Cartesian coordinates; p is calculated by ðV4Þ2 r þ : p ¼ a 344 24*2 42 4*4 k4V2 4 þ 3 2
Static contact angle test
(17)
18 X ui ½4ð! x þ! e i DtÞ 4ð! x Þ V 4¼2 : 2 2 ðDtÞ c i¼1 s 2
(18)
In Eq. (17), m denotes the dimension, m ¼ x, y, z for Cartesian coordinates system. To specify a contact angle in the model, a surface energy jð4s Þ for solid nodes should be added into the free energy calculation and it is assumed to be a simple linear function of 4s [38] (i.e. jð4s Þ ¼ g4s ). A natural boundary constraint for 4s is assumed [38], as ! n $ðV4s Þ ¼ g=k;
(19)
where g is the wetting potential; ! n represents the unit vector with direction pointing to the nearest or a next-nearest fluidside node from the solid node. The 4s and its second order derivatives (V2 4s ) are approximated by [25] P 1
4s ¼
s 4f;1 d v4 vn N
V2 4s y
P h ¼
1
i 4f;1 þ d g=k ; N
P s 2 1 4f;1 4s d v4 vn N
i P h 2 1 4f;1 4s þ d g=k ; ¼ N
61 61 61 lu liquid domain with periodic boundary conditions applied to all the six boundary surfaces, was simulated. We recorded the pressure difference (DP, in lattice unit) across the liquid/gas interface and the radius (R, in lattice unit) of the bubble at the final equilibrium state. Fig. 1 displays the obtained relationship between DP and R. It is obvious that DP is directly proportional to the inverse of R, perfectly amenable to Laplace's Law.
(16)
The first and second derivatives of 4 are calculated by [37] 18 v4 X ui 4ð! x þ! e i DtÞ ¼ eim ; 2 Dt vm c s i¼1
8553
A small volume of liquid on a solid surface tends to spread as film if the surface is hydrophilic or forms a droplet if the surface is hydrophobic. Following the analysis in Refs. [38], the contact angle q of liquid sticking to the solid surface is related to the wetting potential g by pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi 3 3 ~ ~ ; 1g cos q ¼ ð1=2Þ 1þg
pffiffiffiffiffiffiffiffi ~ is a dimensionless parameter, g ~ ¼ g= 2ka. where, g Considering a liquid volume placed on a flat solid surface, we simulated the evolution of the shape of this liquid volume and recorded the shape at steady state. The contact angle q was then easily obtained. Different contact angles were obtained by adjusting the wetting potential related parameters, i.e. k and a. Fig. 2 displays the obtained contact angle q as a function of wetting potential g. The simulated results agree well with the analytical data. In Fig. 2, the simulated liquid droplet shape at final steady state for three representative solid surfaces: hydrophilic (q ¼ 50 ), neutral (q ¼ 90 ) and hydrophobic (q ¼ 130 ) solid surface, are also included.
Numerical reconstruction of carbon paper GDL (20)
Toray TGP-H-060 carbon paper is a typical PEMFC GDL, which is of 0.78 porosity and 7 mm fiber diameter. Its SEM images are displayed in Fig. 3. Morphological anisotropy between the
(21)
where, d is the unit grid space; 4f,l is the order parameter on the nearest or a next-nearest neighboring fluid-side node; N represents the total number of the nearest and next-nearest fluid-side nodes.
Model validation Bubble test For a gas bubble suspended in liquid, the pressure difference (DP) across the liquid/gas interface is related to the radius R of the bubble by Laplace's Law, i.e. DP ¼
2s : R
(23)
(22)
where, s is the surface tension. To verify the model's capability of dealing with the surface tension, the dynamic evolution of a bubble volume (not in spherical shape initially) positioned at the center of a
Fig. 1 e Pressure difference across liquid/gas interface versus bubble radius.
8554
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
direction (i.e. x), the GDL is divided into two regions, one of which has lower PTFE content and the other (neighboring to the gas channel) has higher PTFE content. The sectional images displayed in Fig. 4(c) and (d) give more detailed information about PTFE distribution inside the GDLs: one uniform, the other non-uniform. The reconstructed GDL contains 69 79 79 lattice elements; one lattice element represents a 1.5 1.5 1.5 mm cube, indicating the thickness of the reconstructed GDL is about 103.5 mm.
Results and discussion Influence of PTFE content and distribution on liquid water removal
Fig. 2 e Static contact angle as a function of wetting potential.
through-plane and in-plane direction is clearly seen. Wu and Jiang [39] numerically reconstructed the 3D micro-pore structure of TGP-H-060 carbon paper with the reconstruction method proposed by Schulz et al. [40]. This method was actually adapted from a non-woven structure stochastic generation technique [41]. The present work follows the same method used in Wu and Jiang [39] to reconstruct the Toray TGP-H-060 GDL. The PTFE content in GDL is within 5e40 wt.% and 1.0e13.2% by volume fraction. [43] It is assumed in this work that the PTFE coated on the GDL fibers’ surface has uniform thickness of 0.6 mm, very small compared with the fiber diameter. When distributing the PTFE on the fiber surfaces, we first identified all the numerical elements that represent the surfaces of fibers and randomly specified a certain number of these numerical elements as PTFE. The area integral of all the PTFE elements times the thickness (i.e. 0.6 mm, assumed) of PTFE film gives the total volume of PTFE added to the carbon paper GDL. Fig. 4a shows the micrograph of the reconstructed GDL of 10.0 wt.% PTFE uniformly distributed all through the GDL volume. Fig. 4 (b) shows the reconstructed GDL overall of the same amount of PTFE, 10.0 wt.%, but non-uniformly distributed. Along the GDL thickness (through-plane)
We simulated the liquidegas flow in the computer-generated GDLs (see in Fig. 4 for 3D microstructure) with the LB model described in section (LB Modeling). The PTFE is distinguished from the GDL fibers due to a different contact angle specified. The simulated domain contains 119 79 79 lattice elements, which consists of two sub-domains: one has 69 79 79 lattice elements, and is actually the computer-generated GDL, and the other is an open or empty space. The interface (x ¼ 69 lu) of the two sub-domains mimics the GDL/GC interface. Initially, all the pore space in the GDL and the open space are filled with gas. Liquid is flowing into the GDL from the x ¼ 0 boundary surface (mimics the GDL/CL interface) at a fixed velocity, 2.0 104 lu ts1. The x ¼ 119 lu boundary surface is a fixed pressure boundary with zero pressure value. At all the other boundary surfaces of the simulated domain, the symmetric boundary condition is assumed. Note that all the variables are in the lattice units. In terms of the liquid velocity (2.0 104 lu ts1) applied at the x ¼ 0 boundary surface, the capillary number (Ca) of the fluid flow is calculated to be around 1.08 105, indicating the fluid flow approximates to the capillary fingering flow in PEMFC GDLs [44]. The multiphase LB models generally have difficulty when simulating multiphase flow with a large density ratio between phases, e.g. the air/water system considered in the present work; the simulations are vulnerable to numerical instability [34]. As discussed in Section (LB Modeling), the multiphase flow in PEMFC GDL is normally capillary-fingering flow, dominated by the capillary force. That is to say, the density ratio between phases has relatively little effects on the fluid flow. Therefore, for better numerical stability, the density ratio of the two phases is artificially set to be 1.0 for all the LB
Fig. 3 e SEM images of Toray TGP-H-060 carbon paper GDL [42].
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
8555
Fig. 4 e The computer-generated Toray TGP-H-060 carbon paper GDLs, (a) uniformly distributed PTFE; (b) non-uniformly distributed PTFE, (c) 2D sectional images for the GDL with PTFE uniformly distributed, (d) 2D sectional images for the GDL with PTFE non-uniformly distributed. Red denotes fibers and blue PTFE; in Fig. 4 (d) the gray-colored plane indicates the interface that divides the GDL into two regions of different PTFE contents. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
simulations in the present work. The same treatment to the air/ water density ratio was taken by Hao and Cheng [45] when modeling multiphase transport in PEMFC GDL by LBM. Fig. 5 illustrates the invasion process of liquid water entering into the GDL. The simulated GDL contains 10.0 wt.% PTFE uniformly distributed. The assumed contact angles for the PTFE and carbon fiber are 150 and 92 , respectively. The liquid water is seen to gradually seep into the GDL, As the GDL is hydrophobic and the flow is dominated by the capillary force, the liquid selectively occupies the larger pores firstly. When it arrives to 240,000 time steps (ts), the liquid water is seen to break through the GDL at a few preferential flow paths. Further, Fig. 6 displays the liquid distribution images on three x-sectional planes (with x ¼ 10, 38, and 66 lu, respectively) at three differing time instants: 80,000, 160,000, and
240,000 ts. The liquid water phase is clearly seen to expand toward the y and z directions, occupying local pore space when it penetrates the GDL along the x direction. Fig. 7 presents the liquid saturation profile along the GDL through-plane direction at the final steady state. We justify the onset of steady-state flow according to two criteria: one is the relative deviation of the overall liquid flow rates at the GDL/CL interface and the GDL/GC interface is lower than 0.01, the other is the relative maximum local change of the order parameter 4 in the whole GDL region at two time instants with 200Dt intervals is smaller than 0.01. The displayed liquid saturation values in Fig. 7 are x-section-averaged values at designated x positions. Two cases are considered for comparison. The two cases differ from each other due to the specified PTFE distributions. One has a GDL of 10.0 wt.%
8556
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
Fig. 5 e Illustration of the invasion process of liquid water entering into GDL. Green: liquid water; blue: fibers. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
uniformly distributed PTFE while the other case considers a GDL with non-uniformly distributed PTFE but the overall PTFE content is 10.0 wt.%, the same as the former case. For the latter case, due to the non-uniform distribution of PTFE, the GDL is divided into two regions, indicated by “A” and “B” in Fig. 7. In region A, the PTFE content is 6.5 wt.%, while in region B, it is 21.8 wt.%. For both cases, the contact angle assumed at the GDL fiber surface is 92 while at the PTFE surface is 150 . Seen from Fig. 7, if the GDL contains non-uniformly distributed PTFE, more liquid water is accumulated in region A, harmful to the liquid water removal from CL to GC. For the case with GDL of non-uniformly distributed PTFE, the region B GDL that contains more PTFE acts like a “capillary barrier”, preventing liquid water from penetrating the GDL. This raises the liquid saturation level in region A. The experimental work by Song et al. [46] found that the existence of a less hydrophobic region in the central region of GDL would impede the water removal, qualitatively according with the observations got from Fig. 7. Checking the original data used to draw the two curves in Fig. 7, we find the largest relative difference of sectionaveraged liquid saturation between the two cases is more than 10%. Four insets of liquid water distribution images on typical x-sectional planes are specially depicted in Fig. 7, which show intuitively the differences of the liquid water distribution inside GDL between the two cases.
It is seen from Fig. 7 also that at the final steady state, partial region of the GDL close to the GDL/CL interface is full of liquid water for both cases. When the liquid water flows through the GDL, it will expand toward the in-plane (i.e. y and z) directions and occupy all the pores along its flow path (refer to Fig. 6) if the flow resistance in the through-plane (i.e. x) direction is larger than the capillary force the local smallest hydrophobic pore imposes on the liquid water. We deduce from Fig. 7 that the liquid flow resistance of ~40 mm thick GDL is comparable to the entry pressure of the smallest hydrophobic pore in the GDL. To investigate the influence of PTFE content on the liquid removal, cases with GDL of 10.0 wt.%, 18.1 wt.% and 25.6 wt.% uniformly distributed PTFE are considered for comparison. Fig. 8 gives the liquid saturation profile along the GDL throughplane direction (x) at the final steady state for all the three cases. The displayed liquid saturation values are x-sectionaveraged values. The liquid saturation level is generally seen to decrease if the PTFE content in GDL is higher. Lowering PTFE content means to reduce the hydrophobicity of GDL. A less hydrophobic GDL means the capillary resistance for liquid water enters into a pore is lowered. Hence, if the GDL contains less PTFE, not only larger pores, but also more pores of smaller pore size in the GDL will be occupied by the liquid water; liquid water dwells in more pores after it enters into the GDL, raising
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
8557
Fig. 6 e Images of liquid water distribution on x-sectional planes at three differing time instants. Green: liquid water; blue: fibers. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
the liquid saturation level in the GDL. Qualitatively similar results have been reported by Hao and Cheng [45]. Increasing PTFE content in GDL makes the GDL become more hydrophobic. The liquid water needs to overcome larger capillary force to expand in the in-plane directions and occupy less relatively smaller pores. We see from Fig. 8 that a smaller part of the GDL is full of liquid water if the PTFE content in GDL is increased from 10.0 wt.% to 18.1 wt.%; for the case with GDL of 25.6 wt.% PTFE content, even no GDL x-cross section is full of liquid water except the inlet GDL/CL interface. As described in the Introduction section, there is still lack of a general cognition about the optimum PTFE content in GDL. From the simulation results presented in Fig. 8, the liquid saturation level in GDL diminishes with the increase of PTFE content in GDL. It seems that there is no an optimum PTFE content. We speculate that determining the optimum PTFE content in GDL should comprehensively consider more involved factors, not solely rely on the resultant liquid water saturation level in GDL. As an example, an additional inset in Fig. 8 gives the calculated liquid pressure drops across the GDL for all the three cases. The more hydrophobic GDL significantly elevates the entry pressure of liquid water entering into
the GDL. This may lead to the CL more vulnerable to water flooding. The LB simulation results are dependent on the lattice resolution, which was examined by the Richardson extrapolation technique [47]. For the same computer-generated GDL (microstructure displayed in Fig. 4b), two additional LB simulations with lattice resolutions of 2.5 mm and 1.0 mm respectively were performed. Considering the obtained sectionaverage liquid saturation values on the GDL/GC interface, we do Richardson extrapolation and find the extrapolated accurate value that is from the ideal case with infinite small lattice resolution deviates from the calculated value based on the simulated case with 1.5 mm lattice resolution by only about 4%. This further corroborates the numerical accuracy of the LB simulations carried out in the present work.
Influence of PTFE content and distribution on relative permeability As described in Section (Influence of PTFE content and distribution on liquid water removal), PTFE content and distribution affect the liquidegas flow, which is bound to affect
8558
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
applied in the through-plane direction. The body force applied is well controlled to guarantee the capillary number of fluid flow is within the order of 105. When the fluid flow reaches steady-state, we calculate the average liquid (or gas) velocity in the through-plane direction (ux ) and then derive the relative permeability Kr according to Darcy’ law, as Kr ¼
Fig. 7 e Liquid saturation profile along the GDL thickness direction: effect of PTFE distribution. The insets display the liquid distribution images on x-sectional planes; liquid phase is green-colored. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 8 e Liquid saturation profile along the GDL thickness direction: effect of PTFE content.
the relative permeability of each phase. A series of cases are simulated to derive the relationship between relative permeability and phase saturation. The simulated domain is solely the computer-generated GDLs (illustrated in Fig. 4); periodic boundary condition is considered for the through-plane (x) direction while symmetric boundaries are assumed for the boundary surfaces in the other two directions. Initially, a certain amount of liquid is randomly positioned in the pore space of the GDL. The fluid flow is driven by a body force (Fb)
ux h KFb
(24)
where, Fb is the applied body force; K is the absolute permeability. K is determined to be 4.05 1012 m2 from single phase LB simulations performed in the present work, slightly deviates from the range of K values measured by Mathias et al. [43], 5 1012 m2 e 10 1012 m2. Changing the initially positioned liquid amount in the GDL, a case of different liquid water saturation level in GDL is set. Fig. 9 presents the simulated steady-state liquid water distribution in GDLs of different PTFE contents. The specified contact angles for GDL fiber surface and PTFE surface are 92 and 150 , respectively. The two cases considered have the same phase distribution initially, i.e. the same overall liquid water saturation (sl ¼ 0.42) in the GDL. As evident by Fig. 8, the two cases give significantly different results; the liquid water distribution for the case with GDL of higher PTFE content is more fragmented at the final steady-state, simply due to the more hydrophobic GDL. Fig. 10 presents the derived relationship of relative permeability versus phase saturation. Results from two sets of cases are included in Fig. 10 for comparison. The two sets of cases differ from each other due to the PTFE content. For all the cases, the specified contact angle on GDL fiber surface is 92 , while the contact angle specified on PTFE surface is 150 . It is seen that the relative permeability increases with the increasing phase saturation for the liquid and gas phases both. A natural explanation for this is a higher phase saturation level makes the corresponding fluid phase more easily form connective flow paths all-through the GDL. Inspection on the results of Krl versus sl shown in Fig. 10 finds when the sl is less than ~0.5, the Krl is larger in the GDL of lower PTFE content, whereas when sl is more than ~0.5, the Krl is larger in the GDL of higher PTFE content. There may exist two underlying mechanisms [48]: i) liquid flow in the GDL of higher PTFE content is more vulnerable to fragmentation due to the higher hydrophobicity of GDL; ii) the liquid in more hydrophobic GDL prefers dwelling in larger pores, which have smaller specific surface area interfacing with the solid GDL fibers, leading to smaller liquid flow resistance. It is speculated that the former mechanism dominates the liquid flow when the liquid saturation level in GDL is low. However, when the liquid saturation level in GDL is sufficiently high (say, >0.5), most of the pore space in the GDL is occupied by the liquid phase and there should always have some connective liquid flow paths all-through the GDL. Hence, the latter mechanism takes over the dominating role when the liquid saturation level in GDL is high enough, i.e. sl > 0.5 in the present work. Seen from Fig. 10 also, the PTFE content in GDL has influence on the relationship of Krg versus sg, the Krg in the GDL of lower PTFE content is larger than the one in the GDL of higher PTFE content. This can be explained by the improved
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
8559
Fig. 9 e Simulated liquid water distribution in GDL of (a) 10.0 wt.% and (b) 18.1 wt.% PTFE content. For both cases the overall liquid saturation in GDL is 0.42; on the 2D section cuts, green denotes liquid water. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
connectivity for the gas phase in the GDL of lower PTFE content according to the above-described mechanism i). Both the liquid and gas phase relative permeability values decrease to zero when the phase saturation reduces to a certain low level. For the cases with more hydrophobic GDL, the Krg is seen to decrease to zero when the gas phase saturation is lowered down to about 0.2. Qualitatively similar results were reported in Refs. [26 and 27]. For multi-phase flow in porous medium, if one phase is in very low saturation level, it can be totally fragmented and there is thus no possibility for this phase to form an effective all-through flow path and then gets stuck in the porous medium. This is to some extent similar to the concept of “percolation threshold” of porous media [47]. The inset in Fig. 10 summarizes the relative permeability values predicted by LB models available in open literature. Refs. [25 and 26] used the free energy LB model and the GDL contact angles were considered to be 105 and 145 , respectively; Ref. [27] used the SeC model and 135 GDL contact angle was assumed. The results from different literature vary in similar
trend though differences exist. All the literature data qualitatively agree with the results obtained in the present work. Fig. 11 compares the derived Kr-s relationships for GDL of uniformly distributed PTFE with those for GDL of nonuniformly distributed PTFE. The PTFE content and distribution for both the two sets of cases considered here have been detailed in Section (Influence of PTFE content and distribution on liquid water removal) in relation with Fig. 7. For the GDL of non-uniform PTFE distribution, the region with higher PTFE content acts like a “capillary barrier”, impeding the liquid flow. Seen from Fig. 11, the non-uniform PTFE distribution in GDL does reduce the Krl value when the sl is not more than ~0.55, but with the increase of sl (i.e., when sl is greater than ~0.55), the Krl in the GDL with non-uniform PTFE distribution exceeds the value of Krl in the GDL with uniform PTFE distribution. This can also be explained by the aforementioned two mechanisms. The PTFE distribution shows relatively smaller effect on the Krg; the Krg in the GDL of uniform PTFE distribution is slightly larger than that in the
8560
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
Genuchten function particularly when the phase saturation is low; the van Genuchten function underestimates the values of Krl especially when the non-wetting liquid phase saturation takes an intermediate value.
Summary and conclusions
Fig. 10 e Effects of PTFE content on the liquid and gas relative permeability of GDL.
Fig. 11 e Effects of PTFE distribution on the liquid and gas relative permeability of GDL.
GDL of non-uniform PTFE distribution, indicating the nonuniform distribution of PTFE in GDL may somehow prevents the gas phase transport. Suitable measures taken to avoid non-uniform distribution of PTFE in GDL may relieve the mass transport limitation of PEMFC at high current density operations. In Figs. 10 and 11, the LBM-calculated Kr-s data for both wetting and non-wetting phases are also compared with the van Genuchten function [28], which is formulated as 2ðj1Þ=j ðj1Þ=j ; Kr ¼ s1=3 1 ð1 sÞ
(25)
where the value of j ranges from 2 to 5. It is seen that the Kr-s data for the wetting gas phase generally better fit the van
A 3D multiphase LB model was developed and employed to study the impact of PTFE content and distribution on the liquidegas transport in PEMFC carbon-paper GDLs. The capability and reliability of the developed LB model were first tested by simulating the dynamic evolution of a gas bubble suspended in liquid and of a liquid volume put on a solid surface. The simulation results accord very well with the corresponding analytical results. The 3D micro-pore structure of TGP-H-060 carbon paper was numerically reconstructed and part of the fibers’ surface was randomly specified as PTFE. Two classes of cases were considered. One has PTFE uniformly distributed all through the GDL volume. The other has non-uniformly distributed PTFE: in the GDL sub-region neighboring to the GC, the PTFE content is higher, which mimics the PTFE accumulation toward the GDL/GC interface due to incomplete penetration of PTFE into GDL in practical PTFE treatment processes. The simulated liquid flow in GDL shows evident capillary force dominated characteristics, gradually seeping into the GDL and selectively occupying larger pores first. A GDL with higher PTFE content is more hydrophobic. The GDL sub-region containing more PTFE acts like a “capillary barrier”, which prevents liquid water from passing the GDL and raises the liquid saturation level in the remaining GDL region adjoining to the CL. The simulation results about liquid saturation level in GDL do not support that there exists an optimum PTFE content in GDL. Three cases with GDL of differing PTFE contents, 10.0 wt.%, 18.1 wt.% and 25.6 wt.%, respectively, were simulated. The simulated liquid saturation level in GDL diminishes with the increase of PTFE content in GDL. It is found from simulation results also that the liquid pressure drop across the GDL is higher if the case has more PTFE in GDL. We deduce determining the optimum PTFE content in GDL must consider more involved factors, and should not rely solely on the resultant liquid saturation level in GDL. LB simulations enable the calculation of relative permeability of liquid and gas phase in the carbon paper GDL. We derived the relationship of relative permeability versus phase saturation for both phases from a series of additional simulation cases. The relative permeability is found to increase with the increasing phase saturation for both phases. PTFE content and distribution are found to have important influence on the relationship of liquid phase relative permeability versus phase saturation. If the (overall) liquid saturation is low (say, lower than ~0.5 for the cases considered in the present work), the liquid phase relative permeability is larger in the GDL of lower PTFE content, whereas it is smaller if the liquid saturation is high (say, higher than ~0.5). The liquid phase relative permeability in the GDL with non-uniform PTFE distribution is found to be reduced when the liquid saturation in GDL is low (say, lower than ~0.55), but when the liquid saturation is high (say, higher than ~0.55), it is larger compared with that in the GDL
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
with uniform PTFE distribution. Two underlying mechanisms take effects here: i) liquid flow in the GDL of higher PTFE content is more vulnerable to fragmentation due to the higher hydrophobicity of GDL; ii) the liquid in more hydrophobic GDL prefers dwelling in larger pores, which have smaller specific surface area interfacing with the solid GDL fibers, leading to smaller liquid flow resistance. Mechanism i) dominates the liquid flow when the liquid saturation level in GDL is low, while mechanism ii) takes over the dominating role when the liquid saturation level in GDL is sufficiently high. The PTFE content and distribution in GDL also show influence on the relationship of gas phase relative permeability versus gas saturation. The gas phase relative permeability in the GDL of lower PTFE content is found to be larger than that in the GDL of higher PTFE content and the non-uniform distribution of PTFE in GDL is found to slightly reduce the gas phase relative permeability value. It is found as well that the LBM-predicted relationship of gas phase relative permeability versus phase saturation generally better accord with the van Genuchten function particularly when the phase saturation is low, whereas the van Genuchten function underestimates the liquid phase relative permeability especially when the liquid phase saturation takes an intermediate value.
Acknowledgments Financial support received from the Key Scientific Development Project of Guangdong Province (2015A030308019), the Guangzhou Scientific and Technological Development Plan (2014J4100217), and the CAS “100 talents” Program (FJ) is gratefully acknowledged.
Nomenclature a, k Ai, Bi, Ci, Bo Ca cs d D ! ei fi eq fi Fb g gi eq
gi
j K Kr lm lu M P
parameters related to surface tension Di coefficients in equilibrium distribution function Bond number Capillary number lattice speed of sound, lu ts1 unit grid space, lu characteristic length, m particle speed, lu ts1 density distribution function in the ith direction equilibrium density distribution function body force, lm lu2 ts2 gravity, m/s2 distribution function for the order parameter in the ith direction equilibrium distribution function of the order parameter exponent in the van Genuchten function absolute permeability, m2 relative permeability unit of mass in LBM unit of length in LBM mobility of fluid pressure of fluid, lm lu1 ts2
p0 p Q q R Re s T ts u ! u V ! x
8561
thermodynamic pressure, lm lu1 ts2 bulk pressure, lm lu1 ts2 flow rate, lm lu2 ts1 constant coefficient related to the relaxation time t4 droplet radius, lu Reynolds number saturation of a fluid phase the rate of the order parameter variation; time, ts unit of time in LBM velocity, lu ts1 velocity vector, lu ts1 controlling volume, lu3 spatial position
Greek symbols G parameter used to control the fluid mobility g wetting potential ~ g a dimensionless parameter in relation with Eq. (23) Kronecker operator dxy h viscosity, lm lu2 ts1 m4 chemical potential r density, lm lu3 s surface tension, lm ts2 tr, t4 dimensionless relaxation times, ts y kinematic viscosity, m2/s 4 order parameter 4* order parameter in equilibrium state j bulk free energy density J free energy function weight coefficient ui collision operator in the ith direction Ui Subscripts/superscripts x the direction x in Cartesian coordinates y the direction y in Cartesian coordinates l liquid g gas i the ith discrete distribution direction s solid in inlet out outlet
references
[1] Pasaogullari U, Wang CY. Liquid water transport in gas diffusion layer of polymer electrolyte fuel cells. J Electrochem Soc 2004;151:399e406. [2] Wang CY. Fundamental models for fuel cell engineering. Chem Soc Rev 2004;104:4727e66. [3] Jiang FM, Wang CY. Numerical modeling of liquid water motion in a polymer electrolyte fuel cell. Int J Hydrogen Energy 2014;39:942e50. [4] He W, Yi JS, Nguyen TV. Two-phase flow model of the cathode of PEM fuel cells using interdigitated flow field. Am Inst Chem Eng 2000;46:2053e64. [5] Natarajan D, Nguyen TV. A two-dimensional, two-phase, multicomponent, transient model for the cathode of a proton exchange membrane fuel cell using conventional gas distributors. J Electrochem Soc 2001;148:1324e35.
8562
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 1 ( 2 0 1 6 ) 8 5 5 0 e8 5 6 2
[6] Wu H, Li XG, Berg P. On the modeling of water transport in polymer electrolyte membrane fuel cells. Electrochimica Acta 2009;54:6913e27. [7] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39:201e25. [8] Theodorakakos A, Ous T, Gavaises M. Dynamics of water droplets detached from porous surfaces of relevance to PEM fuel cells. J Colloid Interface Sci 2006;300:673e87. [9] Zhu X, Sui PC, Ned D. Three-dimensional numerical simulations of water droplet dynamics in a PEMFC gas channel. J Power Sources 2008;181:101e15. [10] Park JW, Jiao K, Li XG. Numerical investigations on liquid water removal from the porous gas diffusion layer by reactant flow. Appl Energy 2010;87:2180e6. [11] Aimy B, David S, Ned D. Dynamic water transport and droplet emergence in PEMFC gas diffusion layers. J Power Sources 2008;176:2180e6. [12] Mohamed EA, Ben A, Sass BN. Numerical simulation of droplet dynamics in a proton exchange membrane (PEMFC) fuel cell micro-channel. Int J Hydrogen Energy 2015;40:1333e42. [13] Park J, Li X. Multi-phase micro-scale flow simulation in the electrodes of a PEM fuel cell by lattice Boltzmann method. J Power Sources 2008;178:248e57. [14] Gao Y, Zhang XX, Rama P, Chen R, Ostadi H, Jiang K. Lattice Boltzmann simulation of water and gas flow in porous gas diffusion layers in fuel cells reconstructed from microtomography. Comput Math Appl 2013;65:891e900. [15] Chen L, Luan HB, He YL, Tao WQ. Numerical investigation of liquid water transport and distribution in porous gas diffusion layer of a proton exchange membrane fuel cell using lattice Boltzmann method. Russ J Electrochem 2012;48:712e26. [16] Bevers D, Rogers R, Bradke MV. Examination of the influence of PTFE coating on the properties of carbon paper in polymer electrolyte fuel cells. J Power Source 1998;63:193e201. [17] Lin GY, Nguyen TV. Effect of thickness and hydrophobic polymer content of the Gas diffusion layer on electrode flooding level in a PEMFC. J Electrochem Soc 2005;152:1942e8. [18] Kandlikar SG, Garofalo ML, Lu Z. Water management in a PEMFC: water transport mechanism and material degradation in gas diffusion layers. Fuel Cells 2011;11:814e23. [19] Wua JF, Martina JJ, Orfino FP, Wang HJ, Legzdins C, Yuan XZ, et al. In situ accelerated degradation of gas diffusion layer in proton exchange membrane fuel cell part I: effect of elevated temperature and flow rate. J Power Sources 2010;195:1888e94. [20] Lee C, Merida W. Gas diffusion layer durability under steadystate and freezing conditions. J Power Sources 2007;164:141e53. [21] Paganin VA, Ticianelli EA, Gonzalez ER. Development and electrochemical studies of gas diffusion electrodes for polymer electrolyte fuel cells. J Appl Electrochem 1996;26:297e304. [22] Lim C, Wang CY. Effects of hydrophobic polymer content in GDL on power performance of a PEM fuel cell. Electrochimica Acta 2004;49:4149e56. [23] 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:8401e9. [24] Gostick JT. Random pore network modeling of fibrous PEMFC gas diffusion media using voronoi and delaunay tessellations. J Electrochem Soc 2013;160:731e43. [25] Niu XD, Munekata T, Hyodo SA, Suga K. An investigation of water-gas transport processes in the gas-diffusion-layer. J Power Sources 2007;172:542e52. [26] Hao L, Cheng P. Pore-scale simulations on relative permeabilities of porous media by lattice Boltzmann method. Int J Heat Mass Transf 2010;53:1908e13.
[27] Koido T, Furusawa T, Moriyama K. An approach to modeling two-phase transport in the gas diffusion layer of a proton exchange membrane fuel cell. J Power Sources 2008;175:127e36. [28] van Genuchten MT. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J 1980;44:892e8. [29] Gunstensen AK, Rothman DH, Zaleski S, et al. Lattice Boltzmann models for nonideal gases. Phys Rev E 2000;62:4982. [30] Shan, X.; Chen, H, Lattice Boltzmann model for simulating flows with multiple phases and components. Phys Rev E 1993; 47; 1815. [31] Shan, X.; Chen, H, Simulation of nonideal gases and liquidgas phase transitions by lattice Boltzmann equation. Phys Rev E 1994; 49; 2941. [32] Swift MR, Osborn WR, Yeomans JM. Lattice Boltzmann simulations of nonideal fluids. Phys Rev Lett 1995;75:830. [33] Swift MR, Orlandini E, Osborn WR, Yeomans JM. Lattice Boltzmann simulations of liquid-gas and binary fluid systems. Phys Rev E 1996;54:5041. [34] Huang HB, Wang L, Lu XY. Evaluation of three lattice Boltzmann models for multiphase flows in porous media. Comput Math Appl 2011;61:3606e17. [35] Zheng HW, Shu C, Chew YT. A lattice Boltzmann model for multiphase flows with large density ratio. J Comput Phys 2006;218:353e71. [36] Kendon VM, Gates ME, Pagonabarraga I, Desplat JC, Bladon P. Inertial effects in three-dimensional synodal decomposition of a symmetric binary fluid mixture: a lattice Boltzmann study. J Fluid Mech 2001;440:147e203. [37] Lee T, Lin CL. A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio. J Comput Phys 2005;206:16e47. [38] Briant AJ, Papatzacos P, Yeomans JM. Lattice Boltzmann simulations of contact line motion in a liquidegas system. Philos Trans R Soc Lond A 2002;360:485e95. [39] Wu W, Jiang FM. Microstructure reconstruction and characterization of PEMFC electrodes. Int J Hydrogen Energy 2014;39:15894e906. [40] Schulz VP, Becker J, Wiegmann A, Mukherjee PP, Wang CY. Modeling of two-phase behavior in the gas diffusion medium of PEMFCs via full morphology approach. J Electrochem Soc 2007;154:419e26. [41] Katja S, Stefanie P, Doris RB, Andreas W, Joachim O. Design of acoustic trim based on geometric modeling and flow simulation for non-woven. Comput Mater Sci 2006;38:56e66. [42] Bazylak A, Sinton D, Liu ZS, Njilali N. Effect of compression on liquid water transport and microstructure of PEMFC gas diffusion layers. J Power Sources 2007;163:784e92. [43] Mathias MF, Roth J, Feleming J, Lehnert W. In: Lietsich W, Lamm A, Gasteiger HA, editors. In hand book of fuel cells: fundamentals, technology and applications, Vol. 3. Chichester: John Wiley & Sons; 2003. Part 1. [44] Sinha PK, Wang CY. Pore-network modeling of liquid water transport in gas diffusion layer of a polymer electrolyte fuel cell. Electrochimica Acta 2007;52:7936e45. [45] Hao L, Cheng P. Lattice Boltzmann simulations of water transport in gas diffusion layer of a polymer electrolyte membrane fuel cell. J Power Sources 2010;195:3870e81. [46] Song W, Yu HM, Shao ZG, Yi BL, Lin J, Liu N. Effect of polytetrafluoroethylene distribution in the gas diffusion layer on water flooding in proton exchange membrane fuel cells. Chin J Catal 2014;35:468e73. [47] Jiang FM, Oliveira MSA, Sousa ACM. Mesoscale SPH modeling of fluid flow in isotropic porous media. Comput Phys Commun 2007;176:471e80. [48] Li HN, Pan CX, Miller CT. Pore-scale investigation of viscous coupling effects for two-phase flow in porous media. Phys Rev 2005;72:026705.