On the Effect of Complex Permeability and Thermal Material Properties for 3D-CFD Simulation of PEM Fuel Cells

The more and more stringent regulations for the transportation sector are pushing next-generation Internal Combustion Engines Vehicles (ICEV) to ever growing development costs, motivated by the more complex after-treatment devices needed to comply with current and future legislations. Alongside this, the rapid diffusion of Battery Electric Vehicles (BEV) is still limited by technological constraints in specific transportation sectors (e.g. heavy-duty, off-highway, maritime, aerial). These two factors motivate the renewed interest towards Fuel Cell Electric Vehicles (FCEV), ideally complementing BEV in such sectors thanks to the on-board generation of electric power, allowing to overcome the storage/weight limits of batteries. In the automotive transportation sector, Proton Exchange Membrane Fuel Cells (PEMFC) are considered the most promising type of fuel cells, thanks to the low temperature operation and the absence of corrosive substances. Despite the use of PEMFC in industry has been known for decades, there is still a limited understanding of the complex interplay of physical/electrochemical processes governing the key phenomena for high-efficient PEMFC, as well as the role of realistic material characterization. To fulfill this knowledge gap, multi-dimensional numerical models including all the main processes allow an unprecedented comprehensive study [1-3], shedding a light on limiting phenomena and on system bottlenecks. Among the multiple components present in a PEMFC, gas diffusion layers (GDL) have multiple fundamental roles. They have the structural task of mechanically supporting the Membrane-Electrode Assembly (MEA), but they are also constituted by the materials through which species (reactants/products) flows to/from the active sites of the Catalyst Layers (CL), hence influencing the fuel/oxidizer supply rate, which becomes fundamental at high current density conditions where mass transfer rate limitations are frequent. Furthermore, being composed by carbon-based materials (papers and woven/non-woven cloth) there are optimal thermal and electrical conductors, providing a major transfer pathway for heat and electrons generates by electrochemical reactions at CL. Reviews of GDL characteristics (e.g., porosity, permeability, structure, and transport properties) are presented in [4, 5], providing a clear picture of the concurring design needs, as well as the transport mechanisms developing herein [6].

Detailed studies considering the fibrous nature of GDL are common when the focus is the material characterization [7, 8], where samples of the fiber array are modeled to derive semi-empirical or analytical modeled to be compared with experiments. However, the most common modelling approach for GDL is the so-called macro-homogeneous model, constituting in implying a homogeneous structure of equally-spaced/-sized carbon fibers to reflect their layered composition. This is not directly modelled, but the volume occupation of the solid material and the effect on flow/heat/charge transport are replicated by average (macro-homogeneous) properties, implemented as source terms in the dedicated transport equations. The reason of it is the need to preserve a lower limit to the minimum modeled scale to 1 × 10² µm.

The first objective of multi-dimensional numerical models of porous materials is to replicate their resistance to flow, i.e. the momentum conductance. In this context, the hypotheses of low Reynolds number, constant density, single-phase fluid and steady-state flow, the Darcy equation is used to bundle the absolute material permeability in the i-direction $K_{i}$[m²] to the incoming flow with $u_{i}$  velocity and molecular viscosity $\mu$ , as in Eq. (1). This is implemented in the code via user-coding.

$\frac{d p}{d x_{i}}=\frac{\mu}{K_{i}} u_{i}$     (1)

The mentioned low-Reynolds hypothesis is the foundation of the linear flow rate/pressure differential relationship expressed by Eq. (1), and it is questioned in studies regarding the flow regime in GDL Pharoah [9] for which non-dimensional numbers are proposed Zeng et al. [10] and Feser et al. [11] to weight the relevance of a non-linear (inertial) term to be added to the Darcy equation (Darcy-Forchheimer equation).

The second objective for multi-dimensional flows in GDL is related to their specific fabrication, i.e. their largely non-isotropic (layered) composition. To this aim, models should be able to consider the differentiated direction-dependent flow permeability. In particular, the random orientation of carbon planes induces similar permeability values for the in-plane direction $\left(K_{\|}\right)$, whereas a much-decreased permeability is found for the through-plane direction $\left(K_{\perp}\right)$. The permeability anisotropy in GDL is one of the fundamental aspects of interdigitated flow design [12], where all the flow is forced into the porous material via channel obstruction.  This flow arrangement emphasizes the need of an exact characterization of through- and in-plane permeabilities, whose accurate modelling is a mandatory requirement to accurately predict flow intensity in GDL and pressure losses.

Empirical observations and material-scale models have been used to derive large-scale macroscopic correlations for $K_{\|} \text {and } K_{\perp}$. A periodic geometry-based model was proposed in [7], allowing to relate the porosity to the equally-sized fibers distance and to their diameter. The model is composed of distributed normal and parallel layers and blending techniques were proposed to derive global normal and parallel permeabilities. As for $K_{\perp}$, an analytical solution function of the idealized fiber diameter d was derived in Tamayol et al. [8] and validated against normal flow tests. A similar method was followed for $K_{\|}$  [13], proposing an algebraic correlation for the parallel permeability. The proposed forms for $K_{\perp}$ and $K_{\|}$ are both proportional to $d^{2}$ and they asymptotically enclose the measured permeabilities for uniformly oriented porous materials (normal and parallel to the incoming flow, respectively). In particular, $K_{\perp}$  is systematically lower than $K_{\|}$  and with all weighting techniques lie in-between, with the volume-weighting method granting the best accuracy on a porosity range $\varepsilon=0.5-0.9$  of interest. The model was extended in [8] to account for the porosity reduction due to both GDL compression and ionomer addition, leading to a permeability reduction. Tomadakis and Robertson [14] used an analogy between electrical current and fluid flow and proposed a permeability model again function of $d^{2}$. A further confirmation of the functional relationship with $d^{2}$  is present in the Van Doormaal and Pharoah [15] study and reinforced in Tamayol et al. [7] by the extension of the equally-spaced structure to three-dimensional materials. Gostick et al. [16] measured several GDL samples in both normal and parallel configurations, fitting the measured permeability with the Carman-Kozeny equation and highlighting the ever-present deviation between samples, which could be useful to outline a permeability sensitivity (average deviation 4.7% for normal permeability) applicable to macro-homogeneous models. A similar testing campaign was carried out by Gurau et al. [17], where $K_{\perp}$ and $K_{\|}$ of various GDL samples confirmed the higher in-plane permeability by a factor larger than 2. Feser et al. [18] measured $K_{\|}$  for several GDL (woven, non-woven and carbon type) using a radial flow technique, indicating $K_{\|} \sim 1 \times 10^{-11} m^{2}$  for all GDL, although sensible type-to-type variation is to be considered. Similar conclusions are found in [9, 19], although the non-systematic indication of  $K_{\perp}$ and $K_{\|}$ is still present in this kind of studies.

A similar modelling challenge for GDL regards heat transfer, developing in both solid and fluid phases in a porous material. Their thermal conductivity properties constitute a network of parallel resistances, although the higher conductivity of carbon shifts the dominant role to solid-phase heat conduction. Porous materials are characterized to provide a bundled “effective” thermal conductivity $\lambda$  [W m^-1 K^-1], which can be measured or analytically derived, and which determines the rate of heat dissipation from the active CL surface to cooled BPP.

The layered fibrous nature of GDL has a similar effect on thermal conductivity than that discussed for flow resistance, as in in-plane conductivity ($\lambda_{\|}$) is much higher than the through-plane value ($\lambda_{\perp}$). This is motivated by the analogy between momentum and heat transport. Heat diffuses easily along the carbon fibers direction, whereas the intra-fiber heat transfer is mediated by pressed localized contacts, hence explaining the difference in the preferential heat dissipation pathway. This motivates the importance of using orthotropic thermal conductivity for GDL in multi-dimensional models if the correct heat dissipation contributions are to be predicted, using correlation from experiments or analytical models.

In Sadeghi et al. [20] and Sadeghi et al. [21] a basic cell unit was proposed to analytically determine $\lambda_{\perp}$q of GDL, including multiple heat transfer mechanisms: fiber-to-fiber Hertzian contact for compressed media, bulk solid phase conductivity and gas phase heat conduction. The thermal resistance network model confirmed the dominant role of the inter-fibers contact resistance on the other heat transfer mechanisms and its reduction for high compressive loads. The model was extended in Sadeghi et al. [22] including experimentally measured probability functions of aspect ratio and contact angle distribution for two GDL, confirming that the random fibers orientation can be approximated in a sample unit cell model by an orthogonal arrangement. A practical expression for $\lambda$  is the Dagan equation [23], based on solid/fluid thermal conductivities and material porosity $\mathcal{E}$.

Moreover, an accurate heat transfer modelling must include the Thermal Contact Resistance (TCR [m² K W^-1]) present at the interfaces of GDL with membrane and BPP. The heat transfer loss between separate conductive parts was analytically studied by Sadeghi et al. [24] under various compressive loads. They underlined the superficial nature of TCR, different in concept than the global thermal conductivity of the bulk porous medium. Test acquisition and model prediction based on a parallel network of contact spots agreed in indicating the inter-component contact resistance contributing by 65 - 90% to the overall thermal resistance, mainly due to the limited effective contact area (approx. 1% of the nominal cross-sectional area). Experimental measurements by Khandelwal and Mench [25] confirmed the relevance of TCR in a 1D heat transfer testing arrangement.

In this work the impact of using complex properties for GDL characterization is studied using 3D-CFD models. In the first part the impact of orthotropic K is assessed with respect to standard isotropic models, with focus on inter-channel by-pass flow rate in serpentine-type gas distributors and of the effect on reducing the flow-inactive under-rib regions. In the second part, heat transfer is studied applying orthotropic $\lambda$  for GDL and including TCR at the solid-solid GDL pressed contacts. The effect on global results and on field distribution is critically discussed, and guidelines on modelling practices for these properties are provided.

In this Section the results of the two analyses will be presented, and the practical implications of using complex material properties on the large-scale PEMFC results will be critically discussed.

3.1 Effect of GDL permeability

The analysis of GDL permeability is carried out by imposing three flow regimes (Re=73, 129 and 473, respectively) and by changing the in-plane GDL permeability on three levels ($K_{\|}=5 \times 10^{-12} \mathrm{~m}^{2}, 1 \times 10^{-11} \mathrm{~m}^{2}$  and $5 \times 10^{-11} m^{2}$, respectively), for a total of 9 cases listed in Table 1. These values are typical of GDL with $\varepsilon=0.7$  porosity, as the one considered in the model.

Table 1. List of cases for GDL permeability analyses

In Figure 2 the flow rate distribution for each case is reported, representing the mass flow rate percentage at the outlet channel section. If no flow rate by-pass was present (non-porous GDL) it would read 100%, and for all cases the complement to 100% is the flow rate fraction escaping the channel via the permeable GDL and directly flowing to adjacent gas channels. The practical implication of this is that the under-rib region benefits from an ever more relevant flow convection, facilitating reactants delivery and product water removal. This corresponds to improving the electrochemical activity of a typically under-used catalyst layer area. From a modelling perspective, this better replicates the plane-oriented carbon layers composing GDL, whose random orientation in the cell plane induces a similarly high flow permeability in the plane directions.

A trend of increasing by-pass flow rate for higher in-plane permeability is observed for all the tested flow regimes, indicating a 25-30% of total flow rate escaping through the porous GDL, despite its relevance marginally reduces for high Re conditions. The use of isotropic properties $\left(K_{\|}=K_{\perp}=5 \times 10^{-12} \mathrm{~m}^{2}\right)$ would have largely underpredicted this effect, with more than 90% of the fluid flowing directly to the channel outlet, thus underestimating the entity of the flow in porous diffusion layers.

In Figure 3 the velocity field is represented on GC and GDL sections, visually confirming the entity of higher flow by-pass for cases with higher $K_{\|}$  and highlighting the non-negligible flow in under-rib areas. The well-established literature basis on this aspect indicates that orthotropic permeability with $K_{\|} \gg K_{\perp}$ is to be applied in multi-dimensional models, and that Case 3-9 are considered a more accurate representation of the actual flow behaviour.

2.png

Figure 2. Flow rate percentage at GC outlet, with complement to 100% being the by-pass flow rate through GDL

3.png

Figure 3. Velocity field on GC and GDL sections for Case 1-3 (top) and Case7-9 (bottom)

One of the major aspects to be observed is the degree of flow interaction with the CL surface, being representative of the reactants supply effectiveness of the cell assembly. Despite in this study electrochemical reactions are not included, the efficiency of mass transport and its uniformity on the CL surface are analysed though the Wall Shear Stress (WSS), defined as in Eq. (2). This allows to spatially identify areas of active transport of reactants by near-wall flow convection, as well as inactive regions whose only transport mechanism is wall-normal flow diffusion. The WSS and its spatial distribution are a consequence of the GC and GDL design, hence they can be relevantly altered with different distributors or GDL designs.

$\tau_{w[P a]}=\left.\mu \frac{\partial U_{i}}{\partial x_{j}}\right|_{x_{j}=0}$     (2)

It can be observed in Figure 4 that for increasing in-plane permeability, higher WSS is predicted in under-rib regions, especially in proximity of the inlet section. This indicates them as actively reacting areas despite being far from the reactants-delivering channels. It also underlines how the use of isotropic permeability would largely underestimate this positive effect.

4.png

Figure 4. WSS field on CL for Case 1-3 (top) and Case 7-9 (bottom)

Quantitative measurements of WSS are reported in Figure 5, showing both the mean WSS value on the active surface and its spatial uniformity ($S U_{\tau_{w}}$), defined as in Eq. (3). It is used as a uniformity index for the WSS and considers the surface average WSS on the CL ($\overline{\tau_{w}}$) and the local value $\tau_{w, i}$  on the i-th face having area $A_{i}$. Values of $S U_{\tau_{w}}$  as close as possible to unity are desired, indicating a uniform wall shear stress on the entire CL surface and maximizing the catalyst usage. It is observed that for higher in-plane permeability, both the mean WSS and its uniformity increase, further emphasising the impact of using complex realistic-like properties for GDL on the reactants supply for reactions.

$S U_{\tau_{w}[-]}=1-\frac{\sum_{i}\left|\tau_{w, i}-\overline{\tau_{w}}\right| \cdot A_{i}}{2 \cdot\left|\overline{\tau_{w}}\right| \cdot \sum_{i}^{m} A_{i}}$     (3)

5.png

Figure 5. Surface average WSS at CL (left scale) and its uniformity index (right scale) as a function of in-plane permeability for the 9 cases

3.2 Effect of GDL thermal conductivity and contact resistance

In the second part of the study, the focus moves to the heat dissipation through conductive solid parts. The major heat sources in PEMFC are electrochemical reactions and ohmic resistance, with these last i) considered of reduced entity in highly electrical-conductive materials and ii) not included in this study. Therefore, to mimic electrochemically generated heat, a typical temperature profile is applied on the opposite side of the membrane (i.e. at cathodic catalyst layer, CCL). Thermal conductivity for the membrane is $\lambda_{\text {mem }}=0.5$ W m^-1 K^-1, while for the graphite BPP a $\lambda_{\text {mem }}=20$  W m^-1 K^-1value is used. For both membrane and BPP an isotropic value is justified for $\lambda$, given their different structure than diffusion media.  As for the GDL, the implication of the fiber-oriented nature discussed for flow permeability maintains a similar validity also for heat transfer, as heat diffusion along carbon fibers is much more effective than inter-fibers transport, which is mediated by isolated contact areas. It is to be noted that heat transfer proceeds also along the fluid continuum, although its contribution is marginal compared to the solid-phase process. Therefore, the effect of the fiber-oriented GDL fabrication is considered moving from an equal though-/in-plane conductivity (e.g. isotropic $\lambda$) to a one order of magnitude higher in-plane conductivity $\left(\lambda_{\|} \gg \lambda_{\perp}\right)$.

Moreover, the multi-part composition of the cell assembly has implications on heat transfer routes, i.e. GDL is a separate component than membrane and BPP. Parts are assembled and clamped to a sealing pressure (typically $p_{\text {seal }}=1.0$  Mpa) which facilitates heat conduction between parts, reducing the ever-present TCR between the pressed parts. This is another overlooked aspect in 3D-CFD simulations, whose effect is included and discussed in this study. The list of simulations is reported in Table 2.

A typical temperature profile of CCL section is derived from [27], aiming at obtaining no heating far from the channel region (i.e. T=350 K) and a maximum temperature increase along the channel axis of 100 K, using a modified Gaussian distribution as in Eq. (4). Despite this probably overestimates the reaction heating, it is useful to highlight the different heat dissipation mechanisms to the BPP endplate, which is maintained at 350 K to replicate the presence of a cooling circuit. In Figure 6 the resulting temperature distribution imposed at CCL is illustrated.

Table 2. Simulations list for GDL heat transfer analysis

$\begin{gathered}

T(x)_{[K]}=350+0.05 \frac{1}{2 \times 10^{-4} \sqrt{2 \pi}} \exp [-0.5 \\

\left.\cdot\left(\frac{x}{2 \times 10^{-4}}\right)^{2}\right]

\end{gathered}$     (4)

The mean (volume averaged) part temperature is reported for the three cases in Figure 7. The use of realistic orthotropic $\lambda$ values for the GDL (i.e. from Case 1 to Case 2) leads to a visible decrease in the GDL mean temperature, motivated by the facilitated lateral heat transfer pathway to solid BPP. This is testified also by the increase in the thermal flux measured at the BPP endplate (+227%). This is expected as the heat transfer pathway among components follows the membrane-GDL-BPP path, hence constituting a series network. The introduction of the TCR (i.e. from Case 2 to Case 3) does not significantly alter the mean temperature of any component, although it hampers the heat dissipation to the BPP endplate, which however is higher than in Case 1 (+58%) where both accurate $\lambda$ definition and TCR were omitted. For all the cases the membrane temperature remains similar due to the boundary condition adopted (temperature distribution applied to the membrane-adjacent CCL).

6.png

Figure 6. Temperature field imposed at CCL

Moving to the temperature distribution reported in Figure 8 for the three cases, it further highlights the relevance of both using non-isotropic $\lambda$  values for porous GDL and of the inclusion of TCR at the interfaces of pressed GDL with membrane/BPP.

7.png

Figure 7. Mean (volume average) part temperature (left graph) and surface average heat flux at BPP endplate (right graph) for all the cases

Moving from Case 1 to Case 2 (change to orthotropic $\lambda$), it is noted how the higher $\lambda_{\|}$ value favours the in-plane heat transport, thus reducing thermal gradients on the GDL and replicating the effect of the plane-oriented carbon fibres. When TCR is added (Case 3), the heat transfer from/to the GDL is obstructed: the result is that an even more uniform heat distribution in the GDL is observed.

8.png

Figure 8. Temperature distribution on GC and GDL sections for Cases 1 (top), 2 (middle) and 3 (bottom)

The analysis of the largely different temperature distribution, combined with the effect on mean parts temperature, and of the dissipated heat flux form the conclusion that it is of utmost importance to include accurate and experimentally-derived non-isotropic material properties for GDL, as typical of series-arranged thermal resistances.