Establishing Accuracy of Watershed-Derived Pore Network Extraction for Characterizing In-Plane Effective Diffusivity in Thin Porous Layers

The GDL samples considered in this work (SGL 25AA and SGL 34AA) were examined across a range of compression ratios. Each GDL sample was made of a fibrous carbon paper without an MPL. The detailed, 3D structures were obtained using an X-ray-based microscale CT scanner (Skyscan 1172). The pixel resolution of the radiographs ranged between 2.88 μm/px and 3.23 μm/px. The sample areas were 5 × 5 mm² for the uncompressed samples and 4.0 × 1.5 mm² for the compressed samples. To obtain uncompressed and compressed GDL structures, the GDLs were held with a custom setup as shown in Figures 1a and 1b, respectively. To image the uncompressed GDL structure, the GDL sheet was held on a single edge, while the uncompressed region was imaged. The compressed GDL was imaged using the holder presented in Figure 1b, and the compression ratio was controlled using stacked polyethylene naphthalate (PEN) sheets. The SGL 25AA and SGL 34AA GDLs were imaged for a range of compression ratios. For each compression ratio used, the associated strain was computed as the difference between the uncompressed thickness, l_o, and the compressed thickness, l, divided by the uncompressed thickness as follows: strain = (l₀ − l)/l_o. Table I summarizes all the imaged samples with their associated compression ratios, strains, and compressed thicknesses.

Zoom In Zoom Out Reset image size

Examples of the IP and through-plane greyscale images of the 3D reconstructed GDL (obtained through CT imaging) are presented in Figure 1c. We used an in-house segmentation algorithm similar to Ref. 42 to carefully distinguish between the solid and void phases. The algorithm was used to threshold the greyscale images based on the fiber density computed as:

where σ_w is the areal weight of the sample, and ε₀ is the uncompressed porosity of the sample given by the manufacturer (SGL Sigracet, see Table I). Small variabilities in GDL porosity values have been reported in past literature;⁴³ however, these variabilities do not exceed ±5% (for example Rashapov et al.⁴⁴ measured 0.884 and 0.841 for SGL 25AA and 34AA, respectively). We assume that these small variations in porosity would only impact a few voxels and not affect the global segmentation result.

In order to compute the fiber density, we measured areal weight of the GDL. The areal weight was measured to be 38 g/m² and 77 g/m² for the SGL 25AA and SGL 34AA GDLs, respectively. They were measured by determining the masses of GDL samples with known dimensions. One to six GDL samples were cut using a 1 inch-punch and stacked together to measure their cumulative masses.

Finally, once the images were segmented, a 3D median filter with a radius of 1 voxel was used to clean the binary images.

The GDL porosity was computed based on the binary image obtained after segmentation. The number of solid voxels (SV) – consisting of carbon paper fibers – and void voxels (VV) were summed at each through-plane position, and the porosity distribution was computed as follows:

where ε is the GDL porosity, and z is the through-plane direction.

Similar to the Maximal Ball Method (MBM) presented by Dong and Blunt,⁴⁵ the algorithm employed in this work was initiated with a distance map of a binary pore space, populated with the Euclidean distances between each pore voxel and its nearest solid voxel. Similar to the MBM, the local maxima in this distance map were assumed to be pore centers, as these are the locations and radii of the largest possible spheres that could occupy the space. However, unlike the MBM, the pore space was divided using a customized Watershed⁴⁶ algorithm on the inverse of the distance map. A set of 3D matrices of identical shapes were generated as described below. These matrices, once generated, were accessed at various points in the network extraction algorithm (Procedure 1, below). While the algorithms discussed were applied to 3D datasets, the figures accompanying the following sections were created in 2D for illustrative clarity.

The extraction algorithm presented here is admittedly slower than similar algorithms, which take advantage of highly optimized watershed sub-algorithms (a great example of which was presented by Gostick²⁹). We utilized the following algorithm written from the ground up in efficient languages, such as C++ combined with Matlab for the interface, thereby providing a highly customizable platform to add rules for preventing pore over- or under-segmentation at any point in the process. We believe this flexibility may be of use as the field of CT derived pore networks expands.

The box-shaped binary domain was first oriented such that one pair of opposing faces, labelled "top" and "bottom", of the six sided domain corresponded to the inlet and outlet faces for any intended pore network simulation. The other four walls were then padded with solid voxels in order to simulate a solid enclosure around the sample. This binary matrix was labelled , where , if voxel i is solid (shown as black in Figure 2), and , if voxel i is void (shown as white in Figure 2).

Zoom In Zoom Out Reset image size

In some circumstances, such as with stochastically modelled materials, it may be useful to generate an extracted pore network with periodic side-wall boundary conditions. In these cases, instead of a padding of solid voxels, the four side-walls may be padded with a number of voxel planes from the opposite face of the domain. The pores may then be stitched together at the end of the extraction process.

A new matrix was next formed with the same shape as the binary matrix; however, instead of 0s and 1s, the matrix was filled with the linear index value for each voxel position, where the n voxels were each assigned an index between 0 and n − 1, creating matrix . It is useful to index , such that . In this process sometimes referred to as "flattening", we can now treat an LxWxH matrix as an (L^*W^*H)x1 array.

A Euclidean distance map, , was defined using the bwdist function from the MATLAB Image Processing Toolbox.⁴⁷ In this matrix, , if voxel i was a solid, and , if voxel i was a void, where d was the Euclidean distance between the center of voxel i and the center of the nearest solid voxel.

While the pore and throat radii were ultimately found from , a 3D median filter (using a spherical kernel with radius of 2 voxels) was applied to , resulting in . This filtering process temporarily suppressed the sharp features within , which prevented an unnecessary oversegmentation of the pore space.

A new matrix, , which ranked the voxels according to their value in descending order was generated (lowest number has the highest rank). In the case where multiple voxels had identical values, the rank was decided based on a random permutation of the set of like voxels. This process provided each voxel with a unique, integer value for between 0 and n − 1, and the random permutation prevented a directional bias during extraction.

The matrix that eventually stores the segmented pore space, labelled , was initiated with a duplicate of and modified such that all solid voxels () were replaced with a value of −1. After the pore segmentation process, the void voxels within were populated with pore names.

Pore segmentation algorithm

In this algorithm, each pore name was defined by the index, i, of the ancestor voxel of that pore. The ancestor voxel is the highest-ranking voxel in that pore, meaning that it has the property that for all voxels, j, in the pore. Due to the definition of , this resulted in ancestors with maximal distance values (i.e. for all voxels j in the pore).

The global algorithm followed the logic described hereafter. Initially, all void voxels were defined as their own ancestors (i.e., ). Then, as the focus voxel, i, was incremented from the 0^th rank (i.e. ) to the n − 1^th rank, void neighbors were evaluated. Each neighboring void voxel, j, could potentially change its ancestor value to that of the focus (i.e., can become ), if , (i.e., ).

In some cases, (along a longitudinally symmetric tube, for example) a string of equal values could exist whereby a number of segmented pores could be created, depending on the randomized rank of each of the central voxels. In such an instance, these additional pores would be collapsed into a single pore. This was achieved by checking whether the ancestor's distance value was equal to that of the focus voxel (i.e., ). When this condition was true, all voxels which had the ancestor, , took on the new ancestor, .

Further measures could be taken to prevent over-segmentation of the void phase, and this logic can be either placed within the primary sweep of ranked voxels or as a secondary sweep of identified pores, including throat radius values in the decision logic. However, no additional measures were performed for this study.

This algorithm summarized above is achieved with the following logical steps beginning at voxel i where :

Procedure 1. Algorithm used to reach a segmented state.

This was repeated for i, where , then again for i, where . This process was repeated while .

Figures 2 contains a 2D demonstration of the matrices involved in the segmentation process, including a finalized matrix with the pore voxels outlined. In addition to the matrix, the following pore characteristics were then found at the end of the pore segmentation algorithm:

• Pore radius was defined as the maximum value of within that pore.
• Pore coordinates were defined as the x,y,z values of the position of the maximal value.
• Pore volume was defined as the number of voxels within that pore, multiplied by the voxel volume.

Due to the nature of the pore segmentation sequence, pores interfaced with each other at surfaces containing the three-dimensional saddle points in . These saddle points are the points between the pores that define the size of the largest sphere that could pass from pore to pore. The voxels at and around the saddle point were analyzed for throat radius calculation as follows.

A list of potential throat voxels was created by finding all voxels, i, that have a , that were also neighbors to at least one voxel, j, such that and . For all potential throat voxels, a list of its 26 neighbors in was compiled. If the two most frequent non-negative values present constituted at least 12 of the 26 neighbors, and one value did not constitute more than 21 of the 26 neighbors, that voxel was added to a list of throat voxels connecting the two dominant pores. The above criteria were used to minimize the number of unwanted connections made between pores with irregular shapes. This list was then evaluated on the non-filtered distance map, , to determine the throat diameter. For completeness, the algorithm procedure of this throat extraction is provided below (see Procedure 2).

Procedure 2. Algorithm used to segment the throat voxels.

Particular attention was directed to the throat diameter, since the throat diameter plays a significant role in the mass transport calculation (both using invasion percolation or diffusive process). To prevent the discretized values in from resulting in only a few, discrete throat diameters, the diameter of each throat i, was calculated as follows:

where r₁ and r₂ are the two largest unique numbers in the list of voxels evaluated on for a given throat, n₁ is the number of voxels in the list with a value of r₁, and n₂ is the number of voxels in the list with a value of r₂. This resulted in a radius value slightly larger than r₁; however, this radius value was ranked in a way that made high aspect ratio throats slightly larger than those with aspect ratios closer to 1. This throat diameter calculation is fully detailed in Procedure 3.

Several hours on a regular desktop computer are required to extract the pore network (depending on the domain size). Diffusion calculations are completed within seconds.

Procedure 3. Algorithm used to compute the throat diameter.

Once the pore and throat networks were extracted from the segmented images, the networks were analyzed with respect to the pore and throat size distributions. This analysis was performed by fitting the pore and throat size distributions to the log-normal distribution, similar to Zenyuk et al.²⁸ Unimodal distributions were chosen as suggested by Zenyuk et al.²⁸ for GDL samples compressed using large stains (higher than 0.20).

where the histogram area, A, the mean, μ, and the standard deviation, σ, were numerically fitted onto the extracted pore and throat size distribution. These parameters were then used to compute the average pore and throat diameters, , as well as their standard deviation, d_std, as follows:

Once the pore characteristics (sizes and location) and connectivities were extracted from the CT images of the GDL, the diffusive transport within the GDL pore space was defined. The molar flux of oxygen between two pores connected by a throat (see Figure 3) can be written as,

where N_{1, 2} is the molar rate between pores 1 and 2, g_{1, 2} is the diffusive conductance between pores 1 and 2, and and are the mole fractions of pore 1 and 2, respectively. The diffusive conductance was based on the pore and throat conductances as follows,

where and are the pore conductance of pores 1 and 2, respectively, and g_t is the conductance of the connecting throat. The pore conductances were defined based on their structural properties, which were calculated during the pore network extraction as,

where c_tot is the total air molar concentration, D_b is the bulk diffusivity, is the cross sectional area of the pore i, and is the radius of the pore i. In Equation 9, a small pore has a small cross sectional area which leads to a small conductance. The throat conductance, , was defined in a similar way as,

where A_t is the throat cross sectional area, and L_t is the throat length between the two considered pores. The throat length is defined as the length between two neighboring pores (see Figure 3). The routine, OpenPNM. Geometry.models.throat_length.straight(network = pn,geometry = geom), was used to compute the length between two pore centers less their respective diameters (as illustrated in Figure 3). It is worth noting that in Equation 10 for the case of narrow and long throats, i.e. small A_t and/or large L_t, the throat conductance, g_t, may be small compared to the pore conductance. This relative difference could lead to the dominance of the throat size on the total conductance. Therefore characterizing the throat size becomes critical for accurately modelling the GDL effective diffusivity.

Zoom In Zoom Out Reset image size

Equations 7 to 10 were applied to all the pores in the extracted network using OpenPNM,¹⁶ which led to a linear system of equations where the mole fraction of oxygen within the pore space are the unknowns. Once the linear system was solved, the effective diffusivity of the GDL, D^eff, was computed as follows:

where N_tot is the total molar flow rate entering in the GDL calculated by OpenPNM (in mol/s), L is the thickness of the GDL, A is the cross-sectional area of the GDL, and x_in and x_out are the inlet and outlet molar fractions, respectively, defined as the boundary conditions at the extremities of the GDL. The parameters used in our simulation are given in Table II.

The through-plane porosity distribution of the segmented images was first analyzed for each GDL sample. The reduction in thickness resulted in decreases in porosity for both the SGL 25AA and SGL 34AA GDLs (Figures 4a and 4b). The uncompressed thicknesses reported in Table I correspond to the porosities where ε < 0.95. Each of the porosity profiles presented a wavy-like pattern with a minimum of porosity close to the extremities (e.g., at −100 μm and 75 μm for the uncompressed 34AA GDL). Such wavy-like patterns were also observed by Fishman et al.⁴⁹ for the SGL series GDLs. By increasing the compression of the GDL, the porosity distributions became increasingly flat in the GDL core. This was particularly apparent in the 34AA sample at compressive strains of 0.40, 0.50, and 0.60 (Figure 4b).

Zoom In Zoom Out Reset image size

For a given compressive strain, the measured porosity of SGL 25AA was lower than SGL 34AA. This trend is in agreement with the specifications given by the manufacturer for SGL 25 and 34 series GDLs, as well as with the literature.³⁶ In Figure 4c, the mean through-plane porosity distributions are presented versus the compressive strain of the sample. As previously mentioned, the SGL 25AA GDL samples were less porous than the SGL 34AA GDL samples. This justified the need to go to higher strains for the SGL 25AA GDL in order to significantly decrease the GDL porosity.

In Rashapov et al.,³⁶ it was reported that the porosity in a compressed sample can be obtained with ε = 1 − l₀(1 − ε₀)/l, where ε₀ is the initial thickness. Using the definition of the strain given in CT and segmentation section, this expression leads to the following:

The comparison of Equation 12 to the porosity of the compressed sample versus strain is presented in Figure 4c, and excellent agreement between the model and the experimental data was observed for both GDL samples, thereby validating our methodology used to segment the CT images.

The equivalent pore network was extracted using the algorithm described in Section 2.2. The void space of each GDL binary tomogram (Figure 5a) was divided into pores. A resulting image of the pore network of the uncompressed SGL 34AA GDL is presented in Figure 5b, where each color represents a distinct pore, and the fibers were colored in white. The void space between two fibers in Figure 5a was divided into multiple pores due to the presence of other fibers in the neighboring through-plane positions. In addition to the pore location, the extraction algorithm was also used to compute the pore diameter, volume, and connectivities of the neighboring pores. The throats (not depicted in Figure 5b) were extracted from the binary images, and their diameter and connectivity were also computed, as discussed in Section 2.2. This data was used to define the conductivity between each pore in the network, according to Equations 6 to 9.

Zoom In Zoom Out Reset image size

The pore network was then analyzed with respect to the pore and throat diameter distribution produced via the pore network segmentation. The frequency of pore and throat diameters were computed over a range spanning from 1 to 150 μm every 5 μm for all GDL samples. The resulting distributions are presented in Figures 6a to 6d and fitted to the log-normal law presented in Section 2.2. All parameters obtained from the curve fitting are presented in Table I in the supplementary information section. As the compressive strain increased, the mean and standard deviations of the pores and throats decreased for both GDL series. This effect was expected due to the compression of the GDL samples which reduced the void spaces between fibers, resulting in smaller pore diameters.

Zoom In Zoom Out Reset image size

To characterize pore diameter distributions, a pore network extraction is not necessarily required. Mercury porosimetry invasion (MIP) techniques can also be used to characterize the internal structure of a porous material. However, throat size distributions are not easily measurable, even though throat sizes provide critical data necessary to study fluid transport within the porous media. As such, there is a great advantage to using a pore network extraction algorithm based on CT images to compute the throat size distribution of the porous domain. Using the parameters reported in Table I from the supplementary material section, and the log-normal distribution (Equation 3), the pore and throat diameter distributions can be reconstructed for the SGL 25AA and SGL 34AA GDLs for a range of compressive strains.

The effect of the compressive strain on the mean GDL pore and throat diameters was analyzed in Figures 6e and 6f for both GDLs. The mean pore diameter was observed to linearly decrease with higher compressive strain, an effect already described by Zenyuk et al.²⁸ In the SGL 25AA GDL sample, the mean pore diameter decreased from 53 μm to 28 μm. In the SGL 25AA GDL sample, the mean pore diameter decreased from 53 μm to 28 μm. In comparison with increasing compressive strain, the SGL 34AA GDL mean pore diameters were smaller than the SGL 25AA GDL pores, ranging from 38 μm to 29 μm. In Figure 6c, the slopes of the linear curves of the mean pore diameter were observed to be −22 μm⁻¹ and −17 μm⁻¹ for SGL 34AA and SGL 25AA, respectively, indicating that the SGL 25AA GDL was more sensitive to compression than SGL 34AA. In contrast, the mean throat diameters were found to be relatively constant regardless of the applied compressive strain on both GDL samples. These results may indicate that the compression of the GDL samples had a significant impact on the average pore diameters (but not on the throat diameters).

The inlet (left side) and outlet (right side) oxygen mole fractions were arbitrarily defined as 1 and 0, respectively. Note that the value of these boundary conditions did not change the results of the effective diffusivity calculations. A linear system of N-equations (corresponding to the N-connectivities that exist within the pore network) was solved using the pore network modeling package OpenPNM.¹⁶ The mole fraction in each of the pores of the network was similarly computed. A result of these calculations is given in the supplementary materials, where the map of concentration is presented. This shows that the local information (such as pore molar concentration) can also be obtained using our PNM extraction.

The molar flow rates in all of the throats of the network, which depended on the local pore and throat structure (see Equation 6), were also computed. The total molar flux flowing through the network was then computed at the inlet and used to calculate D^eff using Equation 11. All of these calculations were applied to the nine GDL samples considered in this study.

The normalized diffusivities (effective diffusivity divided by the bulk diffusivity) are presented in Figure 7a for the nine GDL samples considered in this study. These results are compared with the data presented by Rashapov et al.³⁶ For the sake of clarity, the IP effective diffusivity measurements as a function of porosity reported by Rashapov et al. for SGL 25AA and SGL 34AA GDLs were used and fitted to a power law. The effective diffusivities computed from the pore network extractions are in good agreement for both GDLs, with the exception of SGL 25AA GDL samples that underwent a compressive strain of 0.58 and 0.72. These discrepancies are discussed in the next paragraph. The computed diffusivity for the SGL 25AA GDL was higher than that for the SGL 34AA GDL. We attribute this difference to the related structural differences. For example, the SGL 25AA GDL is thinner and more porous (see GDL structural analysis and Pore network analysis sections) than the SGL 34AA GDL, allowing for an increased rate of gas transport through the SGL 25AA GDL. In a similar way, decreasing the GDL porosity resulted in a decreased mean pore diameter, which also led to a significant decrease in the effective diffusivity. Therefore, the computed effective diffusivity presented in this study was in good agreement with the literature and experimental measurements, thereby validating our combined approach of watershed-derived pore networks and computed tomography for characterizing fuel cell GDL effective diffusivity.

Zoom In Zoom Out Reset image size

In Figure 7b, the relative error between the computed effective diffusivity presented in this study, and the fit obtained from the experimental measurements reported in Ref. 36 are compared with respect to the compressive strain of each sample. A relatively low error (below 10%) was found for all the samples compressed with a compressive stain less than 0.55. Above this value, the relative error significantly increased to 60%. Therefore, our results indicate that a limit exists to the applicability of our method. By highly compressing our GDLs (compressive strains greater than 0.55), the sizes of the pore diameters approached the resolution of the CT scanner used in this study (a few microns). At these extremes, the structure of these small pores cannot be adequately captured from microscale CT images.⁴² The structural characteristics of the highly compressed GDLs may not have been fully captured during network extraction process, leading to significant errors in the predicted effective diffusivity of oxygen. Therefore, a sub-micron resolution CT scanner is recommended for studying the transport properties of GDLs that have undergone high compressive strains. Finally, it can be noted that for a given strain of 0.6, for example, this effect was more pronounced for the SGL 25AA GDL (71 μm-thick) compared to SGL 34AA GDL (100 μm-thick). At this high strain, SGL 34AA retained larger pore sizes compared to SGL 25AA.