- Email: [email protected]
Two-phase flow in compressed gas diffusion layer: Finite element and volume of fluid modeling
Journal of Power Sources 437 (2019) 226933
Contents lists available at ScienceDirect
Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour
Two-phase flow in compressed gas diffusion layer: Finite element and volume of fluid modeling Xia Zhou, Zhiqiang Niu, Zhiming Bao, Jingchao Wang, Zhanrui Liu, Yan Yin, Qing Du, Kui Jiao * State Key Laboratory of Engines, Tianjin University, 135 Yaguan Rd, Tianjin, 300350, China
H I G H L I G H T S
� Microstructures of compressed GDL are digitally reconstructed. � Pore-scale two-phase flow in compressed GDL is investigated using VOF model. � Effects of compression on water transport in GDL are discussed. � Quantitative correlations between Pc and s are concluded. A R T I C L E I N F O
A B S T R A C T
Keywords: Gas diffusion layer Compression Capillary pressure Water saturation
In this study, a stochastic model is used to reconstruct the uncompressed gas diffusion layer (GDL) micro structures. Subsequently, the finite element method (FEM) is conducted for assembly pressure simulation to generate the compressed GDL microstructures. The effects of assembly pressure on GDL deformation are investigated. It is found that assembly pressure causes non-uniform deformation of the GDL along the thickness direction. Finally, a volume of fluid (VOF) model is developed to investigate two-phase flow in the compressed GDL. The results show that when the capillary pressure is higher than 4 kPa, the water saturation decreases as the compression ratio increases. But when the capillary pressure is below 3 kPa, compression has little effect on water saturation. Based on the above findings, three regions namely weak deformation region (WDR), moderate deformation region (MDR), and strong deformation region (SDR) are defined. Impacts of compression on water saturation differ in these three regions. Moreover, compression increases the pressure of water breakthrough, but has minor effects on preferential pathways of water breakthrough. Quantitative correlations between water saturation and capillary pressure in the uncompressed and compressed GDL microstructures are also concluded.
1. Introduction Proton exchange membrane fuel cells (PEMFCs) are regarded as one of the potential power sources for vehicles due to their outstanding merits such as high conversion efficiency and zero emission [1]. Improvement of water transport in PEMFCs is essential for fuel cell commercialization [2]. To achieve maximum removal of liquid water, it is important to understand water transport behavior in the gas diffusion layer (GDL) [3]. The GDL is a vital component in PEMFCs and has multiple functions such as providing transport paths for reactant and product liquid water [4]. Deep understanding of liquid water transport in GDL is important to enhance the cell performance. Numerous studies have been conducted to investigate water transport in the gas diffusion layers (GDLs). Liu et al. [5] elaborated the state and transport mechanism of water inside the uncompressed GDLs. It is worth
noting that GDLs in real PEMFCs are usually compressed due to the clamping pressure, which is always ignored in previous studies. Compression decreases the thickness, porosity and electrical resistance, which in turn affects the performance of PEMFCs [6], as shown in Fig. 1a. Several experimental techniques have been adopted to visualize water transport in the compressed GDLs such as fluorescence micro scopy [7], neutron radiography [8,9] and X-ray tomography [10–17]. Bazylak et al. [7] adopted fluorescence microscopy to provide ex situ visualization of water transport through the compressed GDL and found areas under compression coincide with preferential pathways for water transport. Cooper et al. [9] adopted neutron radiography to examine the through-plane liquid water distribution of an operating PEMFC. Ince et al. [10] investigated water transport in two GDL samples by means of synchrotron X-ray tomography. They found pore saturation in the compressed GDL is higher than that in the uncompressed GDL. Zenyuk
* Corresponding author. E-mail address: [email protected] (K. Jiao). https://doi.org/10.1016/j.jpowsour.2019.226933 Received 11 April 2019; Received in revised form 22 July 2019; Accepted 24 July 2019 Available online 31 July 2019 0378-7753/© 2019 Elsevier B.V. All rights reserved.
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Fig. 1. Schematics of FEM model and VOF model based on GDL microstructures in PEMFC: (a) Morphologies of uncompressed and compressed GDL; (b) the computational domain and boundary conditions used for validation of the FEM model; (c) the computational domain used for formal cases of the FEM model; (d) the computational domain and boundary conditions used for validation and formal cases of the VOF model.
et al. [17] found uniform liquid-water front at low compression, while at high compression, the water predominantly advanced at locations under the channel for higher liquid pressure. Although these visualization techniques deliver some fundamental understanding of liquid water transport in the compressed GDL, they are usually time-consuming and have high requirement for facilities. As an alternative, numerical modeling can investigate water trans port inside the GDL efficiently. Several two-phase flow models have been conducted to discuss water dynamics inside uncompressed and compressed GDLs, such as pore network model (PNM), lattice Boltz mann method (LBM), and volume of fluid (VOF) method. Fazeli et al. [18] employed PNM method extracted from synchrotron-based CT im ages of compressed GDL to simulate liquid water transport in the GDL. Satjaritanun et al. [19] conducted the LBM model to predict the breakthrough pressure for GDL samples. Le et al. [20] employed a VOF model to investigate liquid water in PEMFC with serpentine channels. Niu et al. [21] reconstructed a GDL with the experimentally-determined varying porosity, and employed a three-dimensional (3-D) VOF model to investigate two-phase flow inside the GDL. Bao et al. [22] developed a VOF model to investigate water transport behaviors inside a simplified compressed GDL, and found that deformation highly influenced water droplet dynamics. Typical researches on water transport inside un compressed and compressed GDLs using PNM, LBM, and VOF methods are summarized in Table 1. Although many researches have discussed water transport behaviors in GDLs, only a relatively small percentage of the studies [18,20,51–53] have considered compressed GDL microstructures, as shown in Table 1. However, there are still several limitations in these literatures [18,20, 51–53]: 1). Universal correlation between water saturation and
capillary pressure inside compressed GDLs has not been quantitatively concluded yet; 2). Clear effects of compression on two-phase flow in GDLs have not been concluded; 3). The compressed GDL microstructures adopted in Refs. [18,20,51–53] were obtained either through the complicated high spatial resolution X-ray CT or simplified structures. In this study, to efficiently investigate pore-scale liquid water transport behaviors in compressed GDLs, a stochastic algorithm [49,54, 55] is employed for reconstructing seven GDL microstructure samples. A FEM model [64] is developed to study effects of clamping pressure on GDL deformation. The two-phase VOF model [21,22,47–50] is then employed to simulate water transport behaviors in the compressed GDL. Both FEM and VOF models are validated against the experimental lit eratures [56,57]. To the best of the authors’ knowledge, this is the first study which integrates FEM model with VOF model to investigate water transport behaviors in compressed GDL microstructures. The uniqueness of this study shows as follows: (1) The compressed three-dimensional (3D) GDL microstructures are considered; (2) Pore-scale two-phase flow in the compressed GDL is investigated; (3) Correlations between water saturation and capillary pressure in the GDL microstructures with various compression ratios are quantitatively concluded. 2. Model development 2.1. FEM model 2.1.1. Model assumptions In this study, the FEM model is implemented to simulate the defor mation of the GDL microstructures. Then the compressed microstruc tures are captured for the VOF model. The computational domains are 2
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Table 1 Features of LBM, PNM, and VOF simulations conducted to investigate water transport behaviors inside GDLs. Type
Reference and Published year
Uncompressed GDL LBM Molaeimanesh et al. [23] (2014) Molaeimanesh et al. [24] (2014) Chen and Jiang [25] (2016) Molaeimanesh and Akbari [26] (2016) Safi et al. [27] (2017) Satjaritanun et al. [28] (2017) Gao et al. [29] (2017) Sakaida et al. [30] (2017) Shakerinejad et al. [31] (2018) Yu et al. [32] (2018) Kakaee et al. [33] (2018) Jinuntuya et al. [34] (2018) Zhang et al. [35] (2018) Xu et al. [36] (2018) Yu et al. [37] (2018) Fang et al. [38] (2018) Liu et al. [39] (2019) PNM Lee et al. [40] (2014) Shahraeeni and Hoorfar [41] (2014) Fazeli et al. [42] (2015) Agaesse et al. [43] (2016) Qin et al. [44] (2016) Straubhaar et al. [45] (2016) Huang et al. [46] (2019) VOF Chen et al. [47] (2013) Ahmad et al. [48] (2013) Yin et al. [49] (2014) Niu et al. [21] (2018) Niu et al. [50] (2018) Compressed GDL LBM Jeon and Kim [51] (2015) Satjaritanun et al. [52] (2018) Molaeimanesh et al. [20] (2019) PNM Fazeli et al. [18] (2016) VOF Bao et al. [22] (2014) Han and Chen [53]. (2019)
Research Objective
Microstructure reconstruction
Effects of PTFE distribution Effects of wettability and it spanwise and transverse gradient Effects of PTFE content and distribution Effects of wettability and water droplet size Flow and mass diffusion in a binary air/water mixture Effects of MPL existence Effects of porosity and thickness Effects of hydrophobic-hydrophilic pattern Effects of wettability of GDL Dynamic process of liquid water breaking through Effects of 3D PTFE distribution and the binder content Effects of contact angles and pressure differences Coupling two-phase flow and electrochemical reaction Effects of MPL with different thickness and porosity Effects of asymmetric droplet shape Effects of perforation The heterogeneous multi-transport characteristics in GDLs Effects of width of the rib, thickness of the GDL Effects of GDL hydrophobicity and thickness Effects of liquid water inlet boundary conditions Comparation with experiment Effects of MPL Effects of condensation Effects of network structure (uniform and gradient porosity) Effects of GDL surface Effects of PTFE loading Effects of variable contact angle Comparison with Experiment Effects of mixed wettability
Stochastic Stochastic Stochastic Stochastic X- ray CT X-ray CT Stochastic Stochastic Stochastic Stochastic Stochastic X-ray CT Stochastic Stochastic Stochastic Stochastic Stochastic Simplified cubic network Simplified cubic network X- ray CT X-ray CT Simplified cubic network Simplified cubic network X-ray CT Simplified fiber ridges Simplified cylindrical fibers Stochastic Stochastic Stochastic
Effects of compression ratio Water evolution, water saturation, and breakthrough pressure Removal of a droplet from GDL Effects of compression value Effects of position of droplet, contact angles Effects of air velocity and wettability
2D-simplified Stochastic X-ray CT 2D-simplified Stochastic X-ray CT Simplified Plane Simplified fibers (2D)
composed of a bipolar plate (BP), a GDL, and a buffer, as shown in Fig. 1b and c. The deformation of the catalyst layer (CL) or the mem brane layer (ML) is neglected in this study, because the elastic modulus of either CL or ML is much larger than GDL [58]. In addition, it is assumed that both cathode and anode of the PEMFC have same me chanical behaviors. This model only considers the elastic deformation of the GDL microstructures, and GDL damage due to compression is neglected. Herein, a half domain (one bipolar plate and one gas diffusion layer) is chosen to simulate GDL deformation. The normal displacement along z-axis refers to the deformation of the GDL under the rib. A stochastic algorithm [49,54,55] is employed to reconstruct the fibrous GDLs. The process of GDL reconstruction is based on the following assumptions [50]: (a) The fibers are straight cylinders with the constant diameter, and they are placed in the X–Y plane randomly; (b) The fibers orientation is perpendicular to the Z-axis direction; (c) The overlap of the fibers is allowed and PTFE treating is ignored.
where σi , σj and σ j are the normal stresses, τij ,τik and τjk are the shear stresses, Fi ,Fj and Fk are the body forces. 8 > > > > > > > > > 2 > ∂2 γij ∂ εi ∂2 εj > > > þ ¼0 > 2 2 > ∂Xi ∂Xj ∂Xj ∂Xi > > > > > > > > > > ∂2 γjk ∂2 εj ∂2 εk > > þ ¼0 > 2 2 > > ∂Xj ∂Xk ∂Xk ∂Xj > > > > > > > > > ∂2 εk ∂2 εi ∂2 γki > > þ ¼0 > 2 2 > > ∂Xk ∂Xi ∂Xk < ∂Xi (2) > > > � � > 2 > ∂γjk > ∂ ∂γij ∂γ ∂ εi > > þ ik ¼2 > > ∂ X ∂ X ∂ X ∂ X ∂ X > i k j i j ∂Xk > > > > > > � � > > ∂ ∂γji ∂γjk ∂γik ∂ 2 εj > > ¼2 > > ∂Xj ∂Xk þ ∂Xi > ∂Xj ∂Xi ∂Xk > > > > > > � � > > ∂γkj ∂γij ∂2 εk > ∂ ∂γki > > þ ¼2 > > ∂ X ∂ X ∂ X ∂ X ∂ X > k j i k i ∂Xj > > > > > > :
2.1.2. Governing equations This study considers the elastic deformation of the GDL, which is governed by the equilibrium and compatibility equations with the Hooke’s Law. The constitutive relations are as follows [64]: 8 ∂σ ∂τ ∂τ > > > i þ ij þ ik þ Fi ¼ 0 > > ∂Xi ∂Xj ∂Xk > > > > < ∂σ ∂τji ∂τjk j þ þ þ Fj ¼ 0 (1) > ∂Xj ∂Xi ∂Xk > > > > > ∂σk ∂τki ∂τkj > > > : ∂X þ ∂X þ ∂X þ Fk ¼ 0 k i j
where εi , εj and εj are the normal strains, γij ,γik and γjk are the shear strains. Normal stresses are calculated by:
3
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
8 > > σi ¼ λe þ 2Gεi > > < σ ¼ λe þ 2Gε j j σk ¼ λe þ 2Gεk > > > > : γi j ¼ τij =G
(3)
(4)
e ¼ εi þ εj þ εk λ¼
G¼
υE ð1 þ υÞð1
density (ρ) and dynamic viscosity (μ) for the two-phase mixed fluid are determined as follows [50,65]:
E 2ð1 þ υÞ
(6)
ρ ¼ ρl α þ ρg ð1
αÞ
(9) (10)
αÞ
! r⋅ U ¼ 0
(11)
Phase conservation equation (time-dependence of α):
∂α ! ! þ r⋅ðU αÞ þ r⋅½U r αð1 ∂t
2.1.3. Numerical procedures The FEM model is solved in the commercial software ANSYS APDL. The conservation equations are solved by preconditioned conjugate gradient (PCG) method. To validate the FEM model, the GDL with thickness 370 μm and porosity 0.78 is reconstructed. The variables involved in this model must meet the convergence residual less than 10 6. For the formal samples (sample1-sample7), the original thickness of the GDL is 190 μm with a constant porosity 0.84. All the reconstructed GDLs have a constant carbon fiber diameter 8 μm. Hexahedral meshes with grid size 1.5 μm are used to discretize the GDL microstructures, BP and buffers, as presented in Fig. 1b. The number of grids adopted is 3.10 million for the formal cases and 5.17 million for the validation case, as shown in Table 2.
(12)
αÞ� ¼ 0
! where the symbol U represents the effective velocity vector which is ! shared by both two phases throughout the process. U r named ‘compression velocity’ indicates the relative velocity of liquid and gas at ! ! the interface and subscript r denotes “relative velocity”. Both U and U r are calculated respectively by: ! ! U ¼ α U l þ ð1 ! ! Ur ¼ Ul
!
(13)
αÞU g
! Ug
(14)
Momentum equation: ! ∂ðρ U Þ !! þ r⋅ðρ U U Þ ∂t
2.2. VOF model
! r⋅ðμrU Þ
! ðrU Þ⋅rμ ¼
rpd
! g ⋅! x rρ þ σ krα (15)
where α, σ and k denote the phase fraction, surface tension coefficient and mean curvature of the phase interface respectively. pd is the modi fied pressure for simplifying boundary conditions, and it is determined by:
2.2.1. Model assumptions In this study, some assumptions in the VOF model show as follows [21,50]: 1) Two-phase flows inside the GDL are ideal laminar and incompress ible flows; 2) Transport properties of all fluids keep stable; 3) There is no phase change throughout the process.
pd ¼ p
(16)
ρ! g ⋅! x
here ! g indicates the gravity vector and ! x represents the position vector. The subscript d means “dynamic”. The continuum surface force (CSF) model is implemented in this model to take the effects of the surface tension at the gas-liquid interface into account by adding a force fσ to Eq. (15), which is defined as follow
2.2.2. Governing equations VOF model is a common surface tracking method. In VOF, water phase fraction α is introduced as the main variable for the model to solve. The regions full of liquid water are marked as α ¼ 1, while those filled with air are marked as α ¼ 0. Phase fraction between 0 and 1 means the air-liquid water interface. Vl Vp
(8)
where subscripts l and g indicate liquid water and gas respectively. VOF model governing equations [21,67]: Continuity equation:
where λ is the Lame constant, G is the shear elastic modulus, E is the Young’s modulus, υ is the Poisson’s ratio. i,j and k represent the x,y and z coordinates, respectively.
α¼
αÞ
μ ¼ μl α þ μg ð1
(5)
2υÞ
p ¼ pl α þ pg ð1
(17)
fσ ¼ σkrα
where k which is mentioned is the mean curvature of the phase interface k is calculated by: � � rα k ¼ r⋅! n ¼ r⋅ (18) jrαj
(7)
where Vl is the water volume, Vp is the original volume of the pores. The physical properties such as the pressure (p), averaged volume
the surface unit normal ! n is adjusted in the cells close the wall ac
Table 2 Grid size, total number of grids, and different water inlet pressures for the validation and formal cases of FEM model and VOF model. Cases
Compressed By (%)
Model
Grid Size (μm � μm � μm)
Total number of Grids (Million)
Water inlet pressure (kPa)
Validation Cases Formal Cases (samples 1-7)
0 (CPN)
FEM VOF FEM VOF VOF
1.5 � 1.5 � 1.5 2.0 � 2.0 � 1.0 1.5 � 1.5 � 1.5 2.0 � 2.0 � 1.0 2.0 � 2.0 � 1.0
5.17 3.59 3.10 4.59 4.05 3.71 3.37 3.06
1.0 [57] 1.0–7.4
0 (CPN) 8.8 17.7 26.6 35.9
4
1.0–8.8
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
cording to the equation as follow
3. Results and discussion
! ! n ¼! n w cos θ þ t w sin θ
(19)
3.1. Model validation
! where ! n w and t w are the unit vector normal to the wall and unit vector tangential to the wall, the subscript w indicates the wall. Symbol θ is the contact angle, and it is assumed constant here. Further information on the VOF approach can be found in the Refs [50,65,66].
To validate the FEM model, the compression curve for the 370 μm thick GDL sample is compared with the experimental data reported by Mason et al. [56]. Because water saturation is a key parameter to predict water distribution [4], the validation of the VOF model is conducted by comparing porosity ε and water saturation s with that reported in Ref. [57]. Ref. [21] specially introduced the comparison between experimental results and simulated results from the VOF model. The relative thickness “h" and the water saturation “s" were also adopted in Ref. [21] for comparison. Both samples for validation are set the same key geometrical and operating conditions with the experimental sam ples in Refs. [56,57]. Detailed parameters of the samples have been clearly stated in the previous sections. Results of both models have good agreements with the experimental data, as shown in Fig. 2. Grid inde pendence tests are also performed by testing five various meshes for the
2.2.3. Numerical procedures Hexahedral meshes with the grid size 2 μm � 2 μm � 1 μm are used to discretize computational domains. Model validation and grid indepen dency tests have been conducted to guarantee model accuracy and specific results are discussed in the section 3.1. A popular open-source software Open FOAM is utilized to perform the simulation. Governing equations are discretized by a second order scheme. For the coupling solution of the pressure and velocity, a (pressure implicit split operator) PISO scheme coupling the semi-implicit method is employed. The adjustable time step algorithm is employed with the maximum Courant number to 1 and maximum α Courant number to 0.5. In the present work, each case takes 48 h using 112 Intel Xeon @ 2.93 GHz processors. 2.3. Initial and boundary conditions 2.3.1. FEM model Generally, the PEMFC is clamped by the bolts [59]. Fixed support conditions are applied on the bottom surface of the buffer. Displacement constraints are applied to guarantee the movement of GDL and BP in the z-axis direction only, as shown in Fig. 1b. Owing to the existence of the bolts, the contact surfaces between all components are regarded as bonded to avoid slide. Corresponding to compression ratios of 0%, 8.8%, 17.7%, 26.6% and 35.9%, compression stresses of 0 MPa, 1.31 MPa, 1.70 MPa, 2.62 MPa, and 3.40 MPa are applied to the top surface of the bipolar plate respectively. Both compression ratio and compression stress adopted in this study are in an appropriate range reported in Refs. [60,61]. The mechanical properties of BP and carbon fibers of GDLs are gathered in Table 3. The original porosity adopted for the formal cases is constant 0.84 [62]. Initially, there is no pressure on the top surface of the BP. Each case takes about 24 h by using 16 Intel Xeon @3.50 GHz pro cessors in parallel. 2.3.2. VOF model The gravity is set along the negative z-axis direction. The compressed GDL samples used for the VOF model are obtained from the results of the FEM model. The validation sample is 190 μm thick with a cross-section of 300 μm � 300 μm and a varying porosity reported in Ref. [57]. Pres sure inlet of the liquid water with various values and constant pressure outlet are imposed in this model (see Table 2). The outlet pressure is set to fixed value (zero) and all walls are defined as zero gradient. The walls are considered as symmetric planes. Uniform velocity with a given average value 0 m s 1 is defined for the inlet. The zero gradient and no-slip are adopted to the outlet and walls, respectively. The contact angle θ is 109� [63] for the carbon fibers and 90� for the other walls. Initially, there is no liquid water within the GDL. More details of the mesh size and total number of girds of the computational domains under various compression conditions are given in Table 2.
Table 3 Mechanical properties of the BP and the GDL fibers [64]. Component
Young’s modulus (E) (MPa)
Poisson’s ratio (υ)
BP Fibers
13000 3070
0.26 0.256
Fig. 2. Comparisons between the reconstructed GDL of the present study and experimental studies. (a) Compression; (b) porosity; (c) water saturation. 5
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
other pores in the same layer. Similar available phenomenon of pref erential water transport pathways can also be found in the rectangle short-line areas in Fig. 7. Water saturation values are calculated for seven stochastically generated samples with the same porosity distribution respectively. Details of water saturation under various water inlet pressures as well as compression states for sample 3 are presented in Fig. 5a. It is observed that water saturation decreases with increased compression ratio when water inlet pressure is higher than 4 kPa. But compression has minor effects on water saturation when water inlet pressure is below 3 kPa. Both trends are expected, as compression stress will first deform the upper carbon fibers. Because constant compression pressures are applied and the deformation of the upper fibers will increase the effective con tact areas, the lower carbon fibers will bear smaller forces. When the forces transmitted to lower fibers decrease to a certain level, the forces will not have the ability to deform the fibers near the bottom of GDLs. That means constant compression pressure causes uniform deformation of GDL microstructures along the thickness direction, and the upper region experiences greater distortion. Due to the non-uniform deformation, three regions are first defined in compressed GDLs, as shown in Fig. 6. The small region (presented in the red box with a constant blue color) near the bottom of the GDL, named weak deformation region (WDR) here. In WDR, most fibers almost bear the smallest forces and will not be deformed. Relative mechanism of this phenomenon has been elaborated in the previous paragraph. According to Figs. 5a and 6, when the liquid water inlet pressure is low, the liquid water is in WDR. Therefore, for water inlet pressure below 3 kPa, there is little water saturation difference between GDLs with five compression states. As water inlet pressure increases, the liquid phase will arrive in another special region named moderate deformation region (MDR). Compression pressure contributes to the deformation of GDL. There is a small region in the compressed GDL, and the region is next to WDR, as shown in Fig. 6. Similar to WDR, MDR maintains a stable structure when compression pressure changes. MDR is created by compression stress, but it is independent on values of compression pressure in a certain range (in this study, 1.31 MPa–3.4 MPa). Fig. 5a illustrates the range of the capillary pressure for MDR is 3 kPa–4 kPa. The existence of WDR and MDR can also be demonstrated by the red (1.0 kPa) and blue (4.0 kPa) curves in Fig. 5b. The red curve is horizontal, which indicates the water saturation is in dependent on the compression ratio. A slight downward trend can be observed for the blue curve between the uncompressed GDL (with compression ratio 0%) and compressed GDL (with compression ratio 8.8%). Then the curve maintains horizontal, when compression ratio changes from 8.8% to 35.9%. As the water inlet pressure increases, the water passes through MDR and arrives in the region named strong deformation region (SDR). SDR is the main part in compressed GDLs. In SDR, the region is highly influenced by the compression stress, which differs from WDR and MDR. It is observed from Fig. 5 that water satu ration in SDR strongly depends on compression ratio. This phenomenon is in agreement with the result reported in Ref. [18]. A higher compression ratio leads to a lower level of water saturation. The trend is expected because GDL compression leads to a lower porosity, which in turn decreases the water content, as shown in Fig. 5b. Fig. 7 shows the water invasion process in GDLs with different compression ratios, when various water inlet pressures are employed. It is observed in the dotted line areas that water transport behavior in GDLs with five various compression ratios are similar, when water inlet pressure is below 4 kPa. In addition, different water dynamics can be observed as water inlet pressure increases. This phenomenon agrees well with previous studies. It can be seen from elliptical short-line areas (in Fig. 7) that compression has obvious significance on water breakthrough pressure. Water breakthrough pressures in GDLs with five various compression ratios are 6.0 kPa, 6.4 kPa, 7.0 kPa, 7.4 kPa, and 8.4 kPa respectively.
FEM model and three different grids for the VOF model respectively. The results show that mesh sizes adopted in this study are enough to guar antee accuracy, are shown in Fig. 3. 3.2. Effects of compression on liquid-water dynamics Fig. 4 shows compressed GDL microstructures and water break through behaviors in the sample 7. It is observed that clamping pressure deforms the fibrous GDL microstructures, but it has minor effect on preferential pathways for water breakthrough. Water invasion behaviors are similar in these four compressed GDLs. This may be caused by the uniform compression pressures. The uniform compression pressures make the absolute pore size smaller, but have minor effect on the rela tive pore size distribution in the same layer. Although the pore size decreases, the original larger pores are still large pores compared with
Fig. 3. Mesh independence tests for FEM model and VOF model (a) compres sion; (b) porosity; (c) water saturation.
6
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Fig. 4. Results of FEM model and VOF model for sample 7: (a)–(d) Compressed GDL microstructures in FEM model; (e)–(h) Water breakthrough behaviors inside compressed GDL microstructures in VOF model.
3.3. Quantitative correlations between Pc and s Liquid water saturation, which is linked to capillary pressure by empirical relationships, is a key parameter when two-phase flow insides GDLs are considered [4]. There are already some theoretical and nu merical studies on water transport in compressed GDLs, but none of the previous studies specifically clarifies quantitative correlations between water saturation and capillary pressure. To further investigate quanti tative correlations between water saturation, capillary pressure and compression, a series of correlations are first concluded. As shown in Fig. 5b, the original porosity of GDLs is 0.84 and then it decreases to 0.82, 0.80, 0.77, and 0.73, corresponding to compression ratio 8.80%, 17.7%, 26.6%, and 35.9% respectively. Fig. 8 presents the simulated liquid-water saturation inside compressed GDLs with various compres sion ratios. Approximate exponential relationships are observed and the fitting correlations are shown in Eq. (20). 8 > > 1:407 expð0:599Pc Þ ðCompression ¼ 0:00%Þ > > > > < 0:671 expð0:648Pc Þ ðCompression ¼ 8:80%Þ s ¼ 0:460 expð0:671Pc Þ ðCompression ¼ 17:7%Þ (20) > > > 1:518 expð0:455Pc Þ ðCompression ¼ 26:6%Þ > > > : 1:039 expð0:481Pc Þ ðCompression ¼ 35:9%Þ where, s and Pc are the averaged liquid water saturation (%) and capillary pressure with the unit of kPa. A relatively high water-saturation value is observed for all samples when breakthrough pressures are achieved. It is easy to find that the growth rate of water saturation decreases with compression ratio, as shown in Fig. 8. It is worth noting that the same conclusion can be drawn directly by deriving the corresponding Eq. (20). Quantitative correla tions mentioned above may simplify the qualitative analysis process of water transport behavior, and provide a vital reference for quantita tively investigating water content inside the compressed GDLs. 4. Conclusion
Fig. 5. (a) Volume-averaged liquid water saturation as a function of capillary pressure for sample 3 (WDR-weak deformation region; MDR-moderate defor mation region; SDR-strong deformation region; CPN-compression neglected; CP-compression); (b) Volume-averaged liquid water saturation (porosity) as a function of compression for samples 1-7.
In the present study, two-phase flow in uncompressed and com pressed gas diffusion layers (GDLs) is investigated via integrating finite element method (FEM) with volume of fluid (VOF) modeling tech niques. Water transport behaviors and water saturation (s) in the GDLs with various compression ratios are studied. The main conclusions are 7
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Fig. 6. Three regions in the deformed GDL and liquid water emerging behaviors from the bottom of sample 3 at the slice of y ¼ 150. μm
summarized as follows:
uncompressed and compressed GDLs are same. In MDR, Pc -s curves for all compressed GDLs are same. Water saturation in SDR is highly dependent on compression ratios. Corresponding to WDR, MDR and SDR, the range of the water inlet pressure are 0–3 kPa, 3 kPa–4 kPa, and over 4 kPa respectively.
(1) Assembly pressure causes non-uniform deformation of the GDL microstructures along the thickness direction. (2) Three regions (WDR, MDR, and SDR) in compressed GDLs are first defined in this study. In WDR, Pc -s curves for all
Fig. 7. Liquid water dynamics inside sample 3 at five compression states and various water inlet pressures (0 kPa–8.4 kPa). 8
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Fig. 8. Volume-averaged liquid water saturation as a function of capillary pressure for samples 1–7 at five compression states.
(3) Compression blocks the passage of liquid water but has minor effects on the preferential pathways of water breakthrough. (4) Quantitative correlations between water saturation and capillary pressure for the uncompressed and compressed GDL microstruc tures are first concluded.
National Natural Science Foundation of China for International Coop eration and Exchange (Newton Advanced Fellowship) (Grant No. 51861130359, NAF\R1\180146). References [1] H.W. Wu, G.J. Shih, Y.B. Chen, Appl. Energy 220 (2018) 47–58. [2] Y. Wang, K.S. Chen, J. Mishler, S.C. Cho, X.C. Adroher, Appl. Energy 88 (2011) 981–1007. [3] K. Jiao, X. Li, Prog. Energy Combust. Sci. 37 (2011) 221–291.
Acknowledgement This work was supported by the National Key Research and Devel opment Program of China (Grant No. 2017YFB0102703), and the 9
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933 [35] D. Zhang, Q. Cai, S. Gu, Electrochim. Acta 262 (2018) 282–296. [36] P. Xu, S. Xu, Y. Gao, Commun. Comput. Phys. 23 (2018) 1078–1093. [37] J. Yu, D. Froning, U. Reimer, W. Lehnert, Int. J. Hydrogen Energy 43 (2018) 6318–6330. [38] W.Z. Fang, Y.Q. Tang, L. Chen, Q.J. Kang, W.Q. Tao, Int. J. Heat Mass Transf. 126 (2018) 243–255. [39] J. Liu, S. Shin, S. Um, Renew. Energy 139 (2019) 279–291. [40] K.J. Lee, J.H. Kang, J.H. Nam, Int. J. Hydrogen Energy 39 (2014) 6646–6656. [41] M. Shahraeeni, M. Hoorfar, Int. J. Hydrogen Energy 39 (2014) 10697–10709. [42] M. Fazeli, J. Hinebaugh, A. Bazylak, J. Electrochem. Soc. 162 (2015) F661–F668. [43] T. Agaesse, A. Lamibrac, F.N. Büchi, J. Pauchet, M. Prat, J. Power Sources 331 (2016) 462–474. [44] C.Z. Qin, S. Hassanizadeh, L. Van Oosterhout, Computation 4 (2016) 21. [45] B. Straubhaar, J. Pauchet, M. Prat, Int. J. Heat Mass Transf. 102 (2016) 891–901. [46] X. Huang, Y. He, W. Zhou, D. Deng, Y. Zhao, Powder Technol. 343 (2019) 350–361. [47] L. Chen, Y.L. He, W.Q. Tao, Int. J. Heat Mass Transf. 60 (2013) 252–262. [48] Z.Y. Ahmad, S. Didari, J. Moon, T.A.L. Harris, ECS Transact. 45 (2013) 89–100. [49] Y. Yin, T. Wu, P. He, Q. Du, K. Jiao, Int. J. Hydrogen Energy 39 (2014) 15772–15785. [50] Z. Niu, Z. Bao, J. Wu, Y. Wang, K. Jiao, Appl. Energy 232 (2018) 443–450. [51] D.H. Jeon, H. Kim, J. Power Sources 294 (2015) 393–405. [52] P. Satjaritanun, S. Hirano, A.D. Shum, I.V. Zenyuk, A.Z. Weber, J.W. Weidner, S. ShiMPalee, J. Electrochem. Soc. 165 (2018) F1115–F1126. [53] C. Han, Z. Chen, Energy Sources, Part A Recovery, Util. Environ. Eff. 41 (2019) 1253–1271. [54] L. Chen, H. Luan, Y. He, W. Tao, Int. J. Therm. Sci. 51 (2012) 132-144. [55] Z. Nada, X. Li, J. Shen, Appl. Energy 93 (2012) 39-44. [56] T.J. Mason, J. Millichamp, T.P. Neville, A. El-kharouf, B.G. Pollet, D.J. Brett, J. Power Sources 219 (2012) 52–59. [57] R. Flückiger, F. Marone, M. Stampanoni, A. Wokaun, F.N. Büchi, Electrochim. Acta 56 (2011) 2254–2262. [58] Y. Tang, A.M. Karlsson, M.H. Santare, M. Gilbert, S. Cleghorn, W.B. Johnson, Mater. Sci. Eng. A 425 (2006) 297–304. [59] C. Chien, H. Hu, L.Y., T.H. Su, H.T. Liu, C.T. Wang, P.F. Yang, Y.X. Lu, Energy 113 (2016) 1174–1187. [60] J. Wang, J. Yuan, Sund�en Bengt, Int. J. Energy Res. 41 (2017) 985–1003. [61] V. Norouzifard, M. Bahrami, J. Power Sources 264 (2014) 92–99. [62] J.T. Gostick, M.W. Fowler, M.D. Pritzker, M.A. Ioannidis, L.M. Behra, J. Power Sources 162 (2006) 228–238. [63] D.L. Wood, C. Rulison, R.L. Borup, J. Electrochem. Soc. 157 (2010) B195–B206. [64] M. Movahedi, A. Ramiar, A.A. Ranjber, Energy 142 (2018) 617–632. [65] M. Andersson, V. Vuk�cevi�c, S. Zhang, Y. Qi, H. Jasak, S.B. Beale, W. Lehnert, Int. J. Hydrogen Energy 44 (2019) 11088–11096. [66] M. Andersson, S.B. Beale, U. Reimer, W. Lehnert, D., Int. J. Hydrogen Energy 43 (2018) 2961–2976. [67] Z. Niu, J. Wu, Z. Bao, Y. Wang, Y. Yin, K. Jiao, Int. J. Heat Mass Transf. 139 (2019) 58–68.
[4] T.G. Tranter, A.D. Burns, D.B. Ingham, M. Pourkashanian, Int. J. Hydrogen Energy 40 (2015) 652–664. [5] X. Liu, F. Peng, G. Lou, Z. Wen, J. Power Sources 299 (2015) 85–96. [6] J. Ge, A. Higier, H. Liu, J. Power Sources 159 (2006) 922–927. [7] A. Bazylak, D. Sinton, Z.S. Liu, N. Djilali, J. Power Sources 163 (2007) 784–792. [8] M.A. Hickner, N.P. Siegel, K.S. Chen, D.S. Hussey, D.L. Jacobson, M. Arif, J. Electrochem. Soc. 155 (2008) B427–B434. [9] N.J. Cooper, A.D. Santamaria, M.K. Becton, J.W. Park, Int. J. Hydrogen Energy 42 (2017) 16269–16278. [10] U.U. Ince, H. Mark€ otter, M.G. George, H. Liu, N. Ge, J. Lee, A. Bazylak, Int. J. Hydrogen Energy 43 (2018) 391–406. [11] R.W. Atkinson III, Y. Garsany, B.D. Gould, K.E. Swider-Lyons, I.V. Zenyuk, ACS Appl. Energy Mater. 1 (2018) 191–201. [12] M.B. Sassin, Y. Garsany, B.D. Gould, K. Swider-Lyons, J. Electrochem. Soc. 163 (2016) F808–F815. [13] I.V. Zenyuk, D.Y. Parkinson, L.G. Connolly, A.Z. Weber, J. Power Sources 328 (2016) 364–376. [14] F.C. Cetinbas, R.K. Ahluwalia, A.D. Shum, I.V. Zenyuk, J. Electrochem. Soc. 166 (2019) F3001–F3008. [15] T€ otzke Christian, J. Power Sources 324 (2016) 625–636. [16] X.M. Zhang, X.X. Zhang, Fuel Cells 14 (2014) 303–311. [17] I.V. Zenyuk, D.Y. Parkinson, G. Hwang, A.Z. Weber, Electrochem. Commun. 53 (2015) 24–28. [18] M. Fazeli, J. Hinebaugh, Z. Fishman, C. T€ otzke, W. Lehnert, I. Manke, A. Bazylak, J. Power Sources 335 (2016) 162–171. [19] P. Satjaritanun, J.W. Weidner, S. Hirano, Z. Lu, Y. Khunatorn, S. Ogawa, S. ShiMPalee, J. Electrochem. Soc. 164 (2017) E3359–E3371. [20] A.D. Le, B. Zhou, J. Power Sources 182 (2008) 197–222. [21] Z. Niu, Y. Wang, K. Jiao, J. Wu, J. Electrochem. Soc. 165 (2018) F613–F620. [22] N. Bao, Y. Zhou, K. Jiao, Y. Yin, Q. Du, J. Chen, Eng. Appl. Comput. Fluid Mech. 8 (2014) 26–43. [23] G.R. Molaeimanesh, M.H. Akbari, Int. J. Hydrogen Energy 39 (2014) 8401–8409. [24] G.R. Molaeimanesh, M.H. Akbari, Korean J. Chem. Eng. 31 (2014) 598–610. [25] W. Chen, F. Jiang, Int. J. Hydrogen Energy 41 (2016) 8550–8562. [26] G.R. Molaeimanesh, M.H. Akbari, Int. J. Hydrogen Energy 41 (2016) 14872–14884. [27] M.A. Safi, N.I. Prasianakis, J. Mantzaras, A. Lamibrac, F.N. Büchi, Int. J. Heat Mass Transf. 115 (2017) 238–249. [28] P. Satjaritanun, S. Shimpalee, J.W. Weidner, S. Hirano, Z. Lu, A. Shum, I. V. Zenyuk, S. Ogawa, S. Litster, ECS Transact. 80 (2017) 187–195. [29] Y. Gao, A. Montana, F. Chen, J. Power Sources 342 (2017) 252–265. [30] S. Sakaida, Y. Tabe, T. Chikahisa, ECS Transact. 80 (2017) 123–131. [31] E. Shakerinejad, M.H. Kayhani, M. Nazari, A. Tamayol, Int. J. Hydrogen Energy 43 (2018) 2410–2428. [32] J. Yu, D. Froning, U. Reimer, W. Lehnert, J. Power Sources 389 (2018) 56–60. [33] A.H. Kakaee, G.R. Molaeimanesh, M.H.E. Garmaroudi, Int. J. Hydrogen Energy 43 (2018) 15481–15491. [34] F. Jinuntuya, M. Whiteley, R. Chen, A. Fly, J. Power Sources 378 (2018) 53–65.
10
Contents lists available at ScienceDirect
Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour
Two-phase flow in compressed gas diffusion layer: Finite element and volume of fluid modeling Xia Zhou, Zhiqiang Niu, Zhiming Bao, Jingchao Wang, Zhanrui Liu, Yan Yin, Qing Du, Kui Jiao * State Key Laboratory of Engines, Tianjin University, 135 Yaguan Rd, Tianjin, 300350, China
H I G H L I G H T S
� Microstructures of compressed GDL are digitally reconstructed. � Pore-scale two-phase flow in compressed GDL is investigated using VOF model. � Effects of compression on water transport in GDL are discussed. � Quantitative correlations between Pc and s are concluded. A R T I C L E I N F O
A B S T R A C T
Keywords: Gas diffusion layer Compression Capillary pressure Water saturation
In this study, a stochastic model is used to reconstruct the uncompressed gas diffusion layer (GDL) micro structures. Subsequently, the finite element method (FEM) is conducted for assembly pressure simulation to generate the compressed GDL microstructures. The effects of assembly pressure on GDL deformation are investigated. It is found that assembly pressure causes non-uniform deformation of the GDL along the thickness direction. Finally, a volume of fluid (VOF) model is developed to investigate two-phase flow in the compressed GDL. The results show that when the capillary pressure is higher than 4 kPa, the water saturation decreases as the compression ratio increases. But when the capillary pressure is below 3 kPa, compression has little effect on water saturation. Based on the above findings, three regions namely weak deformation region (WDR), moderate deformation region (MDR), and strong deformation region (SDR) are defined. Impacts of compression on water saturation differ in these three regions. Moreover, compression increases the pressure of water breakthrough, but has minor effects on preferential pathways of water breakthrough. Quantitative correlations between water saturation and capillary pressure in the uncompressed and compressed GDL microstructures are also concluded.
1. Introduction Proton exchange membrane fuel cells (PEMFCs) are regarded as one of the potential power sources for vehicles due to their outstanding merits such as high conversion efficiency and zero emission [1]. Improvement of water transport in PEMFCs is essential for fuel cell commercialization [2]. To achieve maximum removal of liquid water, it is important to understand water transport behavior in the gas diffusion layer (GDL) [3]. The GDL is a vital component in PEMFCs and has multiple functions such as providing transport paths for reactant and product liquid water [4]. Deep understanding of liquid water transport in GDL is important to enhance the cell performance. Numerous studies have been conducted to investigate water transport in the gas diffusion layers (GDLs). Liu et al. [5] elaborated the state and transport mechanism of water inside the uncompressed GDLs. It is worth
noting that GDLs in real PEMFCs are usually compressed due to the clamping pressure, which is always ignored in previous studies. Compression decreases the thickness, porosity and electrical resistance, which in turn affects the performance of PEMFCs [6], as shown in Fig. 1a. Several experimental techniques have been adopted to visualize water transport in the compressed GDLs such as fluorescence micro scopy [7], neutron radiography [8,9] and X-ray tomography [10–17]. Bazylak et al. [7] adopted fluorescence microscopy to provide ex situ visualization of water transport through the compressed GDL and found areas under compression coincide with preferential pathways for water transport. Cooper et al. [9] adopted neutron radiography to examine the through-plane liquid water distribution of an operating PEMFC. Ince et al. [10] investigated water transport in two GDL samples by means of synchrotron X-ray tomography. They found pore saturation in the compressed GDL is higher than that in the uncompressed GDL. Zenyuk
* Corresponding author. E-mail address: [email protected] (K. Jiao). https://doi.org/10.1016/j.jpowsour.2019.226933 Received 11 April 2019; Received in revised form 22 July 2019; Accepted 24 July 2019 Available online 31 July 2019 0378-7753/© 2019 Elsevier B.V. All rights reserved.
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Fig. 1. Schematics of FEM model and VOF model based on GDL microstructures in PEMFC: (a) Morphologies of uncompressed and compressed GDL; (b) the computational domain and boundary conditions used for validation of the FEM model; (c) the computational domain used for formal cases of the FEM model; (d) the computational domain and boundary conditions used for validation and formal cases of the VOF model.
et al. [17] found uniform liquid-water front at low compression, while at high compression, the water predominantly advanced at locations under the channel for higher liquid pressure. Although these visualization techniques deliver some fundamental understanding of liquid water transport in the compressed GDL, they are usually time-consuming and have high requirement for facilities. As an alternative, numerical modeling can investigate water trans port inside the GDL efficiently. Several two-phase flow models have been conducted to discuss water dynamics inside uncompressed and compressed GDLs, such as pore network model (PNM), lattice Boltz mann method (LBM), and volume of fluid (VOF) method. Fazeli et al. [18] employed PNM method extracted from synchrotron-based CT im ages of compressed GDL to simulate liquid water transport in the GDL. Satjaritanun et al. [19] conducted the LBM model to predict the breakthrough pressure for GDL samples. Le et al. [20] employed a VOF model to investigate liquid water in PEMFC with serpentine channels. Niu et al. [21] reconstructed a GDL with the experimentally-determined varying porosity, and employed a three-dimensional (3-D) VOF model to investigate two-phase flow inside the GDL. Bao et al. [22] developed a VOF model to investigate water transport behaviors inside a simplified compressed GDL, and found that deformation highly influenced water droplet dynamics. Typical researches on water transport inside un compressed and compressed GDLs using PNM, LBM, and VOF methods are summarized in Table 1. Although many researches have discussed water transport behaviors in GDLs, only a relatively small percentage of the studies [18,20,51–53] have considered compressed GDL microstructures, as shown in Table 1. However, there are still several limitations in these literatures [18,20, 51–53]: 1). Universal correlation between water saturation and
capillary pressure inside compressed GDLs has not been quantitatively concluded yet; 2). Clear effects of compression on two-phase flow in GDLs have not been concluded; 3). The compressed GDL microstructures adopted in Refs. [18,20,51–53] were obtained either through the complicated high spatial resolution X-ray CT or simplified structures. In this study, to efficiently investigate pore-scale liquid water transport behaviors in compressed GDLs, a stochastic algorithm [49,54, 55] is employed for reconstructing seven GDL microstructure samples. A FEM model [64] is developed to study effects of clamping pressure on GDL deformation. The two-phase VOF model [21,22,47–50] is then employed to simulate water transport behaviors in the compressed GDL. Both FEM and VOF models are validated against the experimental lit eratures [56,57]. To the best of the authors’ knowledge, this is the first study which integrates FEM model with VOF model to investigate water transport behaviors in compressed GDL microstructures. The uniqueness of this study shows as follows: (1) The compressed three-dimensional (3D) GDL microstructures are considered; (2) Pore-scale two-phase flow in the compressed GDL is investigated; (3) Correlations between water saturation and capillary pressure in the GDL microstructures with various compression ratios are quantitatively concluded. 2. Model development 2.1. FEM model 2.1.1. Model assumptions In this study, the FEM model is implemented to simulate the defor mation of the GDL microstructures. Then the compressed microstruc tures are captured for the VOF model. The computational domains are 2
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Table 1 Features of LBM, PNM, and VOF simulations conducted to investigate water transport behaviors inside GDLs. Type
Reference and Published year
Uncompressed GDL LBM Molaeimanesh et al. [23] (2014) Molaeimanesh et al. [24] (2014) Chen and Jiang [25] (2016) Molaeimanesh and Akbari [26] (2016) Safi et al. [27] (2017) Satjaritanun et al. [28] (2017) Gao et al. [29] (2017) Sakaida et al. [30] (2017) Shakerinejad et al. [31] (2018) Yu et al. [32] (2018) Kakaee et al. [33] (2018) Jinuntuya et al. [34] (2018) Zhang et al. [35] (2018) Xu et al. [36] (2018) Yu et al. [37] (2018) Fang et al. [38] (2018) Liu et al. [39] (2019) PNM Lee et al. [40] (2014) Shahraeeni and Hoorfar [41] (2014) Fazeli et al. [42] (2015) Agaesse et al. [43] (2016) Qin et al. [44] (2016) Straubhaar et al. [45] (2016) Huang et al. [46] (2019) VOF Chen et al. [47] (2013) Ahmad et al. [48] (2013) Yin et al. [49] (2014) Niu et al. [21] (2018) Niu et al. [50] (2018) Compressed GDL LBM Jeon and Kim [51] (2015) Satjaritanun et al. [52] (2018) Molaeimanesh et al. [20] (2019) PNM Fazeli et al. [18] (2016) VOF Bao et al. [22] (2014) Han and Chen [53]. (2019)
Research Objective
Microstructure reconstruction
Effects of PTFE distribution Effects of wettability and it spanwise and transverse gradient Effects of PTFE content and distribution Effects of wettability and water droplet size Flow and mass diffusion in a binary air/water mixture Effects of MPL existence Effects of porosity and thickness Effects of hydrophobic-hydrophilic pattern Effects of wettability of GDL Dynamic process of liquid water breaking through Effects of 3D PTFE distribution and the binder content Effects of contact angles and pressure differences Coupling two-phase flow and electrochemical reaction Effects of MPL with different thickness and porosity Effects of asymmetric droplet shape Effects of perforation The heterogeneous multi-transport characteristics in GDLs Effects of width of the rib, thickness of the GDL Effects of GDL hydrophobicity and thickness Effects of liquid water inlet boundary conditions Comparation with experiment Effects of MPL Effects of condensation Effects of network structure (uniform and gradient porosity) Effects of GDL surface Effects of PTFE loading Effects of variable contact angle Comparison with Experiment Effects of mixed wettability
Stochastic Stochastic Stochastic Stochastic X- ray CT X-ray CT Stochastic Stochastic Stochastic Stochastic Stochastic X-ray CT Stochastic Stochastic Stochastic Stochastic Stochastic Simplified cubic network Simplified cubic network X- ray CT X-ray CT Simplified cubic network Simplified cubic network X-ray CT Simplified fiber ridges Simplified cylindrical fibers Stochastic Stochastic Stochastic
Effects of compression ratio Water evolution, water saturation, and breakthrough pressure Removal of a droplet from GDL Effects of compression value Effects of position of droplet, contact angles Effects of air velocity and wettability
2D-simplified Stochastic X-ray CT 2D-simplified Stochastic X-ray CT Simplified Plane Simplified fibers (2D)
composed of a bipolar plate (BP), a GDL, and a buffer, as shown in Fig. 1b and c. The deformation of the catalyst layer (CL) or the mem brane layer (ML) is neglected in this study, because the elastic modulus of either CL or ML is much larger than GDL [58]. In addition, it is assumed that both cathode and anode of the PEMFC have same me chanical behaviors. This model only considers the elastic deformation of the GDL microstructures, and GDL damage due to compression is neglected. Herein, a half domain (one bipolar plate and one gas diffusion layer) is chosen to simulate GDL deformation. The normal displacement along z-axis refers to the deformation of the GDL under the rib. A stochastic algorithm [49,54,55] is employed to reconstruct the fibrous GDLs. The process of GDL reconstruction is based on the following assumptions [50]: (a) The fibers are straight cylinders with the constant diameter, and they are placed in the X–Y plane randomly; (b) The fibers orientation is perpendicular to the Z-axis direction; (c) The overlap of the fibers is allowed and PTFE treating is ignored.
where σi , σj and σ j are the normal stresses, τij ,τik and τjk are the shear stresses, Fi ,Fj and Fk are the body forces. 8 > > > > > > > > > 2 > ∂2 γij ∂ εi ∂2 εj > > > þ ¼0 > 2 2 > ∂Xi ∂Xj ∂Xj ∂Xi > > > > > > > > > > ∂2 γjk ∂2 εj ∂2 εk > > þ ¼0 > 2 2 > > ∂Xj ∂Xk ∂Xk ∂Xj > > > > > > > > > ∂2 εk ∂2 εi ∂2 γki > > þ ¼0 > 2 2 > > ∂Xk ∂Xi ∂Xk < ∂Xi (2) > > > � � > 2 > ∂γjk > ∂ ∂γij ∂γ ∂ εi > > þ ik ¼2 > > ∂ X ∂ X ∂ X ∂ X ∂ X > i k j i j ∂Xk > > > > > > � � > > ∂ ∂γji ∂γjk ∂γik ∂ 2 εj > > ¼2 > > ∂Xj ∂Xk þ ∂Xi > ∂Xj ∂Xi ∂Xk > > > > > > � � > > ∂γkj ∂γij ∂2 εk > ∂ ∂γki > > þ ¼2 > > ∂ X ∂ X ∂ X ∂ X ∂ X > k j i k i ∂Xj > > > > > > :
2.1.2. Governing equations This study considers the elastic deformation of the GDL, which is governed by the equilibrium and compatibility equations with the Hooke’s Law. The constitutive relations are as follows [64]: 8 ∂σ ∂τ ∂τ > > > i þ ij þ ik þ Fi ¼ 0 > > ∂Xi ∂Xj ∂Xk > > > > < ∂σ ∂τji ∂τjk j þ þ þ Fj ¼ 0 (1) > ∂Xj ∂Xi ∂Xk > > > > > ∂σk ∂τki ∂τkj > > > : ∂X þ ∂X þ ∂X þ Fk ¼ 0 k i j
where εi , εj and εj are the normal strains, γij ,γik and γjk are the shear strains. Normal stresses are calculated by:
3
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
8 > > σi ¼ λe þ 2Gεi > > < σ ¼ λe þ 2Gε j j σk ¼ λe þ 2Gεk > > > > : γi j ¼ τij =G
(3)
(4)
e ¼ εi þ εj þ εk λ¼
G¼
υE ð1 þ υÞð1
density (ρ) and dynamic viscosity (μ) for the two-phase mixed fluid are determined as follows [50,65]:
E 2ð1 þ υÞ
(6)
ρ ¼ ρl α þ ρg ð1
αÞ
(9) (10)
αÞ
! r⋅ U ¼ 0
(11)
Phase conservation equation (time-dependence of α):
∂α ! ! þ r⋅ðU αÞ þ r⋅½U r αð1 ∂t
2.1.3. Numerical procedures The FEM model is solved in the commercial software ANSYS APDL. The conservation equations are solved by preconditioned conjugate gradient (PCG) method. To validate the FEM model, the GDL with thickness 370 μm and porosity 0.78 is reconstructed. The variables involved in this model must meet the convergence residual less than 10 6. For the formal samples (sample1-sample7), the original thickness of the GDL is 190 μm with a constant porosity 0.84. All the reconstructed GDLs have a constant carbon fiber diameter 8 μm. Hexahedral meshes with grid size 1.5 μm are used to discretize the GDL microstructures, BP and buffers, as presented in Fig. 1b. The number of grids adopted is 3.10 million for the formal cases and 5.17 million for the validation case, as shown in Table 2.
(12)
αÞ� ¼ 0
! where the symbol U represents the effective velocity vector which is ! shared by both two phases throughout the process. U r named ‘compression velocity’ indicates the relative velocity of liquid and gas at ! ! the interface and subscript r denotes “relative velocity”. Both U and U r are calculated respectively by: ! ! U ¼ α U l þ ð1 ! ! Ur ¼ Ul
!
(13)
αÞU g
! Ug
(14)
Momentum equation: ! ∂ðρ U Þ !! þ r⋅ðρ U U Þ ∂t
2.2. VOF model
! r⋅ðμrU Þ
! ðrU Þ⋅rμ ¼
rpd
! g ⋅! x rρ þ σ krα (15)
where α, σ and k denote the phase fraction, surface tension coefficient and mean curvature of the phase interface respectively. pd is the modi fied pressure for simplifying boundary conditions, and it is determined by:
2.2.1. Model assumptions In this study, some assumptions in the VOF model show as follows [21,50]: 1) Two-phase flows inside the GDL are ideal laminar and incompress ible flows; 2) Transport properties of all fluids keep stable; 3) There is no phase change throughout the process.
pd ¼ p
(16)
ρ! g ⋅! x
here ! g indicates the gravity vector and ! x represents the position vector. The subscript d means “dynamic”. The continuum surface force (CSF) model is implemented in this model to take the effects of the surface tension at the gas-liquid interface into account by adding a force fσ to Eq. (15), which is defined as follow
2.2.2. Governing equations VOF model is a common surface tracking method. In VOF, water phase fraction α is introduced as the main variable for the model to solve. The regions full of liquid water are marked as α ¼ 1, while those filled with air are marked as α ¼ 0. Phase fraction between 0 and 1 means the air-liquid water interface. Vl Vp
(8)
where subscripts l and g indicate liquid water and gas respectively. VOF model governing equations [21,67]: Continuity equation:
where λ is the Lame constant, G is the shear elastic modulus, E is the Young’s modulus, υ is the Poisson’s ratio. i,j and k represent the x,y and z coordinates, respectively.
α¼
αÞ
μ ¼ μl α þ μg ð1
(5)
2υÞ
p ¼ pl α þ pg ð1
(17)
fσ ¼ σkrα
where k which is mentioned is the mean curvature of the phase interface k is calculated by: � � rα k ¼ r⋅! n ¼ r⋅ (18) jrαj
(7)
where Vl is the water volume, Vp is the original volume of the pores. The physical properties such as the pressure (p), averaged volume
the surface unit normal ! n is adjusted in the cells close the wall ac
Table 2 Grid size, total number of grids, and different water inlet pressures for the validation and formal cases of FEM model and VOF model. Cases
Compressed By (%)
Model
Grid Size (μm � μm � μm)
Total number of Grids (Million)
Water inlet pressure (kPa)
Validation Cases Formal Cases (samples 1-7)
0 (CPN)
FEM VOF FEM VOF VOF
1.5 � 1.5 � 1.5 2.0 � 2.0 � 1.0 1.5 � 1.5 � 1.5 2.0 � 2.0 � 1.0 2.0 � 2.0 � 1.0
5.17 3.59 3.10 4.59 4.05 3.71 3.37 3.06
1.0 [57] 1.0–7.4
0 (CPN) 8.8 17.7 26.6 35.9
4
1.0–8.8
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
cording to the equation as follow
3. Results and discussion
! ! n ¼! n w cos θ þ t w sin θ
(19)
3.1. Model validation
! where ! n w and t w are the unit vector normal to the wall and unit vector tangential to the wall, the subscript w indicates the wall. Symbol θ is the contact angle, and it is assumed constant here. Further information on the VOF approach can be found in the Refs [50,65,66].
To validate the FEM model, the compression curve for the 370 μm thick GDL sample is compared with the experimental data reported by Mason et al. [56]. Because water saturation is a key parameter to predict water distribution [4], the validation of the VOF model is conducted by comparing porosity ε and water saturation s with that reported in Ref. [57]. Ref. [21] specially introduced the comparison between experimental results and simulated results from the VOF model. The relative thickness “h" and the water saturation “s" were also adopted in Ref. [21] for comparison. Both samples for validation are set the same key geometrical and operating conditions with the experimental sam ples in Refs. [56,57]. Detailed parameters of the samples have been clearly stated in the previous sections. Results of both models have good agreements with the experimental data, as shown in Fig. 2. Grid inde pendence tests are also performed by testing five various meshes for the
2.2.3. Numerical procedures Hexahedral meshes with the grid size 2 μm � 2 μm � 1 μm are used to discretize computational domains. Model validation and grid indepen dency tests have been conducted to guarantee model accuracy and specific results are discussed in the section 3.1. A popular open-source software Open FOAM is utilized to perform the simulation. Governing equations are discretized by a second order scheme. For the coupling solution of the pressure and velocity, a (pressure implicit split operator) PISO scheme coupling the semi-implicit method is employed. The adjustable time step algorithm is employed with the maximum Courant number to 1 and maximum α Courant number to 0.5. In the present work, each case takes 48 h using 112 Intel Xeon @ 2.93 GHz processors. 2.3. Initial and boundary conditions 2.3.1. FEM model Generally, the PEMFC is clamped by the bolts [59]. Fixed support conditions are applied on the bottom surface of the buffer. Displacement constraints are applied to guarantee the movement of GDL and BP in the z-axis direction only, as shown in Fig. 1b. Owing to the existence of the bolts, the contact surfaces between all components are regarded as bonded to avoid slide. Corresponding to compression ratios of 0%, 8.8%, 17.7%, 26.6% and 35.9%, compression stresses of 0 MPa, 1.31 MPa, 1.70 MPa, 2.62 MPa, and 3.40 MPa are applied to the top surface of the bipolar plate respectively. Both compression ratio and compression stress adopted in this study are in an appropriate range reported in Refs. [60,61]. The mechanical properties of BP and carbon fibers of GDLs are gathered in Table 3. The original porosity adopted for the formal cases is constant 0.84 [62]. Initially, there is no pressure on the top surface of the BP. Each case takes about 24 h by using 16 Intel Xeon @3.50 GHz pro cessors in parallel. 2.3.2. VOF model The gravity is set along the negative z-axis direction. The compressed GDL samples used for the VOF model are obtained from the results of the FEM model. The validation sample is 190 μm thick with a cross-section of 300 μm � 300 μm and a varying porosity reported in Ref. [57]. Pres sure inlet of the liquid water with various values and constant pressure outlet are imposed in this model (see Table 2). The outlet pressure is set to fixed value (zero) and all walls are defined as zero gradient. The walls are considered as symmetric planes. Uniform velocity with a given average value 0 m s 1 is defined for the inlet. The zero gradient and no-slip are adopted to the outlet and walls, respectively. The contact angle θ is 109� [63] for the carbon fibers and 90� for the other walls. Initially, there is no liquid water within the GDL. More details of the mesh size and total number of girds of the computational domains under various compression conditions are given in Table 2.
Table 3 Mechanical properties of the BP and the GDL fibers [64]. Component
Young’s modulus (E) (MPa)
Poisson’s ratio (υ)
BP Fibers
13000 3070
0.26 0.256
Fig. 2. Comparisons between the reconstructed GDL of the present study and experimental studies. (a) Compression; (b) porosity; (c) water saturation. 5
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
other pores in the same layer. Similar available phenomenon of pref erential water transport pathways can also be found in the rectangle short-line areas in Fig. 7. Water saturation values are calculated for seven stochastically generated samples with the same porosity distribution respectively. Details of water saturation under various water inlet pressures as well as compression states for sample 3 are presented in Fig. 5a. It is observed that water saturation decreases with increased compression ratio when water inlet pressure is higher than 4 kPa. But compression has minor effects on water saturation when water inlet pressure is below 3 kPa. Both trends are expected, as compression stress will first deform the upper carbon fibers. Because constant compression pressures are applied and the deformation of the upper fibers will increase the effective con tact areas, the lower carbon fibers will bear smaller forces. When the forces transmitted to lower fibers decrease to a certain level, the forces will not have the ability to deform the fibers near the bottom of GDLs. That means constant compression pressure causes uniform deformation of GDL microstructures along the thickness direction, and the upper region experiences greater distortion. Due to the non-uniform deformation, three regions are first defined in compressed GDLs, as shown in Fig. 6. The small region (presented in the red box with a constant blue color) near the bottom of the GDL, named weak deformation region (WDR) here. In WDR, most fibers almost bear the smallest forces and will not be deformed. Relative mechanism of this phenomenon has been elaborated in the previous paragraph. According to Figs. 5a and 6, when the liquid water inlet pressure is low, the liquid water is in WDR. Therefore, for water inlet pressure below 3 kPa, there is little water saturation difference between GDLs with five compression states. As water inlet pressure increases, the liquid phase will arrive in another special region named moderate deformation region (MDR). Compression pressure contributes to the deformation of GDL. There is a small region in the compressed GDL, and the region is next to WDR, as shown in Fig. 6. Similar to WDR, MDR maintains a stable structure when compression pressure changes. MDR is created by compression stress, but it is independent on values of compression pressure in a certain range (in this study, 1.31 MPa–3.4 MPa). Fig. 5a illustrates the range of the capillary pressure for MDR is 3 kPa–4 kPa. The existence of WDR and MDR can also be demonstrated by the red (1.0 kPa) and blue (4.0 kPa) curves in Fig. 5b. The red curve is horizontal, which indicates the water saturation is in dependent on the compression ratio. A slight downward trend can be observed for the blue curve between the uncompressed GDL (with compression ratio 0%) and compressed GDL (with compression ratio 8.8%). Then the curve maintains horizontal, when compression ratio changes from 8.8% to 35.9%. As the water inlet pressure increases, the water passes through MDR and arrives in the region named strong deformation region (SDR). SDR is the main part in compressed GDLs. In SDR, the region is highly influenced by the compression stress, which differs from WDR and MDR. It is observed from Fig. 5 that water satu ration in SDR strongly depends on compression ratio. This phenomenon is in agreement with the result reported in Ref. [18]. A higher compression ratio leads to a lower level of water saturation. The trend is expected because GDL compression leads to a lower porosity, which in turn decreases the water content, as shown in Fig. 5b. Fig. 7 shows the water invasion process in GDLs with different compression ratios, when various water inlet pressures are employed. It is observed in the dotted line areas that water transport behavior in GDLs with five various compression ratios are similar, when water inlet pressure is below 4 kPa. In addition, different water dynamics can be observed as water inlet pressure increases. This phenomenon agrees well with previous studies. It can be seen from elliptical short-line areas (in Fig. 7) that compression has obvious significance on water breakthrough pressure. Water breakthrough pressures in GDLs with five various compression ratios are 6.0 kPa, 6.4 kPa, 7.0 kPa, 7.4 kPa, and 8.4 kPa respectively.
FEM model and three different grids for the VOF model respectively. The results show that mesh sizes adopted in this study are enough to guar antee accuracy, are shown in Fig. 3. 3.2. Effects of compression on liquid-water dynamics Fig. 4 shows compressed GDL microstructures and water break through behaviors in the sample 7. It is observed that clamping pressure deforms the fibrous GDL microstructures, but it has minor effect on preferential pathways for water breakthrough. Water invasion behaviors are similar in these four compressed GDLs. This may be caused by the uniform compression pressures. The uniform compression pressures make the absolute pore size smaller, but have minor effect on the rela tive pore size distribution in the same layer. Although the pore size decreases, the original larger pores are still large pores compared with
Fig. 3. Mesh independence tests for FEM model and VOF model (a) compres sion; (b) porosity; (c) water saturation.
6
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Fig. 4. Results of FEM model and VOF model for sample 7: (a)–(d) Compressed GDL microstructures in FEM model; (e)–(h) Water breakthrough behaviors inside compressed GDL microstructures in VOF model.
3.3. Quantitative correlations between Pc and s Liquid water saturation, which is linked to capillary pressure by empirical relationships, is a key parameter when two-phase flow insides GDLs are considered [4]. There are already some theoretical and nu merical studies on water transport in compressed GDLs, but none of the previous studies specifically clarifies quantitative correlations between water saturation and capillary pressure. To further investigate quanti tative correlations between water saturation, capillary pressure and compression, a series of correlations are first concluded. As shown in Fig. 5b, the original porosity of GDLs is 0.84 and then it decreases to 0.82, 0.80, 0.77, and 0.73, corresponding to compression ratio 8.80%, 17.7%, 26.6%, and 35.9% respectively. Fig. 8 presents the simulated liquid-water saturation inside compressed GDLs with various compres sion ratios. Approximate exponential relationships are observed and the fitting correlations are shown in Eq. (20). 8 > > 1:407 expð0:599Pc Þ ðCompression ¼ 0:00%Þ > > > > < 0:671 expð0:648Pc Þ ðCompression ¼ 8:80%Þ s ¼ 0:460 expð0:671Pc Þ ðCompression ¼ 17:7%Þ (20) > > > 1:518 expð0:455Pc Þ ðCompression ¼ 26:6%Þ > > > : 1:039 expð0:481Pc Þ ðCompression ¼ 35:9%Þ where, s and Pc are the averaged liquid water saturation (%) and capillary pressure with the unit of kPa. A relatively high water-saturation value is observed for all samples when breakthrough pressures are achieved. It is easy to find that the growth rate of water saturation decreases with compression ratio, as shown in Fig. 8. It is worth noting that the same conclusion can be drawn directly by deriving the corresponding Eq. (20). Quantitative correla tions mentioned above may simplify the qualitative analysis process of water transport behavior, and provide a vital reference for quantita tively investigating water content inside the compressed GDLs. 4. Conclusion
Fig. 5. (a) Volume-averaged liquid water saturation as a function of capillary pressure for sample 3 (WDR-weak deformation region; MDR-moderate defor mation region; SDR-strong deformation region; CPN-compression neglected; CP-compression); (b) Volume-averaged liquid water saturation (porosity) as a function of compression for samples 1-7.
In the present study, two-phase flow in uncompressed and com pressed gas diffusion layers (GDLs) is investigated via integrating finite element method (FEM) with volume of fluid (VOF) modeling tech niques. Water transport behaviors and water saturation (s) in the GDLs with various compression ratios are studied. The main conclusions are 7
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Fig. 6. Three regions in the deformed GDL and liquid water emerging behaviors from the bottom of sample 3 at the slice of y ¼ 150. μm
summarized as follows:
uncompressed and compressed GDLs are same. In MDR, Pc -s curves for all compressed GDLs are same. Water saturation in SDR is highly dependent on compression ratios. Corresponding to WDR, MDR and SDR, the range of the water inlet pressure are 0–3 kPa, 3 kPa–4 kPa, and over 4 kPa respectively.
(1) Assembly pressure causes non-uniform deformation of the GDL microstructures along the thickness direction. (2) Three regions (WDR, MDR, and SDR) in compressed GDLs are first defined in this study. In WDR, Pc -s curves for all
Fig. 7. Liquid water dynamics inside sample 3 at five compression states and various water inlet pressures (0 kPa–8.4 kPa). 8
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933
Fig. 8. Volume-averaged liquid water saturation as a function of capillary pressure for samples 1–7 at five compression states.
(3) Compression blocks the passage of liquid water but has minor effects on the preferential pathways of water breakthrough. (4) Quantitative correlations between water saturation and capillary pressure for the uncompressed and compressed GDL microstruc tures are first concluded.
National Natural Science Foundation of China for International Coop eration and Exchange (Newton Advanced Fellowship) (Grant No. 51861130359, NAF\R1\180146). References [1] H.W. Wu, G.J. Shih, Y.B. Chen, Appl. Energy 220 (2018) 47–58. [2] Y. Wang, K.S. Chen, J. Mishler, S.C. Cho, X.C. Adroher, Appl. Energy 88 (2011) 981–1007. [3] K. Jiao, X. Li, Prog. Energy Combust. Sci. 37 (2011) 221–291.
Acknowledgement This work was supported by the National Key Research and Devel opment Program of China (Grant No. 2017YFB0102703), and the 9
X. Zhou et al.
Journal of Power Sources 437 (2019) 226933 [35] D. Zhang, Q. Cai, S. Gu, Electrochim. Acta 262 (2018) 282–296. [36] P. Xu, S. Xu, Y. Gao, Commun. Comput. Phys. 23 (2018) 1078–1093. [37] J. Yu, D. Froning, U. Reimer, W. Lehnert, Int. J. Hydrogen Energy 43 (2018) 6318–6330. [38] W.Z. Fang, Y.Q. Tang, L. Chen, Q.J. Kang, W.Q. Tao, Int. J. Heat Mass Transf. 126 (2018) 243–255. [39] J. Liu, S. Shin, S. Um, Renew. Energy 139 (2019) 279–291. [40] K.J. Lee, J.H. Kang, J.H. Nam, Int. J. Hydrogen Energy 39 (2014) 6646–6656. [41] M. Shahraeeni, M. Hoorfar, Int. J. Hydrogen Energy 39 (2014) 10697–10709. [42] M. Fazeli, J. Hinebaugh, A. Bazylak, J. Electrochem. Soc. 162 (2015) F661–F668. [43] T. Agaesse, A. Lamibrac, F.N. Büchi, J. Pauchet, M. Prat, J. Power Sources 331 (2016) 462–474. [44] C.Z. Qin, S. Hassanizadeh, L. Van Oosterhout, Computation 4 (2016) 21. [45] B. Straubhaar, J. Pauchet, M. Prat, Int. J. Heat Mass Transf. 102 (2016) 891–901. [46] X. Huang, Y. He, W. Zhou, D. Deng, Y. Zhao, Powder Technol. 343 (2019) 350–361. [47] L. Chen, Y.L. He, W.Q. Tao, Int. J. Heat Mass Transf. 60 (2013) 252–262. [48] Z.Y. Ahmad, S. Didari, J. Moon, T.A.L. Harris, ECS Transact. 45 (2013) 89–100. [49] Y. Yin, T. Wu, P. He, Q. Du, K. Jiao, Int. J. Hydrogen Energy 39 (2014) 15772–15785. [50] Z. Niu, Z. Bao, J. Wu, Y. Wang, K. Jiao, Appl. Energy 232 (2018) 443–450. [51] D.H. Jeon, H. Kim, J. Power Sources 294 (2015) 393–405. [52] P. Satjaritanun, S. Hirano, A.D. Shum, I.V. Zenyuk, A.Z. Weber, J.W. Weidner, S. ShiMPalee, J. Electrochem. Soc. 165 (2018) F1115–F1126. [53] C. Han, Z. Chen, Energy Sources, Part A Recovery, Util. Environ. Eff. 41 (2019) 1253–1271. [54] L. Chen, H. Luan, Y. He, W. Tao, Int. J. Therm. Sci. 51 (2012) 132-144. [55] Z. Nada, X. Li, J. Shen, Appl. Energy 93 (2012) 39-44. [56] T.J. Mason, J. Millichamp, T.P. Neville, A. El-kharouf, B.G. Pollet, D.J. Brett, J. Power Sources 219 (2012) 52–59. [57] R. Flückiger, F. Marone, M. Stampanoni, A. Wokaun, F.N. Büchi, Electrochim. Acta 56 (2011) 2254–2262. [58] Y. Tang, A.M. Karlsson, M.H. Santare, M. Gilbert, S. Cleghorn, W.B. Johnson, Mater. Sci. Eng. A 425 (2006) 297–304. [59] C. Chien, H. Hu, L.Y., T.H. Su, H.T. Liu, C.T. Wang, P.F. Yang, Y.X. Lu, Energy 113 (2016) 1174–1187. [60] J. Wang, J. Yuan, Sund�en Bengt, Int. J. Energy Res. 41 (2017) 985–1003. [61] V. Norouzifard, M. Bahrami, J. Power Sources 264 (2014) 92–99. [62] J.T. Gostick, M.W. Fowler, M.D. Pritzker, M.A. Ioannidis, L.M. Behra, J. Power Sources 162 (2006) 228–238. [63] D.L. Wood, C. Rulison, R.L. Borup, J. Electrochem. Soc. 157 (2010) B195–B206. [64] M. Movahedi, A. Ramiar, A.A. Ranjber, Energy 142 (2018) 617–632. [65] M. Andersson, V. Vuk�cevi�c, S. Zhang, Y. Qi, H. Jasak, S.B. Beale, W. Lehnert, Int. J. Hydrogen Energy 44 (2019) 11088–11096. [66] M. Andersson, S.B. Beale, U. Reimer, W. Lehnert, D., Int. J. Hydrogen Energy 43 (2018) 2961–2976. [67] Z. Niu, J. Wu, Z. Bao, Y. Wang, Y. Yin, K. Jiao, Int. J. Heat Mass Transf. 139 (2019) 58–68.
[4] T.G. Tranter, A.D. Burns, D.B. Ingham, M. Pourkashanian, Int. J. Hydrogen Energy 40 (2015) 652–664. [5] X. Liu, F. Peng, G. Lou, Z. Wen, J. Power Sources 299 (2015) 85–96. [6] J. Ge, A. Higier, H. Liu, J. Power Sources 159 (2006) 922–927. [7] A. Bazylak, D. Sinton, Z.S. Liu, N. Djilali, J. Power Sources 163 (2007) 784–792. [8] M.A. Hickner, N.P. Siegel, K.S. Chen, D.S. Hussey, D.L. Jacobson, M. Arif, J. Electrochem. Soc. 155 (2008) B427–B434. [9] N.J. Cooper, A.D. Santamaria, M.K. Becton, J.W. Park, Int. J. Hydrogen Energy 42 (2017) 16269–16278. [10] U.U. Ince, H. Mark€ otter, M.G. George, H. Liu, N. Ge, J. Lee, A. Bazylak, Int. J. Hydrogen Energy 43 (2018) 391–406. [11] R.W. Atkinson III, Y. Garsany, B.D. Gould, K.E. Swider-Lyons, I.V. Zenyuk, ACS Appl. Energy Mater. 1 (2018) 191–201. [12] M.B. Sassin, Y. Garsany, B.D. Gould, K. Swider-Lyons, J. Electrochem. Soc. 163 (2016) F808–F815. [13] I.V. Zenyuk, D.Y. Parkinson, L.G. Connolly, A.Z. Weber, J. Power Sources 328 (2016) 364–376. [14] F.C. Cetinbas, R.K. Ahluwalia, A.D. Shum, I.V. Zenyuk, J. Electrochem. Soc. 166 (2019) F3001–F3008. [15] T€ otzke Christian, J. Power Sources 324 (2016) 625–636. [16] X.M. Zhang, X.X. Zhang, Fuel Cells 14 (2014) 303–311. [17] I.V. Zenyuk, D.Y. Parkinson, G. Hwang, A.Z. Weber, Electrochem. Commun. 53 (2015) 24–28. [18] M. Fazeli, J. Hinebaugh, Z. Fishman, C. T€ otzke, W. Lehnert, I. Manke, A. Bazylak, J. Power Sources 335 (2016) 162–171. [19] P. Satjaritanun, J.W. Weidner, S. Hirano, Z. Lu, Y. Khunatorn, S. Ogawa, S. ShiMPalee, J. Electrochem. Soc. 164 (2017) E3359–E3371. [20] A.D. Le, B. Zhou, J. Power Sources 182 (2008) 197–222. [21] Z. Niu, Y. Wang, K. Jiao, J. Wu, J. Electrochem. Soc. 165 (2018) F613–F620. [22] N. Bao, Y. Zhou, K. Jiao, Y. Yin, Q. Du, J. Chen, Eng. Appl. Comput. Fluid Mech. 8 (2014) 26–43. [23] G.R. Molaeimanesh, M.H. Akbari, Int. J. Hydrogen Energy 39 (2014) 8401–8409. [24] G.R. Molaeimanesh, M.H. Akbari, Korean J. Chem. Eng. 31 (2014) 598–610. [25] W. Chen, F. Jiang, Int. J. Hydrogen Energy 41 (2016) 8550–8562. [26] G.R. Molaeimanesh, M.H. Akbari, Int. J. Hydrogen Energy 41 (2016) 14872–14884. [27] M.A. Safi, N.I. Prasianakis, J. Mantzaras, A. Lamibrac, F.N. Büchi, Int. J. Heat Mass Transf. 115 (2017) 238–249. [28] P. Satjaritanun, S. Shimpalee, J.W. Weidner, S. Hirano, Z. Lu, A. Shum, I. V. Zenyuk, S. Ogawa, S. Litster, ECS Transact. 80 (2017) 187–195. [29] Y. Gao, A. Montana, F. Chen, J. Power Sources 342 (2017) 252–265. [30] S. Sakaida, Y. Tabe, T. Chikahisa, ECS Transact. 80 (2017) 123–131. [31] E. Shakerinejad, M.H. Kayhani, M. Nazari, A. Tamayol, Int. J. Hydrogen Energy 43 (2018) 2410–2428. [32] J. Yu, D. Froning, U. Reimer, W. Lehnert, J. Power Sources 389 (2018) 56–60. [33] A.H. Kakaee, G.R. Molaeimanesh, M.H.E. Garmaroudi, Int. J. Hydrogen Energy 43 (2018) 15481–15491. [34] F. Jinuntuya, M. Whiteley, R. Chen, A. Fly, J. Power Sources 378 (2018) 53–65.
10