- Email: [email protected]
Calculation of effective transport properties of partially saturated gas diffusion layers
Journal of Power Sources 340 (2017) 111e120
Contents lists available at ScienceDirect
Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour
Calculation of effective transport properties of partially saturated gas diffusion layers Tomasz Bednarek*, Georgios Tsotridis European Commission, Joint Research Centre (JRC), Directorate for Energy, Transport and Climate, Westerduinweg 3, NL-1755 LE, Petten, The Netherlands
h i g h l i g h t s Dry and partially saturated gas diffusion layers are investigated. The transport properties as a function of liquid water saturation are calculated. The fine mist approach is used to impregnate GDL without specimen boundary effects. Liquid water distribution depends on exact fibre geometry and capillary forces. Results of Two-Phase Lattice Boltzmann and Pore-Morphology simulations are compared.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 6 August 2016 Received in revised form 20 October 2016 Accepted 27 October 2016
A large number of currently available Computational Fluid Dynamics numerical models of Polymer Electrolyte Membrane Fuel Cells (PEMFC) are based on the assumption that porous structures are mainly considered as thin and homogenous layers, hence the mass transport equations in structures such as Gas Diffusion Layers (GDL) are usually modelled according to the Darcy assumptions. Application of homogenous models implies that the effects of porous structures are taken into consideration via the effective transport properties of porosity, tortuosity, permeability (or flow resistance), diffusivity, electric and thermal conductivity. Therefore, reliable values of those effective properties of GDL play a significant role for PEMFC modelling when employing Computational Fluid Dynamics, since these parameters are required as input values for performing the numerical calculations. The objective of the current study is to calculate the effective transport properties of GDL, namely gas permeability, diffusivity and thermal conductivity, as a function of liquid water saturation by using the Lattice-Boltzmann approach. The study proposes a method of uniform water impregnation of the GDL based on the “Fine-Mist” assumption by taking into account the surface tension of water droplets and the actual shape of GDL pores. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Keywords: Effective transport properties Gas diffusion layer Liquid water saturation Polymer electrolyte fuel cell Lattice-Boltzmann
1. Introduction A common assumption in using Computational Fluid Dynamics (CFD) approaches in modelling of Polymer Electrolyte Fuel Cells (PEMFC) entails the homogenization of the porous structures of the different fuel cell components, such as Gas Diffusion Layers (GDL), Micro-Porous Layers (MPL) and Catalyst Layers (CL). This assumption implies that these porous structures are considered as uniform and continuous, hence the statistical properties (namely porosity, pore size distribution and tortuosity) are used in characterising them. Such an approach for porous PEMFC components was investigated by a number of papers, among others, by Cheddie and
* Corresponding author. E-mail address: [email protected] (T. Bednarek).
Munroe [1], Weber at al. [2], Siegel [3], Bazylak [4], Biyikoglu [5], Wang [6], Squires [7], Djilali [8] and Lomax et al. [9]. In these papers, the porous PEMFC components are treated as homogenous structures by incorporating the Darcy approximation. Since the shape of liquid water droplets and that of the structure of GDL are not taken into account, their contribution to the mass transport within the GDL is neglected, which affects the accuracy of the PEMFC modelling [10,11]. Furthermore, a substantial research effort has been attributed to numerical reproduction of the porous structures of GDLs. Becker et al. [12,13], Zamel et al. [14e16], Mukherjee et al. [17], Daino and Kandlikar [18] propose techniques for reproducing the fibre€tzke et al. [19], Pfrang et al. [20e22], Bosstructure, whereas To omoiu et al. [23] use X-ray Computer Tomography (CT) for reproducing the geometry of the GDL microstructure. Bosomoiu et al.
http://dx.doi.org/10.1016/j.jpowsour.2016.10.098 0378-7753/© 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
112
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
[23] presented the full procedure of CT image acquisition and image manipulation to numerically reproduce real GDL structures. The influence of liquid water on transport properties of GDL was investigated, among others, by Zamel et al. [14], Bosomoiu et al. [23], Schultz et al. [24] using the Pore-Morphology Model (PMM) [25,26]. Stochastic numerical models for creating artificial 3D models of GDL or MPL are considered, among others, by Zamel at al [15,16]. and Becker et al. [12,13]. The reproduction of the GDL structure on the basis of CT imaging and subsequential establishment of effective GDL properties are done, among others, by Pfrang et al. [20], Bosomoiu et al. [23], Gao et al. [27] and García-Salaberri et al. [28,29]. The current study presents a methodology of uniform liquid water impregnating of a GDL sample. Liquid water is introduced as a “fine mist”, by assuming that the water droplets coalesce and subsequently fill the GDL pores. There are no external forces applied to the system, and the flow of liquid water is due to capillary forces only. The “fine mist” approach assumes that the GDL sample is uniformly impregnated with liquid water. The numerical representation of the GDL substrate without the micro-porous layer (MPL) is reproduced from X-ray CT images at a resolution of 2.2 mm, sufficient to determine the 3D anisotropic structure of the GDL fibres. Several liquid water saturation levels have been considered. In contrast to previously published literature [14,23,24,26,28,30] where the capillary forces have no effect on liquid droplets inside the GDL, the current study considers the Lattice-Boltzmann (LB) two-phase approach to simulate water droplet coalescence and movement inside the GDL. The effective transport properties of GDL, namely gas permeability, diffusivity and thermal conductivity, are evaluated based on the Lattice-Boltzmann analysis as a function of liquid water saturation. It is shown that the liquid water has an influence on the effective GDL properties. Results of LB simulations are compared with PoreMorphology Model (PMM) applied in Bosomoiu et al. [23].
Fig. 1. Lattice-Boltzmann fundamentals: a) the movement and collisions of particles, b) the scheme of the D2Q9 lattice.
The most popular Bhatnagar-Gross-Krook (BGK) collision operator [33] provides the standard LB method formula, including the external force term F
1 eq fa ðx þ vea Dt; t þ DtÞ fa ðx; tÞ ¼ fa ðx; tÞ fa ðx; tÞ þ Fa ðx; tÞ t n |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} transport
2. The Lattice Boltzmann method
collision
(3) The Lattice Boltzmann method is based on simulation of particle movement and collisions with other particles [31], shown in Fig. 1a. In particular, the LB method treats the fluid as a set of particle distribution functions, while the motion is described by collisions between particles. The LB simulation algorithm assumes two steps: transport and collision (see Fig. 1)
fa ðx þ vea Dt; tÞ ¼ fa ðx; tÞ fa ðx; t þ DtÞ ¼ fa ðx; tÞ þ Cðx; tÞ
transport collision
(1)
where fa ðx; tÞ is the density distribution function at lattice position x and time t, fa is the density distribution function after transport process, v¼Dx/Dt is the velocity of the particle, a refers to the consecutive neighbours of the central node e0 (see directions in Fig. 1b). C is a collision operator. The discrete velocity depends on the assumed lattice model, i.e. D2Q9 or D3Q19 (it means 2-dimensions 9-nodes or 3-dimensions, 19-nodes respectively). The scheme of D2Q9 lattice is shown in Fig. 1b. From the D2Q9 lattice the discrete velocity vector ea is given as follows [32,33].
8 > > > > > <
ð0 0Þ ða 1Þp ða 1Þp cos sin 2 2 ea ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > ða 5Þp ða 5Þp p > > sin þ : 2 cos 2 2 4
eq
where fa ðx; tÞ is an equilibrium distribution function, tn is the dimensionless relaxation time. The left-hand-side of equation (3) stands for the transport step whereas the right-hand-side of equation (3) corresponds to the collision step (Fig. 1a). The collision term leads to the local Maxwellian equilibrium on a time scale t [32,33]. eq The equilibrium distribution function fa ðx; tÞ for a D2Q9 lattice take the form [32].
3 9 3 eq fa ðx; tÞ ¼ ua r 1 þ 2 ðea $uÞ þ 4 ðea $uÞ2 2 u2 c 2c 2c
The weight factors for D2Q9 lattice are as follows: u0¼4/9, u1and u5-8¼1/36. The fluid density ra and the fluid velocity u are given by the first and the second moments of the density distribution function [32].
4¼1/9
r¼
P
fa aP
ru ¼
a¼0 a ¼ 1; 2; 3; 4 a ¼ 5; 6; 7; 8 (2)
(4)
a
fa ea
(5)
The viscosity of the fluid (in lattice units) is related to the relaxation time by formula
tn ¼ c2s ðtn 0:5ÞDt where cs is the sound speed in lattice units.
(6)
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
following non-ideal equation of state in the lattice units [32,34,35].
3. Two-phase model The gaseous and liquid phases are modelled as two separate lattices with an interface zone between them, Fig. 2a. Originally, Shan and Chen defined the non-local interaction among particles as a force between particles at position x and x' [34,35]
Fðx; x0 Þ ¼ Gðjx x0 jÞjðxÞjðx0 Þðx0 xÞ
(7)
where G is Green's function and j is an effective mass depending on the local density. The force basically arises from cohesion pairwise point coupling between effective mass at location x and neighbouring location x', Fig. 2b. The force given by equation (7) does not conserve the local momentum during the collision. However it has been proved in Ref. [35] that the total momentum of the system is conserved e no net momentum is introduced into the system. The total force acting on the particle at position x is the sum of all interaction forces. In the lattice system, i.e. mentioned D2Q9, constant number of close neighbours N is assumed (Fig. 1b) that interacts with the particle at location x. The interaction force given in equation (7) can be rewritten after [32].
FðxÞ ¼ g jðxÞc2s
N X
a¼1
w jea j2 jðx þ ea Þea
(8)
where g is the cohesion interaction strength between neighbouring particles. The weights wðjea j2 Þ are introduced to account for various distances between particles which maintaining an isotropic reaction force. A positive value of the interaction strength represents a repulsive force between neighbours, whereas negative g corresponds to attractive interaction. Shan and Chen proposed the following form of effective mass [34,35].
j ¼ r0 1 exp
r1 r2
113
(9)
where r1 and r2 are densities of both phases, r0 is a normalization constant usually assumed as 1.0 [32]. The Taylor expansion of the interaction force given in equation (8) leads to the form given by Ref. [36].
c4 FðxÞ ¼ g c2s jVj þ s jVðDjÞ þ … 2
(10)
The interaction force given in equation (10) covers a non-ideal equation of state (the first term) and surface tension (the term corresponding to Dj). The correction of the ideal equation of state (namely p¼r c2s ) due to force given in formula (10) leads to the
g p ¼ rc2s þ c2s j2 2
(11)
Hou at al. and Huang et al. in Refs. [37,38] have validated the presented pseudo-potential model by simulation of a bubble in fully periodic system. It is shown that the presented model correctly simulates the two-phase interaction, i.e. a major improvement has been noticed with comparison to RothmanKeller model [39,40]. For reason of simplicity and computational efficiency the twophase model based on cohesion forces has been widely applied in simulating a variety multiphase flows. However the model developed by Hou et al. and Huang et al. in Refs. [37,38] suffers from some limitations highlighted by Chen et al. [32]. One limitation arises from the relatively low density and kinematic viscosity (relaxation time in particular) ratios. Mukherjee et al. and Kang et al. [30,41] have shown that the density ratio does not play significant role for micro-porous structures since the capillary pressure is the main driving force of liquid flow inside GDL. However to improve the simulation accuracy a dynamic model has been adopted [42], where the particle velocity u is scaled by the density ratio k.
u¼
u1 þ ku2 ; 1þk
k¼
r2 r1
(12)
where subscripts 1 and 2 correspond to separate fluid phases.
4. Calculation of effective transport properties The calculated transport parameters such as permeability, diffusivity and thermal conductivity are considered as scalar values along the direction x through GDL thickness, shown in Fig. 3. The simulations of gas flow, diffusion and heat transfer are performed by using the LB. The transient LB calculations are terminated when the difference between the normalized quantities (velocity, diffusion flux or heat flux in case of permeability, diffusivity and thermal conductivity respectively) from the current and previous time step are smaller than the assumed accuracy ε ¼ 1.0e-4. The calculation of permeability using LB is based on the singlephase flow simulation through the GDL sample. The current study considers dry air as the gas phase. It is assumed that there is no mass exchange between gas and liquid droplets occupying GDL pores e no evaporation/condensation. The effective gas permeability is expressed as follows
keff ¼
um Dp
(13)
where keff is the saturation dependent effective permeability, u is
Fig. 2. The two-phase modelling: a) division of lattices, b) cohesion forces in liquid phase.
114
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
Fig. 3. Location of investigated samples inside the Sigracet®GDL34BC specimen (left hand side) and computational domain (right hand side).
the average out-of-plane velocity component, m is the viscosity and Dp is the pressure gradient between the top and bottom sides of the GDL sample, Fig. 3. In the current study the pressure gradient Dp is assumed equal to 1.0 [kPa]. The effective diffusivity of the GDL partially saturated with liquid water is calculated on the basis of simulation of oxygen diffusion in air. The through plane effective diffusivity Deff is calculated according to Fick's law [16].
Deff ¼
j
Dc
(14)
where j is the out-of-plane component of the average diffusion flux, Dc is the oxygen concentration difference between top and bottom planes of the GDL sample Fig. 3. The effective thermal conductivity of the GDL partially saturated with liquid water is calculated according to Fourier's law
q
leff ¼ DT
active area (2300 2300 pixels). To minimise the statistical error in calculating the effective properties of the GDL, three samples were considered form different locations of the GDL specimen, shown in Fig. 3. The Sigracet® GDL34BC contains the micro-porous layer (MPL), however as it is shown in Ref. [23], the employed CT scan does not have the sufficient resolution to reliably reproduce the MPL since the pore sizes are of the order of nm [16]. Therefore, the MPL was excluded from considerations in this study. The volume of the investigated GDL sample is shown in Fig. 3. The boundary conditions corresponding to gradients of pressure, concentration and temperature (for calculating the permeability, diffusivity and thermal conductivity respectively) are applied at the bottom and the top of the GDL surfaces. All other sides of the sample are assumed as solid walls with no slip (in case of flow simulation) or no flux conditions (for diffusion and heat transfer). 5.1. Water impregnation of GDL
(15)
where q is the out-of-plane component (along x axis, Fig. 3) of the average heat flux at outlet plane, DT is temperature difference between top and bottom planes of the GDL sample, Fig. 3. The liquid droplets occupying the GDL pores are assumed to be rigid hence the model of heat transfer through the GDL takes only the conduction through the carbon fibres and water droplets into account and there is no heat transfer through the empty GDL pores. The heat transfer by radiation is neglected because the normal operating temperature of PEM fuel cell is lower than 100OC.
5. Results The present study considers three effective transport parameters of porous GDLs, namely: gas permeability, diffusivity and thermal conductivity as a function of liquid water saturation. The Sigracet® GDL34BC material is investigated in this study, which was also considered by Bosomoiu et al. [23]. The material was analysed by using the X-ray computer tomography system Phoenix Nanotom S which includes a 180 [kV] nanofocus X-ray tube and a 2D X-ray detector with 120 120 [mm]
During fuel cell operation water vapour diffuses and may condensate inside the diffusion media, hence a number of pores could be occupied by liquid water droplets. However the real conditions of fuel cell operation could not reproduce the contribution of liquid water distribution inside the GDL. During fuel cell operation the flow of reactant gasses will affect the accumulation of liquid water inside GDL due to pressure drop. The convective flow of the humidified gas could possibly remove the liquid water droplets or could force them to move to locations close to boundaries, implying a non-uniform liquid water distribution. However, intrusion of liquid water into a dry GDL results in a non-uniform distribution of liquid water inside the GDL volume. The area close to the uniform inlet is flooded by water, as shown in Ref. [43]. The substitution of the uniform inlet surface by an inlet spot, or a number of inlet spots does not address the problem of a flooded inlet region. The local inlet flooding effects are clearly visible in Fig. 4, where a LB simulation of water injection into GDL is presented and will affect the effective transport properties of the GDL. To avoid the effects of boundary conditions on water distribution inside the GDL an alternative method of water impregnation is proposed. It is assumed that there are no liquid water injectors at
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
115
Fig. 4. Local inlet effects, the lower layers of the GDL are flooded with water.
Fig. 5. Lattice-Boltzmann simulations: a) Two-phase simulation of water droplets creation from fine mist to droplet coalescence. b) Streamlines of gas flow through partially saturated GDL, cross-section through GDL sample.
116
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
the specimen boundaries, therefore all sides of the specimen are impermeable walls. Instead of an inlet boundary condition, the GDL pores are artificially saturated by a fine water mist, i.e. very small water droplets of size corresponding to the size of the voxel (2.2 mm). The water droplets are introduced randomly throughout the whole volume of the sample. To avoid formation of suspended droplets, each artificially introduced water droplet either attach itself to the GDL fibres or would coalese with other water droplets. The LB two-phase model is applied, where both the movement of the liquid phase and droplet coalesce caused by capillary forces due to surface tension and contact angle with GDL fibres, are considered. There are no other external body forces acting on the system. The coalescence and motion of liquid droplets in the GDL sample are shown in Fig. 5a. The comparison of liquid droplets size distribution at two water saturation levels of 25% and 75% is presented in Fig. 6a. Solid lines correspond to Lattice-Boltzmann simulations whereas dotted lines refer to Pore-Morphology Model (PMM) [24,25]. It is shown, that droplets obtained by LB simulation are larger than those obtained by PMM. In the case of 25% water saturation using PMM most of the droplets are 35e65 mm in size, whereas the range of droplet sizes for LB is wider and ranges from 16 up to 85 mm. Such a difference is caused by inclusion of surface tension phenomena and the motion of droplets is due to capillary forces in case of LB simulations, contrary to PMM which does not allow movement and coalescence of neighbouring droplets. This effect is even more visible at higher water saturation level of 75%. In case of PMM only new smaller droplets (form 15e35 mm in size) are introduced, whereas in case of LB simulation, the droplets coalesce surround the fibres and cover a larger volume of GDL sample e the droplets size ranges form 20 up to 140 mm. 5.2. Effective transport properties of the GDL as a function of liquid water saturation The effective permeability of GDL as a function of liquid water saturation is calculated for air at ambient pressure with temperature 65 C (density of 1.050 [kg m3] and kinematic viscosity of 2.15e-5 [Pa s]). The driving force of flow is the pressure difference between the top and the bottom surfaces of the GDL sample and is equal to 1.0 [kPa], Fig. 3. The streamlines of gas flow through the GDL sample are presented in Fig. 5b. The through-plane permeability (x axis in Fig. 3) of dry GDL is
calculated by using formula (13), and it ranges form 3.20e-12 [m2] to 3.55e-12 [m2] corresponding to sample 2 and 3 respectively. These values are similar to those provided by the manufacturer in the GDL “Properties data sheet” [44], namely 3.65e-12 [m2] and to the value of 3.17e-12 [m2] obtained by the authors in a previous investigation [23]. The corresponding distribution of permeability as a function of water saturation is presented in Fig. 7a and Table 1. The relative difference between the permeability values from the samples taken from various locations of the GDL structure does not exceed 10% in case of dry GDL and is less than 3% for the case of 90% of liquid water saturation. As expected, the gas permeability strongly depends on the water saturation of the GDL. Significant decrease of GDL permeability appears when the water saturation level exceeds 50%, Fig. 7a. The reason of such a large decrease in permeability corresponds to liquid water droplets occupying the largest pores of GDL, hence blocking the convective flow through them, as ilustrated in Fig. 6b where the pore size distribution of a dry sample and at two liquid water saturation levels of 25% and 50% respectively are presented. It is seen that in case of saturated GDL the number of the largest open pores decreases, because these pores are occupied by liquid water droplets, hence forcing the gas to pass mainly through smaller pores with a higher flow resistance. The effective diffusivity is calculated by using equation (14). The oxygen concentrations at the top and the bottom surfaces of the samples are fixed. On the top surface oxygen concentration corresponds to fully humidified air, namely 13.5 mol m3, while the oxygen concentration at the bottom surface is equal to 12.0 mol m3. The simulation is performed at temperature T ¼ 65OC. The calculated values of effective diffusivities of dry GDL range from 6.50e-6 m2 s1 to 7.13e-6 m2 s1, and are comparable to the value of 6.6e-6 m2 s1 presented by a previous study [23]. However, they are different to the experimental value of 5.9e-6 m2 s1 obtained by Pant [45] because the experiment by Pant was performed at 28 C whereas both current study and [23] use 65OC, a more typical operating temperature for fuel cell applications. The effect of liquid water saturation on the effective diffusivity is shown in Fig. 7b. Similar to gas permeability, the values for dry GDL differ about 10% and that the difference decreases with increasing saturation level. It is shown that the liquid water saturation has a major impact to the effective diffusivity. It is observed that 20% of water saturation considerably reduces the diffusivity by over 50%, whereas 50% of water saturation decreases diffusivity substantially. The diffusion coefficient of GDL where half of pore space is
Fig. 6. Size distribution of: a) water droplets (comparison between Pore-Morphology and Lattice-Boltzmann methods) and b) open pore size distribution.
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
Fig. 7. Effective properties of GDL as a function of liquid water saturation: a)gas permeability, b) gas diffusivity, c) thermal conductivity.
117
occupied by liquid water (saturation 50%) is about 7 times smaller (1.02e-6 m2 s1) than the diffusion coefficient of the dry GDL. The reason of such significant reduction of the diffusion coefficient is attributed, as in the case of permeability, to the location of liquid water droplets inside the GDL. Since the largest pores are occupied by liquid water (Fig. 6b) the diffusion process is forced to take place through the smaller pores, thereby been significantly reduced. The effective thermal conductivity coefficient is calculated by formula (15). Contrary to the previous cases, thermal conductivity is affected by both the solid fibres and the liquid water droplets. The conductivities of the individual phases, namely carbon fibres and liquid water, at 65 C are assumed to be equal 128.95 W m1K1 and 0.6504 W m1K1 respectively [46,47]. The calculated thermal conductivity coefficient for dry GDL samples varies from 0.83 W m1K1 to 0.92 W m1K1 which is similar to previous investigations [23]. The difference between the most conductive sample 2 and the less conductive sample 3 is about 10%. It is also noticed that with increasing water saturation level the differences between thermal conductivities of samples decreases, since more heat is transported via the liquid water droplets occupying GDL pores. Fig. 7c shows the effect of liquid water on the thermal conductivity of GDL. It is noticed that the effective thermal conductivity of the partially saturated GDL increases with increasing saturation. At 30% water saturation the effective thermal conduction coefficient is nearly doubled as compared to the dry GDL. Further impregnation of GDL does not affect the effective conductivity so strongly, due to the fact that the largest pores are already occupied by liquid water droplets, hence further impregnation fills relatively smaller pores, where most of the heat transfer is performed via the more conductive carbon fibres. The effective conductivity value at increasing liquid water saturation levels tends to the value of approximately 2.0 W m1 K1. However, the current study depicts a different water effect than that presented by Bosomoiu et al. [23]. In an effort to explain the discrepancies of the effective thermal coefficients, the different water saturation levels are re-calculated using the same GeoDict software [48], as in the case of Bosomoiu et al. [23]. The plot of the effective thermal conductivity as a function of water saturation obtained by using Geodict and LB software is also presented in Fig. 7c. The thermal conductivities calculated by both GeoDict and LB are similar for a dry GDL sample. However, the effect of liquid water on thermal conductivity between the current study and [23] (Fig. 7c) arises from different water impregnation methods. In the paper by Bosomoiu et al. [23] the Pore-Morphology impregnation method was used [25,26], whereas in the current study the “ finemist” approach with the two-phase LB simulations have been applied. The differences of liquid water distribution obtained by using the “ fine-mist” approach and the PMM at the pore scale are presented in Fig. 8. The Pore-Morphology Model (Fig. 8a) bases on specifying the spherical water droplet diameter which is then compared with the GDL pores. In cases where the droplet fits to the GDL pore, the droplet remains in the pore, implying that the pore is
Table 1 Effective properties of partially saturated GDL. Liquid water saturation
Gas permeability [m1]
Diffusivity [m2s1]
Thermal conductivity [W m1 K1]
Reference [23]
Current study (average)
Reference [23]
Current study (average)
Reference [23]
Current study (average)
dry 25% 50% 75%
3.14E-12 9.42E-13 3.14E-13 6.28E-14
3.38E-12 8.04E-13 2.30E-13 3.84E-14
6.517E-06 3.375E-06 1.406E-06 1.875E-07
6.83E-06 3.35E-06 1.03E-06 6.83E-07
0.908 1.076 1.380 1.776
0.877 1.513 1.839 1.952
118
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
Fig. 8. Impregnated GDL using: a) Pore-morphology method, b) Lattice-Boltzmann method. In both cases the liquid water saturation is equal ~10%. Red lines indicate the contact area between liquid water (blue) and fibres (grey). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
occupied by liquid water [25,26]. However, from Fig. 8a is shown, that the spherical droplets considered under PMM imply that the contact area between droplet and fibres is considerably smaller than for the LB simulations where capillary interactions are taken into account. Such behaviour of the PMM can lead to underestimating the effects of liquid water saturation on the effective thermal conductivity because the contact area between GDL fibres and water droplets plays a crucial role. In case of the two-phase simulation by LB, the shape of the water droplets is according to real shape of the GDL pores and is influenced by the values of water surface tension. In the presented approach the water droplets can be transported due to capillary forces, even from very inaccessible locations. In particular, moving droplets between GDL pores can coalesce creating flooded zones surrounding some fibres. As a result, the droplets created by “finemist” can lead to formation of droplets larger than the GDL pores themselves. Whereas for the PMM model the formation of droplets larger that the GDL pores is not allowed. Fig. 8b shows that the contact area between droplets and carbon fibres is considerably higher in case of LB two-phase simulation than the case of PMM, hence the effect of water saturation is more significant in the case of thermal conductivity. However, the differences between impregnation methods, employed by PMM and Fine-Mist approach based on LB simulation, do not play a significant role when calculating the effective permeability and diffusivity. This is due to the fact that the mass transport due to pressure or concentration difference takes place in pore space where the influence of contact area between fibres and water droplets is not significant. However, as it was shown in Fig. 7 the effect of water impregnation by the fine mist approach on gas permeability and diffusivity is greater as compared to the pore morphology method. The reason is that for LB simulation the shape of droplets fit to the GDL pores structure, Fig. 8b, closing the pores more tightly as compared to PMM, Fig. 8a. Nevertheless, the differences in the values of the calculated effective properties
between the two methods are not as large as in the case of the thermal conductivity. 6. Conclusions The study shows the effect of liquid water saturation on the effective transport properties of GDL, namely: gas permeability, diffusivity and thermal conductivity, based on the computer tomography image of Sigracet® GDL34BA (without MPL) material. The Fine-Mist-Approach together with two-phase Lattice-Boltzmann simulation was applied to calculate the liquid water distribution inside the GDL sample. The current work, in contrast to Pore-Morphology-Method, considers the liquid water distribution inside the GDL by taking into account the exact geometry of the GDL fibre structure and the capillary forces due to the surface tension of liquid water and contact angle with GDL fibres. In addition, the effective transport properties of GDL for various liquid water saturations were also calculated. The paper suggests a new methodology of GDL impregnation based on the Fine-Mist model. The model assumes that artificial small droplets are introduced in a form of fine mist inside the GDL pores and due to capillary forces acting as driving forces, the water droplets coalesce and try to occupy the void spaces within the fibre GDL structure. As a result, the distribution, size and shape of water droplets within the GDL depend only on the GDL structure and on capillary forces and is not influenced by boundary condition effects, such as flooding areas close to the water injectors. Three GDL samples considered from different locations of GDL specimen are examined. The values of effective transport properties from these locations vary ~10% in case of dry GDL. As expected, increasing water saturation decreases permeability and diffusivity of GDL (due to blocking of open GDL pores), while the impregnating water increases the thermal conductivity (heat transfer occurs via water droplets as well). The results of simulations are compared with Pore-Morphology based methods. It is shown that the effect of
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
liquid water on gas permeability and diffusivity agrees with previous predictions. However the Fine-Mist approach applied to investigate the effects of water saturation on thermal conductivity shows considerably different behaviour as compared to PoreMorphology-Model. It is shown that the size of droplets and contact area between water droplets and solid fibres obtained by using the Pore-Morphology Model is considerably smaller than that obtained by using the Fine-Mist approach with Lattice Boltzmann simulation. Besides the geometrical shape of fibres, the approach applied in the current study takes into account the surface tension of the liquid water droplets. As a result, liquid water can be transported within the GDL volume and water droplets can coalesce following to pore shape. The current study provides the possibility that the water droplets due to surface tension can be entrapped in various shapes provided by the fibre structure, whereas the PMM model considers only “spherical” droplets inside the fibre structure that can not be deformed. Investigations show that the effect of liquid water saturation on thermal conductivity of the partially saturated GDL is more significant based on the Fine-Mist-Approach as compared to PMM. Hence the results obtained using the PMM impregnation method could possibly be under-estimated. References [1] D. Cheddie, N. Munroe, Review and comparison of approaches to proton exchange membrane fuel cell modeling, J. Power Sources 147 (2005) 72e84, http://dx.doi.org/10.1016/j.jpowsour.2005.01.003. [2] A.Z. Weber, R.L. Borup, R.M. Darling, P.K. Das, T.J. Dursch, W. Gu, et al., A critical review of modeling transport phenomena in polymer-electrolyte fuel cells, J. Electrochem. Soc. 161 (2014) F1254eF1299, http://dx.doi.org/ 10.1149/2.0751412jes. [3] C. Siegel, Review of computational heat and mass transfer modeling in polymer-electrolyte-membrane (PEM) fuel cells, Energy 33 (2008) 1331e1352, http://dx.doi.org/10.1016/j.energy.2008.04.015. [4] A. Bazylak, Liquid water visualization in PEM fuel cells: a review, Int. J. Hydrogen Energy 34 (2009) 3845e3857, http://dx.doi.org/10.1016/ j.ijhydene.2009.02.084. [5] A. Biyikoglu, Review of proton exchange membrane fuel cell models, Int. J. Hydrogen Energy 30 (2005) 1181e1212, http://dx.doi.org/10.1016/ j.ijhydene.2005.05.010. [6] C.-Y. Wang, Fundamental models for fuel cell engineering, Chem. Rev. 104 (2004) 4727e4766, http://dx.doi.org/10.1021/cr020718s. [7] T.M. Squires, S.R. Quake, Microfluidics: fluid physics at the nanoliter scale, Rev. Mod. Phys. 77 (2005) 977e1026. [8] N. Djilali, Computational modelling of polymer electrolyte membrane (PEM) fuel cells: challenges and opportunities, Energy 32 (2007) 269e280, http:// dx.doi.org/10.1016/j.energy.2006.08.007. [9] H. Lomax, T. Pulliam, D. Zingg, T. Kowalewski, Fundamentals of computational fluid dynamics, Appl. Mech. Rev. 55 (2002) B61, http://dx.doi.org/10.1115/ 1.1483340. [10] I.V. Zenyuk, P.K. Das, A.Z. Weber, Understanding impacts of catalyst-layer thickness on fuel-cell performance via mathematical modeling, J. Electrochem. Soc. 163 (2016) F691eF703, http://dx.doi.org/10.1149/ 2.1161607jes. [11] I.V. Zenyuk, A.Z. Weber, Understanding liquid-water management in PEFCs using x-ray computed tomography and modeling, ECS Trans. 69 (2015) 1253e1265, http://dx.doi.org/10.1149/06917.1253ecst. [12] J. Becker, C. Wieser, S. Fell, K. Steiner, A multi-scale approach to material modeling of fuel cell diffusion media, Int. J. Heat. Mass Transf. 54 (2011) 1360e1368, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2010.12.003. [13] J. Becker, R. Flückiger, M. Reum, F.N. Büchi, F. Marone, M. Stampanoni, Determination of material properties of gas diffusion layers: experiments and simulations using phase contrast tomographic microscopy, J. Electrochem. Soc. 156 (2009) B1175eB1181, http://dx.doi.org/10.1149/1.3176876. [14] N. Zamel, X. Li, J. Becker, A. Wiegmann, Effect of liquid water on transport properties of the gas diffusion layer of polymer electrolyte membrane fuel cells, Int. J. Hydrog. Energy 36 (2011) 5466e5478, http://dx.doi.org/10.1016/ j.ijhydene.2011.01.146. [15] N. Zamel, X. Li, J. Shen, J. Becker, A. Wiegmann, Estimating effective thermal conductivity in carbon paper diffusion media, Chem. Eng. Sci. 65 (2010) 3994e4006, http://dx.doi.org/10.1016/j.ces.2010.03.047. [16] N. Zamel, J. Becker, A. Wiegmann, Estimating the thermal conductivity and diffusion coefficient of the microporous layer of polymer electrolyte membrane fuel cells, J. Power Sources 207 (2012) 70e80, http://dx.doi.org/ 10.1016/j.jpowsour.2012.02.003. [17] P.P. Mukherjee, V.P. Schulz, J. Becker, E. Glatt, A. Wiegmann, C. Hartnig, et al.,
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]
[39] [40]
[41]
[42]
119
Polymer Electrolyte Membrane and Direct Methanol Fuel Cell Technology, Elsevier, 2012, http://dx.doi.org/10.1533/9780857095473.2.219. M.M. Daino, S.G. Kandlikar, 3D phase-differentiated GDL microstructure generation with binder and PTFE distributions, Int. J. Hydrogen Energy 37 (2012) 5180e5189, http://dx.doi.org/10.1016/j.ijhydene.2011.12.050. € tzke, G. Gaiselmann, M. Osenberg, J. Bohner, T. Arlt, H. Marko € tter, et al., C. To Three-dimensional study of compressed gas diffusion layers using synchrotron X-ray imaging, J. Power Sources 253 (2014) 123e131, http://dx.doi.org/ 10.1016/j.jpowsour.2013.12.062. A. Pfrang, D. Veyret, F. Sieker, G. Tsotridis, X-ray computed tomography of gas diffusion layers of PEM fuel cells: calculation of thermal conductivity, Int. J. Hydrogen Energy 35 (2010) 3751e3757, http://dx.doi.org/10.1016/ j.ijhydene.2010.01.085. A. Pfrang, S. Didas, G. Tsotridis, X-ray computed tomography of gas diffusion layers of {PEM} fuel cells: segmentation of the microporous layer, J. Power Sources 235 (2013) 81e86, http://dx.doi.org/10.1016/j.jpowsour.2013.01.179. A. Pfrang, D. Veyret, G.J.M. Janssen, G. Tsotridis, Imaging of membrane electrode assemblies of proton exchange membrane fuel cells by X-ray computed tomography, J. Power Sources 196 (2011) 5272e5276, http://dx.doi.org/ 10.1016/j.jpowsour.2010.09.020. M. Bosomoiu, G. Tsotridis, T. Bednarek, Study of effective transport properties of fresh and aged gas diffusion layers, J. Power Sources 285 (2015) 568e579, http://dx.doi.org/10.1016/j.jpowsour.2015.03.132. V.P. Schulz, J. Becker, A. Wiegmann, P.P. Mukherjee, C.-Y. Wang, Modeling of two-phase behavior in the gas diffusion medium of PEFCs via full morphology approach, J. Electrochem. Soc. 154 (2007) B419eB426. €lke, V.P. Schulz, M. Krafczyk, K. Roth, Comparison of a latticeH.-J. Vogel, J. To boltzmann model, a full-morphology model, and a pore network model for determining capillary pressureesaturation relationships, Vadose Zo. J. 4 (2005) 380e388. M. Hilpert, C.T. Miller, Pore-morphology-based simulation of drainage in totally wetting porous media, Adv. Water Resour. 24 (2001) 243e255, http:// dx.doi.org/10.1016/S0309-1708(00)00056-7. Y. Gao, X. Zhang, P. Rama, R. Chen, H. Ostadi, K. Jiang, Lattice Boltzmann simulation of water and gas flow in porous gas diffusion layers in fuel cells reconstructed from micro-tomography, Comput. Math. Appl. 65 (2013) 891e900, http://dx.doi.org/10.1016/j.camwa.2012.08.006. P.A. García-Salaberri, J.T. Gostick, G. Hwang, A.Z. Weber, M. Vera, Effective diffusivity in partially-saturated carbon-fiber gas diffusion layers: effect of local saturation and application to macroscopic continuum models, J. Power Sources 296 (2015) 440e453, http://dx.doi.org/10.1016/ j.jpowsour.2015.07.034. P.A. García-Salaberri, J.T. Gostick, G. Hwang, A.Z. Weber, M. Vera, Effective diffusivity in partially-saturated carbon-fiber gas diffusion layers: effect of local saturation and application to macroscopic continuum models, J. Power Sources 296 (2015) 440e453, http://dx.doi.org/10.1016/ j.jpowsour.2015.07.034. P.P. Mukherjee, C.-Y. Wang, Q. Kang, Mesoscopic modeling of two-phase behavior and flooding phenomena in polymer electrolyte fuel cells, Electrochim. Acta 54 (2009) 6861e6875, http://dx.doi.org/10.1016/ j.electacta.2009.06.066. X. He, L.-S. Luo, Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E 56 (1997) 6811e6817, http://dx.doi.org/10.1103/PhysRevE.56.6811. L. Chen, Q. Kang, Y. Mu, Y.-L.L. He, W.-Q.Q. Tao, A critical review of the pseudopotential multiphase lattice Boltzmann model: methods and applications, Int. J. Heat. Mass Transf. 76 (2014) 210e236, http://dx.doi.org/10.1016/ j.ijheatmasstransfer.2014.04.032. P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511e525, http://dx.doi.org/10.1103/PhysRev.94.511. X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815e1819, http://dx.doi.org/ 10.1103/PhysRevE.47.1815. X. Shan, H. Chen, Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E 49 (1994) 2941e2948, http://dx.doi.org/10.1103/PhysRevE.49.2941. X. Shan, Pressure tensor calculation in a class of nonideal gas lattice Boltzmann models, Phys. Rev. E 77 (2008) 66702, http://dx.doi.org/10.1103/ PhysRevE.77.066702. S. Hou, X. Shan, Q. Zou, G.D. Doolen, W.E. Soll, Evaluation of two lattice Boltzmann models for multiphase flows, J. Comput. Phys. 138 (1997) 695e713, http://dx.doi.org/10.1006/jcph.1997.5839. H. Huang, L. Wang, X. Lu, Evaluation of three lattice Boltzmann models for multiphase flows in porous media, Comput. Math. Appl. 61 (2011) 3606e3617, http://dx.doi.org/10.1016/j.camwa.2010.06.034. D.H. Rothman, J.M. Keller, Immiscible cellular-automaton fluids, J. Stat. Phys. 52 (1988) 1119e1127, http://dx.doi.org/10.1007/BF01019743. A.K. Gunstensen, D.H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A 43 (1991) 4320e4327, http://dx.doi.org/ 10.1103/PhysRevA.43.4320. Q. Kang, M. Wang, P.P. Mukherjee, P.C. Lichtner, Mesoscopic modeling of multiphysicochemical transport phenomena in porous media, Adv. Mech. Eng. 2 (2010), http://dx.doi.org/10.1155/2010/142879. K. Kono, T. Ishizuka, H. Tsuda, A. Kurosawa, Application of lattice Boltzmann
120
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
model to multiphase flows with phase transition, Comput. Phys. Commun. 129 (2000) 110e120, http://dx.doi.org/10.1016/S0010-4655(00)00098-9. [43] A. Parmigiani, C. Huber, O. Bachmann, B. Chopard, Pore-scale mass and reactant transport in multiphase porous media flows, J. Fluid Mech. 686 (2011) 40e76, http://dx.doi.org/10.1017/jfm.2011.268. [44] Sigracet®GDL34&35 Series - Datasheet, 2013. [45] L.M. Pant, Experimental and Theoretical Investigation of Mass Transport in
Porous Media of a PEM Fuel Cell, Faculty of Graduate Studies and Research, University of Alberta, 2011. [46] B. Poling, J. Prausnitz, J. O'Connell, Chemical Properties Handbook, McGrawHill, 2001. [47] C. Yaws, Chememical Properties Handbook, McGraw-Hill, 1999. [48] Fraunhofer ITWM, Department Flow and Complex Structures. GeoDict., Kaiserslautern, Germany {available at:, 2014 www.itwm.fraunhofer.de.
Contents lists available at ScienceDirect
Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour
Calculation of effective transport properties of partially saturated gas diffusion layers Tomasz Bednarek*, Georgios Tsotridis European Commission, Joint Research Centre (JRC), Directorate for Energy, Transport and Climate, Westerduinweg 3, NL-1755 LE, Petten, The Netherlands
h i g h l i g h t s Dry and partially saturated gas diffusion layers are investigated. The transport properties as a function of liquid water saturation are calculated. The fine mist approach is used to impregnate GDL without specimen boundary effects. Liquid water distribution depends on exact fibre geometry and capillary forces. Results of Two-Phase Lattice Boltzmann and Pore-Morphology simulations are compared.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 6 August 2016 Received in revised form 20 October 2016 Accepted 27 October 2016
A large number of currently available Computational Fluid Dynamics numerical models of Polymer Electrolyte Membrane Fuel Cells (PEMFC) are based on the assumption that porous structures are mainly considered as thin and homogenous layers, hence the mass transport equations in structures such as Gas Diffusion Layers (GDL) are usually modelled according to the Darcy assumptions. Application of homogenous models implies that the effects of porous structures are taken into consideration via the effective transport properties of porosity, tortuosity, permeability (or flow resistance), diffusivity, electric and thermal conductivity. Therefore, reliable values of those effective properties of GDL play a significant role for PEMFC modelling when employing Computational Fluid Dynamics, since these parameters are required as input values for performing the numerical calculations. The objective of the current study is to calculate the effective transport properties of GDL, namely gas permeability, diffusivity and thermal conductivity, as a function of liquid water saturation by using the Lattice-Boltzmann approach. The study proposes a method of uniform water impregnation of the GDL based on the “Fine-Mist” assumption by taking into account the surface tension of water droplets and the actual shape of GDL pores. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Keywords: Effective transport properties Gas diffusion layer Liquid water saturation Polymer electrolyte fuel cell Lattice-Boltzmann
1. Introduction A common assumption in using Computational Fluid Dynamics (CFD) approaches in modelling of Polymer Electrolyte Fuel Cells (PEMFC) entails the homogenization of the porous structures of the different fuel cell components, such as Gas Diffusion Layers (GDL), Micro-Porous Layers (MPL) and Catalyst Layers (CL). This assumption implies that these porous structures are considered as uniform and continuous, hence the statistical properties (namely porosity, pore size distribution and tortuosity) are used in characterising them. Such an approach for porous PEMFC components was investigated by a number of papers, among others, by Cheddie and
* Corresponding author. E-mail address: [email protected] (T. Bednarek).
Munroe [1], Weber at al. [2], Siegel [3], Bazylak [4], Biyikoglu [5], Wang [6], Squires [7], Djilali [8] and Lomax et al. [9]. In these papers, the porous PEMFC components are treated as homogenous structures by incorporating the Darcy approximation. Since the shape of liquid water droplets and that of the structure of GDL are not taken into account, their contribution to the mass transport within the GDL is neglected, which affects the accuracy of the PEMFC modelling [10,11]. Furthermore, a substantial research effort has been attributed to numerical reproduction of the porous structures of GDLs. Becker et al. [12,13], Zamel et al. [14e16], Mukherjee et al. [17], Daino and Kandlikar [18] propose techniques for reproducing the fibre€tzke et al. [19], Pfrang et al. [20e22], Bosstructure, whereas To omoiu et al. [23] use X-ray Computer Tomography (CT) for reproducing the geometry of the GDL microstructure. Bosomoiu et al.
http://dx.doi.org/10.1016/j.jpowsour.2016.10.098 0378-7753/© 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
112
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
[23] presented the full procedure of CT image acquisition and image manipulation to numerically reproduce real GDL structures. The influence of liquid water on transport properties of GDL was investigated, among others, by Zamel et al. [14], Bosomoiu et al. [23], Schultz et al. [24] using the Pore-Morphology Model (PMM) [25,26]. Stochastic numerical models for creating artificial 3D models of GDL or MPL are considered, among others, by Zamel at al [15,16]. and Becker et al. [12,13]. The reproduction of the GDL structure on the basis of CT imaging and subsequential establishment of effective GDL properties are done, among others, by Pfrang et al. [20], Bosomoiu et al. [23], Gao et al. [27] and García-Salaberri et al. [28,29]. The current study presents a methodology of uniform liquid water impregnating of a GDL sample. Liquid water is introduced as a “fine mist”, by assuming that the water droplets coalesce and subsequently fill the GDL pores. There are no external forces applied to the system, and the flow of liquid water is due to capillary forces only. The “fine mist” approach assumes that the GDL sample is uniformly impregnated with liquid water. The numerical representation of the GDL substrate without the micro-porous layer (MPL) is reproduced from X-ray CT images at a resolution of 2.2 mm, sufficient to determine the 3D anisotropic structure of the GDL fibres. Several liquid water saturation levels have been considered. In contrast to previously published literature [14,23,24,26,28,30] where the capillary forces have no effect on liquid droplets inside the GDL, the current study considers the Lattice-Boltzmann (LB) two-phase approach to simulate water droplet coalescence and movement inside the GDL. The effective transport properties of GDL, namely gas permeability, diffusivity and thermal conductivity, are evaluated based on the Lattice-Boltzmann analysis as a function of liquid water saturation. It is shown that the liquid water has an influence on the effective GDL properties. Results of LB simulations are compared with PoreMorphology Model (PMM) applied in Bosomoiu et al. [23].
Fig. 1. Lattice-Boltzmann fundamentals: a) the movement and collisions of particles, b) the scheme of the D2Q9 lattice.
The most popular Bhatnagar-Gross-Krook (BGK) collision operator [33] provides the standard LB method formula, including the external force term F
1 eq fa ðx þ vea Dt; t þ DtÞ fa ðx; tÞ ¼ fa ðx; tÞ fa ðx; tÞ þ Fa ðx; tÞ t n |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} transport
2. The Lattice Boltzmann method
collision
(3) The Lattice Boltzmann method is based on simulation of particle movement and collisions with other particles [31], shown in Fig. 1a. In particular, the LB method treats the fluid as a set of particle distribution functions, while the motion is described by collisions between particles. The LB simulation algorithm assumes two steps: transport and collision (see Fig. 1)
fa ðx þ vea Dt; tÞ ¼ fa ðx; tÞ fa ðx; t þ DtÞ ¼ fa ðx; tÞ þ Cðx; tÞ
transport collision
(1)
where fa ðx; tÞ is the density distribution function at lattice position x and time t, fa is the density distribution function after transport process, v¼Dx/Dt is the velocity of the particle, a refers to the consecutive neighbours of the central node e0 (see directions in Fig. 1b). C is a collision operator. The discrete velocity depends on the assumed lattice model, i.e. D2Q9 or D3Q19 (it means 2-dimensions 9-nodes or 3-dimensions, 19-nodes respectively). The scheme of D2Q9 lattice is shown in Fig. 1b. From the D2Q9 lattice the discrete velocity vector ea is given as follows [32,33].
8 > > > > > <
ð0 0Þ ða 1Þp ða 1Þp cos sin 2 2 ea ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > ða 5Þp ða 5Þp p > > sin þ : 2 cos 2 2 4
eq
where fa ðx; tÞ is an equilibrium distribution function, tn is the dimensionless relaxation time. The left-hand-side of equation (3) stands for the transport step whereas the right-hand-side of equation (3) corresponds to the collision step (Fig. 1a). The collision term leads to the local Maxwellian equilibrium on a time scale t [32,33]. eq The equilibrium distribution function fa ðx; tÞ for a D2Q9 lattice take the form [32].
3 9 3 eq fa ðx; tÞ ¼ ua r 1 þ 2 ðea $uÞ þ 4 ðea $uÞ2 2 u2 c 2c 2c
The weight factors for D2Q9 lattice are as follows: u0¼4/9, u1and u5-8¼1/36. The fluid density ra and the fluid velocity u are given by the first and the second moments of the density distribution function [32].
4¼1/9
r¼
P
fa aP
ru ¼
a¼0 a ¼ 1; 2; 3; 4 a ¼ 5; 6; 7; 8 (2)
(4)
a
fa ea
(5)
The viscosity of the fluid (in lattice units) is related to the relaxation time by formula
tn ¼ c2s ðtn 0:5ÞDt where cs is the sound speed in lattice units.
(6)
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
following non-ideal equation of state in the lattice units [32,34,35].
3. Two-phase model The gaseous and liquid phases are modelled as two separate lattices with an interface zone between them, Fig. 2a. Originally, Shan and Chen defined the non-local interaction among particles as a force between particles at position x and x' [34,35]
Fðx; x0 Þ ¼ Gðjx x0 jÞjðxÞjðx0 Þðx0 xÞ
(7)
where G is Green's function and j is an effective mass depending on the local density. The force basically arises from cohesion pairwise point coupling between effective mass at location x and neighbouring location x', Fig. 2b. The force given by equation (7) does not conserve the local momentum during the collision. However it has been proved in Ref. [35] that the total momentum of the system is conserved e no net momentum is introduced into the system. The total force acting on the particle at position x is the sum of all interaction forces. In the lattice system, i.e. mentioned D2Q9, constant number of close neighbours N is assumed (Fig. 1b) that interacts with the particle at location x. The interaction force given in equation (7) can be rewritten after [32].
FðxÞ ¼ g jðxÞc2s
N X
a¼1
w jea j2 jðx þ ea Þea
(8)
where g is the cohesion interaction strength between neighbouring particles. The weights wðjea j2 Þ are introduced to account for various distances between particles which maintaining an isotropic reaction force. A positive value of the interaction strength represents a repulsive force between neighbours, whereas negative g corresponds to attractive interaction. Shan and Chen proposed the following form of effective mass [34,35].
j ¼ r0 1 exp
r1 r2
113
(9)
where r1 and r2 are densities of both phases, r0 is a normalization constant usually assumed as 1.0 [32]. The Taylor expansion of the interaction force given in equation (8) leads to the form given by Ref. [36].
c4 FðxÞ ¼ g c2s jVj þ s jVðDjÞ þ … 2
(10)
The interaction force given in equation (10) covers a non-ideal equation of state (the first term) and surface tension (the term corresponding to Dj). The correction of the ideal equation of state (namely p¼r c2s ) due to force given in formula (10) leads to the
g p ¼ rc2s þ c2s j2 2
(11)
Hou at al. and Huang et al. in Refs. [37,38] have validated the presented pseudo-potential model by simulation of a bubble in fully periodic system. It is shown that the presented model correctly simulates the two-phase interaction, i.e. a major improvement has been noticed with comparison to RothmanKeller model [39,40]. For reason of simplicity and computational efficiency the twophase model based on cohesion forces has been widely applied in simulating a variety multiphase flows. However the model developed by Hou et al. and Huang et al. in Refs. [37,38] suffers from some limitations highlighted by Chen et al. [32]. One limitation arises from the relatively low density and kinematic viscosity (relaxation time in particular) ratios. Mukherjee et al. and Kang et al. [30,41] have shown that the density ratio does not play significant role for micro-porous structures since the capillary pressure is the main driving force of liquid flow inside GDL. However to improve the simulation accuracy a dynamic model has been adopted [42], where the particle velocity u is scaled by the density ratio k.
u¼
u1 þ ku2 ; 1þk
k¼
r2 r1
(12)
where subscripts 1 and 2 correspond to separate fluid phases.
4. Calculation of effective transport properties The calculated transport parameters such as permeability, diffusivity and thermal conductivity are considered as scalar values along the direction x through GDL thickness, shown in Fig. 3. The simulations of gas flow, diffusion and heat transfer are performed by using the LB. The transient LB calculations are terminated when the difference between the normalized quantities (velocity, diffusion flux or heat flux in case of permeability, diffusivity and thermal conductivity respectively) from the current and previous time step are smaller than the assumed accuracy ε ¼ 1.0e-4. The calculation of permeability using LB is based on the singlephase flow simulation through the GDL sample. The current study considers dry air as the gas phase. It is assumed that there is no mass exchange between gas and liquid droplets occupying GDL pores e no evaporation/condensation. The effective gas permeability is expressed as follows
keff ¼
um Dp
(13)
where keff is the saturation dependent effective permeability, u is
Fig. 2. The two-phase modelling: a) division of lattices, b) cohesion forces in liquid phase.
114
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
Fig. 3. Location of investigated samples inside the Sigracet®GDL34BC specimen (left hand side) and computational domain (right hand side).
the average out-of-plane velocity component, m is the viscosity and Dp is the pressure gradient between the top and bottom sides of the GDL sample, Fig. 3. In the current study the pressure gradient Dp is assumed equal to 1.0 [kPa]. The effective diffusivity of the GDL partially saturated with liquid water is calculated on the basis of simulation of oxygen diffusion in air. The through plane effective diffusivity Deff is calculated according to Fick's law [16].
Deff ¼
j
Dc
(14)
where j is the out-of-plane component of the average diffusion flux, Dc is the oxygen concentration difference between top and bottom planes of the GDL sample Fig. 3. The effective thermal conductivity of the GDL partially saturated with liquid water is calculated according to Fourier's law
q
leff ¼ DT
active area (2300 2300 pixels). To minimise the statistical error in calculating the effective properties of the GDL, three samples were considered form different locations of the GDL specimen, shown in Fig. 3. The Sigracet® GDL34BC contains the micro-porous layer (MPL), however as it is shown in Ref. [23], the employed CT scan does not have the sufficient resolution to reliably reproduce the MPL since the pore sizes are of the order of nm [16]. Therefore, the MPL was excluded from considerations in this study. The volume of the investigated GDL sample is shown in Fig. 3. The boundary conditions corresponding to gradients of pressure, concentration and temperature (for calculating the permeability, diffusivity and thermal conductivity respectively) are applied at the bottom and the top of the GDL surfaces. All other sides of the sample are assumed as solid walls with no slip (in case of flow simulation) or no flux conditions (for diffusion and heat transfer). 5.1. Water impregnation of GDL
(15)
where q is the out-of-plane component (along x axis, Fig. 3) of the average heat flux at outlet plane, DT is temperature difference between top and bottom planes of the GDL sample, Fig. 3. The liquid droplets occupying the GDL pores are assumed to be rigid hence the model of heat transfer through the GDL takes only the conduction through the carbon fibres and water droplets into account and there is no heat transfer through the empty GDL pores. The heat transfer by radiation is neglected because the normal operating temperature of PEM fuel cell is lower than 100OC.
5. Results The present study considers three effective transport parameters of porous GDLs, namely: gas permeability, diffusivity and thermal conductivity as a function of liquid water saturation. The Sigracet® GDL34BC material is investigated in this study, which was also considered by Bosomoiu et al. [23]. The material was analysed by using the X-ray computer tomography system Phoenix Nanotom S which includes a 180 [kV] nanofocus X-ray tube and a 2D X-ray detector with 120 120 [mm]
During fuel cell operation water vapour diffuses and may condensate inside the diffusion media, hence a number of pores could be occupied by liquid water droplets. However the real conditions of fuel cell operation could not reproduce the contribution of liquid water distribution inside the GDL. During fuel cell operation the flow of reactant gasses will affect the accumulation of liquid water inside GDL due to pressure drop. The convective flow of the humidified gas could possibly remove the liquid water droplets or could force them to move to locations close to boundaries, implying a non-uniform liquid water distribution. However, intrusion of liquid water into a dry GDL results in a non-uniform distribution of liquid water inside the GDL volume. The area close to the uniform inlet is flooded by water, as shown in Ref. [43]. The substitution of the uniform inlet surface by an inlet spot, or a number of inlet spots does not address the problem of a flooded inlet region. The local inlet flooding effects are clearly visible in Fig. 4, where a LB simulation of water injection into GDL is presented and will affect the effective transport properties of the GDL. To avoid the effects of boundary conditions on water distribution inside the GDL an alternative method of water impregnation is proposed. It is assumed that there are no liquid water injectors at
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
115
Fig. 4. Local inlet effects, the lower layers of the GDL are flooded with water.
Fig. 5. Lattice-Boltzmann simulations: a) Two-phase simulation of water droplets creation from fine mist to droplet coalescence. b) Streamlines of gas flow through partially saturated GDL, cross-section through GDL sample.
116
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
the specimen boundaries, therefore all sides of the specimen are impermeable walls. Instead of an inlet boundary condition, the GDL pores are artificially saturated by a fine water mist, i.e. very small water droplets of size corresponding to the size of the voxel (2.2 mm). The water droplets are introduced randomly throughout the whole volume of the sample. To avoid formation of suspended droplets, each artificially introduced water droplet either attach itself to the GDL fibres or would coalese with other water droplets. The LB two-phase model is applied, where both the movement of the liquid phase and droplet coalesce caused by capillary forces due to surface tension and contact angle with GDL fibres, are considered. There are no other external body forces acting on the system. The coalescence and motion of liquid droplets in the GDL sample are shown in Fig. 5a. The comparison of liquid droplets size distribution at two water saturation levels of 25% and 75% is presented in Fig. 6a. Solid lines correspond to Lattice-Boltzmann simulations whereas dotted lines refer to Pore-Morphology Model (PMM) [24,25]. It is shown, that droplets obtained by LB simulation are larger than those obtained by PMM. In the case of 25% water saturation using PMM most of the droplets are 35e65 mm in size, whereas the range of droplet sizes for LB is wider and ranges from 16 up to 85 mm. Such a difference is caused by inclusion of surface tension phenomena and the motion of droplets is due to capillary forces in case of LB simulations, contrary to PMM which does not allow movement and coalescence of neighbouring droplets. This effect is even more visible at higher water saturation level of 75%. In case of PMM only new smaller droplets (form 15e35 mm in size) are introduced, whereas in case of LB simulation, the droplets coalesce surround the fibres and cover a larger volume of GDL sample e the droplets size ranges form 20 up to 140 mm. 5.2. Effective transport properties of the GDL as a function of liquid water saturation The effective permeability of GDL as a function of liquid water saturation is calculated for air at ambient pressure with temperature 65 C (density of 1.050 [kg m3] and kinematic viscosity of 2.15e-5 [Pa s]). The driving force of flow is the pressure difference between the top and the bottom surfaces of the GDL sample and is equal to 1.0 [kPa], Fig. 3. The streamlines of gas flow through the GDL sample are presented in Fig. 5b. The through-plane permeability (x axis in Fig. 3) of dry GDL is
calculated by using formula (13), and it ranges form 3.20e-12 [m2] to 3.55e-12 [m2] corresponding to sample 2 and 3 respectively. These values are similar to those provided by the manufacturer in the GDL “Properties data sheet” [44], namely 3.65e-12 [m2] and to the value of 3.17e-12 [m2] obtained by the authors in a previous investigation [23]. The corresponding distribution of permeability as a function of water saturation is presented in Fig. 7a and Table 1. The relative difference between the permeability values from the samples taken from various locations of the GDL structure does not exceed 10% in case of dry GDL and is less than 3% for the case of 90% of liquid water saturation. As expected, the gas permeability strongly depends on the water saturation of the GDL. Significant decrease of GDL permeability appears when the water saturation level exceeds 50%, Fig. 7a. The reason of such a large decrease in permeability corresponds to liquid water droplets occupying the largest pores of GDL, hence blocking the convective flow through them, as ilustrated in Fig. 6b where the pore size distribution of a dry sample and at two liquid water saturation levels of 25% and 50% respectively are presented. It is seen that in case of saturated GDL the number of the largest open pores decreases, because these pores are occupied by liquid water droplets, hence forcing the gas to pass mainly through smaller pores with a higher flow resistance. The effective diffusivity is calculated by using equation (14). The oxygen concentrations at the top and the bottom surfaces of the samples are fixed. On the top surface oxygen concentration corresponds to fully humidified air, namely 13.5 mol m3, while the oxygen concentration at the bottom surface is equal to 12.0 mol m3. The simulation is performed at temperature T ¼ 65OC. The calculated values of effective diffusivities of dry GDL range from 6.50e-6 m2 s1 to 7.13e-6 m2 s1, and are comparable to the value of 6.6e-6 m2 s1 presented by a previous study [23]. However, they are different to the experimental value of 5.9e-6 m2 s1 obtained by Pant [45] because the experiment by Pant was performed at 28 C whereas both current study and [23] use 65OC, a more typical operating temperature for fuel cell applications. The effect of liquid water saturation on the effective diffusivity is shown in Fig. 7b. Similar to gas permeability, the values for dry GDL differ about 10% and that the difference decreases with increasing saturation level. It is shown that the liquid water saturation has a major impact to the effective diffusivity. It is observed that 20% of water saturation considerably reduces the diffusivity by over 50%, whereas 50% of water saturation decreases diffusivity substantially. The diffusion coefficient of GDL where half of pore space is
Fig. 6. Size distribution of: a) water droplets (comparison between Pore-Morphology and Lattice-Boltzmann methods) and b) open pore size distribution.
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
Fig. 7. Effective properties of GDL as a function of liquid water saturation: a)gas permeability, b) gas diffusivity, c) thermal conductivity.
117
occupied by liquid water (saturation 50%) is about 7 times smaller (1.02e-6 m2 s1) than the diffusion coefficient of the dry GDL. The reason of such significant reduction of the diffusion coefficient is attributed, as in the case of permeability, to the location of liquid water droplets inside the GDL. Since the largest pores are occupied by liquid water (Fig. 6b) the diffusion process is forced to take place through the smaller pores, thereby been significantly reduced. The effective thermal conductivity coefficient is calculated by formula (15). Contrary to the previous cases, thermal conductivity is affected by both the solid fibres and the liquid water droplets. The conductivities of the individual phases, namely carbon fibres and liquid water, at 65 C are assumed to be equal 128.95 W m1K1 and 0.6504 W m1K1 respectively [46,47]. The calculated thermal conductivity coefficient for dry GDL samples varies from 0.83 W m1K1 to 0.92 W m1K1 which is similar to previous investigations [23]. The difference between the most conductive sample 2 and the less conductive sample 3 is about 10%. It is also noticed that with increasing water saturation level the differences between thermal conductivities of samples decreases, since more heat is transported via the liquid water droplets occupying GDL pores. Fig. 7c shows the effect of liquid water on the thermal conductivity of GDL. It is noticed that the effective thermal conductivity of the partially saturated GDL increases with increasing saturation. At 30% water saturation the effective thermal conduction coefficient is nearly doubled as compared to the dry GDL. Further impregnation of GDL does not affect the effective conductivity so strongly, due to the fact that the largest pores are already occupied by liquid water droplets, hence further impregnation fills relatively smaller pores, where most of the heat transfer is performed via the more conductive carbon fibres. The effective conductivity value at increasing liquid water saturation levels tends to the value of approximately 2.0 W m1 K1. However, the current study depicts a different water effect than that presented by Bosomoiu et al. [23]. In an effort to explain the discrepancies of the effective thermal coefficients, the different water saturation levels are re-calculated using the same GeoDict software [48], as in the case of Bosomoiu et al. [23]. The plot of the effective thermal conductivity as a function of water saturation obtained by using Geodict and LB software is also presented in Fig. 7c. The thermal conductivities calculated by both GeoDict and LB are similar for a dry GDL sample. However, the effect of liquid water on thermal conductivity between the current study and [23] (Fig. 7c) arises from different water impregnation methods. In the paper by Bosomoiu et al. [23] the Pore-Morphology impregnation method was used [25,26], whereas in the current study the “ finemist” approach with the two-phase LB simulations have been applied. The differences of liquid water distribution obtained by using the “ fine-mist” approach and the PMM at the pore scale are presented in Fig. 8. The Pore-Morphology Model (Fig. 8a) bases on specifying the spherical water droplet diameter which is then compared with the GDL pores. In cases where the droplet fits to the GDL pore, the droplet remains in the pore, implying that the pore is
Table 1 Effective properties of partially saturated GDL. Liquid water saturation
Gas permeability [m1]
Diffusivity [m2s1]
Thermal conductivity [W m1 K1]
Reference [23]
Current study (average)
Reference [23]
Current study (average)
Reference [23]
Current study (average)
dry 25% 50% 75%
3.14E-12 9.42E-13 3.14E-13 6.28E-14
3.38E-12 8.04E-13 2.30E-13 3.84E-14
6.517E-06 3.375E-06 1.406E-06 1.875E-07
6.83E-06 3.35E-06 1.03E-06 6.83E-07
0.908 1.076 1.380 1.776
0.877 1.513 1.839 1.952
118
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
Fig. 8. Impregnated GDL using: a) Pore-morphology method, b) Lattice-Boltzmann method. In both cases the liquid water saturation is equal ~10%. Red lines indicate the contact area between liquid water (blue) and fibres (grey). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
occupied by liquid water [25,26]. However, from Fig. 8a is shown, that the spherical droplets considered under PMM imply that the contact area between droplet and fibres is considerably smaller than for the LB simulations where capillary interactions are taken into account. Such behaviour of the PMM can lead to underestimating the effects of liquid water saturation on the effective thermal conductivity because the contact area between GDL fibres and water droplets plays a crucial role. In case of the two-phase simulation by LB, the shape of the water droplets is according to real shape of the GDL pores and is influenced by the values of water surface tension. In the presented approach the water droplets can be transported due to capillary forces, even from very inaccessible locations. In particular, moving droplets between GDL pores can coalesce creating flooded zones surrounding some fibres. As a result, the droplets created by “finemist” can lead to formation of droplets larger than the GDL pores themselves. Whereas for the PMM model the formation of droplets larger that the GDL pores is not allowed. Fig. 8b shows that the contact area between droplets and carbon fibres is considerably higher in case of LB two-phase simulation than the case of PMM, hence the effect of water saturation is more significant in the case of thermal conductivity. However, the differences between impregnation methods, employed by PMM and Fine-Mist approach based on LB simulation, do not play a significant role when calculating the effective permeability and diffusivity. This is due to the fact that the mass transport due to pressure or concentration difference takes place in pore space where the influence of contact area between fibres and water droplets is not significant. However, as it was shown in Fig. 7 the effect of water impregnation by the fine mist approach on gas permeability and diffusivity is greater as compared to the pore morphology method. The reason is that for LB simulation the shape of droplets fit to the GDL pores structure, Fig. 8b, closing the pores more tightly as compared to PMM, Fig. 8a. Nevertheless, the differences in the values of the calculated effective properties
between the two methods are not as large as in the case of the thermal conductivity. 6. Conclusions The study shows the effect of liquid water saturation on the effective transport properties of GDL, namely: gas permeability, diffusivity and thermal conductivity, based on the computer tomography image of Sigracet® GDL34BA (without MPL) material. The Fine-Mist-Approach together with two-phase Lattice-Boltzmann simulation was applied to calculate the liquid water distribution inside the GDL sample. The current work, in contrast to Pore-Morphology-Method, considers the liquid water distribution inside the GDL by taking into account the exact geometry of the GDL fibre structure and the capillary forces due to the surface tension of liquid water and contact angle with GDL fibres. In addition, the effective transport properties of GDL for various liquid water saturations were also calculated. The paper suggests a new methodology of GDL impregnation based on the Fine-Mist model. The model assumes that artificial small droplets are introduced in a form of fine mist inside the GDL pores and due to capillary forces acting as driving forces, the water droplets coalesce and try to occupy the void spaces within the fibre GDL structure. As a result, the distribution, size and shape of water droplets within the GDL depend only on the GDL structure and on capillary forces and is not influenced by boundary condition effects, such as flooding areas close to the water injectors. Three GDL samples considered from different locations of GDL specimen are examined. The values of effective transport properties from these locations vary ~10% in case of dry GDL. As expected, increasing water saturation decreases permeability and diffusivity of GDL (due to blocking of open GDL pores), while the impregnating water increases the thermal conductivity (heat transfer occurs via water droplets as well). The results of simulations are compared with Pore-Morphology based methods. It is shown that the effect of
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
liquid water on gas permeability and diffusivity agrees with previous predictions. However the Fine-Mist approach applied to investigate the effects of water saturation on thermal conductivity shows considerably different behaviour as compared to PoreMorphology-Model. It is shown that the size of droplets and contact area between water droplets and solid fibres obtained by using the Pore-Morphology Model is considerably smaller than that obtained by using the Fine-Mist approach with Lattice Boltzmann simulation. Besides the geometrical shape of fibres, the approach applied in the current study takes into account the surface tension of the liquid water droplets. As a result, liquid water can be transported within the GDL volume and water droplets can coalesce following to pore shape. The current study provides the possibility that the water droplets due to surface tension can be entrapped in various shapes provided by the fibre structure, whereas the PMM model considers only “spherical” droplets inside the fibre structure that can not be deformed. Investigations show that the effect of liquid water saturation on thermal conductivity of the partially saturated GDL is more significant based on the Fine-Mist-Approach as compared to PMM. Hence the results obtained using the PMM impregnation method could possibly be under-estimated. References [1] D. Cheddie, N. Munroe, Review and comparison of approaches to proton exchange membrane fuel cell modeling, J. Power Sources 147 (2005) 72e84, http://dx.doi.org/10.1016/j.jpowsour.2005.01.003. [2] A.Z. Weber, R.L. Borup, R.M. Darling, P.K. Das, T.J. Dursch, W. Gu, et al., A critical review of modeling transport phenomena in polymer-electrolyte fuel cells, J. Electrochem. Soc. 161 (2014) F1254eF1299, http://dx.doi.org/ 10.1149/2.0751412jes. [3] C. Siegel, Review of computational heat and mass transfer modeling in polymer-electrolyte-membrane (PEM) fuel cells, Energy 33 (2008) 1331e1352, http://dx.doi.org/10.1016/j.energy.2008.04.015. [4] A. Bazylak, Liquid water visualization in PEM fuel cells: a review, Int. J. Hydrogen Energy 34 (2009) 3845e3857, http://dx.doi.org/10.1016/ j.ijhydene.2009.02.084. [5] A. Biyikoglu, Review of proton exchange membrane fuel cell models, Int. J. Hydrogen Energy 30 (2005) 1181e1212, http://dx.doi.org/10.1016/ j.ijhydene.2005.05.010. [6] C.-Y. Wang, Fundamental models for fuel cell engineering, Chem. Rev. 104 (2004) 4727e4766, http://dx.doi.org/10.1021/cr020718s. [7] T.M. Squires, S.R. Quake, Microfluidics: fluid physics at the nanoliter scale, Rev. Mod. Phys. 77 (2005) 977e1026. [8] N. Djilali, Computational modelling of polymer electrolyte membrane (PEM) fuel cells: challenges and opportunities, Energy 32 (2007) 269e280, http:// dx.doi.org/10.1016/j.energy.2006.08.007. [9] H. Lomax, T. Pulliam, D. Zingg, T. Kowalewski, Fundamentals of computational fluid dynamics, Appl. Mech. Rev. 55 (2002) B61, http://dx.doi.org/10.1115/ 1.1483340. [10] I.V. Zenyuk, P.K. Das, A.Z. Weber, Understanding impacts of catalyst-layer thickness on fuel-cell performance via mathematical modeling, J. Electrochem. Soc. 163 (2016) F691eF703, http://dx.doi.org/10.1149/ 2.1161607jes. [11] I.V. Zenyuk, A.Z. Weber, Understanding liquid-water management in PEFCs using x-ray computed tomography and modeling, ECS Trans. 69 (2015) 1253e1265, http://dx.doi.org/10.1149/06917.1253ecst. [12] J. Becker, C. Wieser, S. Fell, K. Steiner, A multi-scale approach to material modeling of fuel cell diffusion media, Int. J. Heat. Mass Transf. 54 (2011) 1360e1368, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2010.12.003. [13] J. Becker, R. Flückiger, M. Reum, F.N. Büchi, F. Marone, M. Stampanoni, Determination of material properties of gas diffusion layers: experiments and simulations using phase contrast tomographic microscopy, J. Electrochem. Soc. 156 (2009) B1175eB1181, http://dx.doi.org/10.1149/1.3176876. [14] N. Zamel, X. Li, J. Becker, A. Wiegmann, Effect of liquid water on transport properties of the gas diffusion layer of polymer electrolyte membrane fuel cells, Int. J. Hydrog. Energy 36 (2011) 5466e5478, http://dx.doi.org/10.1016/ j.ijhydene.2011.01.146. [15] N. Zamel, X. Li, J. Shen, J. Becker, A. Wiegmann, Estimating effective thermal conductivity in carbon paper diffusion media, Chem. Eng. Sci. 65 (2010) 3994e4006, http://dx.doi.org/10.1016/j.ces.2010.03.047. [16] N. Zamel, J. Becker, A. Wiegmann, Estimating the thermal conductivity and diffusion coefficient of the microporous layer of polymer electrolyte membrane fuel cells, J. Power Sources 207 (2012) 70e80, http://dx.doi.org/ 10.1016/j.jpowsour.2012.02.003. [17] P.P. Mukherjee, V.P. Schulz, J. Becker, E. Glatt, A. Wiegmann, C. Hartnig, et al.,
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]
[39] [40]
[41]
[42]
119
Polymer Electrolyte Membrane and Direct Methanol Fuel Cell Technology, Elsevier, 2012, http://dx.doi.org/10.1533/9780857095473.2.219. M.M. Daino, S.G. Kandlikar, 3D phase-differentiated GDL microstructure generation with binder and PTFE distributions, Int. J. Hydrogen Energy 37 (2012) 5180e5189, http://dx.doi.org/10.1016/j.ijhydene.2011.12.050. € tzke, G. Gaiselmann, M. Osenberg, J. Bohner, T. Arlt, H. Marko € tter, et al., C. To Three-dimensional study of compressed gas diffusion layers using synchrotron X-ray imaging, J. Power Sources 253 (2014) 123e131, http://dx.doi.org/ 10.1016/j.jpowsour.2013.12.062. A. Pfrang, D. Veyret, F. Sieker, G. Tsotridis, X-ray computed tomography of gas diffusion layers of PEM fuel cells: calculation of thermal conductivity, Int. J. Hydrogen Energy 35 (2010) 3751e3757, http://dx.doi.org/10.1016/ j.ijhydene.2010.01.085. A. Pfrang, S. Didas, G. Tsotridis, X-ray computed tomography of gas diffusion layers of {PEM} fuel cells: segmentation of the microporous layer, J. Power Sources 235 (2013) 81e86, http://dx.doi.org/10.1016/j.jpowsour.2013.01.179. A. Pfrang, D. Veyret, G.J.M. Janssen, G. Tsotridis, Imaging of membrane electrode assemblies of proton exchange membrane fuel cells by X-ray computed tomography, J. Power Sources 196 (2011) 5272e5276, http://dx.doi.org/ 10.1016/j.jpowsour.2010.09.020. M. Bosomoiu, G. Tsotridis, T. Bednarek, Study of effective transport properties of fresh and aged gas diffusion layers, J. Power Sources 285 (2015) 568e579, http://dx.doi.org/10.1016/j.jpowsour.2015.03.132. V.P. Schulz, J. Becker, A. Wiegmann, P.P. Mukherjee, C.-Y. Wang, Modeling of two-phase behavior in the gas diffusion medium of PEFCs via full morphology approach, J. Electrochem. Soc. 154 (2007) B419eB426. €lke, V.P. Schulz, M. Krafczyk, K. Roth, Comparison of a latticeH.-J. Vogel, J. To boltzmann model, a full-morphology model, and a pore network model for determining capillary pressureesaturation relationships, Vadose Zo. J. 4 (2005) 380e388. M. Hilpert, C.T. Miller, Pore-morphology-based simulation of drainage in totally wetting porous media, Adv. Water Resour. 24 (2001) 243e255, http:// dx.doi.org/10.1016/S0309-1708(00)00056-7. Y. Gao, X. Zhang, P. Rama, R. Chen, H. Ostadi, K. Jiang, Lattice Boltzmann simulation of water and gas flow in porous gas diffusion layers in fuel cells reconstructed from micro-tomography, Comput. Math. Appl. 65 (2013) 891e900, http://dx.doi.org/10.1016/j.camwa.2012.08.006. P.A. García-Salaberri, J.T. Gostick, G. Hwang, A.Z. Weber, M. Vera, Effective diffusivity in partially-saturated carbon-fiber gas diffusion layers: effect of local saturation and application to macroscopic continuum models, J. Power Sources 296 (2015) 440e453, http://dx.doi.org/10.1016/ j.jpowsour.2015.07.034. P.A. García-Salaberri, J.T. Gostick, G. Hwang, A.Z. Weber, M. Vera, Effective diffusivity in partially-saturated carbon-fiber gas diffusion layers: effect of local saturation and application to macroscopic continuum models, J. Power Sources 296 (2015) 440e453, http://dx.doi.org/10.1016/ j.jpowsour.2015.07.034. P.P. Mukherjee, C.-Y. Wang, Q. Kang, Mesoscopic modeling of two-phase behavior and flooding phenomena in polymer electrolyte fuel cells, Electrochim. Acta 54 (2009) 6861e6875, http://dx.doi.org/10.1016/ j.electacta.2009.06.066. X. He, L.-S. Luo, Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation, Phys. Rev. E 56 (1997) 6811e6817, http://dx.doi.org/10.1103/PhysRevE.56.6811. L. Chen, Q. Kang, Y. Mu, Y.-L.L. He, W.-Q.Q. Tao, A critical review of the pseudopotential multiphase lattice Boltzmann model: methods and applications, Int. J. Heat. Mass Transf. 76 (2014) 210e236, http://dx.doi.org/10.1016/ j.ijheatmasstransfer.2014.04.032. P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511e525, http://dx.doi.org/10.1103/PhysRev.94.511. X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815e1819, http://dx.doi.org/ 10.1103/PhysRevE.47.1815. X. Shan, H. Chen, Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E 49 (1994) 2941e2948, http://dx.doi.org/10.1103/PhysRevE.49.2941. X. Shan, Pressure tensor calculation in a class of nonideal gas lattice Boltzmann models, Phys. Rev. E 77 (2008) 66702, http://dx.doi.org/10.1103/ PhysRevE.77.066702. S. Hou, X. Shan, Q. Zou, G.D. Doolen, W.E. Soll, Evaluation of two lattice Boltzmann models for multiphase flows, J. Comput. Phys. 138 (1997) 695e713, http://dx.doi.org/10.1006/jcph.1997.5839. H. Huang, L. Wang, X. Lu, Evaluation of three lattice Boltzmann models for multiphase flows in porous media, Comput. Math. Appl. 61 (2011) 3606e3617, http://dx.doi.org/10.1016/j.camwa.2010.06.034. D.H. Rothman, J.M. Keller, Immiscible cellular-automaton fluids, J. Stat. Phys. 52 (1988) 1119e1127, http://dx.doi.org/10.1007/BF01019743. A.K. Gunstensen, D.H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A 43 (1991) 4320e4327, http://dx.doi.org/ 10.1103/PhysRevA.43.4320. Q. Kang, M. Wang, P.P. Mukherjee, P.C. Lichtner, Mesoscopic modeling of multiphysicochemical transport phenomena in porous media, Adv. Mech. Eng. 2 (2010), http://dx.doi.org/10.1155/2010/142879. K. Kono, T. Ishizuka, H. Tsuda, A. Kurosawa, Application of lattice Boltzmann
120
T. Bednarek, G. Tsotridis / Journal of Power Sources 340 (2017) 111e120
model to multiphase flows with phase transition, Comput. Phys. Commun. 129 (2000) 110e120, http://dx.doi.org/10.1016/S0010-4655(00)00098-9. [43] A. Parmigiani, C. Huber, O. Bachmann, B. Chopard, Pore-scale mass and reactant transport in multiphase porous media flows, J. Fluid Mech. 686 (2011) 40e76, http://dx.doi.org/10.1017/jfm.2011.268. [44] Sigracet®GDL34&35 Series - Datasheet, 2013. [45] L.M. Pant, Experimental and Theoretical Investigation of Mass Transport in
Porous Media of a PEM Fuel Cell, Faculty of Graduate Studies and Research, University of Alberta, 2011. [46] B. Poling, J. Prausnitz, J. O'Connell, Chemical Properties Handbook, McGrawHill, 2001. [47] C. Yaws, Chememical Properties Handbook, McGraw-Hill, 1999. [48] Fraunhofer ITWM, Department Flow and Complex Structures. GeoDict., Kaiserslautern, Germany {available at:, 2014 www.itwm.fraunhofer.de.