Material and Morphological Heat Transfer Properties of Fuel Cell Porous Transport Layers

The complex structure of the porous transport layer (PTLs), commonly referred to as gas diffusion layers (GDL), makes direct measurement of thermal transport properties difficult. Typically, effective thermal conductivity values are obtained experimentally and computational models are developed to replicate the measurements. This approach is sufficient for generalized one-dimensional models of fuel cell performance, but is not sufficient to predict reactant maldistribution due to the presence of liquid water in the PTL. Further, the anisotropy in the PTL structure results in dissimilar properties for through-plane and in-plane transport. In this work an analytical model is developed to extract material properties from effective thermal conductivity measurements. The model accounts for the change in the effective thermal conductivity as a function of compression. Geometric changes due to compression are separated from morphological changes using a compression modulation function; the shape of which is indicative of how fiber-to-fiber contact morphology changes with respect to strain. Material and morphological properties are determined using effective thermal conductivity data for three commercial PTLs. The properties also enable prediction of in-plane heat transfer using through-plane empirical data and also provide insight into morphological changes due to compression. © The Author(s) 2017. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0341713jes] All rights reserved.

A porous transport layer (PTL), more commonly referred to as a gas diffusion layer (GDL), has multiple purposes in a fuel cell including distribution of reactant gases to the catalyst layer, facilitation of product water removal from the catalyst layer, electronic and thermal conduction to/from the catalyst layer, and redistribution of compressive stresses across the membrane electrode assembly (MEA). A typical PTL consists of carbon fibers, either in fibrous or a non-woven morphology, along with binding agents and non-wetting agents such as polytetrafluoroethylene (PTFE). Examples of a non-woven and a fibrous PTL are shown in Figure 1. PTL designation is used herein because the transport phenomena of interest is not restricted to diffusion of reactant gases, but is focused heat transfer.
Heat transfer across a dry PTL is primarily due to conduction. When liquid water is present, then additional heat may be transferred from the catalyst layer via evaporation and conduction via the liquid water occupying the void space. Either mode of heat transfer, conduction or evaporation, may be the dominant effect depending upon the operating conditions. 1 The location of an evaporation front within a PTL depends upon the temperature and temperature gradient within the PTL. The same can be said for water condensation within the PTL as demonstrated by Caulk and Baker. 2 A highly conductive PTL may increase the conduction of heat, but result in widely distributed condensation that can block reactant transport. A less conductive PTL may decrease the conduction of heat, but result in a sharp temperature gradient that generates a thin, uniform layer of water which also serves as a reactant transport barrier. 3 An accurate prediction of the temperature field within a PTL with and without the presence of liquid water could enable PTL material and morphology to be optimized for specific operating conditions.
Compression of a PTL may significantly change thermal transport and subsequently the temperature field within the PTL. In a proton exchange membrane fuel cell (PEMFC), this compression is non-uniform with high compression between the the 'lands' of the reactant flow field and the catalyst layer and no compression between the reactant 'channels' and the catalyst layer. As a result temperatures and temperature gradients through the PTL thickness (through-plane) will vary from land to channel locations. This non-uniform compression complicates the prediction of internal PTL temperatures and the subsequent prediction of liquid water location. An additional complication arises in that the PTL exhibits severe anisotropy so that heat conduction through the thickness (through-plane) is much less than heat conduction along the span-wise direction (in-plane). Experimental studies have been conducted to determine effective through-plane thermal conductivities of PTLs. [4][5][6][7][8][9][10] The term effective thermal conductivity is used because these values are not strictly material constants. PTLs are porous and chemically inhomogeneous and thermal conductivity is based on temperatures measured at opposite surfaces of the PTL. There is no realistic method for in-situ or ex-situ measurement of temperatures internal to a PTL. Effective thermal conductivity is nearly always based on a linear temperature profile through the PTL.
Attempts to relax this assumption and determine a more realistic temperature distribution have relied upon idealized heat transfer structures with the intent of determining material properties, as opposed effective properties. Radhakrishnan et al. 9 performed an analysis of these models and proposed a fractal model to predict through-plane heat transfer to extract an effective thermal conductivity. Similarly, Sadeghifar et al. 11 proposed a statistical unit cell model for estimating through-plane thermal conductivity. The challenge inherent in these efforts is that the idealized heat transfer structures and the subsequent temperature and property data are only valid for thermal conduction of dry PTLs in the through-plane direction.
Alternatively, heat and mass transport within a PTL has been modeled using geometrically reconstructed PTL computational domains [12][13][14][15] and with geometric abstractions such as pore-network models. [16][17][18][19][20] The complex stochastic structure of the PTL makes it difficult and computationally expensive to directly model the thermal and mass transport using techniques such as computational fluid dynamics applied to a geometrically reconstructed domain. Pore-network models are a relatively simple computational approach, but require accurate pore-level material constants. In principle, thermal transport ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 207.241.231.83 Downloaded on 2018-07-20 to IP and the resulting internal temperature field can be accurately predicted in dry or wet using different approaches provided accurate material constants are available.
The objective of this work is to develop a methodology for extracting accurate material and morphological properties using effective thermal conductivity data. These properties should remain invariant with compression or in the presence of liquid water.

Solid-Cylinder Material Model
The methodology for transforming effective thermal conductivity into material and morphological properties involves building an equivalent thermal resistance to the PTL. To determine these material properties (as opposed to effective properties), the solid portion of the PTL is analyzed in isolation from the overall volume of the PTL sample. This equivalence is referred to as a solid-cylinder material model (SCMM). The solid phase of the PTL is treated as an elastic cylinder with frictionless surfaces so that the effect of geometric changes due to strain (compression) can be separated from material constants and morphological changes (contact resistances).
Nearly all through-plane thermal conductivity measurements on PTLs are based on compression of a stack of 3 to 10 PTLs with insulated peripheral boundaries. Thermal transport property data is extracted using one-dimensional, steady-state conduction through a planar solid with constant thermal conductivity and an assumed linear temperature profile.
In this simplified Fourier's Law, A PTL is the projected area of the PTL stack in the direction of heat conduction, H sample is the thickness of the stack, and T s is the measured surface temperature on either side of the PTL stack. This is often rearranged into an expression of thermal resistance becauseQ and T are the measured quantities. [2] where N is the number of PTLs used and H PTL is the thickness of a single PTL. When considering the thermal resistance of a single PTL under compression, illustrated in Figure 2a, the heat transfer path length is reduced by the change in thickness. [3] where χ = H/H PTL is strain and H PTL is the uncompressed thickness of a single PTL. The projected area of the PTL normal to the direction of heat transfer, A PTL , is considered to remain constant during compression for all data reviewed. Thus, the thermal resistance of a single PTL subject to strain is The exponent is 1 when k eff is strain corrected and 0 when either the strain is not reported or the uncompressed PTL thickness is used in determining k eff . This distinction is necessary because reported values of effective thermal conductivity may be based on the actual thickness of the compressed sample: (1 − χ)H PTL , 9,10 but not always. Effective thermal conductivity may also be reported in relation to the applied stress with an accompanying stress-strain correlation 7,21 Other data is correlated to applied stress, but the original uncompressed sample thickness is used in Equation 2 or the compressed thickness is not reported. 4,5,22 In these instances the primary purpose of the measurements may be to model the overall heat transfer in an operating stack and not examine localized effects in the PTL. Thus, the thickness of the sample used in determining the k eff in Equation 2 may or may not be corrected for strain. In order to separate material and morphological thermal transport properties from PTL geometry and effective thermal conductivity, the PTL thermal resistance R PTL is set equal R cyl which is comprised of two resistances in series. One of the resistances is due to the solid material of the PTL and the second is a contact resistance between fibers. This analog, illustrated in Figure 2b, is referred to as a solid-cylinder material model (SCMM) because only the solid phase of the PTL is considered. The thermal resistance of the SCMM is: where the subscripts s and c indicate solid and contact, respectively. The contact resistance includes a compression modulation function C(χ) that accounts for non-geometric effects of compression on the PTL morphology. In principle, the contact area and length of the two cylinders can be arbitrarily set and the value of the thermal conductivities, k s and k c , adjusted to match effective thermal conductivities. We have adopted the following convention for establishing contact areas and thicknesses. The areas of the solid and contact resistances are set equal to one another and the initial length of the resistances are set to the uncompressed thickness of the original PTL sample. The cylinder representing the solid resistance R s has physical dimensions. The cylinder representing the contact resistance R c has no physical dimension, but is modeled as though it does. With this convention the contact resistance is a virtual cylinder that does not occupy any physical space even though it is modeled as having a length H cyl . Compression of the solid-cylinder material model is treated as elastic with frictionless interfaces. When compressed the contact area for heat transfer increases while the path length decreases. Conservation of volume for the elastically strained solid cylinder requires: As a result the resistance of the solid-cylinder material model when accounting for strain becomes where H cyl and A cyl are the uncompressed length and area, respectively. Geometric changes in thermal resistance due to compression are accounted for in (1 − χ) 2 . The intrinsic thermal conductivity of the solid phase of the PTL is accounted for in k s , which is determined using mass fractions of carbon and PTFE. The intrinsic resistance to heat conduction due to fiber-to-fiber contact is accounted for in 1/k c . And variation in thermal conductivity due to compression-induced morphological changes, as opposed to geometric changes, is accounted for in the compression modulation function C(χ). The contact resistance components k c and C(χ) are determined by equating the experimentally measured resistance R PTL (Equation 4) to the solid cylinder material model resistance R cyl (Equation 7). The geometry of the PTL sample that is experimentally tested is as arbitrary as is the choice of dimensions for the solid cylinder material model. In order to couple the resistances the thicknesses of the PTL and SCMM are set equal, H PTL = H cyl , with the caveat that the contact resistance thickness is not physical. The areas of the PTL and the SCMM, however, are not equal. The relationship between the projected area relative to the direction of heat transfer for the PTL and cylinder is determined by considering only the solid phase of the PTL. The solid volume of the uncompressed PTL and SCCM is [8] where B is the bulk porosity of the PTL. Thus, the relationship between effective heat conduction area of the test sample and SCCM is Equating the PTL and SCCM resistances using this geometric convention provides a relationship between the reported valued of effective thermal conductivity and the material and morphological properties, k c and C(χ), respectively.
Changes in the rate of heat transfer due to strain-induced geometric changes are accounted for in the (1 − χ) 1,2 coefficient. The exponent varies depending on whether reported values of k eff are strain corrected (1) or not strain corrected (2). Non-geometric strain-induced heat transfer effects are accounted for in the compression modulation function C(χ). For convenience, resistivity (inverse of conductivity) is used when fitting of the solid-cylinder material model resistance R cyl to empirical resistance R PTL .

Solid Resistivity
The solid resistivity is determined using thermal conductivity values of the constitutive PTL materials (carbon and PTFE) in proportion to the mass fractions. [11] where M PTFE is the weight percent of PTFE. The thermal conductivity of carbon used in PTLs is taken to be that of dense carbon at 120 W/m K 5 while the value for PTFE is k PTFE = 0.245 W/mK. 23 Addition of PTFE reduces the solid phase thermal conductivity k s , but may increase the effective thermal conductivity k eff of the PTL due to the reduction in porosity per Equation 9. Thermal resistivity is the inverse of the thermal conductivity; ρ s = 1/k s .

Compression Modulation Function
At zero strain the PTL is uncompressed and the compression modulation function is equal to one which provides for a maximum contact resistance in the uncompressed state. At the other extreme the contact resistance component in Equation 10 should diminish to zero as the PTL is compressed to some maximum strain; lim χ→χmax C(χ)ρ c → 0 [12] In practice this maximum strain can never be achieved and so this value represents a limiting condition on the compression modulation function C(χ). In the limit of maximum compression only the solid phase of the PTL will remain so the maximum strain is equivalent to the bulk porosity of the PTL, χ max = B .

Contact Resistivity
The value of ρ c and the compression modulation function C(χ) are determined using empirical data and the limiting conditions for C(χ). Fitting of the contact resistivity to the empirical data begins by rearranging Equation 10 in order to isolate the contact resistivity components. [14] where ρ * = ρ/ρ s . The left side of the Equation 14 is the empirical data with adjustments for porosity and strain. The right side is the magnitude of contact resistivity scaled by the solid resistivity and the compression modulation function. The fitting function used to match the dimensionless effective resistivity that is porosity and strain adjusted (left side of Equation 14) includes a linear component and a non-linear component. The linear component is used to fit to data at small strains (≤ 20%), which has been observed to follow a linear trend for many different types of PTLs. 3 The non-linear component enforces the limiting condition on the contact resistivity C(χ) ρ * c , which decays to zero at maximum strain. For convenience, strain is scaled using the bulk porosity to simplify the fitting of the empirical data to the dimensionless contact resistivity.
The fitting function used herein is of the form g(ξ) = a 0 + a 1 ξ − (a 0 + a 1 )ξ p [16] which is fitted to the the scaled test data, left side of Equation 14. This approach allows for one additional data point not derived from empirical effective thermal conductivities; that is, g(1) = 0. The magnitude of the dimensionless contact resistivity at ξ = 0 is ρ * c = a 0 because C(0) = 1. Therefore, the contact resistivity is ρ c = ρ s a 0 [17] and the compression modulation function becomes where b 1 = a 1 /a 0 . Fitting the scaled test data to g(ξ) yields values for a 0 , a 1 , and p.

Application to Commercial PTLs
The solid-cylinder material model is applied to three commercial PTLs with published effective thermal conductivity data; Toray TGP 060, 10 Freudenberg H2315, 10 and Mitsubishi Rayon Corp. MRC 105. 21 Data for two different PTFE weight percents are available for the Freudenberg H2315 PTL. For Toray TGP 060 and Freudenberg H2315 the effective thermal conductivity is reported as a function of PTL compression thickness, from which the strain can be computed. The effective thermal conductivity for the Mitsubishi MRC 105 is reported as a function of compressive stress along with a correlation between stress and strain. The original effective thermal conductivity data for each PTL are listed in Table I and plotted in Figure 3 against strain. The effective thermal conductivity for Toray TGP 060 is approximately three to four times greater than that for the Freudenberg H2315 and Mitsubishi MRC 105; and also exhibits a greater rate of change with respect to compressive strain.
In order to transform effective thermal conductivity data into the fitting function g(ξ) per Equation 14, the porosity, strain, and value of ρ s are required and are listed in Table II. For these three PTL data sets, the effective thermal conductivity is tabulated against measured or correlated strain. Therefore, in transforming the effective thermal conductivity into g(ξ), the exponent in Equation 14 is equal to one. The porosity of the Freudenberg H2315 with 10% PTFE (wt) was not reported, but can be computed using the density ratio of carbon to PTFE, which can vary depending on the source. The change in porosity due to addition of PTFE based on an unchanged total PTL volume is PTFE [19]  25 Thus, for the Freudenberg H2315 PTL with 10% PTFE (wt), the porosity is reduced to B = 0.747 as compared to the uncoated sample with B = 0.77. The procedure for determining the fitting function constants a 0 , a 1 , and p requires transforming the effective thermal conductivities into the non-dimensional effective resistivity that is strain and porosity adjusted; i.e., the left side of Equation 14. The transformed data is tabulated in Table I as discrete values of g(ξ) and the limiting condition g(1) = 0 is added to the tabulation. A regression analysis using the form in Equation 16 is performed on each PTL data set. The shape of the fitting function, shown in Figure 4, is unique to each PTL. The quality of the fit to data for all three PTLs is very good with correlation coefficients greater than 0.999. Care should be exercised in considering the quality of these fits, however. The effective thermal conductivity data set is sparse with relatively large uncertainty (approximately 10% to 15%) and with large separation in strain from the limiting condition at ξ = 1. As a result, if the effective thermal conductivity exhibits little change with small strain, then the quality of the fit is insensitive to the non-linear component characterized by the exponent p. This can be observed in the fit to the two Freudenberg PTLs. For the 0% PTFE (wt) data, the value of p was selected to be 3.5. Whereas for the 10% PTFE (wt) data, the value of p was selected to be 10. Both fits match both sets of data, indicating that   the choice of p for the Freudenberg PTLs is somewhat arbitrary. The insensitivity in the non-linear component of the Freudenberg PTLs is due to an absence of morphological changes at small strains. The magnitude of the material constants indicate that the throughplane thermal transport in PTLs primarily depends largely on the contact resistivity, ρ c . The relative effects of solid vs contact on throughplane thermal transport can be more easily seen in the respective thermal conductivities k s and k c . The value of k s is 165 to 217 times greater than k c depending on the PTL. Thus, thermal transport is limited by the contact resistance.
The shape of the fitting function and subsequent compression modulation function provides indication of what happens to the contact morphology due to compression. The non-geometric compression effects for all three PTLs are compared using the compression modulation functions, C(ξ), shown in Figure 4d. For the Toray TGP 060 PTL data there is a significant change in effective thermal conductivity that cannot be accounted for by geometric changes due to strain. The steep decrease in C(ξ) indicates substantial changes in contact morphology has occurred at small strain. Whereas for the Freudenberg H2315 data there is negligible change in C(ξ) small strain; the variation is within the uncertainty of the original test data. This indicates that the change in effective thermal conductivity due to compression is only due to geometric changes associated with strain. A moderate change in contact morphology occurs for the Mitsubishi MRC 105 and this change appears to occur over a large range of strain.
Using the material and morphological constants listed in Table  II the through-plane effective thermal conductivity can be predicted. Figure 5 illustrates the predicted k eff vs normalized strain for all three PTLs. The inset in each plot is a magnified perspective for strains less than 0.25 for which the trend is approximately linear as observed by Roos. 3 The prediction matches the empirical data as expected. The utility in this methodology, as opposed to linear fits to effective thermal conductivity with strain, is in application to multiunidimensional heat and mass transfer in a PTL with the presence of liquid water. The role of thermal conduction through the PTL material can be isolated from that through liquid water in the pore phase. In addition, this methodology allows for prediction of in-plane conduction heat with the PTL under compression using only through-plane data.

In-Plane Effective Thermal Conductivity
Manufacturing of the PTLs results in the orientation of fibers along the span-wise (in-plane) direction, which results in fewer fiber-to-fiber contacts per unit length. Subsequently, effective thermal conductivities in the in-plane direction should be greater than for through-plane. Measurements of in-plane effective thermal conductivity are limited, but values of in-plane k eff range from 3 to 21 W/m K, 22,27-29 which is 10 to 40 times greater than typical values for through-plane k eff .
The ratio of the fiber diameter and the projected in-plane length of the fiber provides a scaling parameter proportional to the number of fiber-to-fiber contacts in each direction. Multiplying the contact resistance in the SCMM by this ratio reduces the overall thermal resistance in the in-plane direction. [20] In terms of effective thermal conductivity, the change includes the ratio of fiber diameter to in-plane length as a multiplier on the contact resistivity.
This expression is the same as Equation 9, but with the contact resistivity reduced by the ratio of fiber diameter to length. The strain correction term, which has an exponent of 1 or 2 in Equation 9 depending on how empirical data is reported, is simplified to an exponent of 1 because Equation 21 is predictive. The mean fiber diameter, d f , is determined by analysis of SEM images and optical micrographs. The values of d f for each of the PTLs is listed in Table II. To identify the mean fiber lengths in the inplane direction, images of PTL surface were analyzed to identify the fibers and their orientation with respect to one of the orthogonal axis. First, the distribution of fiber lengths and orientation relative to an orthogonal direction was established using the procedure similar to that described in Parikh et al. 30 The projected fiber length in the in-plane direction, f , is found by multiplying the mean fiber length by the cosine of the mean fiber orientation. This projected length characterizes the mean, in-plane length over which heat transfer occurs without a fiber-to-fiber contact. For a non-woven PTL such as the Freudenberg H2315, the in-plane fiber length distribution is determined using the longest straight segments of individual fibers. This results in wider distributions of in-plane fiber lengths and fiber orientation as compared to the fibrous TGP 060 and MRC 105 PTLs; though the mean is still easily identified. The ratio f /d f is 175, 143.6, and 158.6 for TGP 060, H2315, and MRC 105, respectively.
The effective in-plane thermal conductivity for the three PTLs with compression is shown in Figure 6. The same PTL-specific values of k s , k c , b 1 and p were used to calculate through-plane k eff (Equation 9) and in-plane k eff (Equation 21). At strains less than 0.25 the inplane effective thermal conductivity for these three PTLs is between 10 and 25 W/m K, which is 40 to 100 times greater than the throughplane effective thermal conductivity. At maximum strain, in-plane and through-plane conductivities become equal. Comprehensive experimental studies for in-plane effective thermal conductivity are lacking, especially for the three commercial PTLs considered. Nonetheless, the predicted values of in-plane effective thermal conductivity for these PTLs are comparable to reported values for similar PTLs, which vary from 3 to 21 W/m K for similar PTLs. 22,27-29 The effect of PTFE on the in-plane heat transfer remains uncertain. Sadeghi et al. 27 observed a small increase in in-plane thermal conductivity from 17.39 W/m K to 17.81 W/m K with an increase in PTFE content from 5% to 30%, respectively. In contrast, Zamel et al. 31 observed a sharp decrease of approximately 11 W/m K at 60 to 70 • C with addition of 5% PTFE to Toray TGP 120. Teertstra et al. 28 reported a decrease in in-plane thermal conductivity with increasing PTFE content for SpectraCarb, while Alhazmi et al. 29 observed large increase for SGL with increasing PTFE content. These discrepancies may be due to changes in fiber-to-fiber contact morphology in PTLs with, geometric changes with strain, changes in porosity, or some combination of each. The shape of the compression modulation function for these PTLs with varying PTFE content may help elucidate these relative importance of these various effects.

Conclusions
A solid-cylinder model has been developed extracting material and morphological thermal transport properties of PTLs based on throughplane effective thermal conductivity data. The effect of porosity and geometric changes due to strain are isolated from the intrinsic thermal conductivity of the solid material and changed in contact morphology due to strain. As a result, the material properties along with the compression modulation function can be used to predict in-plane and through-plane heat transfer and temperature distributions in a PTL that is non-uniformly compressed.
Through-plane heat transfer is primarily limited by the fiber-tofiber contact resistance. When compressed, both the geometry of the PTL and the morphology of the fiber-to-fiber contact may change. The shape of the compression modulation function C(χ) can be used to assess the role of geometry and morphology on heat transfer changed due to compression. Data from the two Freudenberg PTLs with differing PTFE content show no measurable changes in contact morphology at strains below 0.25. The measured differences in k eff with strain are most likely due to geometric changes in the PTL. Changes in contact morphology do not occur until the strain reaches values that cannot be realistically obtained. In contrast, the Toray TGP 060 PTL data indicates that even small application of strain results in rapid changes in contact morphology, with the rate of change decaying at large strains. In between these two extremes, the Mitsubishi MRC 105 exhibits a relatively linear change in contact morphology across the full range of strain.