- Email: [email protected]
Modeling of composite fibrous porous diffusion media
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Modeling of composite fibrous porous diffusion media Sima Didari a, Arash Asadi b, Yan Wang a, Tequila A.L. Harris a,* a
George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, 813 Ferst Dr., Atlanta, GA 30332, USA b School of Mathematics, Georgia Institute of Technology, 686 Cherry Street, Atlanta, GA 30332, USA
article info
abstract
Article history:
To engineer the desired properties of fibrous porous media, a parametric modeling
Received 7 October 2013
approach is needed to support the rational design of the materials before the fabrication. In
Received in revised form
this study, we propose a methodology that enables the accurate representation of three-
27 March 2014
dimensional (3D) microstructures of fibrous porous media and prediction of their trans-
Accepted 2 April 2014
port properties. Toray TGP-H-060 gas diffusion layer (GDL) is selected as an example to
Available online 29 April 2014
demonstrate the feasibility of the suggested design methodology. The detailed microstructure of the GDL with the inclusion of locally distributed binder is constructed using an
Keywords:
extended periodic surface (PS) modeling technique. A 3D morphological approach is taken
Modeling
to create the binder distribution within the fibrous microstructure. Transport properties
Composite
including permeability, relative diffusivity, and tortuosity and local structure characteris-
Gas diffusion layer
tics of the generated microstructure, under different binder loading are calculated. It is
Porous media
shown that the detailed model of the fiber-binder composite has a strong influence on the
Transport properties
predicted properties.
Design
Copyright ª 2014, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Introduction Fibrous porous media have been used for many decades in a plethora of industries such as bioengineering, electronics, energy and environmental processes [1e6]. As an example, gas diffusion layers (GDL) are used in polymer electrolyte membrane fuel cells (PEMFC), a burgeoning renewable energy resource. GDLs serve as electrodes in PEMFCs, which permeate gases and help remove water droplets that are formed during fuel cell operation. Furthermore, these carbon based GDLs have other properties such as high porosity, strength-to-weight ratio, and electrical and thermal
conductivities that make them attractive for use in PEMFCs and other areas [7]. It has been shown that the performance of PEMFCs depends upon the morphological, transport, and mechanical properties of the GDLs [8e12]. However, no systematic and parametric modeling approach that help predict and optimize the key characteristics of porous materials such as GDLs for the purpose of materials design exists to date [13]. Rather, studies have been conducted to determine the morphological and transport properties of the fibrous diffusion media that have been fabricated [14]. In these studies, the microstructure of the porous media was generated by image-based reconstruction or geometrical construction approaches. In the
* Corresponding author. E-mail addresses: [email protected], [email protected] (T.A.L. Harris). http://dx.doi.org/10.1016/j.ijhydene.2014.04.011 0360-3199/Copyright ª 2014, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
9376
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
image-based approaches, 2D images or 3D voxel representation of the internal pore structures are obtained by scanning electron microscopy (SEM), synchrotron-based tomography, confocal microscopy, or computed microtomography (mCT) [15e24]. Although imaging techniques accurately capture the geometry of porous structures, the model of the microstructures reconstructed from the images cannot be modified. Thus, imaged-based approaches are not suitable for designing microstructures of porous media. The geometrical construction approaches that utilize conventional computer aided design (CAD) tools are not efficient in modeling complex structures. Because in these approaches, a large number of discrete control points are required to accurately represent the topology of the porous media. Conversely, implicit surface modeling approaches that construct the microstructures based on closed-form mathematical functions are much more promising [25]. The geometrical construction of fibrous diffusion media has been primarily based on statistical methods. Fibers are assumed to have simple shapes such as long cylinders and randomly distributed in a representative volume element (RVE). In the work of Schulz et al. [26], the fiber skeleton of GDLs was modeled as a collection of intersecting cylinders oriented along randomly distributed axes. The directional distribution of the fibers was assumed as normal or determined from the SEM images. Thiedmann et al. [24] constructed GDL microstructures by implementing a Poisson planar line tessellations (PLT) process to represent the fibers’ axes. Gaiselmann et al. [27] represented the pore space of fibrous media with random geometric graphs. The parameters required for modeling were specified such that the distributions of vertex degrees and edge lengths conform with X-ray synchrotron tomography images. GDLs have also been represented with simplified geometries such as arrays of parallel and perpendicular cylinders or channels [28,29]. Fibrous porous media are often more complex than a collection of fibers. They also include some binding material(s). To model the binder in the microstructure of the GDL, the pore space of the generated fiber skeleton is filled using various approaches such as stochastic methods with predefined probability functions or morphological operations [24,30]. To quantitatively characterize the pore space of a porous medium and to predict its bulk material properties, the morphological properties of the generated microstructure (e.g., pore size and shortest path distributions) have been also studied. The pore size has been calculated as the spherical distance of the points inside the void space to the surrounding solid boundaries [31e33]. However, in some studies, the pores were defined as 2D closed areas and the pore size is defined as the diameter of a circle with an area equal to the closed area [34,35]. To determine the inhomogeneity and directionality of the distribution of pores in the microstructure, the chord length distribution of a porous medium was also found by evaluating angularly resolved lineal path functions [36]. The distribution of the shortest path within a porous microstructure has been found from the voxel mesh or the graph representation of the porous domain [24,36], using Dijkstra’s algorithm [37]. The existing numerical schemes offer useful platforms and tools for the evaluation
of some morphological properties. However, they are computationally expensive. In the current study, these schemes are modified for higher accuracy and lower computational time. Transport properties of microstructures can be predicted using different approaches such as pore network (PN), lattice Boltzmann models (LBM) and computational fluid dynamics (CFD). The effects of change in porosity due to the compression forces and fiber size and spatial orientation have been investigated through permeability, tortuosity, diffusion and relative permeability of fibrous porous media [14]. The CFD approach is more time efficient than LBM and can simulate larger domains. The PN modeling approach uses a simplified representation of the pore structure and is efficient for large domains. However, it requires the information of length correlations that need to be extracted from an existing structure. The process itself is time consuming. The objective of this work is to introduce a platform for designing composite porous media, with the characteristics that are similar to those of GDLs. A new methodology is developed to generate the fiber skeleton of composite porous diffusion media. In this study, a periodic surface (PS) modeling technique developed by Wang [25] is used to generate the fibrous component of the microstructures, and the binder distribution is modeled by the Minkowski sum of fibers and a structure element (SE). The microstructure generation process is directly integrated with the morphology and transport characterization analyses.
Material Toray TGP-H-060 carbon paper GDL, shown in Fig. 1, is chosen as the fibrous composite microstructure to examine the feasibility of the microstructure generation method. Toray TGP 060 is made of a matrix of carbonaceous binder that is a thermoset resin and graphitized rigid carbon fibers [22]. Typical parameters of the GDL include, a fiber diameter of 7 mm, a thickness of 200 mm, porosity of 0.78 and 0.27 solid fraction of carbonaceous thermoset binder [30].
Implicit surface modeling of fibrous structures An implicit surface S ¼ {r ˛ R3:j(r) ¼ 0} is defined as a set of points that are the roots of a given function j:R3 / R with realvalued coordinates. Recently a periodic surface (PS) model was developed by Wang [25] to implicitly model structures with complex topology, which is define as jðrÞ ¼
L X M X
mlm cos 2pkl PTm $r ¼ 0
(1)
l¼1 m¼1
where j(r) is the implicit function, mlm is periodic moment, kl is the scale parameter, Pm ¼ [am, bm, cm, am]T is a basis vector for the major phase, [am, bm, cm]T represents the normal vector of the basis plane, am is the phase of the plane, r ¼ [x, y, z, 1]T is the location vector in homogenous coordinates. The PS model of a rod surface with the periodicity of 1 along the z-axis direction in an RVE is:
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
9377
Fig. 1 e SEM images of Toray GDL, (a) Front view, (b) Cross sectional view.
jðx; y; zÞ ¼ 4 cosð2pxÞ þ 4 cosð2pyÞ þ 3 ¼ j0
(2)
where j0 is the adjustable isovalue. To construct many fibers in a section of a fibrous medium, the PS model of each fiber has to be constructed and inserted in an RVE. The whole microstructure with n fibers can be generated in an RVE with the dimensionless lengths using the union operation jW ðrÞ ¼ minðj1 ðrÞ; j2 ðrÞ; .; jn ðrÞÞ
(3)
over the implicit functions of all individual fibers in the domain. Due to the scalable nature of PS model, the generated RVE can be readily scaled to the desired size.
Geometrical modeling of Toray TGP-H-060 GDL The composite structure of the GDL is constructed in two steps. First, the fiber skeleton is generated using the PS model. Then the binder is integrated into the skeleton of the microstructure. To generate the structure of fiber skeleton, the diameter of the fibers (D) and the solid volume fraction or the porosity of fiber skeleton have to be known. The diameter of the fibers is determined by the isovalue (j0) of the rod PS model. As the isovalue changes, the diameter of the rod surface changes. The diameter of the fiber can be approximated by the distance between the positions with the maximum ycoordinate Ymax and minimum y-coordinate Ymin in the cross section of the rod as D ¼ Ymax Ymin [4]. Since the period is set to be 1, the cross section of the rod in the xey plane of the RVE is symmetric about x ¼ 1/2. Both Ymax and Ymin are located at x ¼ 1/2. Therefore, 1 þ 4 cosð2pyÞ þ 3 ¼ j0 4 cos 2p 2
(4)
The simplified form of Equation (4) is cos(2py) ¼ (1 þ j0)/4, which has two solutions for y ˛ [0,1], which are Ymax and Ymin. Thus, the diameter can be found as: D ¼ Ymax Ymin ¼ 1
arccosðð1 þ j0 Þ=4Þ p
(5)
The porosity of the GDL is related to the number of fibers in the RVE. To decide the orientation of a fiber in the RVE, two
angles are needed. One is the in-plane angle in the xey plane, q, and the other is the side angle in the zex plane, 4. These two angles determine the orientation of major phase vector Pm of each of the implicit functions. In a semi-layered structure, the ranges in which q and 4 can change are set to be [0, p] and [p/ 36, p/36] respectively. The positions and orientations of the fibers follow the normal distributions in the RVE. The randomly generated fibers may intersect with each other in the RVE. The overall intersection volume of the generated structure is quantified by U ¼ jX ðrÞ ¼ min max ji ðrÞ; jj ðrÞ 1ijn
(6)
To avoid the unrealistic phenomenon of intersection, U in Equation (6) is minimized as min U (qi, 4i, xi, yi, zi)i¼1:n where n is the number of the fibers, qi, 4i, xi, yi, zi represent the orientation and position of the ith fiber. The optimization problem is solved by the genetic algorithm (GA) and implement in MATLAB. The optimization procedure is described schematically in Fig. 2. Assuming negligible deformation for carbon fibers due to their rigidity, the microstructure of the fiber skeleton is generated as shown in Fig. 3(a).
Binder impregnation in Toray TGP-H-060 GDL Because of the low viscosity and high wettability of the resin and the rapid impregnation process, in this work only the final composite structure is modeled. The resin-fiber bond that forms the composite nature of a GDL is modeled in two steps. In the first step, a Minkowski sum between the fiber skeleton (denoted by A) and an SE that represents the resin is performed, which is defined as B ¼ A 4 SE ¼ {a þ seja ˛ A, se ˛ SE}. The result is the thicken fiber bundle as shown in Fig. 3(b). In the second step, an erosion operation is conducted on set B. The erosion, defined as B . (SE), is the Minkowski subtraction (shown by (.)) of B with the reflected set of SE. As shown in Fig. 3(c), the erosion process leads to a reduction of set B [38,39], whereby the pores that are smaller than the SE are filled. The PS model and binder impregnation model are developed in MATLAB software and can be easily coupled with
9378
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
of the binder in the generated composite structure is calculated. If the target weight fraction is not met, the filling procedure is conducted with a larger SE. Using a disk shape SE, the microstructure of Toray TGP-H-060 carbon paper is generated, as shown in Fig. 4, in comparison with scanning electron microscope (SEM) images. As depicted in the figure, the structure has a semi-layer-by-layer configuration that resembles an actual GDL.
Morphology characterization Calculating porosity
Fig. 2 e Block diagram of generating the fibers skeleton.
microstructure characterization. The fiber skeleton generated using the PS model is converted into a 3D binary voxel representation that serves as an input for the binder impregnation model. The Minkowski summation and erosion procedures are implemented with the MATLAB image processing toolbox [39]. In the binary voxel representation of the fibrous domain coded as a binary matrix, the elements of the matrix that correspond to the points inside the solid phase are set to be 1 and the elements in the void or pore space are set to be 0. Given the PS model defined in the RVE, the value of the implicit function at each node in the discretized grid space is calculated. This resulting value j is compared to the given isovalue j0. If j j0, the node is in the solid phase, otherwise it is located in the pore space. Then, a value of 1 or 0 is assigned to the node. After generating the fiber skeleton using the smallest size for the SE, i.e., the mesh size, the two-step resin impregnation procedure described above is applied and the weight fraction
The binary matrix representation of the porous structure is used to calculate the morphological properties of interest. For a composite porous structure, three separate binary matrices are used to calculate the porosity of the whole structure, porosity of the fibrous skeleton, and the volume fraction of the binder. To find the solid volume fraction of each phase used in the composite structure, i.e., the resin and fibers, the number of elements in the binary matrix that have the value of 1 are counted. Then the number is divided by the matrix size. The porosity of a composite is determined by counting “0” elements in the binary matrix, and the number is divided by the matrix size.
Determining pore size distribution The pore size distribution is determined by finding a set of the largest fitting spheres in the pore space. Given that each node in the pore space, i.e., 0’s of the binary matrix, can be a potential pore center, the implemented search algorithm starts with placing a sphere with the smallest radius (mesh size) at each of these nodes. The radius is increased gradually until that the sphere fits within the pore space. A sphere is considered as fit if it has at least three points of contacts with the surrounding solid surfaces and the three points are not in the same direction. As such, multiple spheres can be contained inside a pore space. Those spheres with centers located inside a larger sphere are deleted. The remaining spheres then form a set of mutually overlapped pores, where none of their centers are located inside another. The radii of the spheres are the pore sizes.
Fig. 3 e (a) Fibrous bundle, (b) Minkowski summation with structure element, (c) Minkowski erosion with the same structure element.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
9379
Fig. 4 e (a) Isometric view of constructed Toray TGP-H-060 GDL, (b) and (c) cross sectional views, virtual GDL (left side) and SEM images (right side).
To compute the pore size distribution that is closer to the one found experimentally, the closed and dead end pores are discarded and not counted in the statistics. The connectivity of the pores is found by determining if there exists a pathway that connects one pore center to the neighboring pore centers. The resulting pore size distribution is shown in Fig. 5. The average pore size for the model generated for Toray TGP-H060 GDLs is 21 mm, which is close to the average pore size of 18 mm measured by Flu¨ckiger et al. [12].
Determining chord length distribution The pore size distribution provides an insight on the void space within the porous structure, but it cannot provide any information regarding the directionality of the pore space. Hence, the chord length (Lch) defined as the distance between the internal solid boundaries in the pore space is used to represent the degree of anisotropy of the microstructure in a desired direction. The chord length is a measure of the available free space in which the particles can move inside a porous structure. Thus, using the chord length distribution, the dominant mode of diffusion inside porous media can also be found. To find the chord length distribution, the cross-
Fig. 5 e Numerical pore size distribution of virtual Toray TGP-H-060 GDL.
sectional 2D binary matrices of the porous medium are extracted from the 3D binary matrix. The chord length is the number of consecutive 0’s in a principal direction multiplied by the size of discrete elements. The chord length distributions for the through-plane and in-plane directions of the model for Toray 060 GDL are shown in Fig. 6, where the chord lengths are normalized as the ratio to the thickness of the GDL. As depicted in the figure, in the through-plane direction the chord lengths are smaller due to the semi-layered structure of the GDL, while in the in-plane direction the chord lengths have a wider distribution.
Determining the shortest path distribution To find the distribution of the shortest path inside the porous structure, Dijkstra’s algorithm [37] is implemented and modified to calculate all of the shortest paths starting from each pore. To implement Dijkstra’s algorithm, a graph representation of the porous domain is extracted. To construct the graph representation of the porous medium, first it is assumed that there is no solid phase in the RVE. A graph G with vertices corresponding to the 3D grid points with the size of nx ny nz is generated to show the connectivity of nodes inside the RVE. Any point in the 3D grid with the position index (i, j, k), where 0 i nx, 0 j ny, and 0 k nz, is a vertex of G. To construct the edges of G, the connectivity of the vertices is investigated. Two vertices (i, j, k) and (a, b, c) are assumed to be connected if and only if ji-aj þ jj-bj þ jk-cj ¼ 1. Graph G is then modified to account for the effect of the solid phase in the RVE. The solid phase acts as an edge eliminator, vertex eliminator, or both, to determine whether a vertex belongs to the void space or the solid phase. Then the vertex and the associated edges in G are deleted. The graph G is a weighted graph where the weight of each edge is 1. It is assumed that all the vertices in the inlet and outlet are in the pore space and that the sides of the RVE are in the solid space. For each vertex at the inlet, the modified version of Dijkstra’s algorithm is run and the length of the shortest path is found. The shortest path from each vertex at the inlet to each vertex
9380
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
Fig. 6 e Chord length distribution for virtual Toray TGP-H-060 GDL, (a) through-plane direction, (b) in-plane direction.
at the outlet is found. The algorithm is implemented in Cþþ and the MEX format is used for data exchange with MATLAB. The distributions of the shortest paths inside the models for Toray 060 GDL along the through-plane and in-plane directions are shown in Fig. 7. The computed average length of the shortest path in the through-plane direction is 1.64, which agrees well with the experimentally measured value of 1.51 [24]. As expected, due to the semi-layer structure of the GDL in the through-plane direction, the GDL has longer pathways or higher tortuosity.
Transport property characterization A high-quality volumetric mesh of the computational domain is needed to determine the transport properties of interest, i.e., permeability, tortuosity, and diffusivity. For geometries that can be represented by CAD software, common mesh generation packages can be used. However, for the complex geometries of the composite structures presented in this work, conventional mesh generation methods are not applicable. Thus, a voxel mesh representation of the computational domain is used.
conservation equations that govern the flow behavior through the porous media are V$u ¼ 0 and ru$Vu ¼ VP þ Vbs
(8)
respectively, where u is the velocity, r is density, P is the pressure field, and b t is the stress tensor. By importing the volumetric mesh into ANSYS FLUENT and solving the mass and momentum equations given by Equations (7) and (8), the permeability of the microstructures is found from Darcy’s law: K V ¼ VP h
(9)
where V is the superficial velocity (i.e., volumetric flow rate per unit cross-sectional area), h is fluid viscosity, and VP is the pressure drop per unit length through the porous material. A mesh sensitivity study was conducted to find appropriate mesh sizes. It is also found that a computational domain of 200 mm 200 mm 200 mm is sufficient to accurately compute the transport properties. To determine the permeability in a direction, the inlet and outlet boundary conditions are set for the specific flow direction, as shown in Fig. 8. At the inlet a
Volumetric mesh representation of the computational domain Using the MATLAB image processing toolbox, the crosssectional images of the porous structure are generated based on 2D arrays, which are extracted from the 3D binary matrix, as illustrated in Fig. 8. These images are saved as a sequence of images and are inserted into Geo-Dict software to generate a volumetric voxel mesh. The whole procedure of mesh generation, including reading the 3D binary matrix, creating the 2D cross-sectional images and importing them into the Geo-Dict software, is fully automated using custom-built MATLAB scripts. This procedure is an efficient and flexible approach for generating the voxel mesh representation of the domain; hence, the microstructures can be easily altered or changed.
Determining permeability of GDLs The flow field inside the pore space is solved to find the permeability (K). The steady state momentum and mass
(7)
Fig. 7 e Distribution of the shortest path inside the pore space of TGP-H-060 GDL.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
9381
Determining the relative diffusivity and tortuosity factor in porous structures
Fig. 8 e Schematic illustrating the generation of a volumetric mesh.
pressure boundary is defined, while at the outlet a pressure outlet boundary condition is imposed for a gauge pressure of zero. The no slip boundary condition is applied to the internal solid fluid interfaces, which are the outer surfaces of the fibers. The momentum equation is discretized by the second order upwind scheme with the SIMPLE algorithm for coupling pressure and velocity [40]. The under relaxation factors of 0.25 and 0.55 for pressure and momentum are used, respectively. After the simulations have converged, the superficial velocity is calculated by volume averaging the velocity fields in the GDL. Since the GDL thickness is known, the pressure gradient is then calculated. Applying the calculated values for the superficial velocity and the pressure gradient to Equation (9), the permeability is computed.
Determining the tortuosity in GDLs
Diffusivity in porous media (Dbulk) is defined as Dbulk ¼ Drel D0, where Drel is relative diffusivity, which represent the reduction of the diffusivity in the porous media due to the presence of solid phase and D0 is the binary diffusion coefficient. D0 ¼ 1/ 3ly, where l is the mean free path of the gas that passes through the porous media and y is the thermal velocity. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y ¼ 8kB T=pm, where kB is the Boltzmann constant, T is temperature, and m is the molecular mass. To determine the mechanism by which diffusion takes place inside the pore space, the Knudson number (Kn) is calculated, which is the ratio of the mean free path of a specific fluid (l) to the available length scale inside the porous structure [41]. This length scale, which is set as the average chord length ðLch Þ, indicates the freedom of fluids to move inside the pore space. In the current study, the diffusivity of O2, H2, and H2O, which are the fuel and byproduct in PEMFCs, are computed. The mean free path for H2, O2, and H2O are presented in Table 1 [36]. Using l and the average value of the chord length found from the morphology characterization, the average Kn numbers are calculated, as given in Table 1. The Knudsen numbers are less than one, indicating that the pore diameters are large compared to the mean free path length of the fluid molecules. Thus, the diffusion in the pores is governed by the concentration Laplace equation given by DC ¼ 0
where C is the concentration of the gas in the pore space. To find Drel, Equation (11) is solved using Geo-Dict software by applying zero flux at the solid surfaces and a concentration gradient in the direction of the desired Drel. For a known diffusive flux (J), diffusivity is then calculated using Fick’s law J ¼ Dbulk VC
1
s¼N
PN i¼1
L
lðri Þ
(12)
Then, Drel can be directly computed and used to calculate the tortuosity factor (sf) by sf ¼
The other important transport property of porous media is tortuosity (s), which represents the complexity of the microscopic flow paths in the porous medium. Tortuosity is defined as the ratio of average length of the flow path to the characteristic length of the porous medium. The thickness of porous medium in the flow direction is usually assumed as the characteristic length. For inhomogeneous porous structures, depending on the direction of flow (in-plane or throughplane), tortuosity can be different. Tortuosity, in a desired direction, is calculated from
(11)
ε Drel
(13)
where ε is the porosity of the porous domain. The through-plane and in-plane permeability, relative diffusivity, tortuosity, and tortuosity factor of Toray 060 GDL are calculated, using the methods discussed, and are given in Table 2. There exists very good agreement between the experimentally measured data and the numerically computed data based on the microstructure(s) generated using the PS model and binder impregnation model developed. As expected, due to the semi-layered structure of the GDL, the transport properties in the through-plane directions are lower
(10)
where N is the number of the particles that are tracked, l(ri) is the length of the path line that belongs to the ith particle located at ri, and L is the GDL thickness. After convergence of the numerical simulation, for 2500 particles, the path line length is computed and the tortuosity of the GDL is calculated for both the through-plane and in-plane directions.
Table 1 e Average Knudson number in Toray 060 GDL at T [ 100 C and P [ 100 kPa. Fluid type H2O O2 H2
l (m)
Kn number 7
7.5 10 6.1 107 1.3 107
1.2 102 9.7 103 2.1 103
9382
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
Table 2 e Transport properties of Toray 060 GDL. Numerical Experimental
Direction
Permeability (m2)
Relative diffusivity
Tortuosity factor
Tortuosity
Through-plane In-plane Through-plane In-plane
6.8 1012 1.19 1011 5.2 1012 [22] 2.1 1011 [43]
0.42 0.67 0.38 [42] 0.57 [42]
1.85 1.17 2.05 [42] 1.36 [42]
1.06 1.19 NA NA
than the in-plane direction. However, larger discrepancy is observed in the through-plane direction compared to the inplane direction.
Effect of binder impregnation on the characteristics of GDL As mentioned in Binder Impregnation in Toray TGP-H-060 GDL Section, to model composite structures impregnated with binder, a two-step procedure is implemented. The procedure includes Minkowski summation and erosion of the fiber skeleton model by a filler SE. To determine the SE that can accurately represent the composite fibrous structure, a disk SE (2D) and a cube SE (3D) are used to construct the microstructure of Toray 060 GDL, and the morphological and transport properties of both microstructures are found. As shown in Fig. 9(a), the pore size distributions for the two microstructures are different. However, the pore size distribution corresponding to the 2D disk SE shows higher resemblance to the pore size distribution found experimentally [42]. Distributions of the shortest pathways in the through-plane and in-plane directions are found and presented in Fig. 9(b). As illustrated, the GDL generated using the disk SE has an inhomogeneous topology as expected, due to the semi-layered structure of the porous medium of interest. However, the distributions of the pathways in the through-plane and inplane directions of the microstructure generated using the 3D cube SE, are approximately equal. Using the 3D SE results in a uniform distribution of the binder in both directions, which may lead to microstructures that are not semi-layered, and in this case, not desired. Transport properties of microstructures generated using the 2D and 3D SEs, including permeability, relative diffusivity, tortuosity factor, and tortuosity, are computed and given in Table 3. The transport properties of both structures in the in-
plane direction are close. However, the microstructure generated using the 3D SE has higher through-plane permeability compared to the microstructure generated using a 2D SE. The average chord lengths in the through-plane and inplane directions for the microstructures generated by the 2D and 3D SEs are 20 mm/46 mm and 26 mm/40 mm, respectively. The average chord lengths of both microstructures are close in the in-plane direction compared to the through-plane direction (13% difference in the in-plane direction, while 30% variation for the through-plane direction). Thus, the transport properties of both microstructures are close in the in-plane direction. Due to the larger chord lengths in the throughplane direction, the microstructure generated by 3D SE exerts lower resistance to the fluid resulting in larger throughplane permeability. The through-plane permeability of the microstructures generated using a 2D SE is closer to the experimentally measured data for Toray 060, which is 5.2 1012 m2 (23% difference for microstructure generated by 2D SE versus 35% differences for microstructure generated by 3D SE). The relative diffusivities in the in-plane and through-plane directions of the microstructure generated by 3D SE are close compared to the 2D microstructure and the experimental data. This may be due to topological similarities of the generated microstructure in the in-plane and through-plane directions when using a 3D cubical SE, as shown in Fig. 9(b). Also, the microstructure generated by 3D SE has fewer tortuous pathways in the through-plane direction resulting in the higher through-plane relative diffusivity compared to the 2D microstructure (25% variation form experimental data compared to 9% variation for microstructure generated by 2D SE). Negligible changes are observed in tortuosity for both microstructures. Thus, it is concluded that the tortuosity is not sensitive to topological changes in the porous medium. Since the microstructure generated by 2D SE resembles the semi-layered structure of GDL and results in more accurate
Fig. 9 e (a) Pore size (b) shortest path distributions for Toray 060 GDL modeled with cube and disk structure elements.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
9383
Table 3 e Transport properties of Toray 060 GDL modeled with different structure elements. SE type
Direction
Disk, 2D
Through-plane In-plane Through-plane In-plane
Cube, 3D
Permeability (m2) 6.8 1.19 8.1 1.14
1012 1011 1012 1011
prediction for the transport properties, the 2D disk SE is selected as a better filler for modeling the impregnated binder in the GDL.
Effect of binder volume fraction on the properties of the GDL To study the effects of binder quantity on the morphological and transport properties of composite structures, the volume fraction of binder in the GDL is changed from 20% to 47% while keeping the solid volume fraction of the fibrous bundle constant and equal to 0.16. The relationship between the volume fraction of binder and porosity of the composite structure is given in Table 4. As shown in Fig. 10(a), utilizing a higher amount of binder in the microstructure results in smaller pores. To study the inhomogeneity of the composite structures under various binder loading, the chord length distributions for the principal directions are found and the change in the average value of chord lengths are plotted in Fig. 10(b). As shown, the average chord length in the in-plane direction slightly increases with the increasing amount of binder, while the through-plane direction chord lengths sharply decrease. Due
Table 4 e Porosity of the composite structure and its binder volume fraction. Volume fraction of binder 20% 27% 36% 43% 47%
Volume fraction of fiber skeleton
Porosity of composite structure
0.16 0.16 0.16 0.16 0.16
0.80 0.78 0.75 0.72 0.70
Relative diffusivity
Tortuosity factor
Tortuosity
0.42 0.67 0.51 0.61
1.85 1.17 1.53 1.28
1.06 1.19 1.06 1.18
to the semi-layered structure of the GDL and high wettability of the binder, the small pores that exist between the fibers in the in-plane direction are filled. This results in an increase of chord lengths in the in-plane direction and a decrease of chord lengths in the through-plane direction. Although the reduction of the free space in the porous microstructure with the increase of binder loading can be determined from the pore size distribution, the directionality of the topology of the porous structure cannot. Thus, the pore size distribution provides little knowledge regarding the complex and inhomogeneous microstructure of the GDL. The presented results indicate the importance of conducting directional morphological analysis for composite fibrous microstructures. The dimensionless permeability (defined as K/r2, where r is the fiber radius) and relative diffusivity of the microstructures are calculated and graphed in Fig. 11. In addition, the numerically predicted dimensionless permeability in the through-plane direction is compared to the analytical results found by Spielman and Goren [44]. As shown in Fig. 11, the numerical predictions and the analytical results of permeability are in good agreement. Furthermore, as the binder volume fraction increases, both properties decrease, but larger changes are observed in the through-plane direction. By increasing the binder volume fraction, permeability decreases by 37% and 67%, and relative diffusivity reduces by 22% and 69%, in the in-plane and through-plane directions respectively. It is shown that although increasing binder volume fraction reduces both properties, the permeability in the inplane direction is more sensitive to the volume of binder compared to the relative diffusivity. The variations of tortuosity and tortuosity factor versus the amount of binder in the GDLs are shown in Fig. 12(a) and (b). The tortuosity factor in the through-plane direction depends heavily on the volume fraction of binder. This is due to the morphological properties of the GDLs. As depicted in Fig. 10(b),
Fig. 10 e Effect of volume fraction of binder in the GDL’s structure on (a) Pore size distribution, (b) The average chord length.
9384
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
Fig. 11 e Effect of volume fraction of binder in the GDL’s structure on permeability and relative diffusivity.
the average chord length in the in-plane direction changes by 17%, while in through-plane direction it changes by 55%. Because the tortuosity factor is a function of relative diffusivity (Drel) and the free space available in the porous structure, it is sensitive to the volume fraction of binder in GDL. However, the change in tortuosity versus increasing the volume fraction of binder is relatively insignificant.
Conclusion In this paper, a systematic approach has been presented to model and design 3D structures of composite fibrous porous media. The methodology includes fiber skeleton and binder impregnation modeling, and integrated characterization of the morphological and transport properties. The feasibility of the proposed methodology is verified by comparing the results to an existing porous microstructure used in PEMFCs, Toray 060 GDLs. It has been demonstrated that the PS model and binder impregnation model can be utilized to construct 3D geometries that resemble real or physical microstructures, e.g., GDLs, with good accuracy. In this work, the rod PS model was used to represent the fibrous semi-layered microstructure. To
construct a composite resembling a GDL, a 2D disk SE was selected for the impregnation model to incorporate the binder. Using the constructed microstructure, the interplay between the morphological and transport properties of the composite fibrous microstructures was studied by building the GDL models with various binder volume fractions. It was shown that the volume fraction of binder in the microstructure not only changes its morphological characteristic, but also the transport properties of the GDL. However, larger changes are observed in the through-plane direction compared to the inplane direction. It is found that binder volume fraction has a larger impact on permeability, relative diffusivity, and tortuosity factor than on the tortuosity. The high sensitivity of these properties to the volume fraction of binder can be attributed to the large reduction of the chord length in the through-plane direction. A major advantage of the proposed modeling approach for composite fibrous porous media over the microscopy approaches is that the ability to construct microstructures with tailored properties, because the fiber size and density, the composition and orientation, the shape and size of the binder, etc. can be altered. Thus, this framework is essential for designing and engineering such complex materials.
Fig. 12 e Effect of volume fraction of binder in the GDL’s structure on, (a) tortuosity, (b) tortuosity factor.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
Acknowledgments [19]
This work is supported by financial support of National Science Foundation under the grant numbers CMMI-0928093 and CMMI-1001040.
references
[20]
[21]
[22] [1] Yang F, Murugan R, Wang S, Ramakrishna S. Electrospinning of nano/micro scale poly(l-lactic acid) aligned fibers and their potential in neural tissue engineering. Biomaterials 2005;26(15):2603e10. [2] Foitzik RC, Kaynak A, Pfeffer FM. Conductive poly(a,u-bis(3pyrrolyl)alkanes)-coated wool fabrics. Synth Met 2007;157(13e15):534e9. [3] Ding X, Didari S, Fuller TF, Harris TA. Membrane electrode assembly fabrication process for directly coating catalyzed gas diffusion layers. J Electrochem Soc 2012;159(6):B746e53. [4] Didari S, Harris TAL, Huang W, Tessier SM, Wang Y. Feasibility of periodic surface models to develop gas diffusion layers: a gas permeability study. Int J Hydrogen Energy 2012;37(19):14427e38. [5] Dobrego KV, Shmelev ES, Koznacheev IA, Suvorov AV. Water purification of organic inclusions by the method of combustion within an inert porous media. Int J Heat Mass Transf 2010;53(11e12):2484e90. [6] Wang Y, Tang Y, Dong A, Wang X, Ren N, Shan W, et al. Selfsupporting porous zeolite membranes with sponge-like architecture and zeolitic microtubes. Adv Mater 2002;14(13e14):994e7. [7] Vielstich W, Gasteiger HA, Lamm A. Handbook of fuel cells: fuel cell technology and applications, vol. 4. Weinheim: Wiley-VCH Verlag GmbH; 2009. [8] Hoogers G. Fuel cell technology handbook. CRC; 2002. [9] Jordan L, Shukla A, Behrsing T, Avery N, Muddle B, Forsyth M. Diffusion layer parameters influencing optimal fuel cell performance. J Power Sources 2000;86(1):250e4. [10] de Frank Bruijn A, Janssen GJ. PEM Fuel cell materials: costs, performance and durability. Fuel cells. Springer; 2013. pp. 249e303. [11] Mortazavi M, Tajiri K. Effect of the PTFE content in the gas diffusion layer on water transport in polymer electrolyte fuel cells (PEFCs). J Power Sources 2014;245:236e44. [12] Nguyen T, Lin G, Ohn H, Wang X. Measurement of capillary pressure property of gas diffusion media used in proton exchange membrane fuel cells. Electrochem Solid-State Lett 2008;11(8):B127e31. [13] Advani SG, Sozer EM. Process modeling in composites manufacturing. CRC; 2002. [14] Mukherjee PP, Kang Q, Wang C-Y. Pore-scale modeling of two-phase transport in polymer electrolyte fuel cells e progress and perspective. Energy Environ Sci 2010;4(2):346e69. [15] Jaganathan S, Vahedi Tafreshi H, Pourdeyhimi B. A realistic approach for modeling permeability of fibrous media: 3-D imaging coupled with CFD simulation. Chem Eng Sci 2008;63(1):244e52. [16] Kinney JH, Nichols MC. X-ray tomographic microscopy (XTM) using synchrotron radiation. Annu Rev Mater Sci 1992;22(1):121e52. [17] Fredrich J, Menendez B, Wong T. Imaging the pore structure of geomaterials. Science (New York, NY) 1995;268(5208):276. [18] Lee J, Hinebaugh J, Bazylak A. Synchrotron X-ray radiographic investigations of liquid water transport
[23]
[24]
[25] [26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
9385
behaviour in a PEMFC with MPL-coated GDLs. J Power Sources 2013;227:123e30. Lee H-K, Park J-H, Kim D-Y, Lee T-H. A study on the characteristics of the diffusion layer thickness and porosity of the PEMFC. J Power Sources 2004;131(1):200e6. Fishman Z, Bazylak A. Heterogeneous through-plane distributions of tortuosity, effective diffusivity, and permeability for PEMFC GDLs. J Electrochem Soc 2011;158(2):B247e52. Fishman Z, Hinebaugh J, Bazylak A. Microscale tomography investigations of heterogeneous porosity distributions of PEMFC GDLs. J Electrochem Soc 2010;157(11):B1643e50. Becker J, Flu¨ckiger R, Reum M, Bu¨chi FN, Marone F, Stampanoni M. Determination of material properties of gas diffusion layers: experiments and simulations using phase contrast tomographic microscopy. J Electrochem Soc 2009;156(10):B1175e81. Becker J, Schulz V, Wiegmann A. Numerical determination of two-phase material parameters of a gas diffusion layer using tomography images. J Fuel Cell Sci Technol 2008;5(2):21006. Thiedmann R, Hartnig C, Manke I, Schmidt V, Lehnert W. Local structural characteristics of pore space in GDLs of PEM fuel cells based on geometric 3D graphs. J Electrochem Soc 2009;156(11):B1339e47. Wang Y. Periodic surface modeling for computer aided nano design. Computer-Aided Des 2007;39(3):179e89. Schulz VP, Becker J, Wiegmann A, Mukherjee PP, Wang C-Y. Modeling of two-phase behavior in the gas diffusion medium of PEFCs via full morphology approach. J Electrochem Soc 2007;154(4):B419e26. Gaiselmann G, Thiedmann R, Manke I, Lehnert W, Schmidt V. Stochastic 3D modeling of fiber-based materials. Comput Mater Sci 2012;59:75e86. Wang X, Zhang K, Yang Y, Wang L, Zhou Z, Zhu M, et al. Development of hydrophilic barrier layer on nanofibrous substrate as composite membrane via a facile route. J Membr Sci 2010;356(1e2):110e6. Tamayol A, McGregor F, Bahrami M. Single phase throughplane permeability of carbon paper gas diffusion layers. J Power Sources 2012;204:94e9. Daino MM, Kandlikar SG. 3D phase-differentiated GDL microstructure generation with binder and PTFE distributions. Int J Hydrogen Energy 2012;37(6):5180e9. Thiedmann R, Manke I, Lehnert W, Schmidt V. Random geometric graphs for modelling the pore space of fibre-based materials. J Mater Sci 2011;46(24):7745e59. Gostick JT, Fowler MW, Pritzker MD, Ioannidis MA, Behra LM. In-plane and through-plane gas permeability of carbon fiber electrode backing layers. J Power Sources 2006;162(1):228e38. Thiedmann R, Fleischer F, Hartnig C, Lehnert W, Schmidt V. Stochastic 3D modeling of the GDL structure in PEMFCs based on thin section detection. J Electrochem Soc 2008;155(4):B391e9. Gostick JT, Ioannidis MA, Fowler MW, Pritzker MD. Pore network modeling of fibrous gas diffusion layers for polymer electrolyte membrane fuel cells. J Power Sources 2007;173(1):277e90. Inoue G, Yoshimoto T, Matsukuma Y, Minemoto M. Development of simulated gas diffusion layer of polymer electrolyte fuel cells and evaluation of its structure. J Power Sources 2008;175(1):145e58. C¸ec¸en A, Wargo E, Hanna A, Turner D, Kalidindi S, Kumbur E. 3-D microstructure analysis of fuel cell materials: spatial distributions of tortuosity, void size and diffusivity. J Electrochem Soc 2012;159(3):B299e307. Dijkstra EW. A note on two problems in connexion with graphs. Numer Math 1959;1(1):269e71.
9386
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
[38] Gonzalez RC, Woods RE, Eddins SL. Digital image processing using MATLAB. Pearson Education India; 2004. [39] Image processing toolbox documentation. Available from: www.mathworks.com/help/images/index.html; 2010. [40] ANSYS FLUENT user’s guide, release 13.0. ANSYS; 2010. [41] Becker J, Wieser C, Fell S, Steiner K. A multi-scale approach to material modeling of fuel cell diffusion media. Int J Heat Mass Transf 2011;54(7):1360e8. [42] Flu¨ckiger R, Freunberger SA, Kramer D, Wokaun A, Scherer GG, Bu¨chi FN. Anisotropic, effective diffusivity of
porous gas diffusion layer materials for PEFC. Electrochim Acta 2008;54(2):551e9. [43] Feser J, Prasad A, Advani S. Experimental characterization of in-plane permeability of gas diffusion layers. J Power Sources 2006;162(2):1226e31. [44] Spielman L, Goren SL. Model for predicting pressure drop and filtration efficiency in fibrous media. Environ Sci Technol 1968;2(4):279e87.
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Modeling of composite fibrous porous diffusion media Sima Didari a, Arash Asadi b, Yan Wang a, Tequila A.L. Harris a,* a
George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, 813 Ferst Dr., Atlanta, GA 30332, USA b School of Mathematics, Georgia Institute of Technology, 686 Cherry Street, Atlanta, GA 30332, USA
article info
abstract
Article history:
To engineer the desired properties of fibrous porous media, a parametric modeling
Received 7 October 2013
approach is needed to support the rational design of the materials before the fabrication. In
Received in revised form
this study, we propose a methodology that enables the accurate representation of three-
27 March 2014
dimensional (3D) microstructures of fibrous porous media and prediction of their trans-
Accepted 2 April 2014
port properties. Toray TGP-H-060 gas diffusion layer (GDL) is selected as an example to
Available online 29 April 2014
demonstrate the feasibility of the suggested design methodology. The detailed microstructure of the GDL with the inclusion of locally distributed binder is constructed using an
Keywords:
extended periodic surface (PS) modeling technique. A 3D morphological approach is taken
Modeling
to create the binder distribution within the fibrous microstructure. Transport properties
Composite
including permeability, relative diffusivity, and tortuosity and local structure characteris-
Gas diffusion layer
tics of the generated microstructure, under different binder loading are calculated. It is
Porous media
shown that the detailed model of the fiber-binder composite has a strong influence on the
Transport properties
predicted properties.
Design
Copyright ª 2014, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Introduction Fibrous porous media have been used for many decades in a plethora of industries such as bioengineering, electronics, energy and environmental processes [1e6]. As an example, gas diffusion layers (GDL) are used in polymer electrolyte membrane fuel cells (PEMFC), a burgeoning renewable energy resource. GDLs serve as electrodes in PEMFCs, which permeate gases and help remove water droplets that are formed during fuel cell operation. Furthermore, these carbon based GDLs have other properties such as high porosity, strength-to-weight ratio, and electrical and thermal
conductivities that make them attractive for use in PEMFCs and other areas [7]. It has been shown that the performance of PEMFCs depends upon the morphological, transport, and mechanical properties of the GDLs [8e12]. However, no systematic and parametric modeling approach that help predict and optimize the key characteristics of porous materials such as GDLs for the purpose of materials design exists to date [13]. Rather, studies have been conducted to determine the morphological and transport properties of the fibrous diffusion media that have been fabricated [14]. In these studies, the microstructure of the porous media was generated by image-based reconstruction or geometrical construction approaches. In the
* Corresponding author. E-mail addresses: [email protected], [email protected] (T.A.L. Harris). http://dx.doi.org/10.1016/j.ijhydene.2014.04.011 0360-3199/Copyright ª 2014, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
9376
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
image-based approaches, 2D images or 3D voxel representation of the internal pore structures are obtained by scanning electron microscopy (SEM), synchrotron-based tomography, confocal microscopy, or computed microtomography (mCT) [15e24]. Although imaging techniques accurately capture the geometry of porous structures, the model of the microstructures reconstructed from the images cannot be modified. Thus, imaged-based approaches are not suitable for designing microstructures of porous media. The geometrical construction approaches that utilize conventional computer aided design (CAD) tools are not efficient in modeling complex structures. Because in these approaches, a large number of discrete control points are required to accurately represent the topology of the porous media. Conversely, implicit surface modeling approaches that construct the microstructures based on closed-form mathematical functions are much more promising [25]. The geometrical construction of fibrous diffusion media has been primarily based on statistical methods. Fibers are assumed to have simple shapes such as long cylinders and randomly distributed in a representative volume element (RVE). In the work of Schulz et al. [26], the fiber skeleton of GDLs was modeled as a collection of intersecting cylinders oriented along randomly distributed axes. The directional distribution of the fibers was assumed as normal or determined from the SEM images. Thiedmann et al. [24] constructed GDL microstructures by implementing a Poisson planar line tessellations (PLT) process to represent the fibers’ axes. Gaiselmann et al. [27] represented the pore space of fibrous media with random geometric graphs. The parameters required for modeling were specified such that the distributions of vertex degrees and edge lengths conform with X-ray synchrotron tomography images. GDLs have also been represented with simplified geometries such as arrays of parallel and perpendicular cylinders or channels [28,29]. Fibrous porous media are often more complex than a collection of fibers. They also include some binding material(s). To model the binder in the microstructure of the GDL, the pore space of the generated fiber skeleton is filled using various approaches such as stochastic methods with predefined probability functions or morphological operations [24,30]. To quantitatively characterize the pore space of a porous medium and to predict its bulk material properties, the morphological properties of the generated microstructure (e.g., pore size and shortest path distributions) have been also studied. The pore size has been calculated as the spherical distance of the points inside the void space to the surrounding solid boundaries [31e33]. However, in some studies, the pores were defined as 2D closed areas and the pore size is defined as the diameter of a circle with an area equal to the closed area [34,35]. To determine the inhomogeneity and directionality of the distribution of pores in the microstructure, the chord length distribution of a porous medium was also found by evaluating angularly resolved lineal path functions [36]. The distribution of the shortest path within a porous microstructure has been found from the voxel mesh or the graph representation of the porous domain [24,36], using Dijkstra’s algorithm [37]. The existing numerical schemes offer useful platforms and tools for the evaluation
of some morphological properties. However, they are computationally expensive. In the current study, these schemes are modified for higher accuracy and lower computational time. Transport properties of microstructures can be predicted using different approaches such as pore network (PN), lattice Boltzmann models (LBM) and computational fluid dynamics (CFD). The effects of change in porosity due to the compression forces and fiber size and spatial orientation have been investigated through permeability, tortuosity, diffusion and relative permeability of fibrous porous media [14]. The CFD approach is more time efficient than LBM and can simulate larger domains. The PN modeling approach uses a simplified representation of the pore structure and is efficient for large domains. However, it requires the information of length correlations that need to be extracted from an existing structure. The process itself is time consuming. The objective of this work is to introduce a platform for designing composite porous media, with the characteristics that are similar to those of GDLs. A new methodology is developed to generate the fiber skeleton of composite porous diffusion media. In this study, a periodic surface (PS) modeling technique developed by Wang [25] is used to generate the fibrous component of the microstructures, and the binder distribution is modeled by the Minkowski sum of fibers and a structure element (SE). The microstructure generation process is directly integrated with the morphology and transport characterization analyses.
Material Toray TGP-H-060 carbon paper GDL, shown in Fig. 1, is chosen as the fibrous composite microstructure to examine the feasibility of the microstructure generation method. Toray TGP 060 is made of a matrix of carbonaceous binder that is a thermoset resin and graphitized rigid carbon fibers [22]. Typical parameters of the GDL include, a fiber diameter of 7 mm, a thickness of 200 mm, porosity of 0.78 and 0.27 solid fraction of carbonaceous thermoset binder [30].
Implicit surface modeling of fibrous structures An implicit surface S ¼ {r ˛ R3:j(r) ¼ 0} is defined as a set of points that are the roots of a given function j:R3 / R with realvalued coordinates. Recently a periodic surface (PS) model was developed by Wang [25] to implicitly model structures with complex topology, which is define as jðrÞ ¼
L X M X
mlm cos 2pkl PTm $r ¼ 0
(1)
l¼1 m¼1
where j(r) is the implicit function, mlm is periodic moment, kl is the scale parameter, Pm ¼ [am, bm, cm, am]T is a basis vector for the major phase, [am, bm, cm]T represents the normal vector of the basis plane, am is the phase of the plane, r ¼ [x, y, z, 1]T is the location vector in homogenous coordinates. The PS model of a rod surface with the periodicity of 1 along the z-axis direction in an RVE is:
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
9377
Fig. 1 e SEM images of Toray GDL, (a) Front view, (b) Cross sectional view.
jðx; y; zÞ ¼ 4 cosð2pxÞ þ 4 cosð2pyÞ þ 3 ¼ j0
(2)
where j0 is the adjustable isovalue. To construct many fibers in a section of a fibrous medium, the PS model of each fiber has to be constructed and inserted in an RVE. The whole microstructure with n fibers can be generated in an RVE with the dimensionless lengths using the union operation jW ðrÞ ¼ minðj1 ðrÞ; j2 ðrÞ; .; jn ðrÞÞ
(3)
over the implicit functions of all individual fibers in the domain. Due to the scalable nature of PS model, the generated RVE can be readily scaled to the desired size.
Geometrical modeling of Toray TGP-H-060 GDL The composite structure of the GDL is constructed in two steps. First, the fiber skeleton is generated using the PS model. Then the binder is integrated into the skeleton of the microstructure. To generate the structure of fiber skeleton, the diameter of the fibers (D) and the solid volume fraction or the porosity of fiber skeleton have to be known. The diameter of the fibers is determined by the isovalue (j0) of the rod PS model. As the isovalue changes, the diameter of the rod surface changes. The diameter of the fiber can be approximated by the distance between the positions with the maximum ycoordinate Ymax and minimum y-coordinate Ymin in the cross section of the rod as D ¼ Ymax Ymin [4]. Since the period is set to be 1, the cross section of the rod in the xey plane of the RVE is symmetric about x ¼ 1/2. Both Ymax and Ymin are located at x ¼ 1/2. Therefore, 1 þ 4 cosð2pyÞ þ 3 ¼ j0 4 cos 2p 2
(4)
The simplified form of Equation (4) is cos(2py) ¼ (1 þ j0)/4, which has two solutions for y ˛ [0,1], which are Ymax and Ymin. Thus, the diameter can be found as: D ¼ Ymax Ymin ¼ 1
arccosðð1 þ j0 Þ=4Þ p
(5)
The porosity of the GDL is related to the number of fibers in the RVE. To decide the orientation of a fiber in the RVE, two
angles are needed. One is the in-plane angle in the xey plane, q, and the other is the side angle in the zex plane, 4. These two angles determine the orientation of major phase vector Pm of each of the implicit functions. In a semi-layered structure, the ranges in which q and 4 can change are set to be [0, p] and [p/ 36, p/36] respectively. The positions and orientations of the fibers follow the normal distributions in the RVE. The randomly generated fibers may intersect with each other in the RVE. The overall intersection volume of the generated structure is quantified by U ¼ jX ðrÞ ¼ min max ji ðrÞ; jj ðrÞ 1ijn
(6)
To avoid the unrealistic phenomenon of intersection, U in Equation (6) is minimized as min U (qi, 4i, xi, yi, zi)i¼1:n where n is the number of the fibers, qi, 4i, xi, yi, zi represent the orientation and position of the ith fiber. The optimization problem is solved by the genetic algorithm (GA) and implement in MATLAB. The optimization procedure is described schematically in Fig. 2. Assuming negligible deformation for carbon fibers due to their rigidity, the microstructure of the fiber skeleton is generated as shown in Fig. 3(a).
Binder impregnation in Toray TGP-H-060 GDL Because of the low viscosity and high wettability of the resin and the rapid impregnation process, in this work only the final composite structure is modeled. The resin-fiber bond that forms the composite nature of a GDL is modeled in two steps. In the first step, a Minkowski sum between the fiber skeleton (denoted by A) and an SE that represents the resin is performed, which is defined as B ¼ A 4 SE ¼ {a þ seja ˛ A, se ˛ SE}. The result is the thicken fiber bundle as shown in Fig. 3(b). In the second step, an erosion operation is conducted on set B. The erosion, defined as B . (SE), is the Minkowski subtraction (shown by (.)) of B with the reflected set of SE. As shown in Fig. 3(c), the erosion process leads to a reduction of set B [38,39], whereby the pores that are smaller than the SE are filled. The PS model and binder impregnation model are developed in MATLAB software and can be easily coupled with
9378
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
of the binder in the generated composite structure is calculated. If the target weight fraction is not met, the filling procedure is conducted with a larger SE. Using a disk shape SE, the microstructure of Toray TGP-H-060 carbon paper is generated, as shown in Fig. 4, in comparison with scanning electron microscope (SEM) images. As depicted in the figure, the structure has a semi-layer-by-layer configuration that resembles an actual GDL.
Morphology characterization Calculating porosity
Fig. 2 e Block diagram of generating the fibers skeleton.
microstructure characterization. The fiber skeleton generated using the PS model is converted into a 3D binary voxel representation that serves as an input for the binder impregnation model. The Minkowski summation and erosion procedures are implemented with the MATLAB image processing toolbox [39]. In the binary voxel representation of the fibrous domain coded as a binary matrix, the elements of the matrix that correspond to the points inside the solid phase are set to be 1 and the elements in the void or pore space are set to be 0. Given the PS model defined in the RVE, the value of the implicit function at each node in the discretized grid space is calculated. This resulting value j is compared to the given isovalue j0. If j j0, the node is in the solid phase, otherwise it is located in the pore space. Then, a value of 1 or 0 is assigned to the node. After generating the fiber skeleton using the smallest size for the SE, i.e., the mesh size, the two-step resin impregnation procedure described above is applied and the weight fraction
The binary matrix representation of the porous structure is used to calculate the morphological properties of interest. For a composite porous structure, three separate binary matrices are used to calculate the porosity of the whole structure, porosity of the fibrous skeleton, and the volume fraction of the binder. To find the solid volume fraction of each phase used in the composite structure, i.e., the resin and fibers, the number of elements in the binary matrix that have the value of 1 are counted. Then the number is divided by the matrix size. The porosity of a composite is determined by counting “0” elements in the binary matrix, and the number is divided by the matrix size.
Determining pore size distribution The pore size distribution is determined by finding a set of the largest fitting spheres in the pore space. Given that each node in the pore space, i.e., 0’s of the binary matrix, can be a potential pore center, the implemented search algorithm starts with placing a sphere with the smallest radius (mesh size) at each of these nodes. The radius is increased gradually until that the sphere fits within the pore space. A sphere is considered as fit if it has at least three points of contacts with the surrounding solid surfaces and the three points are not in the same direction. As such, multiple spheres can be contained inside a pore space. Those spheres with centers located inside a larger sphere are deleted. The remaining spheres then form a set of mutually overlapped pores, where none of their centers are located inside another. The radii of the spheres are the pore sizes.
Fig. 3 e (a) Fibrous bundle, (b) Minkowski summation with structure element, (c) Minkowski erosion with the same structure element.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
9379
Fig. 4 e (a) Isometric view of constructed Toray TGP-H-060 GDL, (b) and (c) cross sectional views, virtual GDL (left side) and SEM images (right side).
To compute the pore size distribution that is closer to the one found experimentally, the closed and dead end pores are discarded and not counted in the statistics. The connectivity of the pores is found by determining if there exists a pathway that connects one pore center to the neighboring pore centers. The resulting pore size distribution is shown in Fig. 5. The average pore size for the model generated for Toray TGP-H060 GDLs is 21 mm, which is close to the average pore size of 18 mm measured by Flu¨ckiger et al. [12].
Determining chord length distribution The pore size distribution provides an insight on the void space within the porous structure, but it cannot provide any information regarding the directionality of the pore space. Hence, the chord length (Lch) defined as the distance between the internal solid boundaries in the pore space is used to represent the degree of anisotropy of the microstructure in a desired direction. The chord length is a measure of the available free space in which the particles can move inside a porous structure. Thus, using the chord length distribution, the dominant mode of diffusion inside porous media can also be found. To find the chord length distribution, the cross-
Fig. 5 e Numerical pore size distribution of virtual Toray TGP-H-060 GDL.
sectional 2D binary matrices of the porous medium are extracted from the 3D binary matrix. The chord length is the number of consecutive 0’s in a principal direction multiplied by the size of discrete elements. The chord length distributions for the through-plane and in-plane directions of the model for Toray 060 GDL are shown in Fig. 6, where the chord lengths are normalized as the ratio to the thickness of the GDL. As depicted in the figure, in the through-plane direction the chord lengths are smaller due to the semi-layered structure of the GDL, while in the in-plane direction the chord lengths have a wider distribution.
Determining the shortest path distribution To find the distribution of the shortest path inside the porous structure, Dijkstra’s algorithm [37] is implemented and modified to calculate all of the shortest paths starting from each pore. To implement Dijkstra’s algorithm, a graph representation of the porous domain is extracted. To construct the graph representation of the porous medium, first it is assumed that there is no solid phase in the RVE. A graph G with vertices corresponding to the 3D grid points with the size of nx ny nz is generated to show the connectivity of nodes inside the RVE. Any point in the 3D grid with the position index (i, j, k), where 0 i nx, 0 j ny, and 0 k nz, is a vertex of G. To construct the edges of G, the connectivity of the vertices is investigated. Two vertices (i, j, k) and (a, b, c) are assumed to be connected if and only if ji-aj þ jj-bj þ jk-cj ¼ 1. Graph G is then modified to account for the effect of the solid phase in the RVE. The solid phase acts as an edge eliminator, vertex eliminator, or both, to determine whether a vertex belongs to the void space or the solid phase. Then the vertex and the associated edges in G are deleted. The graph G is a weighted graph where the weight of each edge is 1. It is assumed that all the vertices in the inlet and outlet are in the pore space and that the sides of the RVE are in the solid space. For each vertex at the inlet, the modified version of Dijkstra’s algorithm is run and the length of the shortest path is found. The shortest path from each vertex at the inlet to each vertex
9380
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
Fig. 6 e Chord length distribution for virtual Toray TGP-H-060 GDL, (a) through-plane direction, (b) in-plane direction.
at the outlet is found. The algorithm is implemented in Cþþ and the MEX format is used for data exchange with MATLAB. The distributions of the shortest paths inside the models for Toray 060 GDL along the through-plane and in-plane directions are shown in Fig. 7. The computed average length of the shortest path in the through-plane direction is 1.64, which agrees well with the experimentally measured value of 1.51 [24]. As expected, due to the semi-layer structure of the GDL in the through-plane direction, the GDL has longer pathways or higher tortuosity.
Transport property characterization A high-quality volumetric mesh of the computational domain is needed to determine the transport properties of interest, i.e., permeability, tortuosity, and diffusivity. For geometries that can be represented by CAD software, common mesh generation packages can be used. However, for the complex geometries of the composite structures presented in this work, conventional mesh generation methods are not applicable. Thus, a voxel mesh representation of the computational domain is used.
conservation equations that govern the flow behavior through the porous media are V$u ¼ 0 and ru$Vu ¼ VP þ Vbs
(8)
respectively, where u is the velocity, r is density, P is the pressure field, and b t is the stress tensor. By importing the volumetric mesh into ANSYS FLUENT and solving the mass and momentum equations given by Equations (7) and (8), the permeability of the microstructures is found from Darcy’s law: K V ¼ VP h
(9)
where V is the superficial velocity (i.e., volumetric flow rate per unit cross-sectional area), h is fluid viscosity, and VP is the pressure drop per unit length through the porous material. A mesh sensitivity study was conducted to find appropriate mesh sizes. It is also found that a computational domain of 200 mm 200 mm 200 mm is sufficient to accurately compute the transport properties. To determine the permeability in a direction, the inlet and outlet boundary conditions are set for the specific flow direction, as shown in Fig. 8. At the inlet a
Volumetric mesh representation of the computational domain Using the MATLAB image processing toolbox, the crosssectional images of the porous structure are generated based on 2D arrays, which are extracted from the 3D binary matrix, as illustrated in Fig. 8. These images are saved as a sequence of images and are inserted into Geo-Dict software to generate a volumetric voxel mesh. The whole procedure of mesh generation, including reading the 3D binary matrix, creating the 2D cross-sectional images and importing them into the Geo-Dict software, is fully automated using custom-built MATLAB scripts. This procedure is an efficient and flexible approach for generating the voxel mesh representation of the domain; hence, the microstructures can be easily altered or changed.
Determining permeability of GDLs The flow field inside the pore space is solved to find the permeability (K). The steady state momentum and mass
(7)
Fig. 7 e Distribution of the shortest path inside the pore space of TGP-H-060 GDL.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
9381
Determining the relative diffusivity and tortuosity factor in porous structures
Fig. 8 e Schematic illustrating the generation of a volumetric mesh.
pressure boundary is defined, while at the outlet a pressure outlet boundary condition is imposed for a gauge pressure of zero. The no slip boundary condition is applied to the internal solid fluid interfaces, which are the outer surfaces of the fibers. The momentum equation is discretized by the second order upwind scheme with the SIMPLE algorithm for coupling pressure and velocity [40]. The under relaxation factors of 0.25 and 0.55 for pressure and momentum are used, respectively. After the simulations have converged, the superficial velocity is calculated by volume averaging the velocity fields in the GDL. Since the GDL thickness is known, the pressure gradient is then calculated. Applying the calculated values for the superficial velocity and the pressure gradient to Equation (9), the permeability is computed.
Determining the tortuosity in GDLs
Diffusivity in porous media (Dbulk) is defined as Dbulk ¼ Drel D0, where Drel is relative diffusivity, which represent the reduction of the diffusivity in the porous media due to the presence of solid phase and D0 is the binary diffusion coefficient. D0 ¼ 1/ 3ly, where l is the mean free path of the gas that passes through the porous media and y is the thermal velocity. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y ¼ 8kB T=pm, where kB is the Boltzmann constant, T is temperature, and m is the molecular mass. To determine the mechanism by which diffusion takes place inside the pore space, the Knudson number (Kn) is calculated, which is the ratio of the mean free path of a specific fluid (l) to the available length scale inside the porous structure [41]. This length scale, which is set as the average chord length ðLch Þ, indicates the freedom of fluids to move inside the pore space. In the current study, the diffusivity of O2, H2, and H2O, which are the fuel and byproduct in PEMFCs, are computed. The mean free path for H2, O2, and H2O are presented in Table 1 [36]. Using l and the average value of the chord length found from the morphology characterization, the average Kn numbers are calculated, as given in Table 1. The Knudsen numbers are less than one, indicating that the pore diameters are large compared to the mean free path length of the fluid molecules. Thus, the diffusion in the pores is governed by the concentration Laplace equation given by DC ¼ 0
where C is the concentration of the gas in the pore space. To find Drel, Equation (11) is solved using Geo-Dict software by applying zero flux at the solid surfaces and a concentration gradient in the direction of the desired Drel. For a known diffusive flux (J), diffusivity is then calculated using Fick’s law J ¼ Dbulk VC
1
s¼N
PN i¼1
L
lðri Þ
(12)
Then, Drel can be directly computed and used to calculate the tortuosity factor (sf) by sf ¼
The other important transport property of porous media is tortuosity (s), which represents the complexity of the microscopic flow paths in the porous medium. Tortuosity is defined as the ratio of average length of the flow path to the characteristic length of the porous medium. The thickness of porous medium in the flow direction is usually assumed as the characteristic length. For inhomogeneous porous structures, depending on the direction of flow (in-plane or throughplane), tortuosity can be different. Tortuosity, in a desired direction, is calculated from
(11)
ε Drel
(13)
where ε is the porosity of the porous domain. The through-plane and in-plane permeability, relative diffusivity, tortuosity, and tortuosity factor of Toray 060 GDL are calculated, using the methods discussed, and are given in Table 2. There exists very good agreement between the experimentally measured data and the numerically computed data based on the microstructure(s) generated using the PS model and binder impregnation model developed. As expected, due to the semi-layered structure of the GDL, the transport properties in the through-plane directions are lower
(10)
where N is the number of the particles that are tracked, l(ri) is the length of the path line that belongs to the ith particle located at ri, and L is the GDL thickness. After convergence of the numerical simulation, for 2500 particles, the path line length is computed and the tortuosity of the GDL is calculated for both the through-plane and in-plane directions.
Table 1 e Average Knudson number in Toray 060 GDL at T [ 100 C and P [ 100 kPa. Fluid type H2O O2 H2
l (m)
Kn number 7
7.5 10 6.1 107 1.3 107
1.2 102 9.7 103 2.1 103
9382
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
Table 2 e Transport properties of Toray 060 GDL. Numerical Experimental
Direction
Permeability (m2)
Relative diffusivity
Tortuosity factor
Tortuosity
Through-plane In-plane Through-plane In-plane
6.8 1012 1.19 1011 5.2 1012 [22] 2.1 1011 [43]
0.42 0.67 0.38 [42] 0.57 [42]
1.85 1.17 2.05 [42] 1.36 [42]
1.06 1.19 NA NA
than the in-plane direction. However, larger discrepancy is observed in the through-plane direction compared to the inplane direction.
Effect of binder impregnation on the characteristics of GDL As mentioned in Binder Impregnation in Toray TGP-H-060 GDL Section, to model composite structures impregnated with binder, a two-step procedure is implemented. The procedure includes Minkowski summation and erosion of the fiber skeleton model by a filler SE. To determine the SE that can accurately represent the composite fibrous structure, a disk SE (2D) and a cube SE (3D) are used to construct the microstructure of Toray 060 GDL, and the morphological and transport properties of both microstructures are found. As shown in Fig. 9(a), the pore size distributions for the two microstructures are different. However, the pore size distribution corresponding to the 2D disk SE shows higher resemblance to the pore size distribution found experimentally [42]. Distributions of the shortest pathways in the through-plane and in-plane directions are found and presented in Fig. 9(b). As illustrated, the GDL generated using the disk SE has an inhomogeneous topology as expected, due to the semi-layered structure of the porous medium of interest. However, the distributions of the pathways in the through-plane and inplane directions of the microstructure generated using the 3D cube SE, are approximately equal. Using the 3D SE results in a uniform distribution of the binder in both directions, which may lead to microstructures that are not semi-layered, and in this case, not desired. Transport properties of microstructures generated using the 2D and 3D SEs, including permeability, relative diffusivity, tortuosity factor, and tortuosity, are computed and given in Table 3. The transport properties of both structures in the in-
plane direction are close. However, the microstructure generated using the 3D SE has higher through-plane permeability compared to the microstructure generated using a 2D SE. The average chord lengths in the through-plane and inplane directions for the microstructures generated by the 2D and 3D SEs are 20 mm/46 mm and 26 mm/40 mm, respectively. The average chord lengths of both microstructures are close in the in-plane direction compared to the through-plane direction (13% difference in the in-plane direction, while 30% variation for the through-plane direction). Thus, the transport properties of both microstructures are close in the in-plane direction. Due to the larger chord lengths in the throughplane direction, the microstructure generated by 3D SE exerts lower resistance to the fluid resulting in larger throughplane permeability. The through-plane permeability of the microstructures generated using a 2D SE is closer to the experimentally measured data for Toray 060, which is 5.2 1012 m2 (23% difference for microstructure generated by 2D SE versus 35% differences for microstructure generated by 3D SE). The relative diffusivities in the in-plane and through-plane directions of the microstructure generated by 3D SE are close compared to the 2D microstructure and the experimental data. This may be due to topological similarities of the generated microstructure in the in-plane and through-plane directions when using a 3D cubical SE, as shown in Fig. 9(b). Also, the microstructure generated by 3D SE has fewer tortuous pathways in the through-plane direction resulting in the higher through-plane relative diffusivity compared to the 2D microstructure (25% variation form experimental data compared to 9% variation for microstructure generated by 2D SE). Negligible changes are observed in tortuosity for both microstructures. Thus, it is concluded that the tortuosity is not sensitive to topological changes in the porous medium. Since the microstructure generated by 2D SE resembles the semi-layered structure of GDL and results in more accurate
Fig. 9 e (a) Pore size (b) shortest path distributions for Toray 060 GDL modeled with cube and disk structure elements.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
9383
Table 3 e Transport properties of Toray 060 GDL modeled with different structure elements. SE type
Direction
Disk, 2D
Through-plane In-plane Through-plane In-plane
Cube, 3D
Permeability (m2) 6.8 1.19 8.1 1.14
1012 1011 1012 1011
prediction for the transport properties, the 2D disk SE is selected as a better filler for modeling the impregnated binder in the GDL.
Effect of binder volume fraction on the properties of the GDL To study the effects of binder quantity on the morphological and transport properties of composite structures, the volume fraction of binder in the GDL is changed from 20% to 47% while keeping the solid volume fraction of the fibrous bundle constant and equal to 0.16. The relationship between the volume fraction of binder and porosity of the composite structure is given in Table 4. As shown in Fig. 10(a), utilizing a higher amount of binder in the microstructure results in smaller pores. To study the inhomogeneity of the composite structures under various binder loading, the chord length distributions for the principal directions are found and the change in the average value of chord lengths are plotted in Fig. 10(b). As shown, the average chord length in the in-plane direction slightly increases with the increasing amount of binder, while the through-plane direction chord lengths sharply decrease. Due
Table 4 e Porosity of the composite structure and its binder volume fraction. Volume fraction of binder 20% 27% 36% 43% 47%
Volume fraction of fiber skeleton
Porosity of composite structure
0.16 0.16 0.16 0.16 0.16
0.80 0.78 0.75 0.72 0.70
Relative diffusivity
Tortuosity factor
Tortuosity
0.42 0.67 0.51 0.61
1.85 1.17 1.53 1.28
1.06 1.19 1.06 1.18
to the semi-layered structure of the GDL and high wettability of the binder, the small pores that exist between the fibers in the in-plane direction are filled. This results in an increase of chord lengths in the in-plane direction and a decrease of chord lengths in the through-plane direction. Although the reduction of the free space in the porous microstructure with the increase of binder loading can be determined from the pore size distribution, the directionality of the topology of the porous structure cannot. Thus, the pore size distribution provides little knowledge regarding the complex and inhomogeneous microstructure of the GDL. The presented results indicate the importance of conducting directional morphological analysis for composite fibrous microstructures. The dimensionless permeability (defined as K/r2, where r is the fiber radius) and relative diffusivity of the microstructures are calculated and graphed in Fig. 11. In addition, the numerically predicted dimensionless permeability in the through-plane direction is compared to the analytical results found by Spielman and Goren [44]. As shown in Fig. 11, the numerical predictions and the analytical results of permeability are in good agreement. Furthermore, as the binder volume fraction increases, both properties decrease, but larger changes are observed in the through-plane direction. By increasing the binder volume fraction, permeability decreases by 37% and 67%, and relative diffusivity reduces by 22% and 69%, in the in-plane and through-plane directions respectively. It is shown that although increasing binder volume fraction reduces both properties, the permeability in the inplane direction is more sensitive to the volume of binder compared to the relative diffusivity. The variations of tortuosity and tortuosity factor versus the amount of binder in the GDLs are shown in Fig. 12(a) and (b). The tortuosity factor in the through-plane direction depends heavily on the volume fraction of binder. This is due to the morphological properties of the GDLs. As depicted in Fig. 10(b),
Fig. 10 e Effect of volume fraction of binder in the GDL’s structure on (a) Pore size distribution, (b) The average chord length.
9384
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
Fig. 11 e Effect of volume fraction of binder in the GDL’s structure on permeability and relative diffusivity.
the average chord length in the in-plane direction changes by 17%, while in through-plane direction it changes by 55%. Because the tortuosity factor is a function of relative diffusivity (Drel) and the free space available in the porous structure, it is sensitive to the volume fraction of binder in GDL. However, the change in tortuosity versus increasing the volume fraction of binder is relatively insignificant.
Conclusion In this paper, a systematic approach has been presented to model and design 3D structures of composite fibrous porous media. The methodology includes fiber skeleton and binder impregnation modeling, and integrated characterization of the morphological and transport properties. The feasibility of the proposed methodology is verified by comparing the results to an existing porous microstructure used in PEMFCs, Toray 060 GDLs. It has been demonstrated that the PS model and binder impregnation model can be utilized to construct 3D geometries that resemble real or physical microstructures, e.g., GDLs, with good accuracy. In this work, the rod PS model was used to represent the fibrous semi-layered microstructure. To
construct a composite resembling a GDL, a 2D disk SE was selected for the impregnation model to incorporate the binder. Using the constructed microstructure, the interplay between the morphological and transport properties of the composite fibrous microstructures was studied by building the GDL models with various binder volume fractions. It was shown that the volume fraction of binder in the microstructure not only changes its morphological characteristic, but also the transport properties of the GDL. However, larger changes are observed in the through-plane direction compared to the inplane direction. It is found that binder volume fraction has a larger impact on permeability, relative diffusivity, and tortuosity factor than on the tortuosity. The high sensitivity of these properties to the volume fraction of binder can be attributed to the large reduction of the chord length in the through-plane direction. A major advantage of the proposed modeling approach for composite fibrous porous media over the microscopy approaches is that the ability to construct microstructures with tailored properties, because the fiber size and density, the composition and orientation, the shape and size of the binder, etc. can be altered. Thus, this framework is essential for designing and engineering such complex materials.
Fig. 12 e Effect of volume fraction of binder in the GDL’s structure on, (a) tortuosity, (b) tortuosity factor.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
Acknowledgments [19]
This work is supported by financial support of National Science Foundation under the grant numbers CMMI-0928093 and CMMI-1001040.
references
[20]
[21]
[22] [1] Yang F, Murugan R, Wang S, Ramakrishna S. Electrospinning of nano/micro scale poly(l-lactic acid) aligned fibers and their potential in neural tissue engineering. Biomaterials 2005;26(15):2603e10. [2] Foitzik RC, Kaynak A, Pfeffer FM. Conductive poly(a,u-bis(3pyrrolyl)alkanes)-coated wool fabrics. Synth Met 2007;157(13e15):534e9. [3] Ding X, Didari S, Fuller TF, Harris TA. Membrane electrode assembly fabrication process for directly coating catalyzed gas diffusion layers. J Electrochem Soc 2012;159(6):B746e53. [4] Didari S, Harris TAL, Huang W, Tessier SM, Wang Y. Feasibility of periodic surface models to develop gas diffusion layers: a gas permeability study. Int J Hydrogen Energy 2012;37(19):14427e38. [5] Dobrego KV, Shmelev ES, Koznacheev IA, Suvorov AV. Water purification of organic inclusions by the method of combustion within an inert porous media. Int J Heat Mass Transf 2010;53(11e12):2484e90. [6] Wang Y, Tang Y, Dong A, Wang X, Ren N, Shan W, et al. Selfsupporting porous zeolite membranes with sponge-like architecture and zeolitic microtubes. Adv Mater 2002;14(13e14):994e7. [7] Vielstich W, Gasteiger HA, Lamm A. Handbook of fuel cells: fuel cell technology and applications, vol. 4. Weinheim: Wiley-VCH Verlag GmbH; 2009. [8] Hoogers G. Fuel cell technology handbook. CRC; 2002. [9] Jordan L, Shukla A, Behrsing T, Avery N, Muddle B, Forsyth M. Diffusion layer parameters influencing optimal fuel cell performance. J Power Sources 2000;86(1):250e4. [10] de Frank Bruijn A, Janssen GJ. PEM Fuel cell materials: costs, performance and durability. Fuel cells. Springer; 2013. pp. 249e303. [11] Mortazavi M, Tajiri K. Effect of the PTFE content in the gas diffusion layer on water transport in polymer electrolyte fuel cells (PEFCs). J Power Sources 2014;245:236e44. [12] Nguyen T, Lin G, Ohn H, Wang X. Measurement of capillary pressure property of gas diffusion media used in proton exchange membrane fuel cells. Electrochem Solid-State Lett 2008;11(8):B127e31. [13] Advani SG, Sozer EM. Process modeling in composites manufacturing. CRC; 2002. [14] Mukherjee PP, Kang Q, Wang C-Y. Pore-scale modeling of two-phase transport in polymer electrolyte fuel cells e progress and perspective. Energy Environ Sci 2010;4(2):346e69. [15] Jaganathan S, Vahedi Tafreshi H, Pourdeyhimi B. A realistic approach for modeling permeability of fibrous media: 3-D imaging coupled with CFD simulation. Chem Eng Sci 2008;63(1):244e52. [16] Kinney JH, Nichols MC. X-ray tomographic microscopy (XTM) using synchrotron radiation. Annu Rev Mater Sci 1992;22(1):121e52. [17] Fredrich J, Menendez B, Wong T. Imaging the pore structure of geomaterials. Science (New York, NY) 1995;268(5208):276. [18] Lee J, Hinebaugh J, Bazylak A. Synchrotron X-ray radiographic investigations of liquid water transport
[23]
[24]
[25] [26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
9385
behaviour in a PEMFC with MPL-coated GDLs. J Power Sources 2013;227:123e30. Lee H-K, Park J-H, Kim D-Y, Lee T-H. A study on the characteristics of the diffusion layer thickness and porosity of the PEMFC. J Power Sources 2004;131(1):200e6. Fishman Z, Bazylak A. Heterogeneous through-plane distributions of tortuosity, effective diffusivity, and permeability for PEMFC GDLs. J Electrochem Soc 2011;158(2):B247e52. Fishman Z, Hinebaugh J, Bazylak A. Microscale tomography investigations of heterogeneous porosity distributions of PEMFC GDLs. J Electrochem Soc 2010;157(11):B1643e50. Becker J, Flu¨ckiger R, Reum M, Bu¨chi FN, Marone F, Stampanoni M. Determination of material properties of gas diffusion layers: experiments and simulations using phase contrast tomographic microscopy. J Electrochem Soc 2009;156(10):B1175e81. Becker J, Schulz V, Wiegmann A. Numerical determination of two-phase material parameters of a gas diffusion layer using tomography images. J Fuel Cell Sci Technol 2008;5(2):21006. Thiedmann R, Hartnig C, Manke I, Schmidt V, Lehnert W. Local structural characteristics of pore space in GDLs of PEM fuel cells based on geometric 3D graphs. J Electrochem Soc 2009;156(11):B1339e47. Wang Y. Periodic surface modeling for computer aided nano design. Computer-Aided Des 2007;39(3):179e89. Schulz VP, Becker J, Wiegmann A, Mukherjee PP, Wang C-Y. Modeling of two-phase behavior in the gas diffusion medium of PEFCs via full morphology approach. J Electrochem Soc 2007;154(4):B419e26. Gaiselmann G, Thiedmann R, Manke I, Lehnert W, Schmidt V. Stochastic 3D modeling of fiber-based materials. Comput Mater Sci 2012;59:75e86. Wang X, Zhang K, Yang Y, Wang L, Zhou Z, Zhu M, et al. Development of hydrophilic barrier layer on nanofibrous substrate as composite membrane via a facile route. J Membr Sci 2010;356(1e2):110e6. Tamayol A, McGregor F, Bahrami M. Single phase throughplane permeability of carbon paper gas diffusion layers. J Power Sources 2012;204:94e9. Daino MM, Kandlikar SG. 3D phase-differentiated GDL microstructure generation with binder and PTFE distributions. Int J Hydrogen Energy 2012;37(6):5180e9. Thiedmann R, Manke I, Lehnert W, Schmidt V. Random geometric graphs for modelling the pore space of fibre-based materials. J Mater Sci 2011;46(24):7745e59. Gostick JT, Fowler MW, Pritzker MD, Ioannidis MA, Behra LM. In-plane and through-plane gas permeability of carbon fiber electrode backing layers. J Power Sources 2006;162(1):228e38. Thiedmann R, Fleischer F, Hartnig C, Lehnert W, Schmidt V. Stochastic 3D modeling of the GDL structure in PEMFCs based on thin section detection. J Electrochem Soc 2008;155(4):B391e9. Gostick JT, Ioannidis MA, Fowler MW, Pritzker MD. Pore network modeling of fibrous gas diffusion layers for polymer electrolyte membrane fuel cells. J Power Sources 2007;173(1):277e90. Inoue G, Yoshimoto T, Matsukuma Y, Minemoto M. Development of simulated gas diffusion layer of polymer electrolyte fuel cells and evaluation of its structure. J Power Sources 2008;175(1):145e58. C¸ec¸en A, Wargo E, Hanna A, Turner D, Kalidindi S, Kumbur E. 3-D microstructure analysis of fuel cell materials: spatial distributions of tortuosity, void size and diffusivity. J Electrochem Soc 2012;159(3):B299e307. Dijkstra EW. A note on two problems in connexion with graphs. Numer Math 1959;1(1):269e71.
9386
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 9 3 7 5 e9 3 8 6
[38] Gonzalez RC, Woods RE, Eddins SL. Digital image processing using MATLAB. Pearson Education India; 2004. [39] Image processing toolbox documentation. Available from: www.mathworks.com/help/images/index.html; 2010. [40] ANSYS FLUENT user’s guide, release 13.0. ANSYS; 2010. [41] Becker J, Wieser C, Fell S, Steiner K. A multi-scale approach to material modeling of fuel cell diffusion media. Int J Heat Mass Transf 2011;54(7):1360e8. [42] Flu¨ckiger R, Freunberger SA, Kramer D, Wokaun A, Scherer GG, Bu¨chi FN. Anisotropic, effective diffusivity of
porous gas diffusion layer materials for PEFC. Electrochim Acta 2008;54(2):551e9. [43] Feser J, Prasad A, Advani S. Experimental characterization of in-plane permeability of gas diffusion layers. J Power Sources 2006;162(2):1226e31. [44] Spielman L, Goren SL. Model for predicting pressure drop and filtration efficiency in fibrous media. Environ Sci Technol 1968;2(4):279e87.