Stochastic Analysis of the Gas Flow at the Gas Diffusion Layer/Channel Interface of a High-Temperature Polymer Electrolyte Fuel Cell

Featured Application: Featured applications are channel or cell simulations that can consider information about the distribution of steam coming out of the GDL. Abstract: Gas diffusion layers (GDLs) play a signiﬁcant role in the efﬁcient operation of high-temperature polymer electrolyte fuel cells. They connect the electrodes to the gas channels of the bipolar plate by porous material with a meso-scale geometric structure. The electrodes must be sufﬁciently supplied by gases from the channels to operate fuel cells efﬁciently. Furthermore, reaction products must be transported in the other direction. The gas transport is simulated in the through-plane direction of the GDL, and its microstructure created by a stochastic model is equivalent to the structure of real GDL material. Continuum approaches in cell-scale simulations have model parameters for porous regions that can be taken from effective properties calculated from the meso-scale simulation results, as one feature of multi-scale simulations. Another signiﬁcant issue in multi-scale simulations is the interface between two regions. The focus is on the gas ﬂow at the interface between GDL and the gas channel, which is analyzed using statistical methods. Quantitative relationships between functionality and microstructure can be detected. With this approach, virtual GDL materials can possibly be designed with improved transport properties. The evaluation of the surface ﬂow with stochastic methods offers substantiated beneﬁts that are suitable for connecting the meso-scale to larger spatial scales.


Introduction
In all types of polymer electrolyte fuel cells (PEFC), also called proton exchange membrane (PEM) fuel cells, the gas diffusion layers (GDLs) are components with high relevance for efficient operation of fuel cells. Two major requirements can be identified. One is the electric contact to be provided between the bipolar plates and the catalyst layer (CL). This determines the choice of the material. The second requirement is to facilitate efficient mass transport. For this purpose, an appropriate microstructure is required. To meet the requirements, carbon fibers are typically used for the fabrication of different types of GDLs, e.g., paper, woven, or non-woven textiles. The lattice Boltzmann method (LBM) is water droplets in a channel from the surface of different types of GDL. Wang et al. [29], meanwhile, identified the GDL surface near the channel as a topic of high interest for the simulation of gas flow on the channels of the flow field of a fuel cell. Niu et al. [30] simulated two-phase flow in fuel cell channels, assuming a static contact angle at the GDL surface. Kim et al. [31], in turn, simulated the hydrodynamics of water droplets in the gas channels. They investigated droplets leaving the GDL surface at two distinct positions. Koz and Kandlikar [32] simulated the inhibition of oxygen transport in flow channels in the presence of liquid water. They worked with regular patterns at the GDL interface where the water was entering the channel. The position of liquid water transported from the GDL into the gas channel can be inherently transferred by coupling the simulation domains of both spatial scales, as was done by Chen et al. [33]. They found that such tight coupled simulations can require enormous computational resources.
Two adjacent regions of a fuel cell are sometimes simulated on different spatial scales. From a macroscopic view, the interface is a 2D element. The role of the GDL/channel interface in this scenario is illustrated in Figure 1. The GDL/channel interface was investigated by Yu et al. [34], who analyzed irregular contact angles of water droplets at the GDL surface and their pattern when they passed the GDL/channel interface at several positions. The relevance of the interface for multi-scale simulations was shown by Qin et al. [35,36] and Aghihi et al. [37], who used pore network modeling (PNM) to bridge the scales (still on water transport in low temperature PEFCs). Niu et al. [38] coupled the LBM in the porous GDL structure with OpenFOAM simulations in the air channel.  The relevance of the characterization of the interface for multi-scale simulations still holds for other types of fuel cells, e.g., HT-PEFCs. Yang et al. [39] analyzed the gas flow in a channel over a regular porous structure for different Reynolds numbers. The interface was included in the comprehensive analytical studies of Kulikovsky [40]. The studies of Chevalier et al. [41] showed that the Peclet number at the GDL/channel interface is relevant for the current density profile along the channel. They focused their studies on oxygen transport in the GDL and channel, as well as on charge transport in the membrane.
In this manuscript, through-plane transport in GDLs is simulated with the LBM. The stochastic geometry model of Thiedmann et al. [8] is used to create 25 representations of the microstructure; Froning et al. [4,5] sed the same geometries. The statistical variation of the microstructure was completed with various compression levels. For this manuscript, virtual microstructures were selected from the studies mentioned above, both uncompressed and compressed. The focus of the new investigations is the analysis of the two-dimensional region between the GDL and gas channels; this is called the GDL/channel interface. Areas are classified according to the total amount of gas leaving the GDL with the highest velocity. The location of areas at the GDL surface where the most gas is flowing is the resulting 2D information on the GDL surface. The knowledge can possibly assist the development of new methods in the field of channel/cell-level simulations.

Methods
The LBM was applied in transport simulations of stochastic microstructures, with both methods being the same as was presented earlier [4,5,26].

Geometric Data
Thiedmann et al. [8] developed a stochastic geometry model that describes layers of paper-type GDL-in particular, Toray 090-as intersecting lines, with every layer being independent of the others. Three-dimensional realizations of a series of such layers were validated against the real 3D structure of the GDL using synchrotron data from the Helmholtz Center in Berlin [42]. Binder was covered on random polygons built by the fibers in a layer. It was shown that an 18-µm thickness of this binder coating is reasonable for transport simulations in this kind of geometry [4].
Three realizations of the geometry model are shown in Figure 2. Because of its stochastic nature, the microstructure shows large differences at local positions. Statistical differences between adjacent fiber layers are illustrated by the black-colored top layer at the GDL exit and the gray-colored layer below. The different looking images in Figure 2a-c show local differences between the realizations, all of which represent the same material. The images of size 512 × 512 represent a section of 768 µm × 768 µm with a resolution of 1.5 µm per pixel. In this way, five images form a layer of fibers with a thickness of 7.5 µm each, while 130 images representing 26 fiber layers of a GDL. Compressed GDL structures were generated by merging images of adjacent fiber layers. Froning et al. [4] presented the merging algorithm in detail and how compression levels in steps of 10% can be accomplished. It was also shown that the algorithm led to sufficient accuracy of flow calculations for compression levels of up to 30%.

Lattice Boltzmann Method
The through-plane transport of water vapor is simulated with the LBM. The method statistically describes ensembles of many molecules according to the principles of gas kinetics [43]. On a regular lattice-in this work, the D3Q19 scheme is used-in this three-dimensional scheme at each grid location 19 neighbors are connected, i.e., functions f i ( x, p, t), i = 0, ..., 18, are defined, specifying the probability that a molecule can be found at the place x and time t that moves to a neighboring node in the lattice with momentum p [44]. With the Bhatnagar-Gross-Krook (BGK) scheme [45], the equilibrium state can be numerically calculated by using Equation (1): Here, f eq is the Maxwellian distribution, f (n) and f (n+1) are the fields of the f i . The superscripts (n) and (n+1) specify the time t and the next time step t + ∆t.
This variant is feasible for single-phase, single-component transport simulations to be applied on the given kind of microstructures [4,5]. From the resulting velocity field, the volume-based characteristics of permeability κ = −q · µ ∇P and tortuosity τ = |v| v x were calculated [46][47][48].

Simulation Frame
The transport simulations were applied in the small section shown below the channel in Figure 1. The gas is transported from the bottom to the top, towards the air channel. For the purpose of characterization of the GDL material, the gas flow in the channel of an operating fuel cell is not considered, avoiding the composition of different effects that could influence the results.
The microstructure of the GDL is specified by a series of binary images, as shown in Figure 2. In this way, 130 images of size 512 × 512 define a 512 × 512 × 130 lattice. With a resolution of 1.5 µm, this defines a section of 768 µm × 768 µm of a GDL with a 195 µm thickness, consisting of 26 fiber layers, in accordance with the stochastic geometry model of Thiedmann et al. [8]. A series of images specifies an irregular porous structure. A bounce back condition is applied at the fiber surface, representing no-slip conditions. The definition of relevant boundary conditions is enabled by free space added upstream and downstream, as was already used in previous studies [4,5]. At the inlet, a Dirichlet boundary condition is specified, implemented as a fixed velocity that corresponds to a total mass flow of steam produced by an average current density of 1 A/cm 2 . At the outlet, a Neumann boundary condition, also known as constant pressure, is specified. At the side walls, a slip condition was applied. The goal is the detection of the relationships between well-known volumetric properties like permeability or tortuosity and the characteristics of the surface flow. Single-phase, single-component flow is simulated to obtain characteristics related to the pure GDL material without being affected by physical effects from outside the GDL. In Table 1, the values are completed with the amount of hydrogen and oxygen needed for the electrochemical conversion of 1 A/cm 2 as a reminder of Faraday's law.
Compared to earlier work [4,5,26], the binary images representing the microstructure are used in reverse order. For the calculation of average volumetric values and statistical evaluation, this is unnecessary, and the only reason is to avoid local inconsistencies with the evaluation of the other interface (GDL/electrode) presented by Froning et al. [26].
Geometries with 30% compression were created by merging adjacent fiber layers as mentioned in Section 2.1 and described in more detail by Froning et al. [4]. In this case, 92 images represent a GDL of 138 µm in thickness.

Analysis of the Interface
Similar to as discussed by Froning et al. [26] for the analysis of the GDL/electrode interface, the areas where most of the steam leaves the GDL shall be identified. For this purpose, the x component of the velocity at the GDL exit, which is the top layer of the microstructure, is evaluated. For incompressible flow, the velocity is related to the mass flow. Mass-related quantiles are defined at the GDL exit, which is the GDL/channel interface according to the transport simulation of steam through the GDL. Based on the total mass flow: where A is the area of the GDL beneath the flow field (related to the cell-scale) and u x is the through-plane velocity. The remaining constant c is not needed in the subsequent evaluation, because only relative values are used for the definition of q ∈ [0, 1] and the total mass fraction M q : The mass-related quantile z is implicitly defined by the following: Finally, the root of the function U(q, z) is the quantile z needed in Equation (4): In the final step, contour levels u q (z) were used to visualize the mass-related quantiles of the through-plane velocity at the GDL surface using paraview [49]. While the contour levels u q (z) are useful for visualization, their integral functions U(q, z) can quantify the area.

Results
Through-plane transport was simulated in 25 realizations created by the geometry model, oriented in reverse order to distinguish the GDL/channel surface from the GDL/electrode surface presented by Froning et al. [26]. The operating conditions, according to Section 2.3, are summarized in Table 1.
The conditions led to a Reynolds number of 2.4 × 10 −4 . This in turn led to velocity vectors at the GDL exit that were almost parallel to the through-plane direction. The free space downstream of the GDL mentioned in Section 2.3 was required to arrange the velocity vectors properly behind the porous structure.
The variation of the volumetric characteristics-permeability κ and tortuosity τ-is shown in Figure 3a. A comprehensive study of these was presented by Froning et al. [4,5]. The quantile levels z were implicitly defined by u q (z), as introduced by Equation (4). They define the total quantile areas: that can be summarized to the values illustrated in Figure 3b. Like the volumetric characteristics, they showed statistical variation, because the simulation results were based on statistical microstructures. In addition to the total quantile area defined by Equation (6), the size of the largest of these regions is shown in Figure 3c. For this purpose, the function U(q, z) from Equation (5) needs to be applied on the largest area identified by the contour levels u q (z) from Equation (4), which was done via the visualization tool paraview and an external R [50] script. In contrast to the display of the key properties in Figure 3, the detail of the two realizations, Nos. 5 and 15, is depicted in Figure 4. With the decrease of the quantile level q, the areas resulting from the contour lines became smaller, while the number of areas decreased, as well. The contour levels addressed in Section 2.4 look completely different in detail, which depicts the large variation of the areas summarized in Figure 3c. Mass-related quantiles for q = 70%, 50%, and 20%. Labels "a"-"d" are described in Section 3.1.

Analysis of the GDL Surface
The GDL/channel interface is at the location where the through-plane flow leaves the GDL. This interface was analyzed by evaluating the mass-related quantiles defined by Equation (4). For two realizations, the 70%, 50%, and 20% quantiles are shown in Figure 4. With lower quantile levels, the sizes of the areas illustrated by the contour lines in the post-processing step decreased, and some characteristic situations could be identified. Areas that are concentric for the three quantile levels with similar shapes are marked with labels, e.g., "a" in Figure 4. In a similar manner, the areas can still be concentric, but with a modified shape of the inner area (for q = 20%), labeled with "b". Another situation is marked with "c". The inner area was split into two, identifying two small pores in the GDL that were hidden, when higher quantile levels were shown. Positions where only one or two areas were nested-labeled as "d"-showed pores, where less gas left the GDL. Furthermore, positions where "b" and "c" are combined are labeled accordingly.
The gas within the mass-related quantiles was flowing through a certain fraction of the total quantile area of the GDL, which is shown in Figure 3b. The sizes of the largest area show large variation in Figure 3c.
Through-plane transport was simulated in microstructures as specified in Sections 2.1 and 2.3. They were applied to 25 realizations of the geometry model. Table 2 summarizes the detailed results of the analysis of the mass-related quantiles obtained from the velocity field. The total quantile area s tot is presented for q = 70%, 50%, and 20% and, also, the number of areas (number of regions, #reg.) belonging to the quantiles and the average s ∅ and maximum size s max of the areas. The maximum size is also presented in Figure 3c. The largest area can be of interest for fuel cell and stack modeling when the location of the pores emphasized in the GDL surface is required. The large spread of the values in Figure 3c is in line with the variation coefficient of s max in Table 2.

Statistical Evaluation
The transport simulations provided velocity fields that enabled the calculation of volume-based properties-κ and τ-and surface properties, as discussed in Section 3.1. All of these showed statistical variations and were different in their ranges of values. For the volume-based properties of porous material, many approaches have been developed that describe the relationships between several properties of the material. The Kozeny-Carman relation is only one of them [47]. It is desired to have also a simple approach to calculate the surface properties from volume-based characteristics. As a first step, the statistical correlation between the results can be calculated. Hedderich and Sachs [51] provided not only the fundamentals of various statistical tests, but also criteria for the proper choice of the right test for the desired evaluation. Kendall's correlation test [51] allows for correlation coefficients to be calculated under the conditions mentioned above. Table 3 shows the coefficients r(a, b), a, b ∈ {s tot (q), κ, τ}, r ∈ [−1, 1] that were obtained using the R software [50,52]. The total quantile areas s tot (q) for q ∈ {70%, 50%, 20%} were chosen as candidates for the correlation test with the volumetric properties: permeability κ and tortuosity τ. Correlation coefficients can indicate a positive (r > 0) or negative (r < 0) relationship. Every value other than an extreme one ∈ {−10, 1} can only be used for comparisons. For this reason, the absolute value |r(κ, τ)| was chosen as an indicator of whether a correlation coefficient indicates a possible physical relation. It was already shown by Froning et al. [4] that κ and τ are non-linear with respect to the Kozeny-Carman relation: with only the geometric properties of the microstructure on the right side: the porosity ε, total volume V p , and inner surface S p . In Table 3, the correlation between κ and τ is r(κ, τ) = −0.593. The absolute value of every correlation coefficient between the surface-related properties s tot and κ or τ was much lower than this value. This indicates that there was no statistical correlation between surface-based and volume-based properties, although all of them were calculated from the same velocity field. Therefore, it is still necessary to perform transport simulations on the 3D structure to obtain the 2D properties presented in Section 3.1. The same conclusion was found in the analysis of the other interface (GDL/electrode) by Froning et al. [26].

Impact of the Compression
The through-plane transport was simulated under 30% compression. Figure 5 shows the velocity extracted at the GDL/channel interface with contour levels according to the mass-related quantiles z. The velocity is colored in the same range as is used in Figure 4. Although the transport properties-permeability and tortuosity-changed under compression according to the Kozeny-Carman equation, the distribution of the relative total quantile areas at the interface did not change. The permeability and tortuosity are shown in Figure 6a. Compared to the uncompressed properties in Figure 3a, the tortuosity increased and the permeability decreased under compression, which is consistent with earlier studies [4]. Related to the quantiles of the total mass flow, Figure 6b,c shows the relative total quantile areas of the gas flow and the largest of these. The overall picture of the diagrams in Figure 4 looks very similar to that in Figure 6, but with reduced absolute values. It is noticeable that also the permeability and tortuosity in Figure 6 changed their absolute values compared to the numbers in Figure 3, which was already discussed by Froning et al. [4,5]. Detailed surface characteristics of the compressed GDL are shown in Table 4. The GDL was compressed here by 30%.
The size of the total area s tot was smaller under compression, i.e., the average value was reduced from 167,000 µm 2 (Table 2) to 150,000 µm 2 for the 70% quantile level. The percentage of s tot related to the total surface of 590,000 µm 2 is shown in Table 5. The upper limit of the total quantile area was roughly the porosity of the GDL, which was 78% in the representations of the geometric model [4,8]. In particular, s tot was limited by the local porosity of the upper fiber layer, which may vary slightly from the average. A similar trend can be observed on the other average values in the tables. The reduced porosity of a compressed porous structure leads (on average) to a higher absolute velocity under the condition of a fixed total amount of gas to be transported through the GDL. The higher velocity is caused by smaller pore sizes. The correlation between the surface characteristics of the compressed and uncompressed material is presented in Table 6. Because the range of the values changed under compression, Kendall's test was again chosen as the test method. The entries in the table are the correlation coefficients of the given property, i.e., s tot for the 70% quantile in the first line, first column, obtained from simulation results on uncompressed and compressed microstructures. As before, the value of 0.593-the correlation between κ and τ from Table 3-was taken as a lower limit for judging two characteristics as being related to each other or not. On the basis of this value, the total quantile area s tot was evaluated as being correlated under compression for all quantile levels, as well as the average sizes s ∅ for the quantile levels of 70% and 50%. The average size s ∅ was uncorrelated for the 20% quantile level, and the largest area s max was not correlated in any case. The transfer of the detailed information from transport simulations in the GDL to larger scales of cell and stack simulations requires additional investigations of the channel and flow field modeling. Simulation domains of GDLs were typically much smaller than any flow field of real fuel cells, which is illustrated in Figure 1. Furthermore, the GDL/channel interface connected a region of deterministic geometry with an irregular microstructure. As a consequence, many simulations on stochastic equivalent representations of the microstructure are required to consider the results from GDL transport simulations in the operating or boundary conditions of the cell-level models. Transport simulations on the microstructure can be evaluated statistically to characterize material or surface properties.

Discussion
The use of the volumetric effective properties of porous materials is the state of the art [53]. This is reflected by simulations that use continuum-based approaches for porous regions. These methods are used, e.g., in fuel cell simulations based on commercial CFD software [21,22,54]. Other investigations, not only fuel cell related, used open source CFD software [55][56][57].
Domain sizes for the pore-scale simulations of GDLs are often in the mm range [4][5][6], while cell and stack simulations require domain sizes of many cm, the size of real fuel cell stacks [21]. The tight coupling of both scales is therefore impossible because of computational resources. The GDL/channel interface was already addressed by PEFC modeling reviews focusing on microfluidics [58] and cell-scale modeling [2]. Some kinds of simulation tasks use boundary conditions at the GDL/channel interface. Weber et al. [1] identified the GDL/channel interface as being of high significance when they discussed two-phase phenomena in PEFCs. The path of further evaluation of the GDL/channel interface depends on its use for upscaling the results to cell-/stack-level simulations. In PNM approaches, the GDL can be represented by regularly-located pores and their characterization by randomly-distributed radii and flow behavior at the interface [35,37]. This increases the relevance of the knowledge about the sizes and positions of such pores. Cai et al. [59] studied meander-type flow field channels by placing inlet regions on the GDL/channel interface of their PEFC model. Further investigations in the analysis of the exit surface from transport simulations in GDL can improve such cell-level simulations. The combination of this kind of transport simulations leads to multi-scale approaches in fuel cell modeling.
The approach of analyzing surfaces can potentially also be applied to different simulation techniques. Modeling approaches based on PNM have the potential to cover multiple scales [36,37], and the interfaces between domains of different spatial scales are of high interest for such investigations. In the field of multi-scale simulations, the presented methods can be a vehicle for combining the simulation domains, especially when domains of different spatial scales are connected.

Conclusions
The through-plane transport of water vapor was simulated by means of the lattice Boltzmann method in 25 realizations of a stochastic geometry model representing Toray 090 GDL of an HT-PEFC. The results at the GDL/channel interface were statistically analyzed. For this purpose, mass-related quantiles were specified. Based on the total mass flow the 70%, 50%, and 20% quantiles were evaluated for uncompressed and compressed GDL material. The surface-related results were not correlated with the volume characteristics of the GDL. This was shown by Kendall's correlation test, which was applied to the quantiles of the surface flow and the permeability and tortuosity of the gas flow through the GDL (all data obtained from the same transport simulations). The surface-related analysis of 3D transport simulations in microstructures can possibly support multi-scale investigations in fuel cell modeling.