- Email: [email protected]
An integrated simulator of structure and anisotropic flow in gas diffusion layers with hydrophobic additives
Journal of Power Sources 365 (2017) 179e189
Contents lists available at ScienceDirect
Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour
An integrated simulator of structure and anisotropic flow in gas diffusion layers with hydrophobic additives Vasilis N. Burganos a, *, Eugene D. Skouras a, b, Alexandros N. Kalarakis a, b a b
Institute of Chemical Engineering Sciences, Foundation for Research and Technology e Hellas, Patras 26504, Greece Department of Mechanical Engineering, Technological-Educational Institute of Western Greece, Patras 26334, Greece
h i g h l i g h t s Three dimensional reconstruction of gas diffusion layers from SEM is achieved. Addition of binder and hydrophobic agents is simulated using a two-phase flow model. Flow anisotropy index is predicted for different PTFE loadings and compression levels.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 14 March 2017 Received in revised form 13 June 2017 Accepted 18 August 2017
The lattice-Boltzmann (LB) method is used in this work to reproduce the controlled addition of binder and hydrophobicity-promoting agents, like polytetrafluoroethylene (PTFE), into gas diffusion layers (GDLs) and to predict flow permeabilities in the through- and in-plane directions. The present simulator manages to reproduce spreading of binder and hydrophobic additives, sequentially, into the neat fibrous layer using a two-phase flow model. Gas flow simulation is achieved by the same code, sidestepping the need for a post-processing flow code and avoiding the usual input/output and data interface problems that arise in other techniques. Compression effects on flow anisotropy of the impregnated GDL are also studied. The permeability predictions for different compression levels and for different binder or PTFE loadings are found to compare well with experimental data for commercial GDL products and with computational fluid dynamics (CFD) predictions. Alternatively, the PTFE-impregnated structure is reproduced from Scanning Electron Microscopy (SEM) images using an independent, purely geometrical approach. A comparison of the two approaches is made regarding their adequacy to reproduce correctly the main structural features of the GDL and to predict anisotropic flow permeabilities at different volume fractions of binder and hydrophobic additives. © 2017 Elsevier B.V. All rights reserved.
Keywords: Gas diffusion layer Simulation Reconstruction Flow computation Teflon addition
1. Introduction The role of gas diffusion layers in the operation of proton exchange membrane fuel cells (PEMFCs) is to provide a multitude of supporting functions in order to improve the performance of the cell. More specifically, GDLs aim to offer mechanical support to the catalyst layer and the membrane, fast conducting pathways for the gases and for the water that is produced, homogeneous distribution of the gas entrance front, thermally and electrically conductive pathways [1]. In fact, it is known that increased mass transfer rates in the porous electrodes allow higher current densities and provide
* Corresponding author. E-mail address: [email protected] (V.N. Burganos). http://dx.doi.org/10.1016/j.jpowsour.2017.08.070 0378-7753/© 2017 Elsevier B.V. All rights reserved.
higher concentrations of reactants in the catalyst layer, thus reducing the catalyst quantities needed. On the other hand, given that the catalyst layer is usually made of small particles, GDL pores that face the catalyst layer must have suitable size for a functional interface to develop. To facilitate the analysis and guide GDL design, several theoretical models and simulators have been developed, which, in combination with experimental studies, aim to improve our understanding of the role of the internal structure of GDLs in achieving the aforementioned functions [2e9]. The models that address these issues can be classified in two major categories: (a) Macroscopic models that involve phenomenological transport equations. They can provide a quick parametric guidance to the electrode design from the perspective of its interface with the rest of the membrane electrode assembly. (b) Simulators that pay attention to the internal structure of the porous electrode and
180
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
attempt to offer quantitative relationships between structural details and transport properties. The latter offer the additional advantage that they can offer more detailed guidance to the design and tailored construction of GDLs with improved function during cell operation. The main structural parameters that are known to affect the GDL performance are the underlying carbon fiber configuration, the compression level as the result of the applied clamping pressure, the addition of binder to improve mechanical stability of the GDL, the addition of PTFE or other hydrophobicity-inducing agent, the eventual porosity, as well as the tortuosity of the GDL structure that affects both the mass transport coefficients and the gas feed uniformity into the porous electrodes. Several GDL structure models have appeared in the literature that refer to the reconstruction of the fiber system alone and not to the subsequent stages of binder and PTFE addition (for instance [10e12]). Such models assume typically a random placement of fibers in space, allowing or prohibiting mutual overlapping, of uniform diameter, and in random or in layered arrangements parallel to the macroscopic plane of the GDL. Analysis of SEM images is quite useful for the extraction of pore structure data that can be used in the stochastic reconstruction of GDL samples [13]. The number density of fibers and the fiber diameter can be selected in a way that allows reproduction of the GDL porosity in agreement with literature values or with supplier's data. Although useful in understanding the effect of neat fiber arrangement on the GDL transport properties, the role of binder and PTFE addition may be dominant in this direction, depending of course on the respective loading degree. In fact, visual inspection of SEM images of commercial GDLs that are treated with binder and, subsequently, undergo bulk treatment with PTFE reveals the development of extensive fiber coating as well as extensive web formation between adjacent fibers. Binding of fibers that belong to different layers is mainly achieved through local coating around the fiber-touching or fiber-crossing regions. Depending on the type of GDL production, PTFE may be added as a binder itself (e.g., wet-laid filled papers) or during post-processing as a bulk treatment of the GDL. The morphology of the PTFE addition in the interior of the GDL depends, in turn, on the processing route and conditions, for instance, dipping or spraying, PTFE particle sintering etc [14]. Usually, PTFE forms coatings and extended webs, the precise geometry and size of which depend on a range of factors, including the degree of loading, the fiber wettability (quantified typically by the contact angle with PTFE), processing type and conditions, the viscosity of the additive, etc. In any case, it is noted that web planes are almost parallel to the fiber layers and to the macroscopic GDL plane. This particular orientation of webs causes a dramatic reduction of the GDL permeability in the through-plane direction and is also expected to affect the in-plane permeability as well. This raises the main issue of the dependence of the well-known GDL anisotropy on the binder or PTFE loading and the concomitant effect on the gas and water distribution as well as the water percolation through the porous electrode [9]. In fact, it was found that although diffusion is the main transport mechanism across the GDL and towards the catalyst and membrane layers, convection may also be important, driven by the pressure drop that develops; in addition, pressure drop in the in-plane direction is also known to develop, owing to the generation of pressure differences between neighboring gas channels. Such a lateral flow is expected to affect gas flow configuration through the employment of a number of lateral transport pathways and, eventually, to delay the outset of mass transfer limitations. To address the aforementioned issues, different attempts have been made in the literature to include binder and/or PTFE presence in GDL structure modeling. Typically, these attempts involve the
placement of certain PTFE quantities in selected regions within the fibrous assembly [5], using also sequential spherical structuring elements [2] or the approximation of the void space by a network of pore bodies and pore throats [15] without paying particular attention to the local physicochemistry of the impregnation or to the conditions that must be met for the particular web geometries. Probabilistic models have also been developed to include the presence of binder within the fibrous assembly [8]. Other attempts include the addition of binder or PTFE in the form of fiber coating, as well as within the intersection regions that form at fiber crossings [16e18]. The employment of X-ray microtomography is another option to explore the internal structure of GDLs, offering serial sections of the GDL structure, with or without PTFE addition, that can be subsequently used for the 3D reconstruction of the material [3,19e21]. However, although quite direct, this approach is certainly far more expensive than a statistical or stochastic reconstruction approach that is based on few SEM images only. In addition, it cannot be adopted to run a parametric study of the effects of the structural parameters on the GDL performance regarding gas distribution and conduction that would facilitate the design of improved GDL and guide preparation processing details. In the present work, a new simulator is developed that generates GDL structures considering sequential spreading of binder and PTFE around and between touching fibers. The simulator is based on a two-phase lattice-Boltzmann model and, in addition to fiber coating, allows the formation of webs between fibers, thus reproducing in a natural and realistic way the corresponding web formation that is noted on SEM images. The extent of these webs is dependent upon the loadings of these additives, which are treated as input variables. The GDL fibers can be randomly distributed in three dimensions or be restricted within layers that are parallel to the macroscopic GDL surface. An additional advantage of the new simulator is that the generation of the GDL structure can be directly followed by the computation of the flow permeability using the same code in both the through- and in-plane directions. Freezing of the additives at their equilibration positions is assumed whereas the solid fiber volume of the GDL is expanded by the volume of the additives that were introduced in the sample. In this way, an integrated structure and flow simulator is developed to study structural effects on GDL permeability and to predict binder and PTFE effects on GDL flow anisotropy. The effect of GDL compression is also investigated here and the model predictions are compared to literature data. Finally, a comparison is made between the performance of the new structure simulator and that of an alternative structure model, also developed here, that regenerates PTFEimpregnated GDL samples in three dimensions using geometrical aspects only, drawing all needed information from SEM images. 2. Methodology of GDL reconstruction from SEM images e geometrical approach The computational reconstruction of a GDL proceeds in two stages: (a) reconstruction of the fibrous part of the layer; (b) application of binder and PTFE coatings around and between the fibers. The procedure allows the three-dimensional reconstruction of the GDL layer with sequential three-dimensional binder and PTFE coatings. That is, the methodology is free of the commonly employed approximation of single-layer (essentially planar) fiber arrangements, placed sequentially on top of each other. It also allows the addition of binder and PTFE in a fashion that is compatible with wettability requirements and respects fundamental capillarity conditions. Although computationally and conceptually more complicated, the methodology offers two significant advantages, namely, it offers a closer representation of the actual structure and
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
it allows for a closer study of compression effects on transport properties of the GDL layer thanks to the fully three-dimensional (as compared to planar) reconstruction process. A detailed description of the two reconstruction stages follows below. 2.1. Reconstruction of the fibrous part of GDL A working domain that is shaped as rectangular parallelepiped is used with size L1, L2, L3 in the Cartesian directions x, y, z, respectively (Fig. 1a). A prescribed number of randomly oriented fibers are placed at random positions within the working domain. The diameter of each fiber is considered constant but can be easily also sampled from a prescribed size distribution. Both structural properties, namely, the number density of fibers and the fiber diameter distribution can be extracted from SEM images following appropriate image analysis or can be obtained directly from the supplier. The traces of the fiber axes on the planes x ¼ 0, x ¼ L1, y ¼ 0, y ¼ L2 are distributed randomly over domains that ensure complete positioning of the elliptic traces of the fibers within the thickness L3 of the GDL. In the case that planar arrangement of fibers is desired, the traces can be simply confined to the same z-level. Alternatively, positioning of the traces on opposite x and y faces can be limited within a fixed interval of z positions. The fibers can be either allowed to overlap freely or prohibited from mutual overlapping. The latter is checked with the help of the distance, d (i, j), of the axes of fibers i and j, which is compared to the sum of their radii, ri þ rj. The Appendix lists the main working expressions for the distance of a point from a line in the 3D space and for the distance between non-crossing lines in 3D space. If d (i, j)< ri þ rj, the fiber that was placed last is placed at a new position with a new random orientation. The above check is made against all previously placed fibers and so on, until the desired number of fibers is inserted within the working domain. 2.2. Application of binder and PTFE coating Depending on the binder or PTFE addition process, such a coating may extend around some or all of the fibers and wet the angular areas that form between touching or crossing fibers. Coating around the fibers is straightforward to implement: a coating of desired thickness is added forming a coaxial cylindrical shell outside the fiber. The simulation of a web between crossing fibers in three dimensions, having a finite, prescribed thickness, will reproduce realistic representations of actual GDL structures but is far more cumbersome than fiber coating. The procedure that is suggested in this work is outlined below. First, a thorough check is made to register potential positions for binder or PTFE coating throughout the working domain. These positions fall in two categories: (a) fiber-to-fiber crossing positions, which is meaningful if the fibers are allowed to overlap; (b) fiber distance smaller than twice the prescribed thickness of the fiber coating. A certain fiber may cross more than one fiber along its length, within the working domain. Formation of a web at the aforementioned positions can be either imposed on each and every “crossing” position or decided randomly subject to a prescribed probability. A web can develop only in the interior of an acute angle; the interior of obtuse angles is assumed to be free of webs, in accord with visual observation of several microscopy images of actual GDL samples. In fact, discernible webbing of obtuse angles sometimes appears in microscopy images but is limited to a very thin coating near the position of fiber crossing. Construction of a web proceeds as follows. The overlapped region between two fibers is located around the shortest line
181
segment, AB, that is normal to both fiber axes. The length of this segment is defined as the distance of the fiber axes. A plane is defined by the midpoint, M, of segment AB and two points C and D that lie on the two fiber axes surrounding the web and define the extension of the particular web. In order to simulate the curvature of the web edges that results from the wetting procedure, a void sphere is generated that is tangential to the coatings of the two fibers at points C and D, with center at point N, as shown in Fig. 2. The volume of the web that is cut out by this sphere is considered as part of the void volume of the GDL. The distances of C and D from M are defined in a way that, eventually, the desired binder or PTFE content is reproduced. In general, MCsMD, but here it is assumed that MC ¼ MD to simplify the analysis. The thickness of the web, w, can be also obtained from microscopy images and can be assumed equal to the thickness of the coating shells around the fibers in order to reduce the number of coating parameters. The identification of voxels that lie within the web thickness proceeds as follows. The distance of the center of the voxel from the plane of the web is calculated first using the corresponding expression that is mentioned in the Appendix. The projection on the plane of the web is found and a check is made as to whether this point lies in the interior of the web triangle or is external to it. The criteria for this check are summarized in the Appendix. If the projection of the voxel center is internal to the web triangle, the voxel is considered as part of the web if the distance of its center from the web triangle is less than half the thickness of the web, w/2. Qualitatively, the particular reconstruction procedure is in line with the requirement set by the Young Laplace equation [22], namely, that the sum of the inverse values of the local radii of curvature is constant. On several occasions, depending on the preparation route and conditions, the coating around the fibers and/or the webs are not completely filled with the binding material but are made up of particulate matter with solid volume fraction 4i<1, where i indicates the phase of interest, that is, i ¼ c for the coating layer, i ¼ w for the web area. This is simulated as follows. The working domain is discretized into voxels of uniform size. If the center of a voxel is located within a distance w/2 from the plane MCD, the voxel is registered as part of the web, provided that R < 4i, where R is a randomly selected number in the interval (0,1). The same procedure can be applied to the construction of the coating shell around a fiber. A typical example of the reconstructed GDL is shown in Fig. 1b. The porosity of the sample is calculated as follows. A check is made whether the center of each and every voxel of the discretized sample lies within either the fibrous part of the sample or some coating or web region. The number fraction of these voxels defines the solid volume fraction 4T of the sample, that is, 4T ¼ 1-εT, where εT is the total porosity. The solid volume fraction of the uncoated fibers is defined as 40 ¼ 1-ε0, where ε0 is the porosity of the fibrous structure prior to application of coating or webs. Obviously, 4T ¼ 40 þ 4w and ε0 ¼ 4w þ εT. The PTFE content, u, that is usually mentioned in commercial GDL specifications is the mass fraction of PTFE and can be calculated from
u¼
ε0 εT
að1 ε0 Þ þ ε0 εT
(1)
where a is the ratio of the fiber density to the PTFE density. Wrapping of fibers by binder or PTFE can also be simulated in a fashion that reproduces the particulate nature of these additives as detected from SEM images. Specifically, a certain voxel may or may not be occupied by additive even if it belongs to the coating shell of predefined radius, as mentioned above. One may prescribe a probability of additive attachment on the fiber surface that defines the desired coverage of fibers (Fig. 1c). The corresponding
182
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
Fig. 1. (a) Simulated GDL structure. The fibrous part of the structure is shown. (b) Geometrically simulated GDL structure following PTFE coating and web formation. (c) Simulated GDL structure with webs developing between selected fibers. The particulate nature of webs and coatings is shown to simulate limited PTFE sintering.
technique along with the source of their values. 2.3. Simulation of GDL compression
Fig. 2. Definition of web triangle in three-dimensional representation of fiber crossing.
probability value can be set to reproduce the desired solid fraction of the shells and webs, following observations of SEM images. Compact PTFE films are obtained for the value of 1, which is used in the present simulations to keep the number of parameters to a minimum. Although several types of structural data may be requested for the maximum degree of reliability and accuracy in this reconstruction, in the present work it is suggested that the fiber coating and the webs share the same geometrical properties, namely, thickness and threshold probability for solid voxels. The size of the webs is more easily extracted from SEM images; one may use an average value for it or sample random sizes from the size distribution that is extracted from SEM images. Table 1 lists the structural parameters that are associated with the reconstruction
Depending on the PEM type and configuration, the GDL may undergo different levels of compression during assembling of the MEA. The compression of the GDL is known to have significant effects on critical GDL properties, like permeability, diffusivity, conductivity, water condensation, and water transport. In addition, the effects of the transport properties on the anisotropy, which is always present in GDLs, are also expected to be quite significant. In this work, compression is simulated in a simple fashion, namely, by translating the z-coordinates of each point in the reconstructed GDL by a prescribed compression factor. Although this procedure may slightly alter the degree of overlapping of certain fibers, it is expected that the concomitant changes in the GDL configuration will be negligible. Images of simulated uncompressed and compressed variants of the same GDL structure are shown in Fig. 3. 3. Reconstruction of GDL using the lattice-Boltzmann technique The two-phase lattice Boltzmann technique has been developed with the aim to simplify the simulation of flow of two phases under either displacement or steady state flow conditions, and in addition, to model the equilibrium configuration of immiscible phases both in the bulk and also in the vicinity of solid boundaries. Wetting phenomena and effects on the behavior of two-phase systems can also be studied taking into account the degree of surface wettability
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
183
Table 1 Structural parameters in GDL reconstruction. Parameter
Source
Base case
Porosity of neat fiber layer Porosity after addition of binder/PTFE PTFE content Thickness of fiber coating Thickness of webs Solid fraction of fiber coating Solid fraction of webs Web length Number fraction of coated fibers Number of fraction of webbed fibers Fiber size Number density of fibers Thickness of GDL
supplier, porosimetry supplier, porosimetry Supplier SEM SEM SEM, extracted SEM, extracted SEM SEM SEM SEM, supplier SEM, supplier supplier
0.84 0.80 20% 4 mm 4 mm 0.1 0.7 100 mm 0.5 0.5 8 mm (mean) 4.7 105mm3 190 mm
Fig. 3. (a) Uncompressed GDL layer, reconstructed using the method presented in this work. Porosity: 0.86. (b) Compressed GDL layer in the z-direction. Porosity: 0.72. Rate of compression 2:1. Only the fibrous part of the GDL sample is shown.
as quantified through the contact angle. In this work, the free energy-based two-phase lattice Boltzmann model [23] is employed, as it was incorporated in the lattice Boltzmann simulator [24,25]. The mathematical treatment that was developed in this simulator allowed restoration of the Galilean invariance and the realistic simulation of moving droplets and interfaces in two and three dimensions. Briefly, the evolution of particle distribution functions, discrete in both time and space that obey the lattice-Boltzmann equation is monitored. The particle populations move along unit vectors
connecting neighboring nodes. For the present calculations, the three-dimensional Faced Centered Cubic (FCC) lattice is used with 14 different directions of particle propagation at each time step (D3Q15 lattice, Fig. 4a). Using droplet simulation and calculations, it has been shown [25] that the simulator is consistent with the Laplace law and Gibbs-Thomson equations for two-phase systems. The wettability of solid walls can be adjusted by prescribing a fluid density profile along the interfacial sites. The simulator accounts also for phase transition, which allows the simulation of evaporation or condensation phenomena. The reconstruction of binder or PTFE coatings and webs between the GDL fibers proceeds as follows. Once the fiber system is constructed, the desired amount of additive (binder or hydrophobicity promoter) is set and the contact angle of the additive on the fiber surfaces is fixed. Next, a fraction of sites that corresponds to the additive volume fraction is assigned a fixed “liquid” density value whereas the rest are assigned a “gas” density value. The terms liquid and gas here are coined arbitrarily and are understood as distinct, immiscible phase density values. The site densities are then allowed to redistribute themselves subject to the latticeBoltzmann equations, which eventually lead to phase separation thanks to the action of surface tension: for small contact angles, the additive adheres to the fibers and creates connecting areas (webs) at fiber intersections, the extent of which depends on the additive fraction and the contact angle. In fact, the same method can be extended to simulate several stages of GDL processing, including addition of binder to stabilize the layer, and subsequent addition of PTFE to promote further hydrophobicity. The simulation of sequential processing stages involves successive addition of the material of interest and intermediate solidification. More specifically, binder material can be added to the original fiber assembly using the lattice-Boltzmann technique that was presented above. Once the binder phase has organized itself around the fiber connections, the binding “liquid” phase of the two-phase lattice-Boltzmann simulator is assumed to solidify at that position, ready to admit the addition of PTFE in a similar way, using again the two-phase lattice-Boltzmann simulator. The PTFE amount that was introduced during this stage is then allowed to solidify after completion of PTFE-organization around and between the fibers and binder (Fig. 4). The procedure can be used for as many layers of additive as desired, depending on the application. Clearly, one can define different wettability conditions at each step, for instance through the definition of different contact angle, and different volume fraction of the additive. The main advantage of the present simulator is that all these stages can be simulated in a sequential yet automated fashion, reproducing GDL fiber coating and additive spreading between the fibers in a spontaneous manner that respects wettability requirements. The
184
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
Fig. 4. Reconstruction of GDL and hydrophobic additives using the lattice-Boltzmann method. (a) Lattice-Boltzmann unit cell; (b) untreated fiber system; (c) simulated addition of binder using the LB technique; (d) simulated subsequent addition of PTFE using the LB technique. Thickness of binder layer: 1 mm. PTFE content: 30%.
method allows incorporation of practically every structural parameter that appears on Table 1. 4. Flow permeability calculation The reconstructed samples of GDL can be directly used for the numerical calculation of transport properties of gas species along the principal directions of the layer (through-plane and in-plane properties). Based on visual observations of SEM images of GDL samples but also on available experimental data, anisotropy in GDL samples is quite pronounced and clearly increases upon impregnation of the samples with binder and PTFE, ought to the formation of extended webs. The latter develop in direction parallel to the macroscopic surface of the gas diffusion layer and are expected to reduce more drastically the transport properties in the throughplane direction than in the in-plane direction. Two different numerical techniques are used in this work for the calculation of the flow permeability, namely, the single-phase flow lattice-Boltzmann technique and the Finite Element Method (FEM). Both techniques can be applied to the purely geometrical
reconstruction domains as well as to the ones that are produced by the two-phase lattice-Boltzmann method. The advantage of employing the lattice-Boltzmann option in the latter case is that it is free of input/output (I/O) issues as it is applied directly to the LB reconstructed domain. In this way, an integrated lattice-Boltzmann technique is formulated that uses a single code to handle both the reconstruction stage and the permeability computation stage. A brief description of the LB flow simulator is given next. The simulator uses the three-dimensional, 15 speed (D3Q15) unit cell model and involves the evolution of population distribution functions in space and time. The evolution of the system follows the LB equation [26].
fi ðr þ ei ; t þ 1Þ ¼ fi ðr; tÞ
fi fieq
t
;
i ¼ 0; …; 14
(2)
where r is the position vector, ei is the unit velocity vector in direction i, fi is the corresponding population distribution function, t is the single relaxation parameter according to the Bhatnagar, Gross, and Krook (BGK) approach [27] and feq i is the equilibrium
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
185
distribution function at position i, given by the following expansion in terms of the local velocity, y [25,28]:
1 1 1 fieq ¼ wi r 1 þ 2 ðei ,yÞ þ 4 ðei ,yÞ2 2 y2 ; cs 2cs 2cs w0 ¼
2 9
;
1 wi ¼ ; i ¼ 1; …; 6 9
;
wi ¼
1 ; i ¼ 7; …; 14 72
(3)
where r is the fluid density. The speed of sound in this system is
1 cs ¼ pffiffiffi 3
(4)
and the corresponding kinematic viscosity for the D3Q15 model is:
n0 ¼
2t 1 : 6
(5)
An orthogonal computational domain is used for the present study. Flow is induced by prescribing a pressure difference between the inlet and outlet faces of the working domain, whereas the rest of the faces are treated as impermeable boundaries (walls) or as symmetric boundaries. The latter choice is adopted for the sake of comparison with permeability results from other computational techniques, like the FEM method that is mentioned below. The no slip boundary condition is applied at wall sites and is implemented by assuming particle bounce-back at a position that is located halfway between the solid sites and the adjacent fluid sites. An open source LB-BGK software using the D3Q19 unit cell [http://palabos. org] has also been used in the literature [29,30] to calculate flow permeability in GDL domains, adopting the probabilistic geometric approach of [8] to fill the fibrous structure with the desired binder quantity. To simplify the computation, the single relaxation time option is used, which is found to converge sufficiently rapidly for domain sizes of 4 106 sites. In order to avoid problems with flow calculations that depend on the value of the relaxation time, a parametric study was carried out using as varying quantity the kinematic viscosity or, equivalently, the relaxation time. It was found that the flow permeability result remained nearly constant provided that the viscosity was varied within a sufficiently wide interval, namely, between 0.01 and 0.1. The second numerical technique that is used here makes use of commercial software COMSOL and allows various levels of meshing of the reconstructed domain. An appropriate software was required to allow interfacing of the flow solver with the structure regeneration process, including the main fiber assembly and the PTFE coating and web structures.
Fig. 5. Simulated variation of mean pore diameter and specific area with PTFE content.
also plotted in this figure, as extracted from the expression ¼ 4εT/Sv. Both quantities are rendered dimensionless using as reference values the values that correspond to the neat fiber domain, that is, prior to impregnation with PTFE. As expected, the
5. Results and discussion Input data for the simulation studies were taken from the literature and from SEM image analysis. Typical input values are listed in Table 1 for a base-case study. Both of the algorithms that were presented above for the reconstruction of the PTFE coating and web formation, namely, the purely geometrical approach and the two-phase lattice-Boltzmann model, were used to generate GDL samples in three dimensions. In the latter case, strong wettability conditions were implemented to reproduce a small contact angle, equal to 10 . Subsequently, the flow permeabilities were calculated using both the single-phase lattice-Boltzmann approach and the FEM method. In Fig. 5, the variation of the specific surface area, Sv, with the PTFE content is shown, using the lattice-Boltzmann method for the impregnation of GDL with PTFE. The mean pore diameter is
Fig. 6. Typical flow field in the interior of a reconstructed GDL sample. Results are shown on the same cross section of the sample, perpendicular to the GDL orientation, before (a) and after (b) the addition of hydrophobic agent (u ¼ 0.3). Block arrows indicate main flow direction (through-plane).
186
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
internal surface area clearly increases with the addition of PTFE. Flow and permeability calculations are presented next. Typical flow fields that are obtained before and after the addition of hydrophobic agents are shown in Fig. 6. As expected, addition of PTFE (u ¼ 0.3)leads to narrower openings between fibers, which increase flow resistance in the through-plane direction and reduce the velocity component in this direction, diverting flow into the in-plane direction, as will be discussed below. Fig. 7a shows the reduction of the through-plane permeability upon compression. The case of a neat fiber structure is considered here, which is free of binder or PTFE presence. At the uncompressed stage, the fibers are placed within the working domain exercising care to avoid overlapping, as explained in the reconstruction section above. Evidently, some limited degree of overlapping will develop during compression as the result of the simplified thickness reduction procedure that was adopted in this work. The compression level is expressed in this figure through the reduction of the porosity of the original, uncompressed structure and, alternatively, through the reduction of the thickness of the GDL. The relationship between the GDL thickness and porosity is given by
Lz 1 ε0 ¼ Lzi 1ε
(6)
where Lz and ε indicate the thickness and the porosity of the GDL sample, respectively, whereas superscript 0 is used to indicate values at the original, uncompressed stage. The simulation results were also compared with experimental data for Toray paper [7] using as reference value the permeability of the uncompressed structure. The agreement is excellent and indicates very satisfactory representation of the GDL for every compressed level employed here. The dimensional permeability of the uncompressed sample was calculated to 8.5 1012 m2, which also compares very well with the experimental value of 8.7 1012 m2 in Ref. [7] and 8.99 1012 m2 in Ref. [4]. This reflects the success in the reconstruction of the uncompressed GDL sample. It is also important to note that the two numerical solvers of the flow equations, namely, the LB and the FEM solvers, yield permeability predictions that are in very good agreement with each other, despite the fact that the LB method uses the centers of the reconstruction voxels as sites for the density function propagation, and, hence, for the flow “meshing” part, without further refinement. Based on this observation, and in order to save on computational cost and effort, flow calculations are carried out using the LB technique unless otherwise noted. The pressure drop between neighboring flow channels gives rise to flow in direction normal to the main transport direction, which is known to delay onset of mass transfer limitations during cell operation. To elucidate structural effects on the in-plane flow phenomena, flow computations were carried out in this direction as well, leading to the estimation of the in-plane permeability. The results are shown in Fig. 7b and compared to experimental data from Ref. [4]. Again, very good agreement is noticed over the entire compression range (0e30%) with respect to the uncompressed GDL sample. The simulation results in the two principal in-plane directions, x- and y-, were averaged in order to produce orientationally averaged results that can be compared to experiments. The degree of anisotropy of the GDL structure is then expressed by the ratio of inplane to through-plane permeabilities and is found to range between 1.7 and 2.4 depending on the degree of compression (Fig. 7c) using the LB technique. This is of strong interest in the design and performance of the electrode assembly: increased transport in the direction normal to the main transport direction is known to affect gas distribution at the entrance of the catalyst layer and the
Fig. 7. (a) Dependence of reduced permeability on compression level expressed in the form of reduced porosity (bottom axis) and thickness reduction (top axis) for a reconstructed GDL sample, free of binder or PTFE addition. Comparison of the predictions of the two numerical techniques that are employed here of (a) through-plane values with experimental ones from Ref. [7], (b) in-plane values with experimental ones from Ref. [4], and (c) compression effect on the GDL anisotropy, expressed as ratio of in-plane to through-plane permeabilities.
concomitant demands for catalyst performance. Results on the effect of GDL impregnation with binder and
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
187
Fig. 8. (a) Effect of PTFE loading on through-plane permeability for a GDL reconstruction. PTFE addition was implemented through the LB technique (filled symbols) and, alternatively, through the geometry-based technique (open symbols). Comparison of numerical predictions for a GDL reconstruction with experimental values from Ref. [7]. (b) Prediction of the effect of PTFE loading on the GDL anisotropy, expressed as the ratio of in-plane to through-plane permeability. PTFE addition was implemented through the LB technique.
hydrophobicity-enhancing materials, like PTFE, on the flow permeability are presented next. Both techniques that were developed in the present work are employed, namely, the latticeBoltzmann technique using the two-phase simulator, and the reconstruction from SEM images using a purely geometrical
approach. In the lattice-Boltzmann case, the volume fraction of the impregnating additive was varied keeping all other properties constant, that is, surface tension, viscosity, and density. For each value of the volume fraction, the two-phase lattice-Boltzmann simulator was employed to allow spreading of the liquid around
188
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
and between the GDL fibers, as described above. Once equilibrated to a configuration that is no longer propagating, the liquid part is assumed to “solidify” and the flow permeability is calculated using the lattice-Boltzmann simulator for single phase flow. The results for the through-plane permeability are shown in Fig. 8a and compared to experimental data from Ref. [7]. The agreement to the experimental values is quite satisfactory if one considers the stochastic nature of the GDL reconstruction and the fact that the impregnation process was simulated using a two-phase simulator neglecting dynamic invasion phenomena. The computation process can be repeated for the in-plane permeability and its ratio to the corresponding through-plane permeability is plotted against the PTFE content in Fig. 8b. It is noteworthy that anisotropy is enhanced by PTFE addition, due to the unequal effect of PTFE addition on the permeability in the various directions of the GDL sample. Specifically, the flow permeability is more drastically reduced upon increase of PTFE loading in the through-plane direction than in the in-plane direction. In the case of PTFE spreading reproduction by the geometrical approach, SEM images are directly used in the way that was described above. The corresponding flow permeability results are also shown in Fig. 8a in the form of open symbols, for the sake of comparison with the lattice-Boltzmann reconstruction procedure. Although some deviation between the two sets of permeability predictions is apparent, the general trend in the behavior of the permeability is the same. This deviation is obviously attributed to the different approach that was used to simulate PTFE spreading and to the simplifying assumptions that were made in the geometrical reconstruction case regarding the shape and thickness of the webs. 6. Conclusions An integrated simulator is developed in this work that is capable of reproducing in three dimensions the structure of gas diffusion layers and simulating flow and diffusion in the reconstructed samples. The simulator is based on the LB technique, which is appropriately adjusted to simulate addition of binder to the carbon fiber system and the subsequent addition of PTFE or any other additive to promote further the hydrophobic characteristics of the layer. Among the advantages of the present simulator compared to other approaches is the fact that fiber coating by binder and further coating by PTFE accompanied by web formation between adjacent fibers can be handled in a physicochemically consistent manner, as it takes into account wettability conditions that shape the web configuration across the structure. A second advantage is that all stages of the process, that is, reconstruction in three dimensions, simulation of flow, and calculation of permeabilities, can be tackled by the same simulator without having to resort to different numerical solvers for the different stages of simulation. An alternative reconstruction procedure was also developed in this work that is based merely on geometrical features, as these are extracted from SEM images of GDL samples. Comparison with experimental data for the GDL permeability reveals very satisfactory agreement in both main directions, that is, for the through- and the in-plane permeability. Compression effects are also simulated and the results are shown to be in quantitative agreement with experiments. Addition of PTFE leads to significant reduction of permeability in both main directions. The degree of flow anisotropy is also estimated and is found to be dependent on the compression level and on the PTFE loading. This anisotropy is known to have a significant role in the design and performance of the fuel cell as it affects the gas distribution in both anode and cathode sides, which in turn influences the onset of mass transfer limitations during cell operation. The simulator that is developed in this work can be further
utilized for a comprehensive simulation of the phenomena that take place inside a GDL, namely, gas transport, water condensation, and water transport in both the anode and the cathode of the cell. Breakthrough conditions can also be identified using as parameters the water content, the hydrophobicity of the layer, and the applied pressure drop. In fact, these conditions can be investigated using the same LB simulator given that it inherently accounts for phase transition and for dynamic two-phase flow, provided of course certain appropriate adjustments and fine-tuning are carried out. This work is already in progress by the authors. Acknowledgements Financial support to this work was provided by the Institute of Chemical Engineering Sciences of the Foundation for Research and Technology, Hellas, and by the European Union's Seventh Framework Programme (FP7/2007e2013) under FCH-JU Project Grant number 303024, EURECA. APPENDIX
A. Distance of point re0 from line r ¼ rei þ t a: e e
dðr ;r Þ¼ðr r Þa
e0 e
e0 e1 e a e
B. Distance of two lines in 3D Cartesian space. yy1 yy2 zz1 xx2 zz2 1 The distance of lines xx [1 ¼ m1 ¼ n1 and [2 ¼ m2 ¼ n2 is given by
x2 x1 y2 y1 z2 z1 [1 m1 n1 [2 m2 n2 jð[1 ; m1 ; n1 Þxð[2 ; m2 ; n2 Þj C. Conditions for point (x,y,z) to lie in the interior of the triangle that is defined by the corners (xi, yi, zi), i ¼ 1,2,3: (a) a, b > 0 (b) aþb < 1where
x x1 x3 x1 a ¼ x2 x1 x3 x1
x x1 y y1 x x1 y3 y1 ; b ¼ 2 y2 y1 x2 x1 x3 x1 y3 y1
y y1 y2 y1 y2 y1 y3 y1
References [1] F. Barbir, Sol. Energy 78 (2005) 661e669. [2] M.M. Daino, S.G. Kandlikar, Int. J. Hydrogen Energy 37 (2012) 5180e5189. [3] Y. Gao, X.X. Zhang, P. Rama, Y. Liu, R. Chen, H. Ostadi, K. Jiang, Fuel Cells 12 (2012) 365e381. [4] J.T. Gostick, M.W. Fowler, M.D. Pritzker, M.A. Ioannidis, L.M. Behra, J. Power Sources 162 (2006) 228e238. [5] L. Hao, P. Cheng, J. Power Sources 186 (2009) 104e114. [6] J.H. Nam, M. Kaviany, Int. J. Heat. Mass Transf. 46 (2003) 4595e4611. [7] A. Tamayol, F. McGregor, M. Bahrami, J. Power Sources 204 (2012) 94e99. [8] R. Thiedmann, F. Fleischer, C. Hartnig, W. Lehnert, V. Schmidt, J. Electrochem. Soc. 155 (2008) B391eB399. [9] N. Zamel, X.G. Li, Prog. Energy Combust. Sci. 39 (2013) 111e146. [10] G. Inoue, T. Yoshimoto, Y. Matsukuma, M. Minemoto, J. Power Sources 175 (2008) 145e158. [11] X.B. Nie, G.D. Doolen, S.Y. Chen, J. Stat. Phys. 107 (2002) 279e289. [12] V.P. Schulz, J. Becker, A. Wiegmann, P.P. Mukherjee, C.Y. Wang, J. Electrochem. Soc. 154 (2007) B419eB426. [13] N. Parikh, J.S. Allen, R.S. Yassar, Fuel Cells 12 (2012) 382e390. [14] M.F. Mathias, J. Roth, J. Fleming, W. Lehnert, Diffusion media materials and characterisation, in: Handbook of Fuel Cells, John Wiley & Sons, Ltd, 2010. [15] J.T. Gostick, M.A. Ioannidis, M.W. Fowler, M.D. Pritzker, J. Power Sources 173 (2007) 277e290.
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189 [16] A. Nabovati, J. Hinebaugh, A. Bazylak, C.H. Amon, J. Power Sources 248 (2014) 83e90. [17] H. Sadeghifar, N. Djilali, M. Bahrami, J. Power Sources 248 (2014) 632e641. [18] J. Yablecki, A. Bazylak, J. Power Sources 217 (2012) 470e478. [19] H. Ostadi, P. Rama, Y. Liu, R. Chen, X.X. Zhang, K. Jiang, J. Membr. Sci. 351 (2010) 69e74. [20] A. Pfrang, S. Didas, G. Tsotridis, J. Power Sources 235 (2013) 81e86. [21] P.K. Sinha, P. Halleck, C.Y. Wang, Electrochem. Solid State Lett. 9 (2006) A344eA348. [22] F.A.L. Dullien, Porous Media: Fluid Transport and Pore Structure, first ed., Academic Press, New York, 1979.
[23] [24] [25] [26] [27] [28] [29]
189
M.R. Swift, W.R. Osborn, J.M. Yeomans, Phys. Rev. Lett. 75 (1995) 830e833. A.N. Kalarakis, V.N. Burganos, A.C. Payatakes, Phys. Rev. E 65 (2002). A.N. Kalarakis, V.N. Burganos, A.C. Payatakes, Phys. Rev. E 67 (2003). Y.H. Qian, D. Dhumieres, P. Lallemand, Europhys. Lett. 17 (1992) 479e484. P.L. Bhatnagar, E.P. Gross, M. Krook, Phys. Rev. 94 (1954) 511e525. S. Chen, Z. Wang, X. Shan, G.D. Doolen, J. Stat. Phys. 68 (1992) 379e400. D. Froning, J. Brinkmann, U. Reimer, V. Schmidt, W. Lehnert, D. Stolten, Electrochim. Acta 110 (2013) 325e334. [30] D. Froning, G. Gaiselmann, U. Reimer, J. Brinkmann, V. Schmidt, W. Lehnert, Transp. Porous Media 103 (2014) 469e495.
Contents lists available at ScienceDirect
Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour
An integrated simulator of structure and anisotropic flow in gas diffusion layers with hydrophobic additives Vasilis N. Burganos a, *, Eugene D. Skouras a, b, Alexandros N. Kalarakis a, b a b
Institute of Chemical Engineering Sciences, Foundation for Research and Technology e Hellas, Patras 26504, Greece Department of Mechanical Engineering, Technological-Educational Institute of Western Greece, Patras 26334, Greece
h i g h l i g h t s Three dimensional reconstruction of gas diffusion layers from SEM is achieved. Addition of binder and hydrophobic agents is simulated using a two-phase flow model. Flow anisotropy index is predicted for different PTFE loadings and compression levels.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 14 March 2017 Received in revised form 13 June 2017 Accepted 18 August 2017
The lattice-Boltzmann (LB) method is used in this work to reproduce the controlled addition of binder and hydrophobicity-promoting agents, like polytetrafluoroethylene (PTFE), into gas diffusion layers (GDLs) and to predict flow permeabilities in the through- and in-plane directions. The present simulator manages to reproduce spreading of binder and hydrophobic additives, sequentially, into the neat fibrous layer using a two-phase flow model. Gas flow simulation is achieved by the same code, sidestepping the need for a post-processing flow code and avoiding the usual input/output and data interface problems that arise in other techniques. Compression effects on flow anisotropy of the impregnated GDL are also studied. The permeability predictions for different compression levels and for different binder or PTFE loadings are found to compare well with experimental data for commercial GDL products and with computational fluid dynamics (CFD) predictions. Alternatively, the PTFE-impregnated structure is reproduced from Scanning Electron Microscopy (SEM) images using an independent, purely geometrical approach. A comparison of the two approaches is made regarding their adequacy to reproduce correctly the main structural features of the GDL and to predict anisotropic flow permeabilities at different volume fractions of binder and hydrophobic additives. © 2017 Elsevier B.V. All rights reserved.
Keywords: Gas diffusion layer Simulation Reconstruction Flow computation Teflon addition
1. Introduction The role of gas diffusion layers in the operation of proton exchange membrane fuel cells (PEMFCs) is to provide a multitude of supporting functions in order to improve the performance of the cell. More specifically, GDLs aim to offer mechanical support to the catalyst layer and the membrane, fast conducting pathways for the gases and for the water that is produced, homogeneous distribution of the gas entrance front, thermally and electrically conductive pathways [1]. In fact, it is known that increased mass transfer rates in the porous electrodes allow higher current densities and provide
* Corresponding author. E-mail address: [email protected] (V.N. Burganos). http://dx.doi.org/10.1016/j.jpowsour.2017.08.070 0378-7753/© 2017 Elsevier B.V. All rights reserved.
higher concentrations of reactants in the catalyst layer, thus reducing the catalyst quantities needed. On the other hand, given that the catalyst layer is usually made of small particles, GDL pores that face the catalyst layer must have suitable size for a functional interface to develop. To facilitate the analysis and guide GDL design, several theoretical models and simulators have been developed, which, in combination with experimental studies, aim to improve our understanding of the role of the internal structure of GDLs in achieving the aforementioned functions [2e9]. The models that address these issues can be classified in two major categories: (a) Macroscopic models that involve phenomenological transport equations. They can provide a quick parametric guidance to the electrode design from the perspective of its interface with the rest of the membrane electrode assembly. (b) Simulators that pay attention to the internal structure of the porous electrode and
180
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
attempt to offer quantitative relationships between structural details and transport properties. The latter offer the additional advantage that they can offer more detailed guidance to the design and tailored construction of GDLs with improved function during cell operation. The main structural parameters that are known to affect the GDL performance are the underlying carbon fiber configuration, the compression level as the result of the applied clamping pressure, the addition of binder to improve mechanical stability of the GDL, the addition of PTFE or other hydrophobicity-inducing agent, the eventual porosity, as well as the tortuosity of the GDL structure that affects both the mass transport coefficients and the gas feed uniformity into the porous electrodes. Several GDL structure models have appeared in the literature that refer to the reconstruction of the fiber system alone and not to the subsequent stages of binder and PTFE addition (for instance [10e12]). Such models assume typically a random placement of fibers in space, allowing or prohibiting mutual overlapping, of uniform diameter, and in random or in layered arrangements parallel to the macroscopic plane of the GDL. Analysis of SEM images is quite useful for the extraction of pore structure data that can be used in the stochastic reconstruction of GDL samples [13]. The number density of fibers and the fiber diameter can be selected in a way that allows reproduction of the GDL porosity in agreement with literature values or with supplier's data. Although useful in understanding the effect of neat fiber arrangement on the GDL transport properties, the role of binder and PTFE addition may be dominant in this direction, depending of course on the respective loading degree. In fact, visual inspection of SEM images of commercial GDLs that are treated with binder and, subsequently, undergo bulk treatment with PTFE reveals the development of extensive fiber coating as well as extensive web formation between adjacent fibers. Binding of fibers that belong to different layers is mainly achieved through local coating around the fiber-touching or fiber-crossing regions. Depending on the type of GDL production, PTFE may be added as a binder itself (e.g., wet-laid filled papers) or during post-processing as a bulk treatment of the GDL. The morphology of the PTFE addition in the interior of the GDL depends, in turn, on the processing route and conditions, for instance, dipping or spraying, PTFE particle sintering etc [14]. Usually, PTFE forms coatings and extended webs, the precise geometry and size of which depend on a range of factors, including the degree of loading, the fiber wettability (quantified typically by the contact angle with PTFE), processing type and conditions, the viscosity of the additive, etc. In any case, it is noted that web planes are almost parallel to the fiber layers and to the macroscopic GDL plane. This particular orientation of webs causes a dramatic reduction of the GDL permeability in the through-plane direction and is also expected to affect the in-plane permeability as well. This raises the main issue of the dependence of the well-known GDL anisotropy on the binder or PTFE loading and the concomitant effect on the gas and water distribution as well as the water percolation through the porous electrode [9]. In fact, it was found that although diffusion is the main transport mechanism across the GDL and towards the catalyst and membrane layers, convection may also be important, driven by the pressure drop that develops; in addition, pressure drop in the in-plane direction is also known to develop, owing to the generation of pressure differences between neighboring gas channels. Such a lateral flow is expected to affect gas flow configuration through the employment of a number of lateral transport pathways and, eventually, to delay the outset of mass transfer limitations. To address the aforementioned issues, different attempts have been made in the literature to include binder and/or PTFE presence in GDL structure modeling. Typically, these attempts involve the
placement of certain PTFE quantities in selected regions within the fibrous assembly [5], using also sequential spherical structuring elements [2] or the approximation of the void space by a network of pore bodies and pore throats [15] without paying particular attention to the local physicochemistry of the impregnation or to the conditions that must be met for the particular web geometries. Probabilistic models have also been developed to include the presence of binder within the fibrous assembly [8]. Other attempts include the addition of binder or PTFE in the form of fiber coating, as well as within the intersection regions that form at fiber crossings [16e18]. The employment of X-ray microtomography is another option to explore the internal structure of GDLs, offering serial sections of the GDL structure, with or without PTFE addition, that can be subsequently used for the 3D reconstruction of the material [3,19e21]. However, although quite direct, this approach is certainly far more expensive than a statistical or stochastic reconstruction approach that is based on few SEM images only. In addition, it cannot be adopted to run a parametric study of the effects of the structural parameters on the GDL performance regarding gas distribution and conduction that would facilitate the design of improved GDL and guide preparation processing details. In the present work, a new simulator is developed that generates GDL structures considering sequential spreading of binder and PTFE around and between touching fibers. The simulator is based on a two-phase lattice-Boltzmann model and, in addition to fiber coating, allows the formation of webs between fibers, thus reproducing in a natural and realistic way the corresponding web formation that is noted on SEM images. The extent of these webs is dependent upon the loadings of these additives, which are treated as input variables. The GDL fibers can be randomly distributed in three dimensions or be restricted within layers that are parallel to the macroscopic GDL surface. An additional advantage of the new simulator is that the generation of the GDL structure can be directly followed by the computation of the flow permeability using the same code in both the through- and in-plane directions. Freezing of the additives at their equilibration positions is assumed whereas the solid fiber volume of the GDL is expanded by the volume of the additives that were introduced in the sample. In this way, an integrated structure and flow simulator is developed to study structural effects on GDL permeability and to predict binder and PTFE effects on GDL flow anisotropy. The effect of GDL compression is also investigated here and the model predictions are compared to literature data. Finally, a comparison is made between the performance of the new structure simulator and that of an alternative structure model, also developed here, that regenerates PTFEimpregnated GDL samples in three dimensions using geometrical aspects only, drawing all needed information from SEM images. 2. Methodology of GDL reconstruction from SEM images e geometrical approach The computational reconstruction of a GDL proceeds in two stages: (a) reconstruction of the fibrous part of the layer; (b) application of binder and PTFE coatings around and between the fibers. The procedure allows the three-dimensional reconstruction of the GDL layer with sequential three-dimensional binder and PTFE coatings. That is, the methodology is free of the commonly employed approximation of single-layer (essentially planar) fiber arrangements, placed sequentially on top of each other. It also allows the addition of binder and PTFE in a fashion that is compatible with wettability requirements and respects fundamental capillarity conditions. Although computationally and conceptually more complicated, the methodology offers two significant advantages, namely, it offers a closer representation of the actual structure and
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
it allows for a closer study of compression effects on transport properties of the GDL layer thanks to the fully three-dimensional (as compared to planar) reconstruction process. A detailed description of the two reconstruction stages follows below. 2.1. Reconstruction of the fibrous part of GDL A working domain that is shaped as rectangular parallelepiped is used with size L1, L2, L3 in the Cartesian directions x, y, z, respectively (Fig. 1a). A prescribed number of randomly oriented fibers are placed at random positions within the working domain. The diameter of each fiber is considered constant but can be easily also sampled from a prescribed size distribution. Both structural properties, namely, the number density of fibers and the fiber diameter distribution can be extracted from SEM images following appropriate image analysis or can be obtained directly from the supplier. The traces of the fiber axes on the planes x ¼ 0, x ¼ L1, y ¼ 0, y ¼ L2 are distributed randomly over domains that ensure complete positioning of the elliptic traces of the fibers within the thickness L3 of the GDL. In the case that planar arrangement of fibers is desired, the traces can be simply confined to the same z-level. Alternatively, positioning of the traces on opposite x and y faces can be limited within a fixed interval of z positions. The fibers can be either allowed to overlap freely or prohibited from mutual overlapping. The latter is checked with the help of the distance, d (i, j), of the axes of fibers i and j, which is compared to the sum of their radii, ri þ rj. The Appendix lists the main working expressions for the distance of a point from a line in the 3D space and for the distance between non-crossing lines in 3D space. If d (i, j)< ri þ rj, the fiber that was placed last is placed at a new position with a new random orientation. The above check is made against all previously placed fibers and so on, until the desired number of fibers is inserted within the working domain. 2.2. Application of binder and PTFE coating Depending on the binder or PTFE addition process, such a coating may extend around some or all of the fibers and wet the angular areas that form between touching or crossing fibers. Coating around the fibers is straightforward to implement: a coating of desired thickness is added forming a coaxial cylindrical shell outside the fiber. The simulation of a web between crossing fibers in three dimensions, having a finite, prescribed thickness, will reproduce realistic representations of actual GDL structures but is far more cumbersome than fiber coating. The procedure that is suggested in this work is outlined below. First, a thorough check is made to register potential positions for binder or PTFE coating throughout the working domain. These positions fall in two categories: (a) fiber-to-fiber crossing positions, which is meaningful if the fibers are allowed to overlap; (b) fiber distance smaller than twice the prescribed thickness of the fiber coating. A certain fiber may cross more than one fiber along its length, within the working domain. Formation of a web at the aforementioned positions can be either imposed on each and every “crossing” position or decided randomly subject to a prescribed probability. A web can develop only in the interior of an acute angle; the interior of obtuse angles is assumed to be free of webs, in accord with visual observation of several microscopy images of actual GDL samples. In fact, discernible webbing of obtuse angles sometimes appears in microscopy images but is limited to a very thin coating near the position of fiber crossing. Construction of a web proceeds as follows. The overlapped region between two fibers is located around the shortest line
181
segment, AB, that is normal to both fiber axes. The length of this segment is defined as the distance of the fiber axes. A plane is defined by the midpoint, M, of segment AB and two points C and D that lie on the two fiber axes surrounding the web and define the extension of the particular web. In order to simulate the curvature of the web edges that results from the wetting procedure, a void sphere is generated that is tangential to the coatings of the two fibers at points C and D, with center at point N, as shown in Fig. 2. The volume of the web that is cut out by this sphere is considered as part of the void volume of the GDL. The distances of C and D from M are defined in a way that, eventually, the desired binder or PTFE content is reproduced. In general, MCsMD, but here it is assumed that MC ¼ MD to simplify the analysis. The thickness of the web, w, can be also obtained from microscopy images and can be assumed equal to the thickness of the coating shells around the fibers in order to reduce the number of coating parameters. The identification of voxels that lie within the web thickness proceeds as follows. The distance of the center of the voxel from the plane of the web is calculated first using the corresponding expression that is mentioned in the Appendix. The projection on the plane of the web is found and a check is made as to whether this point lies in the interior of the web triangle or is external to it. The criteria for this check are summarized in the Appendix. If the projection of the voxel center is internal to the web triangle, the voxel is considered as part of the web if the distance of its center from the web triangle is less than half the thickness of the web, w/2. Qualitatively, the particular reconstruction procedure is in line with the requirement set by the Young Laplace equation [22], namely, that the sum of the inverse values of the local radii of curvature is constant. On several occasions, depending on the preparation route and conditions, the coating around the fibers and/or the webs are not completely filled with the binding material but are made up of particulate matter with solid volume fraction 4i<1, where i indicates the phase of interest, that is, i ¼ c for the coating layer, i ¼ w for the web area. This is simulated as follows. The working domain is discretized into voxels of uniform size. If the center of a voxel is located within a distance w/2 from the plane MCD, the voxel is registered as part of the web, provided that R < 4i, where R is a randomly selected number in the interval (0,1). The same procedure can be applied to the construction of the coating shell around a fiber. A typical example of the reconstructed GDL is shown in Fig. 1b. The porosity of the sample is calculated as follows. A check is made whether the center of each and every voxel of the discretized sample lies within either the fibrous part of the sample or some coating or web region. The number fraction of these voxels defines the solid volume fraction 4T of the sample, that is, 4T ¼ 1-εT, where εT is the total porosity. The solid volume fraction of the uncoated fibers is defined as 40 ¼ 1-ε0, where ε0 is the porosity of the fibrous structure prior to application of coating or webs. Obviously, 4T ¼ 40 þ 4w and ε0 ¼ 4w þ εT. The PTFE content, u, that is usually mentioned in commercial GDL specifications is the mass fraction of PTFE and can be calculated from
u¼
ε0 εT
að1 ε0 Þ þ ε0 εT
(1)
where a is the ratio of the fiber density to the PTFE density. Wrapping of fibers by binder or PTFE can also be simulated in a fashion that reproduces the particulate nature of these additives as detected from SEM images. Specifically, a certain voxel may or may not be occupied by additive even if it belongs to the coating shell of predefined radius, as mentioned above. One may prescribe a probability of additive attachment on the fiber surface that defines the desired coverage of fibers (Fig. 1c). The corresponding
182
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
Fig. 1. (a) Simulated GDL structure. The fibrous part of the structure is shown. (b) Geometrically simulated GDL structure following PTFE coating and web formation. (c) Simulated GDL structure with webs developing between selected fibers. The particulate nature of webs and coatings is shown to simulate limited PTFE sintering.
technique along with the source of their values. 2.3. Simulation of GDL compression
Fig. 2. Definition of web triangle in three-dimensional representation of fiber crossing.
probability value can be set to reproduce the desired solid fraction of the shells and webs, following observations of SEM images. Compact PTFE films are obtained for the value of 1, which is used in the present simulations to keep the number of parameters to a minimum. Although several types of structural data may be requested for the maximum degree of reliability and accuracy in this reconstruction, in the present work it is suggested that the fiber coating and the webs share the same geometrical properties, namely, thickness and threshold probability for solid voxels. The size of the webs is more easily extracted from SEM images; one may use an average value for it or sample random sizes from the size distribution that is extracted from SEM images. Table 1 lists the structural parameters that are associated with the reconstruction
Depending on the PEM type and configuration, the GDL may undergo different levels of compression during assembling of the MEA. The compression of the GDL is known to have significant effects on critical GDL properties, like permeability, diffusivity, conductivity, water condensation, and water transport. In addition, the effects of the transport properties on the anisotropy, which is always present in GDLs, are also expected to be quite significant. In this work, compression is simulated in a simple fashion, namely, by translating the z-coordinates of each point in the reconstructed GDL by a prescribed compression factor. Although this procedure may slightly alter the degree of overlapping of certain fibers, it is expected that the concomitant changes in the GDL configuration will be negligible. Images of simulated uncompressed and compressed variants of the same GDL structure are shown in Fig. 3. 3. Reconstruction of GDL using the lattice-Boltzmann technique The two-phase lattice Boltzmann technique has been developed with the aim to simplify the simulation of flow of two phases under either displacement or steady state flow conditions, and in addition, to model the equilibrium configuration of immiscible phases both in the bulk and also in the vicinity of solid boundaries. Wetting phenomena and effects on the behavior of two-phase systems can also be studied taking into account the degree of surface wettability
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
183
Table 1 Structural parameters in GDL reconstruction. Parameter
Source
Base case
Porosity of neat fiber layer Porosity after addition of binder/PTFE PTFE content Thickness of fiber coating Thickness of webs Solid fraction of fiber coating Solid fraction of webs Web length Number fraction of coated fibers Number of fraction of webbed fibers Fiber size Number density of fibers Thickness of GDL
supplier, porosimetry supplier, porosimetry Supplier SEM SEM SEM, extracted SEM, extracted SEM SEM SEM SEM, supplier SEM, supplier supplier
0.84 0.80 20% 4 mm 4 mm 0.1 0.7 100 mm 0.5 0.5 8 mm (mean) 4.7 105mm3 190 mm
Fig. 3. (a) Uncompressed GDL layer, reconstructed using the method presented in this work. Porosity: 0.86. (b) Compressed GDL layer in the z-direction. Porosity: 0.72. Rate of compression 2:1. Only the fibrous part of the GDL sample is shown.
as quantified through the contact angle. In this work, the free energy-based two-phase lattice Boltzmann model [23] is employed, as it was incorporated in the lattice Boltzmann simulator [24,25]. The mathematical treatment that was developed in this simulator allowed restoration of the Galilean invariance and the realistic simulation of moving droplets and interfaces in two and three dimensions. Briefly, the evolution of particle distribution functions, discrete in both time and space that obey the lattice-Boltzmann equation is monitored. The particle populations move along unit vectors
connecting neighboring nodes. For the present calculations, the three-dimensional Faced Centered Cubic (FCC) lattice is used with 14 different directions of particle propagation at each time step (D3Q15 lattice, Fig. 4a). Using droplet simulation and calculations, it has been shown [25] that the simulator is consistent with the Laplace law and Gibbs-Thomson equations for two-phase systems. The wettability of solid walls can be adjusted by prescribing a fluid density profile along the interfacial sites. The simulator accounts also for phase transition, which allows the simulation of evaporation or condensation phenomena. The reconstruction of binder or PTFE coatings and webs between the GDL fibers proceeds as follows. Once the fiber system is constructed, the desired amount of additive (binder or hydrophobicity promoter) is set and the contact angle of the additive on the fiber surfaces is fixed. Next, a fraction of sites that corresponds to the additive volume fraction is assigned a fixed “liquid” density value whereas the rest are assigned a “gas” density value. The terms liquid and gas here are coined arbitrarily and are understood as distinct, immiscible phase density values. The site densities are then allowed to redistribute themselves subject to the latticeBoltzmann equations, which eventually lead to phase separation thanks to the action of surface tension: for small contact angles, the additive adheres to the fibers and creates connecting areas (webs) at fiber intersections, the extent of which depends on the additive fraction and the contact angle. In fact, the same method can be extended to simulate several stages of GDL processing, including addition of binder to stabilize the layer, and subsequent addition of PTFE to promote further hydrophobicity. The simulation of sequential processing stages involves successive addition of the material of interest and intermediate solidification. More specifically, binder material can be added to the original fiber assembly using the lattice-Boltzmann technique that was presented above. Once the binder phase has organized itself around the fiber connections, the binding “liquid” phase of the two-phase lattice-Boltzmann simulator is assumed to solidify at that position, ready to admit the addition of PTFE in a similar way, using again the two-phase lattice-Boltzmann simulator. The PTFE amount that was introduced during this stage is then allowed to solidify after completion of PTFE-organization around and between the fibers and binder (Fig. 4). The procedure can be used for as many layers of additive as desired, depending on the application. Clearly, one can define different wettability conditions at each step, for instance through the definition of different contact angle, and different volume fraction of the additive. The main advantage of the present simulator is that all these stages can be simulated in a sequential yet automated fashion, reproducing GDL fiber coating and additive spreading between the fibers in a spontaneous manner that respects wettability requirements. The
184
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
Fig. 4. Reconstruction of GDL and hydrophobic additives using the lattice-Boltzmann method. (a) Lattice-Boltzmann unit cell; (b) untreated fiber system; (c) simulated addition of binder using the LB technique; (d) simulated subsequent addition of PTFE using the LB technique. Thickness of binder layer: 1 mm. PTFE content: 30%.
method allows incorporation of practically every structural parameter that appears on Table 1. 4. Flow permeability calculation The reconstructed samples of GDL can be directly used for the numerical calculation of transport properties of gas species along the principal directions of the layer (through-plane and in-plane properties). Based on visual observations of SEM images of GDL samples but also on available experimental data, anisotropy in GDL samples is quite pronounced and clearly increases upon impregnation of the samples with binder and PTFE, ought to the formation of extended webs. The latter develop in direction parallel to the macroscopic surface of the gas diffusion layer and are expected to reduce more drastically the transport properties in the throughplane direction than in the in-plane direction. Two different numerical techniques are used in this work for the calculation of the flow permeability, namely, the single-phase flow lattice-Boltzmann technique and the Finite Element Method (FEM). Both techniques can be applied to the purely geometrical
reconstruction domains as well as to the ones that are produced by the two-phase lattice-Boltzmann method. The advantage of employing the lattice-Boltzmann option in the latter case is that it is free of input/output (I/O) issues as it is applied directly to the LB reconstructed domain. In this way, an integrated lattice-Boltzmann technique is formulated that uses a single code to handle both the reconstruction stage and the permeability computation stage. A brief description of the LB flow simulator is given next. The simulator uses the three-dimensional, 15 speed (D3Q15) unit cell model and involves the evolution of population distribution functions in space and time. The evolution of the system follows the LB equation [26].
fi ðr þ ei ; t þ 1Þ ¼ fi ðr; tÞ
fi fieq
t
;
i ¼ 0; …; 14
(2)
where r is the position vector, ei is the unit velocity vector in direction i, fi is the corresponding population distribution function, t is the single relaxation parameter according to the Bhatnagar, Gross, and Krook (BGK) approach [27] and feq i is the equilibrium
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
185
distribution function at position i, given by the following expansion in terms of the local velocity, y [25,28]:
1 1 1 fieq ¼ wi r 1 þ 2 ðei ,yÞ þ 4 ðei ,yÞ2 2 y2 ; cs 2cs 2cs w0 ¼
2 9
;
1 wi ¼ ; i ¼ 1; …; 6 9
;
wi ¼
1 ; i ¼ 7; …; 14 72
(3)
where r is the fluid density. The speed of sound in this system is
1 cs ¼ pffiffiffi 3
(4)
and the corresponding kinematic viscosity for the D3Q15 model is:
n0 ¼
2t 1 : 6
(5)
An orthogonal computational domain is used for the present study. Flow is induced by prescribing a pressure difference between the inlet and outlet faces of the working domain, whereas the rest of the faces are treated as impermeable boundaries (walls) or as symmetric boundaries. The latter choice is adopted for the sake of comparison with permeability results from other computational techniques, like the FEM method that is mentioned below. The no slip boundary condition is applied at wall sites and is implemented by assuming particle bounce-back at a position that is located halfway between the solid sites and the adjacent fluid sites. An open source LB-BGK software using the D3Q19 unit cell [http://palabos. org] has also been used in the literature [29,30] to calculate flow permeability in GDL domains, adopting the probabilistic geometric approach of [8] to fill the fibrous structure with the desired binder quantity. To simplify the computation, the single relaxation time option is used, which is found to converge sufficiently rapidly for domain sizes of 4 106 sites. In order to avoid problems with flow calculations that depend on the value of the relaxation time, a parametric study was carried out using as varying quantity the kinematic viscosity or, equivalently, the relaxation time. It was found that the flow permeability result remained nearly constant provided that the viscosity was varied within a sufficiently wide interval, namely, between 0.01 and 0.1. The second numerical technique that is used here makes use of commercial software COMSOL and allows various levels of meshing of the reconstructed domain. An appropriate software was required to allow interfacing of the flow solver with the structure regeneration process, including the main fiber assembly and the PTFE coating and web structures.
Fig. 5. Simulated variation of mean pore diameter and specific area with PTFE content.
also plotted in this figure, as extracted from the expression
5. Results and discussion Input data for the simulation studies were taken from the literature and from SEM image analysis. Typical input values are listed in Table 1 for a base-case study. Both of the algorithms that were presented above for the reconstruction of the PTFE coating and web formation, namely, the purely geometrical approach and the two-phase lattice-Boltzmann model, were used to generate GDL samples in three dimensions. In the latter case, strong wettability conditions were implemented to reproduce a small contact angle, equal to 10 . Subsequently, the flow permeabilities were calculated using both the single-phase lattice-Boltzmann approach and the FEM method. In Fig. 5, the variation of the specific surface area, Sv, with the PTFE content is shown, using the lattice-Boltzmann method for the impregnation of GDL with PTFE. The mean pore diameter
Fig. 6. Typical flow field in the interior of a reconstructed GDL sample. Results are shown on the same cross section of the sample, perpendicular to the GDL orientation, before (a) and after (b) the addition of hydrophobic agent (u ¼ 0.3). Block arrows indicate main flow direction (through-plane).
186
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
internal surface area clearly increases with the addition of PTFE. Flow and permeability calculations are presented next. Typical flow fields that are obtained before and after the addition of hydrophobic agents are shown in Fig. 6. As expected, addition of PTFE (u ¼ 0.3)leads to narrower openings between fibers, which increase flow resistance in the through-plane direction and reduce the velocity component in this direction, diverting flow into the in-plane direction, as will be discussed below. Fig. 7a shows the reduction of the through-plane permeability upon compression. The case of a neat fiber structure is considered here, which is free of binder or PTFE presence. At the uncompressed stage, the fibers are placed within the working domain exercising care to avoid overlapping, as explained in the reconstruction section above. Evidently, some limited degree of overlapping will develop during compression as the result of the simplified thickness reduction procedure that was adopted in this work. The compression level is expressed in this figure through the reduction of the porosity of the original, uncompressed structure and, alternatively, through the reduction of the thickness of the GDL. The relationship between the GDL thickness and porosity is given by
Lz 1 ε0 ¼ Lzi 1ε
(6)
where Lz and ε indicate the thickness and the porosity of the GDL sample, respectively, whereas superscript 0 is used to indicate values at the original, uncompressed stage. The simulation results were also compared with experimental data for Toray paper [7] using as reference value the permeability of the uncompressed structure. The agreement is excellent and indicates very satisfactory representation of the GDL for every compressed level employed here. The dimensional permeability of the uncompressed sample was calculated to 8.5 1012 m2, which also compares very well with the experimental value of 8.7 1012 m2 in Ref. [7] and 8.99 1012 m2 in Ref. [4]. This reflects the success in the reconstruction of the uncompressed GDL sample. It is also important to note that the two numerical solvers of the flow equations, namely, the LB and the FEM solvers, yield permeability predictions that are in very good agreement with each other, despite the fact that the LB method uses the centers of the reconstruction voxels as sites for the density function propagation, and, hence, for the flow “meshing” part, without further refinement. Based on this observation, and in order to save on computational cost and effort, flow calculations are carried out using the LB technique unless otherwise noted. The pressure drop between neighboring flow channels gives rise to flow in direction normal to the main transport direction, which is known to delay onset of mass transfer limitations during cell operation. To elucidate structural effects on the in-plane flow phenomena, flow computations were carried out in this direction as well, leading to the estimation of the in-plane permeability. The results are shown in Fig. 7b and compared to experimental data from Ref. [4]. Again, very good agreement is noticed over the entire compression range (0e30%) with respect to the uncompressed GDL sample. The simulation results in the two principal in-plane directions, x- and y-, were averaged in order to produce orientationally averaged results that can be compared to experiments. The degree of anisotropy of the GDL structure is then expressed by the ratio of inplane to through-plane permeabilities and is found to range between 1.7 and 2.4 depending on the degree of compression (Fig. 7c) using the LB technique. This is of strong interest in the design and performance of the electrode assembly: increased transport in the direction normal to the main transport direction is known to affect gas distribution at the entrance of the catalyst layer and the
Fig. 7. (a) Dependence of reduced permeability on compression level expressed in the form of reduced porosity (bottom axis) and thickness reduction (top axis) for a reconstructed GDL sample, free of binder or PTFE addition. Comparison of the predictions of the two numerical techniques that are employed here of (a) through-plane values with experimental ones from Ref. [7], (b) in-plane values with experimental ones from Ref. [4], and (c) compression effect on the GDL anisotropy, expressed as ratio of in-plane to through-plane permeabilities.
concomitant demands for catalyst performance. Results on the effect of GDL impregnation with binder and
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
187
Fig. 8. (a) Effect of PTFE loading on through-plane permeability for a GDL reconstruction. PTFE addition was implemented through the LB technique (filled symbols) and, alternatively, through the geometry-based technique (open symbols). Comparison of numerical predictions for a GDL reconstruction with experimental values from Ref. [7]. (b) Prediction of the effect of PTFE loading on the GDL anisotropy, expressed as the ratio of in-plane to through-plane permeability. PTFE addition was implemented through the LB technique.
hydrophobicity-enhancing materials, like PTFE, on the flow permeability are presented next. Both techniques that were developed in the present work are employed, namely, the latticeBoltzmann technique using the two-phase simulator, and the reconstruction from SEM images using a purely geometrical
approach. In the lattice-Boltzmann case, the volume fraction of the impregnating additive was varied keeping all other properties constant, that is, surface tension, viscosity, and density. For each value of the volume fraction, the two-phase lattice-Boltzmann simulator was employed to allow spreading of the liquid around
188
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189
and between the GDL fibers, as described above. Once equilibrated to a configuration that is no longer propagating, the liquid part is assumed to “solidify” and the flow permeability is calculated using the lattice-Boltzmann simulator for single phase flow. The results for the through-plane permeability are shown in Fig. 8a and compared to experimental data from Ref. [7]. The agreement to the experimental values is quite satisfactory if one considers the stochastic nature of the GDL reconstruction and the fact that the impregnation process was simulated using a two-phase simulator neglecting dynamic invasion phenomena. The computation process can be repeated for the in-plane permeability and its ratio to the corresponding through-plane permeability is plotted against the PTFE content in Fig. 8b. It is noteworthy that anisotropy is enhanced by PTFE addition, due to the unequal effect of PTFE addition on the permeability in the various directions of the GDL sample. Specifically, the flow permeability is more drastically reduced upon increase of PTFE loading in the through-plane direction than in the in-plane direction. In the case of PTFE spreading reproduction by the geometrical approach, SEM images are directly used in the way that was described above. The corresponding flow permeability results are also shown in Fig. 8a in the form of open symbols, for the sake of comparison with the lattice-Boltzmann reconstruction procedure. Although some deviation between the two sets of permeability predictions is apparent, the general trend in the behavior of the permeability is the same. This deviation is obviously attributed to the different approach that was used to simulate PTFE spreading and to the simplifying assumptions that were made in the geometrical reconstruction case regarding the shape and thickness of the webs. 6. Conclusions An integrated simulator is developed in this work that is capable of reproducing in three dimensions the structure of gas diffusion layers and simulating flow and diffusion in the reconstructed samples. The simulator is based on the LB technique, which is appropriately adjusted to simulate addition of binder to the carbon fiber system and the subsequent addition of PTFE or any other additive to promote further the hydrophobic characteristics of the layer. Among the advantages of the present simulator compared to other approaches is the fact that fiber coating by binder and further coating by PTFE accompanied by web formation between adjacent fibers can be handled in a physicochemically consistent manner, as it takes into account wettability conditions that shape the web configuration across the structure. A second advantage is that all stages of the process, that is, reconstruction in three dimensions, simulation of flow, and calculation of permeabilities, can be tackled by the same simulator without having to resort to different numerical solvers for the different stages of simulation. An alternative reconstruction procedure was also developed in this work that is based merely on geometrical features, as these are extracted from SEM images of GDL samples. Comparison with experimental data for the GDL permeability reveals very satisfactory agreement in both main directions, that is, for the through- and the in-plane permeability. Compression effects are also simulated and the results are shown to be in quantitative agreement with experiments. Addition of PTFE leads to significant reduction of permeability in both main directions. The degree of flow anisotropy is also estimated and is found to be dependent on the compression level and on the PTFE loading. This anisotropy is known to have a significant role in the design and performance of the fuel cell as it affects the gas distribution in both anode and cathode sides, which in turn influences the onset of mass transfer limitations during cell operation. The simulator that is developed in this work can be further
utilized for a comprehensive simulation of the phenomena that take place inside a GDL, namely, gas transport, water condensation, and water transport in both the anode and the cathode of the cell. Breakthrough conditions can also be identified using as parameters the water content, the hydrophobicity of the layer, and the applied pressure drop. In fact, these conditions can be investigated using the same LB simulator given that it inherently accounts for phase transition and for dynamic two-phase flow, provided of course certain appropriate adjustments and fine-tuning are carried out. This work is already in progress by the authors. Acknowledgements Financial support to this work was provided by the Institute of Chemical Engineering Sciences of the Foundation for Research and Technology, Hellas, and by the European Union's Seventh Framework Programme (FP7/2007e2013) under FCH-JU Project Grant number 303024, EURECA. APPENDIX
A. Distance of point re0 from line r ¼ rei þ t a: e e
dðr ;r Þ¼ðr r Þa
e0 e
e0 e1 e a e
B. Distance of two lines in 3D Cartesian space. yy1 yy2 zz1 xx2 zz2 1 The distance of lines xx [1 ¼ m1 ¼ n1 and [2 ¼ m2 ¼ n2 is given by
x2 x1 y2 y1 z2 z1 [1 m1 n1 [2 m2 n2 jð[1 ; m1 ; n1 Þxð[2 ; m2 ; n2 Þj C. Conditions for point (x,y,z) to lie in the interior of the triangle that is defined by the corners (xi, yi, zi), i ¼ 1,2,3: (a) a, b > 0 (b) aþb < 1where
x x1 x3 x1 a ¼ x2 x1 x3 x1
x x1 y y1 x x1 y3 y1 ; b ¼ 2 y2 y1 x2 x1 x3 x1 y3 y1
y y1 y2 y1 y2 y1 y3 y1
References [1] F. Barbir, Sol. Energy 78 (2005) 661e669. [2] M.M. Daino, S.G. Kandlikar, Int. J. Hydrogen Energy 37 (2012) 5180e5189. [3] Y. Gao, X.X. Zhang, P. Rama, Y. Liu, R. Chen, H. Ostadi, K. Jiang, Fuel Cells 12 (2012) 365e381. [4] J.T. Gostick, M.W. Fowler, M.D. Pritzker, M.A. Ioannidis, L.M. Behra, J. Power Sources 162 (2006) 228e238. [5] L. Hao, P. Cheng, J. Power Sources 186 (2009) 104e114. [6] J.H. Nam, M. Kaviany, Int. J. Heat. Mass Transf. 46 (2003) 4595e4611. [7] A. Tamayol, F. McGregor, M. Bahrami, J. Power Sources 204 (2012) 94e99. [8] R. Thiedmann, F. Fleischer, C. Hartnig, W. Lehnert, V. Schmidt, J. Electrochem. Soc. 155 (2008) B391eB399. [9] N. Zamel, X.G. Li, Prog. Energy Combust. Sci. 39 (2013) 111e146. [10] G. Inoue, T. Yoshimoto, Y. Matsukuma, M. Minemoto, J. Power Sources 175 (2008) 145e158. [11] X.B. Nie, G.D. Doolen, S.Y. Chen, J. Stat. Phys. 107 (2002) 279e289. [12] V.P. Schulz, J. Becker, A. Wiegmann, P.P. Mukherjee, C.Y. Wang, J. Electrochem. Soc. 154 (2007) B419eB426. [13] N. Parikh, J.S. Allen, R.S. Yassar, Fuel Cells 12 (2012) 382e390. [14] M.F. Mathias, J. Roth, J. Fleming, W. Lehnert, Diffusion media materials and characterisation, in: Handbook of Fuel Cells, John Wiley & Sons, Ltd, 2010. [15] J.T. Gostick, M.A. Ioannidis, M.W. Fowler, M.D. Pritzker, J. Power Sources 173 (2007) 277e290.
V.N. Burganos et al. / Journal of Power Sources 365 (2017) 179e189 [16] A. Nabovati, J. Hinebaugh, A. Bazylak, C.H. Amon, J. Power Sources 248 (2014) 83e90. [17] H. Sadeghifar, N. Djilali, M. Bahrami, J. Power Sources 248 (2014) 632e641. [18] J. Yablecki, A. Bazylak, J. Power Sources 217 (2012) 470e478. [19] H. Ostadi, P. Rama, Y. Liu, R. Chen, X.X. Zhang, K. Jiang, J. Membr. Sci. 351 (2010) 69e74. [20] A. Pfrang, S. Didas, G. Tsotridis, J. Power Sources 235 (2013) 81e86. [21] P.K. Sinha, P. Halleck, C.Y. Wang, Electrochem. Solid State Lett. 9 (2006) A344eA348. [22] F.A.L. Dullien, Porous Media: Fluid Transport and Pore Structure, first ed., Academic Press, New York, 1979.
[23] [24] [25] [26] [27] [28] [29]
189
M.R. Swift, W.R. Osborn, J.M. Yeomans, Phys. Rev. Lett. 75 (1995) 830e833. A.N. Kalarakis, V.N. Burganos, A.C. Payatakes, Phys. Rev. E 65 (2002). A.N. Kalarakis, V.N. Burganos, A.C. Payatakes, Phys. Rev. E 67 (2003). Y.H. Qian, D. Dhumieres, P. Lallemand, Europhys. Lett. 17 (1992) 479e484. P.L. Bhatnagar, E.P. Gross, M. Krook, Phys. Rev. 94 (1954) 511e525. S. Chen, Z. Wang, X. Shan, G.D. Doolen, J. Stat. Phys. 68 (1992) 379e400. D. Froning, J. Brinkmann, U. Reimer, V. Schmidt, W. Lehnert, D. Stolten, Electrochim. Acta 110 (2013) 325e334. [30] D. Froning, G. Gaiselmann, U. Reimer, J. Brinkmann, V. Schmidt, W. Lehnert, Transp. Porous Media 103 (2014) 469e495.