• No results found

A NUMERICAL MODEL OF HEAT- AND MASS TRANSFER IN POLYMER ELECTROLYTE FUEL CELLS : A two-dimensional 1+1D approach to solve the steady-state temperature- and mass- distributions

N/A
N/A
Protected

Academic year: 2021

Share "A NUMERICAL MODEL OF HEAT- AND MASS TRANSFER IN POLYMER ELECTROLYTE FUEL CELLS : A two-dimensional 1+1D approach to solve the steady-state temperature- and mass- distributions"

Copied!
51
0
0

Loading.... (view fulltext now)

Full text

(1)

A NUMERICAL MODEL OF

HEAT- AND MASS TRANSFER IN

POLYMER ELECTROLYTE FUEL CELLS

A two-dimensional 1+1D approach to solve the steady-state

temperature- and mass- distributions.

EMIL SKOGLUND

(2)

ABSTRACT

Methods of solving the steady state characteristics of a node matrix equation system over a polymer electrolyte fuel cell (PEFC) were evaluated. The most suitable method, referred to as the semi-implicit method, was set up in a MATLAB program. The model covers heat transfer due to thermal diffusion throughout the layers and due to thermal advection+diffusion in the gas channels. Included mass transport processes cover only transport of water vapor and consist of the same diffusion/advection schematics as the heat transfer processes. The mass transport processes are hence Fickian diffusion throughout all the layers and diffusion+advection in the gas channels. Data regarding all the relevant properties of the layer materials were gathered to simulate these heat- and mass transfer processes.

Comparing the simulated temperature profiles obtained with the model to the temperature profiles of a previous work’s model, showed that the characteristics and behavior of the temperature profile are realistic. There were however differences between the results, but due to the number of unknown parameters in the previous work’s model it was not possible to draw conclusions regarding the accuracy of the model by comparing the results.

Comparing the simulated water concentration profiles of the model and measured values, showed that the model produced concentration characteristics that for the most part aligned well with the measurement data. The part of the fuel cell where the concentration profile did not match the measured data was the cathode side gas diffusion layer (GDL). This comparison was however performed with the assumption that relative humidity corresponds to liquid water concentration, and that this liquid water concentration is in the same range as the measured data. Because of this assumption it was not possible to determine the accuracy of the model.

Keywords: PEFC, PEM, convection, mass transport, heat transport, co-flow, temperature profile, heat transfer model, transient, steady state, water concentration profile, water vapor profile, relative humidity, through-plane cross section

(3)

PREFACE

This work was made possible thanks to Michael Striednig, at the Paul Scherrer Institute, Switzerland, and his contributions. Many thanks to you, Michael.

Andreas Michalski was also a key contributor to the project. Thank you for the help, Andi. During my time at the university of Mälardalen, I have made many friends that truly have given me the support I needed. Thank you all.

___________

Till min familj, jag älskar er väldigt mycket.

Västerås in June of 2021

(4)

CONTENT

1 INTRODUCTION ...1

1.1 Background ... 1

1.1.1 The generic PEFC ... 1

1.1.2 PSI’s evaporative cooling approach ... 2

ADVANTAGES OF THE COCHET ET. AL. EC PEFC CONCEPT ...3

1.2 Purpose/Aim ... 5

1.3 Research questions ... 5

1.4 Delimitation and model assumptions ... 5

2 METHOD ...6

3 EVALUATION GROUNDWORK ...7

3.1 Calculation method descriptions ... 7

3.1.1 Implicit scheme ... 7

3.1.2 Semi-implicit scheme ... 7

3.1.3 Steady state equation system ... 8

3.2 Evaluation of solvers and calculation methods ... 8

4 MODEL DESCRIPTION ... 11 4.1 Thermal diffusion ...11 4.2 Mass diffusion ...12 4.3 Thermal advection ...12 4.4 Mass advection ...13 5 PREVIOUS WORK ... 14

(5)

6.1.5 Species diffusivity ...19

6.2 Gas diffusion layer morphology ...20

6.2.1 Tortuosity and porosity ...21

6.2.2 Thermal conductivity ...22

6.2.3 Species diffusivity ...23

6.3 Gas channel morphology ...23

6.4 Bipolar plate morphology ...23

6.5 Catalyst layers morphology ...24

6.5.1 Density ...24

6.5.2 Specific heat capacity ...25

6.5.3 Thermal conductivity ...26

6.5.4 Mass diffusivity ...26

7 RESULTS ... 27

7.1 Inputs used in simulations...27

7.2 Comparison with simulation by Weber & Hickner ...28

7.3 Temperature distribution ...31

7.4 Water concentration distribution...32

8 DISCUSSION... 33

8.1 Temperature profiles ...33

8.2 Water concentration profiles ...34

9 CONCLUSIONS ... 35

10 SUGGESTIONS FOR FURTHER WORK ... 36

(6)

LIST OF FIGURES

Figure 1 Hydrophilic lines and water supply channel in a test cell not in operation ... 3

Figure 2 Voltage and resistivity for a test cell with and without evaporative cooling (left). Water distribution inside a test cell for the different current densities marked (1), (2) and (3) ... 4

Figure 3 EC PEFC concept sketch (left) and generic PEFC sketch (right) with gas- and coolant flows ... 4

Figure 4 Computation time for equation systems of different size for different solvers. ... 9

Figure 5 Mass- and heat transfer mechanisms present in each layer in the EC PEFC model. . 11

Figure 6 measured water content (left) and simulated temperature profiles (right) of the MEA in a PEFC at various temperatures and current densities at 60 °C cell temperature ...14

Figure 7 Molecular structure of membrane materials Nafion (a) and 3M ionomer (b) ... 15

Figure 8 Density of Nafion and 3M PFSA ionomers with various EW ...16

Figure 9 Specific heat capacity of Nafion 1100 as a function of water content... 17

Figure 10 Thermal conductivity of dry Nafion as a function of temperature at atmospheric pressure ... 18

Figure 11 Thermal conductivity of Nafion as a function of water content at atmospheric pressure and room temperature ... 18

Figure 12 Thermal conductivity of water as a function of temperature at atmospheric pressure ... 18

Figure 13 Water diffusivity as functions of average water content (right) and water content as function of relative humidity (left) in a Nafion 1110 membrane ...19

Figure 14 Through plane (left) and in-plane (right) cross section imaging of a Toray GDL. .. 20

Figure 15 Through-plane thermal resistance of Toray fiber-paper ... 22

Figure 16 In-plane conduction coefficient of TGP-H-120 GDL with various PTFE content ... 22

Figure 17 Density of a CL with 20% Pt/C loading as a function of water content ... 25

Figure 18 Specific heat capacity of a CL with 20% Pt/C loading as a function of water content ………..26

Figure 19 Velocity profiles ... 27

Figure 20 Temperature profile for the GC ribs cross section. ... 29

Figure 21 Averaged temperature profile over the BPP/GC ribs. ... 29

Figure 22 Concentration profile for the GC ribs cross section ... 30

Figure 23 Averaged concentration profile over the BPP/GC ribs. ... 30

Figure 24 Temperature distribution at various velocity profiles ... 31

(7)

LIST OF TABLES

Table 1 Generic PEFC layers and their function. ... 2

Table 2 Computation time and accuracy, steady state equation system. ... 9

Table 3 Semi-implicit scheme. ... 9

Table 4 Implicit scheme (Broyden method). ... 10

Table 5 Specific heat capacity of Nafion 1100 at lambda = 14. ... 17

Table 6 Toray carbon fiber paper TGP-H material data ...21

Table 7 Porosity and tortuosity of ELAT carbon fiber paper GDLs ...21

Table 8 Sigracell® bipolar plates material data ... 23

Table 9 Catalyst layer material properties ... 24

Table 10 Densities, weight fractions and total density of the CL. 20% Pt/C loading ... 24

Table 11 Specific heat capacity of CL at λ = 14. ... 25

Table 12 Parameters used in FC simulation ... 27

LIST OF EQUATIONS

Equation 1 Energy balance over node located in point (x). ... 7

Equation 2 Energy balance over node located in point (x). ... 8

Equation 3 Energy balance over node x at steady state. ... 8

Equation 4 Concentration of water in membrane ... 15

Equation 5 Rule-of-mixture, membrane density ...16

Equation 6 Coefficient of effective diffusivity in porous material ...19

Equation 7 Water content in Nafion 1110 membrane as function of relative humidity and temperature ...19

Equation 8 Coefficient of effective diffusivity in Nafion 1110 membrane material as a function of water content ... 20

Equation 9 Coefficient of thermal conduction for TGP-H-60 GDL material (through plane) 22 Equation 10 Coefficient of thermal conduction for TGP-H-60 GDL material (in-plane) ... 23

Equation 12 Coefficient of thermal conductivity for a 20% Pt/C CL26 Equation 13 Coefficient of effective diffusivity in porous media with spherical particle structure ... 26

(8)

NOMENCLATURE

Symbol Description Unit

𝑇 Temperature 𝐾

𝜌 Density 𝑘𝑔

𝑚3

𝑐𝑝 Specific heat capacity 𝐽

𝑘𝑔, 𝐾 𝑘 Coefficient of thermal conductivity 𝑊

𝑚, 𝐾

𝑉 Volume 𝑚3

𝐴 Area 𝑚2

𝐷 Coefficient of mass diffusivity 𝑚2

𝑠 𝜆 Water content in membrane material 𝑚𝑜𝑙 𝐻2𝑂

𝑚𝑜𝑙 𝑆𝑂3

𝐿 Distance between nodes 𝑚

𝑢 Flow velocity 𝑚

𝑠

𝑚 Mass flow 𝑘𝑔

𝑠 𝐶 Concentration of the species in question

(water) 𝑚𝑔/𝑐𝑚

3

𝜀 Porosity −

𝜏 Tortuosity −

𝑅𝐻 Relative humidity %

Subscript Description Unit

𝑎𝑚𝑏 Ambient −

𝑥 Node number / node index 𝑁𝑢𝑚𝑏𝑒𝑟

(9)

ABBREVIATIONS

Abbreviation Description BPP Bipolar plate

GC Gas channel

GDL Gas diffusion layer CCL Cathode catalyst layer ACL Anode catalyst layer

MEM Membrane (polymer electrolyte) MEA Membrane-electrode-assembly

MPL Micro-porous layer (GDL with finer pores)

FC Fuel cell

PEFC Polymer electrolyte fuel cell

EC PEFC Evaporation cooled polymer electrolyte fuel cell DTC Down-the-channel (Gas channel)

EC Evaporative cooling

PFSA perfluoro-sulfonic-acid (polymer electrolyte material) PTFE Polytetrafluoroethylene (polymer electrolyte material) ORR Oxygen reduction reaction

HOR Hydrogen oxidation reaction

DEFINITIONS

Definition Description

Ambient Refers to the surroundings of the fuel cell (e.g. room’s conditions).

Down-the-channel Direction along the gas channel.

(10)

1

INTRODUCTION

The world’s population is steadily increasing at the same time as the overall wealth is going up, eventually enabling world-wide adoption of technologies and their related power consumption. At the same time global warming is on a rise and will be further catalyzed by this increased power consumption. It can not be stressed enough that renewable energy sources are a key factor in mitigating the environmental impacts that global warming may bring and the conflicts of interest between nations that could emerge.

According to data from the International Energy Agency, 25 percent of the global CO2 emissions from fuel combustion originated in the transport sector (International Energy Agency, 2018). Electrification of the transport sector is therefore of interest for mitigating the environmental impacts of the transport sector. Batteries are agreed upon to be a suitable source of energy for electric vehicles but come with a capacity limit which is restricted by the size of the battery. Fuel cells (FC’s) are an alternative to batteries, as the energy source in electric vehicles, that utilizes hydrogen gas to produce electricity. The travel range of a FC-powered electric vehicle is essentially limited by the capacity of the hydrogen storage, which can reach a power density much greater than batteries and can therefore achieve travel ranges many times that of battery-driven electric vehicles.

1.1 Background

Fuel cells are a technology which utilizes hydrogen gas and oxygen, usually contained in air, to produce electricity. The most common type of fuel cell is the polymer electrolyte fuel cell (PEFC) which features an ion-conducting polymer membrane. The membrane is thereby the electrolyte, which separates the anode- and cathode electrodes. As hydrogen and oxygen gets introduced at the anode and cathode, respectively, a water formation reaction occurs. The membrane only conducts the hydrogen ions, so the electrons are forced through an external circuit to produce electrical power. However, heat is also produced in the reaction and cooling is required. Typically, PEFC’s are cooled by convective flows in designated channels, but by new gas diffusion layer (GDL) technology the cooling can be performed by evaporating water within the layer. The evaporative cooling (EC) approach developed by the Paul Scherrer Institute (PSI) reduces the total size of the fuel cell stack substantially, making it more suitable

(11)

(MEA). The MEA is the middle part of the cell and is enclosed on both sides by first a gas diffusion layer (GDL), and then the gas channels (GC) and the bipolar plate (BPP). The fuel cell sandwich is essentially symmetrical, except that each layer on the anode- and cathode side of the MEA have different properties according to their respective function. Table 1 shows an overview of the layer structure and function.

Table 1 Generic PEFC layers and their function.

Layer name Layer function Layer

structure An o d e sid e

BPP Conduct electrons and contain the gas channels Solid

GC Hydrogen gas flow Hollow

GDL Transports hydrogen towards the catalyst layer Porous Catalyst layer (ACL) Hydrogen oxidation reaction (HOR) Porous Membrane Conducts hydrogen ions and block electrons Porous/Solid

Catho

d

e

sid

e

Catalyst layer (CCL) Oxygen reduction reaction (ORR) Porous GDL Transports oxygen towards the catalyst layer and

product water towards gas channel.

Porous

GC Oxygen/air flow Hollow

BPP + Coolant channels Conduct electrons and contains the gas- and coolant channels

Solid

It is important to note that the GDL rarely is one single layer. The GDL usually consists of several layers with decreasing pore size closer to the MEA. Actually, not any layer except perhaps the membrane and BPP should be considered to have uniform structure throughout the layer. Another important note is that cooling of the cell is done by designated channels in the BPP, contrary to the evaporation cooled PEFC concept explained in section 1.1.2.

1.1.2 PSI’s evaporative cooling approach

The evaporation cooled PEFC (EC PEFC), features the same sandwich layout as the generic PEFC. However, the EC PEFC concept features the possibility of introducing water to some layer of the cell in order to cool it through evaporation.

A promising EC PEFC concept is described by Cochet et. al (2020) in the work “Enabling High

Power Density Fuel Cells by Evaporative Cooling with Advanced Porous Media”. The core of

this cooling concept lies in the GDL, which contains numerous 500 µm wide waterfilled rows referred to as hydrophilic lines. Water is supplied to the hydrophilic lines through water supply channels, in which water flows parallel to the hydrogen gas flow of the anode-side gas channel. The hydrophilic lines are capable of wicking the water from the supply. Very numerous and

(12)

diffusive resistance. This EC PEFC concept can thereby be optimized for different operation conditions by adjusting the number of water supply channels to achieve the desired evaporation rate (Cochet, et al., 2020).

Figure 1 Hydrophilic lines and water supply channel in a test cell not in operation. Pure nitrogen at cathode side and hydrogen at anode side gas channels.

(Cochet, et al., 2020)

Advantages of the Cochet et. al. EC PEFC concept

The membrane (electrolyte) of the PEFC consists of a polymer material which requires to be humidified in order to conduct ions effectively. In other words, the membrane needs a high relative humidity (RH). While keeping the membrane humid, flooding of the membrane must be prevented. Flooding is when saturated water clogs the porous materials of the cell, reducing its effectiveness by preventing gas diffusion and ion-conduction. Since water is already produced at the cathode as a bi-product of the energy conversion reaction and flooding is an issue, the hydrophilic lines are located in the anode-side GDL (Cochet, et al., 2020).

This set up, with water injected at the anode-side GDL has several advantages over the generic type PEFC. Firstly, it humidifies the hydrogen gas as it comes in contact with the hydrophilic lines, decreasing the need of the gas pre-humidification typically done in generic PEFCs. Secondly, the gas channels dedicated to the coolant can be dispersed which allows for up to 33 % reduction of the cell volume, as illustrated in Figure 3. However, as shown in Figure 2 there are some voltage losses present at high current densities with this EC PEFC concept due to flooding. These voltage losses, at a current density of 1.5 𝐴/𝑐𝑚2, combined with the volume

reduction increases the power density of the overall cell by up to 35 %. At lower current densities, the increase in power density can be up to 50 % (Cochet, et al., 2020). The power density increase of the complete system is however debatable as the EC PEFC system requires the addition of a water condenser, and the generic PEFC system features gas pre-humidifiers.

(13)

Figure 2 Voltage and resistivity for a test cell with and without evaporative cooling (left). Water distribution inside a test cell for the different current densities marked (1), (2) and (3). Pure nitrogen at cathode side- and pure hydrogen at anode side gas channels.

Figure 3 EC PEFC concept sketch (left) and generic PEFC sketch (right) with gas- and coolant flows. (Cochet, et al., 2020)

(14)

1.2 Purpose/Aim

The aim of this work is to produce a two-dimensional heat- and mass transfer model of a PEFC to determine the steady state temperature- and gas species distribution.

Such a model could provide insight in how the system may react to changes of the internal- and external inputs. Furthermore, it could be the base of more advanced models which allows to simulate the EC PEFC concept by the Paul Scherrer Institute.

1.3 Research questions

What is the steady state temperature distribution within the PEFC? What is the steady state concentration distribution within the PEFC? Can the simulated distributions be considered accurate?

1.4 Delimitation and model assumptions

PEFC’s are subject to degradation, due to e.g., contaminants contained in gases and other causes, which effects have not been included in the model produced for this study. Water mass transfer within the solid cell layers is modeled as pure Fickian diffusion in air, using the effective diffusion coefficients provided by empirical data. The mass- and heat transfer in the gas channels are modeled as pure convection with the assumptions of perfect laminar flow and constant velocity profile. No electrochemical processes concerning the electrical properties of the cell are included in the model. Water adsorption is not included, and the presence of liquid water is neglected. Uniform current density throughout the CCL is assumed. Thermal contact resistances were not included.

(15)

2

METHOD

Different calculation methods and different solvers were evaluated before a final version of the model was set up. The evaluation of the calculation methods showed that the methods utilizing non-linear solvers were orders of magnitude slower than the method referred to as the semi-implicit scheme which utilizes a linear solver. The semi-semi-implicit scheme did not produce any resulting errors compared to the non-linear methods that were evaluated unless very few time step iterations were made.

The semi-implicit scheme was determined to be the most suitable method of calculating the steady state temperature- and species distribution. Hence, a model of the fuel cell was set up to cover the through-plane cross section. Temperature and species concentration were represented by 𝐾 and 𝑘𝑔/𝑚3, respectively, in each node.

To get an as dynamic and easily used model as possible, a UI was set up to allow to quickly change the layers, node count and property relations included in the simulation. Densities, specific heat capacities, heat conductivities and diffusion coefficients were represented by equations relating the properties to temperature and concentration in the best scope possible. If empirical data of these relationships could not be obtained, which was especially true for the thermal conductivity, density and specific heat capacity of the different materials, constant values were used. For the membrane material, which absorbs water, a rule-of-mixture solution to calculating the density a specific heat capacity was applied. The effect of contact resistances was included in some of the empirical data obtained but remained unknown for most of the layer boundaries. Hence, the effect of thermal contact resistances was not included through the addition of separate sections between the FC’s layers.

The model consists of a node matrix spanning the through-plane cross section of the FC. Each layer has the same number of nodes in both through-plane and down-the-channel directions. Between each layer there is a single node column making up the intra-layer boundary. In these nodes the properties density, specific heat capacity, heat conductivity and mass diffusion coefficients are averaged for the two merging layers to approximate the conditions at the layer-boundary nodes. For the layer-boundary nodes between the gas channels and solid/porous layers, there is no averaging and the properties are given by the solid/porous layer material. The model calculates the heat- and mass transfer time differentials of each node for a given time step size and uses this to estimate the temperature and concentration of the node in the next time step. All the variable properties used in this calculation are based on the properties of the node in question, and not the average value between the node pair that is interacting. Each node is interacting with its adjacent nodes, making up to four node pairs that are included in the calculations of the differentials.

The results were compared to the simulated temperature profiles and measured water concentration profile by Weber & Hickner (2008).

(16)

3

EVALUATION GROUNDWORK

There are several methods that can be used to solve the steady state temperature- and species distribution. In this work, three different methods were evaluated before the most suitable was implemented. These methods are as follows, and described in section 3.1.

• Transient implicit scheme (non-linear equation system) • Transient semi-implicit scheme (linear equation system) • Steady state equation system (non-linear equation system)

Since the linear equation systems handled in this work make up singular matrices, the Moore-Penrose pseudoinverse is used as the linear solver. Two non-linear solvers were compared for solving the non-linear equation systems, the interior-point method and the Broyden method. The computation time of these solvers and their respective methods are evaluated in section 3.2.

3.1 Calculation method descriptions

See this section for a description of the calculation methods: Implicit scheme, Semi-implicit scheme, and Steady state equation system.

3.1.1 Implicit scheme

The implicit scheme (a.k.a. Euler’s backward scheme) gives a transient behavior of the system variables. The steady state distribution is solved when there is no more change to the system’s variables.

Since the heat conductivity 𝑘, specific heat 𝑐𝑝 and density 𝜌 are obtained through non-linear expressions of e.g. temperature at the current timestep, the resulting equation system is of linear nature. The equation system that is set up for all nodes can then be solved using a non-linear solver like the interior-point method or the Broyden method. The value in each node at the next time step is then calculated until there is minimal change to the system and steady state is reached.

Equation 1 Energy balance over node located in node x for the transient Implicit scheme. 𝑑𝐸𝑥 𝑡+1

(17)

using values from the previous time step, making the equation system of all nodes a linear equation system. The equation system can thereby be solved using linear solvers like the Moore-Penrose pseudoinverse. The value in each node at the next time step is then calculated until there is minimal change to the system and steady state is reached.

Equation 2 Energy balance over node located in node x for the transient Semi-implicit scheme. 𝐸𝑥𝑡+1= 𝐸𝑥𝑡+

𝑑𝐸𝑥 𝑡

𝑑𝑡

3.1.3 Steady state equation system

The steady state solution can also be solved for directly, instead of computing the transient change until convergence occurs like in the Implicit- and Semi-implicit schemes. This is done through the steady state equation system where heat balances, which are true only for the steady state condition, are set up for each node. Since many properties are calculated through non-linear functions of other properties, a non-linear solver like the interior-point method or the Broyden method is used.

Equation 3 Energy balance over node x at steady state. 𝑑𝐸𝑥

𝑑𝑡 = 0

3.2 Evaluation of solvers and calculation methods

The different solvers. Moore-Penrose pseudoinverse, interior point method and Broyden method were evaluated in regards of computation time. This was done by recording the time to complete one iteration for a 1D node matrix/equation system, ranging from a total of 3 nodes/equations to a total of 1000 nodes/equations. The results, in Figure 4, showed that the linear Moore-Penrose pseudoinverse was an order of magnitude faster than both of the non-linear interior point and Broyden solvers. Comparing the two non-non-linear solvers, the Broyden method is an order of magnitude faster than the interior point method.

(18)

Figure 4 Computation time for equation systems of different size (number of nodes = number of equations) for different solvers.

In order to determine the most suitable approach to calculating the steady state temperature- and species distribution, a comparison of the different calculation methods was made in regards of computation time and steady state solution accuracy. This comparison was made through numerous simulations of a sandwiched cell structure without flow fields (purely conductive heat transfer). The steady state criteria for the transient schemes were defined as when the summed squared change between the two last time steps is less than 10−6 K. It can

be noted that the semi-implicit method produce a much more inaccurate result when the number of iterations are few, as seen in Table 3 Computation and accuracy for various time steps, semi-implicit scheme..

Table 2 Computation time and accuracy, steady state equation system.

Solver Computation time (s) Squared error sum

Interior point 0.392 0

Broyden 0. 122 0

Table 3 Computation and accuracy for various time steps, semi-implicit scheme.

Time step (s) Computation time (s) Squared error sum Simulation time (s)

0.0001 160.1 0 20

2E-10x3+ 1E-07x2- 8E-06x + 0.0023

1E-12x4- 2E-09x3+ 3E-06x2+ 0.0003x - 0.0102

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 200 400 600 800 1000 Tim e to comp le te o n e it era tio n (s ) Number of nodes

Computation time per iteration of number of nodes for different solvers

Moore Penrose pseudoinverse

Broyden method

Interior point method

(19)

Table 4 Computation time and accuracy for various time steps, Implicit scheme (Broyden method). Time step (s) Computation time (s) Squared error sum Simulation time (s)

- - - - 0.001 446.7 0 20 0.01 63.14 0 20 0.1 5.022 0 20 1 0.757 0 20 10 0.289 0.155 20

(20)

4

MODEL DESCRIPTION

The model is set up to feature thermal diffusion, thermal advection, mass diffusion and mass advection. All layers of the FC are subject to the diffusion mechanisms, but advection mechanisms are only present in the gas channel. Processes marked in red are not included in this work.

Figure 5 Mass- and heat transfer mechanisms present in each layer in the EC PEFC model.

4.1 Thermal diffusion

The thermal diffusion, or conduction, is modeled using empirical data of the effective thermal conductivity coefficient 𝑘, specific heat capacity 𝑐𝑝, and density 𝜌.

𝑑𝐸𝑥 𝑑𝑡 = ∑ 𝐸̇𝑜𝑢𝑡𝑖𝑛 [𝐽 𝑠] → 𝑑(𝜌𝑥 𝑐𝑝𝑥 𝑇𝑥 𝑉𝑥) 𝑑𝑡 = ∑ 𝑘𝑥 𝐿 ∆𝑇 [ 𝐽 𝑠] → 𝑑𝑇𝑥 𝑑𝑡 = ∑ 𝑘𝑥 𝐿 𝑑𝑡 𝜌𝑥 𝑐𝑝𝑥∗ 𝑉𝑥 ∆𝑇 [𝑇 𝑠]

(21)

4.2 Mass diffusion

The mass diffusion is modeled as Fickian diffusion using empirical data of the effective mass diffusivity coefficient 𝐷. 𝑑𝐶𝑥 𝑑𝑡 = ∑ 𝐶̇𝑜𝑢𝑡𝑖𝑛 [ 𝑘𝑔 𝑚3, 𝑠] → 𝑑𝐶𝑥 𝑑𝑡 = ∑ 𝐷𝑥 𝐴 ∆𝐶 [ 𝑘𝑔 𝑚3, 𝑠] → 𝑑𝐶𝑥 𝑑𝑡 = ∑ 𝐷𝑥 2𝐿∆𝐶 [ 𝑘𝑔 𝑚3, 𝑠]

𝐴 is the through-plane area that the distance and width of the pathway between the node pair covers and is commonly replaced by 2𝐿 for 1D+1D diffusion systems such as the one in this work. ∆𝐶 is the difference in concentration between the interacting node pair.

4.3 Thermal advection

The thermal advection is present only in the gas channels and is modeled through the use of empirical data of the gas’s densitiy and specific heat capacity. A velocity profile is imported for each gas channel and the velocity is interpolated for each node in the through plane direction of the gas channel. Down the channel, velocities are constant.

𝑑𝐸𝑥 𝑑𝑡 = ∑ 𝐸̇𝑜𝑢𝑡𝑖𝑛 [𝐽 𝑠] → 𝑑(𝜌𝑥 𝑐𝑝𝑥 𝑇𝑥 𝑉𝑥) 𝑑𝑡 = 𝑐𝑝𝑖𝑛 𝜌𝑖𝑛 𝑉𝑥 𝑢 𝐿 ∗ 𝑇𝑖𝑛− 𝑐𝑝𝑜𝑢𝑡 𝜌𝑜𝑢𝑡 𝑉𝑥 𝑢 𝐿 ∗ 𝑇𝑜𝑢𝑡 [ 𝐽 𝑠] → 𝑑𝑇𝑥 𝑑𝑡 = 𝑐𝑝𝑖𝑛 𝜌𝑖𝑛 𝑢 𝐿 𝜌𝑥 𝑐𝑝𝑥 ∗ 𝑇𝑖𝑛− 𝑐𝑝𝑜𝑢𝑡 𝜌𝑜𝑢𝑡 𝑢 𝐿 𝜌𝑥 𝑐𝑝𝑥 ∗ 𝑇𝑜𝑢𝑡 [ 𝑇 𝑠] 𝑇𝑖𝑛 is the temperature that is flowing in to node 𝑥, which is temperature of the adjacent node

up-the-channel. 𝑇𝑜𝑢𝑡 is the temperature of the flow leaving the node, which is the temperature

(22)

4.4 Mass advection

The mass advection is present only in the gas channels and is modelled through the principles of mass transfer based on the velocity profile and concentrations.

𝑑𝐶𝑥 𝑑𝑡 = ∑ 𝐶̇𝑜𝑢𝑡𝑖𝑛 [ 𝑘𝑔 𝑚3, 𝑠] → 𝑑𝐶𝑥 𝑑𝑡 = 𝑢 𝐿𝐶𝑖𝑛 − 𝑢 𝐿𝐶𝑜𝑢𝑡 [ 𝑘𝑔 𝑚3, 𝑠] → 𝑑𝐶𝑥 𝑑𝑡 = 𝑢 𝑑𝑡 𝐿 (𝐶𝑖𝑛 − 𝐶𝑜𝑢𝑡) [ 𝑘𝑔 𝑚3, 𝑠]

𝐶𝑖𝑛 is the concentration that is flowing into node 𝑥, i.e. the concentration in the adjacent

up-stream node. 𝐶𝑜𝑢𝑡 is the concentration flowing out of node 𝑥, which is the concentration in

(23)

5

PREVIOUS WORK

Weber & Hickner (2008) performed a series of experiments to measure the liquid water content of the MEA and GDL sandwich in a PEFC. The experiments featured a test-PEFC with active area of 2.1 cm x 7.7 cm of which the middle length of 2 cm was imaged with neutron imaging. The cell was oriented with co-flow gas flows, where gas inlets were located at the top of the cell and followed a single channel serpentine arrangement downwards (along gravity). The gas GC-BPP surface had a 1mm land- and 1mm channel width. The anode flow rate was maintained at 480 𝑠𝑡𝑑. 𝑐𝑚3/𝑚𝑖𝑛. Both the anode and cathode flow were pre-humidified to dew

point. The catalyst layers consisted of Nafion 112 membranes with catalyst contents of 1 𝑚𝑔 /𝑐𝑚2. The catalyst used was E-TEK 40% Pt/C. The GDLs used were carbon cloth with

integrated microporous layer (MPL) E-TEK LT-1400W. The resolution of the neutron imaging was 16.4 µ𝑚 𝑥 16.4 µ𝑚 per pixel. The results of their experiments are shown in Figure 6 (left). Additionally, they also set a model up over the test cell to simulate the temperature distribution between the GDL/GC boundaries for various current densities at 60°C cell temperature. The result of the simulation is shown in Figure 6 (right). No further information regarding these experiments and the simulation were disclosed (Weber & Hickner, 2008).

(24)

6

THEORETICAL FRAMEWORK/ LITERATURE STUDY

This section presents all the empirical data that was used in the model, categorized per layer of the FC.

6.1 Membrane morphology

The polymer material commonly used as the membrane of PEFCs, and work as the electrolyte in the MEA, are referred to as perfluoro-sulfonic-acid (PFSA) ionomers. The molecules of these ionomers feature a single sulfonic acid (−𝑆𝑂3𝐻) pendant group on the side chain (Wang, et al., 2012).

Figure 7 Molecular structure of membrane materials Nafion (a) and 3M ionomer (b). (Wang, et al., 2012)

6.1.1 Water uptake

The water uptake of PFSA membranes decreases with increased equivalent weight 𝐸𝑊 = 𝑔𝑝𝑜𝑙𝑦𝑚𝑒𝑟/𝑚𝑜𝑙𝑆𝑂3− (Kusoglu & Weber, 2017). The water content of the membrane, defined as the

number of water molecules per ionic group (𝜆 = 𝑚𝑜𝑙𝐻2𝑂/𝑚𝑜𝑙𝑆𝑂3−) (Springer, et al., 1991). The average concentration of water absorbed in the membrane, denoted 𝑐𝑤 (𝑚𝑜𝑙𝐻2𝑂/𝑚3), is then

described by Equation 4, where 𝑉̅𝑤 and 𝑉̅𝑝 is the molar volume of water and dry polymer

(25)

6.1.2 Density changes

Water uptake of the membrane causes swelling. This swelling of the membrane causes many changes to the morphology of the membrane’s material. One of these changes regards density decrease of the PFSA ionomer, which is shown in Figure 8. The dry density of the polymer ionomer is in the range 𝜌𝑃𝑑𝑟𝑦= 2.05 ± 0.05 𝑔

𝑐𝑚3 for most PFSA membranes, but changes with EW (Kusoglu & Weber, 2017).

Figure 8 Density of Nafion and 3M PFSA ionomers with various EW. Plotted as function of absorbed-water content @25℃. Nafion data from (Takamatsu & Eisenberg, 1979; Morris & Sun, 2005) and 3M data from (Zawodzinski, et al., 2013).

The dotted line in Figure 8 represents an approximation of the Nafion 1100 membrane’s density, derived from the rule of mixture Equation 5. The symbol 𝜙 denotes the weight fraction of water and PFSA of the membrane (Kusoglu & Weber, 2017).

Equation 5 Rule-of-mixture, membrane density

𝜌𝑀𝐸𝑀 = 𝜙𝑃∗ 𝜌𝑃 𝑑𝑟𝑦

+ 𝜙𝑊∗ 𝜌𝑊

𝜙𝑊= 1 − 𝜙𝑃

𝜌𝑀𝐸𝑀, 𝑁𝑎𝑓𝑖𝑜𝑛 1100= 0.5679𝜆2− 31.311𝜆 + 2098.8

6.1.3 Specific heat capacity

The specific heat capacity also changes as the membrane absorbs water. A rule of mixture solution was applied to approximate the effects of water the uptake on the specific heat capacity, based on values taken from Liu & Zhong (2002).

(26)

Table 5 Specific heat capacity of Nafion 1100 at λ = 14. Molar weight (g/mol) EW g polymer / mol Fraction (g / g polymer) Weight fraction

Specific heat capacity (J/kg,K)

Nafion 1100 (dry) - ~1 1 71% 4188 (Liu & Zhong, 2002)

Water 18.02 44 2.44 29% 4200

SO3- 80.06 ~1100 1100 0.06% 3221

Total: - - - 100% 4189

Figure 9 Specific heat capacity of Nafion 1100 as a function of water content. Based on weight fractions and properties in Table 5.

6.1.4 Thermal conductivity

One study on the thermal conductivity in Nafion membranes at constant room temperature show that the thermal conductivity increases linearly with water uptake, from 0.18 𝑊/𝑚, 𝐾 to 0.27 𝑊/𝑚, 𝐾 . This study covered water contents of near 0 𝑚𝑜𝑙𝐻2𝑂/𝑚𝑜𝑙𝑆𝑂3− to 22 𝑚𝑜𝑙𝐻2𝑂/

𝑚𝑜𝑙𝑆𝑂3− (Burheim, et al., 2010). A second study evaluated the thermal conductivity of dry

Nafion at various temperatures, and found the through-plane thermal conductivity to decrease linearly with temperature, from 0.17 𝑊/𝑚, 𝐾 at room temperature to 0.14 𝑊/𝑚, 𝐾 at 65 ℃ (Khandelwal & Mench, 2006). The in-plane thermal conductivity of dry Nafion was measured by Alhazmi et al. (2013) and they concluded values of 0.19 𝑊/𝑚, 𝐾 at 35 ℃ and 0.14 𝑊/𝑚, 𝐾 at 65 ℃, which suggests that the difference in in-plane and through-plane thermal conductivity of Nafion is not significant at higher operation temperatures (Burheim, et al., 2014). The

4187 4188 4189 4190 4191 4192 0 5 10 15 20 25 30 Sp ecif ic h eat cap acity (J /kg, K)

Water content λ (mol H2O/mol SO3)

(27)

Figure 10 Thermal conductivity of dry Nafion as a function of temperature at atmospheric pressure. (Khandelwal & Mench, 2006; Alhazmi, et al., 2013)

Figure 11 Thermal conductivity of Nafion as a function of water content at atmospheric pressure and room temperature. (Burheim, et al., 2010)

Figure 12 Thermal conductivity of water as a function of temperature at atmospheric pressure. (The Engineering ToolBox, u.d.) y = -0.0008x + 0.3936 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 298.15 308.15 318.15 328.15 338.15 348.15 358.15 Th erm al con d u ctiv ity (W/m,K) Temperature (K)

Thermal conductivity of Dry Nafion as a function of temperature

y = 0.0041x + 0.18 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0 5 10 15 20 25 Th erm al con d u ctiv ity (W, m ,K)

Water content (mol H2O / mol SO3)

Thermal conductivity of Nafion as a function of water content

y = -8E-06x2+ 0.0019x + 0.5643 0.575 0.6 0.625 0.65 0.675 0.7 20 30 40 50 60 70 80 90 100 110 Th erm al con d u ctiv ity (W/m,K) Temperature (℃)

(28)

6.1.5 Species diffusivity

The effective diffusivity of a species inside a porous media generally considers the pathway characteristics such as porosity 𝜀 and tortuosity 𝜏 (Gu, et al., 2010).

Equation 6 Coefficient of effective diffusivity in porous material. (Gu, et al., 2010) 𝐷𝐻 2𝑂 𝑒𝑓𝑓 = 𝐷𝐻2𝑂∗ 𝜀 𝜏

However, measurements of the effective diffusivity are favorable for the preciseness of models. Gu. et. al. (2010) performed a series of experiments to determine the through-plane effective diffusivity of water vapor in air of various Nafion membranes, which result is presented by the solid lines in Figure 13 (right). The dotted lines in Figure 13 represents another experiment, conducted by Springer et. al (Gu, et al., 2010).

These relationships, concluded by Gu. et. al. (2010), are well represented by the expressions in Equation 8 and

Equation 7. The last term in Equation 8 (a) was added by Caulk, et al. (2012) as an adjustment to cover water transport at very low water contents (Caulk, et al., 2012). In these equations the water content 𝜆 is in 𝑚𝑜𝑙𝐻2𝑂/𝑚𝑜𝑙𝑆𝑂3𝐻, temperature 𝑇 in °C, effective water diffusivity 𝐷 in 𝑚2/𝑠.

Figure 13 Water diffusivity as functions of average water content (right) and water content as function of relative humidity (left) in a Nafion 1110 membrane. Experiments performed at 80 °C and by Gu. et. al. and Springer et. al. (Gu, et al., 2010). Illustration by Andreas Michalsky.

(29)

Equation 8 Coefficient of effective diffusivity in Nafion 1110 membrane material as a function of water content. (Gu, et al., 2010; Caulk, et al., 2012)

𝐷𝜆<3.5= 10−6∗ 1.3𝜆(−1 + 𝑒0.21𝜆) + 0.3 ∗ 10−6 (𝑎) 𝐷𝜆≥ 3.5= 10−6∗ 0.05𝜆(1 + 247𝑒−0.66𝜆) (𝑏) 𝐷𝑀𝐸𝑀,𝑒𝑓𝑓 = 𝜌𝑚𝑒𝑚,𝑑𝑟𝑦 𝑀𝐻2𝑂 𝐸𝑊 𝐷𝜆 𝜆 𝑐𝑤,𝑠𝑎𝑡

6.2 Gas diffusion layer morphology

The GDL material is not hydrophilic, so there is no water uptake in the GDL. Some important changes to the morphology of the GDL considered for this work is however regarding the density, specific heat capacity, heat conductivity coefficient and mass diffusion coefficient. The GDL can have many different designs, with stacked carbon fibers being one of the most common ones. This stacking of fiber threads gives a different cross section depending on which angle it is viewed from, which is what fulfills the main purpose of the GDL. The GDL is diffusing reactant gases towards the CL while distributing them more uniformly to achieve a coating of the CL. This is achieved by the through plane transportation resistance that the GDL incurs to the flow, forcing the gases to distribute in the in-plane directions. Figure 1Figure 14 shows imaging of the through-plane and in-plane sections.

Figure 14 Through plane (left) and in-plane (right) cross section imaging of a Toray GDL.

The GDL manufacturer Toray released the following product specs for their TGP-H line of carbon fiber paper GDLs.

(30)

Table 6 Toray carbon fiber paper TGP-H material data. (Toray, u.d.) Product ID TGP-H-030 TGP-H-060 TGP-H-090 TGP-H-120 Thickness 0.11 0.19 0.28 0.37 mm Bulk density 400 440 440 450 kg/m^3 Porosity 80 (0.8) 78 (0.78) 78 (0.78) 78 (0.78) % (-) Surface roughness 8 8 8 8 µm Gas permeability 2500 1900 1700 1500 ml mm/cm^2, hr, mmAq Thermal conductivity (through plane @25C) - 1.7 1.7 1.7 W/m, k Thermal conductivity (in plane @25C) - 21 21 21 W/m, k Thermal conductivity (in plane @100C) - 23 23 23 W/m, k

Coefficient of thermal expansion (in plane @25C~100C)

-0.8 -0.8 -0.8 -0.8 10^-6/K Available dimensions: 15" x 15", 400mm x 400mm, 500mm x 500mm, 800mm x 800mm

6.2.1 Tortuosity and porosity

Gurau, et. al (2000) suggests an effective porosity of 0.4 − 0.7 for ELAT carbon fiber papers of different fineness, and a through-plane tortuosity of 1.5. The effective porosity is considering the presence of liquid water that is partially filling the pores of the GDL structure at thermodynamic and hydrodynamic equilibrium (Gurau, et al., 2000).

Table 7 Porosity and tortuosity of ELAT carbon fiber paper GDLs. (Gurau, et al., 2000)

GDL position Thickness (mm) Effective porosity Through-plane tortuosity

First (GC boundary) 0.1 0.7 1.5

Second 0.15 0.5 1.5

Third (CL boundary) 0.1 0.4 1.5

Garcia-Salaberri et al. (2011) fitted the reported effective in-plane and though-plane diffusivities for Toray’s TGH-H-060 of various PTFE content and derived the following expressions for the in-plane and through-plane tortuosity (Garcia-Salaberri, et al., 2011).

(31)

6.2.2 Thermal conductivity

Gu. et. al. (2010) conducted a series of experiments to determine the through-plane thermal conduction coefficient of a Toray carbon fiber paper. They concluded a thermal conduction coefficient of 0.0176 𝑊/𝑐𝑚, 𝐾 from the obtained slope in Figure 15. This value is valid for a compressive load of 2.2 MPa (Gu, et al., 2010).

Figure 15 Through-plane thermal resistance of Toray fiber-paper (5 wt. % PTFE content and no MPL). (Gu, et al., 2010)

Equation 9 Coefficient of thermal conduction for TGP-H-60 GDL material (through plane). (Gu, et al., 2010) 𝑘𝐺𝐷𝐿(𝑇𝐺𝑃−𝐻, 𝑡ℎ𝑟𝑜𝑢𝑔ℎ−𝑝𝑙𝑎𝑛𝑒) = 1.76 𝑊/𝑚, 𝐾

In a paper by Sadeghi et al., the in-plane thermal conductivity is evaluated through experimental and analytical methods for TGP-H-120 GDLs with varying PTFE content. They state that the effective thermal conduction coefficient of the in-plane direction to be about 12 times higher than the through-plane conductivity and concluded it to be approximately 17.5 𝑊/𝑚, 𝐾 (Sadeghi, et al., 2011).

(32)

Equation 10 Coefficient of thermal conduction for TGP-H-60 GDL material (in-plane). (Gu, et al., 2010) 𝑘𝐺𝐷𝐿(𝑇𝐺𝑃−𝐻, 𝑖𝑛−𝑝𝑙𝑎𝑛𝑒) = 17.5 𝑊/𝑚, 𝐾

6.2.3 Species diffusivity

The effective diffusivity of the GDLs can be approximated by the diffusion coefficient and the general characteristics such as porosity 𝜀 and tortuosity 𝜏 through Equation 6 (Gu, et al., 2010).

Equation 6 Coefficient of effective diffusivity in porous material. (Gu, et al., 2010) 𝐷𝐻 2𝑂 𝑒𝑓𝑓 = 𝐷𝐻2𝑂∗ 𝜀 𝜏

6.3 Gas channel morphology

The gas channels (GCs) of the PEFC are enclosed by the BPP and the GDL. The GCs are located in the ribs of the BPP (see Figure 3), meaning that the gas flow field does not make contact with the GDL throughout all of the GDL surface. At the instances between the ribs/gas channels the BPP is in direct contact with the GDL (Cochet, et al., 2020).

6.4 Bipolar plate morphology

The bipolar plate of the PEFC is commonly made from graphite and is considered to have a solid structure. Therefore, mass diffusion is barely present in the BPP. The important morphological changes in the BPP concerns density, specific heat capacity and the heat conductivity coefficient.

There are many types of bipolar plates which properties vary depending on the manufacturer. Sigracell® is the producer of the product sheet is used in this work.

Table 8 Sigracell® bipolar plates material data. (SGL Carbon, 2020)

(33)

6.5 Catalyst layers morphology

The CLs of the PEM FC is where the ORR and HOR reactions mainly occur. The catalyst layer of PEM FC can be described as a mixture of four ingredients. There is the base of electronically conductive carbon grains (black carbon) making up 200 − 400𝜇𝑚 agglomerates, Pt that is supported on this carbon, and a hydrophobic PTFE+Teflon binder. The last ingredient making up the CL is perfluoro sulfonate ionomer (PFSI), soaked with water (Gurau, et al., 2000). The hydrophobic binder and PFSI ionomer are commonly represented by the Nafion ionomer material (Wang, et al., 2005). The CL has micropores ranging between 20 − 40𝜇𝑚 and macropores ranging between 40 − 200𝜇𝑚. It has been shown that increasing the PTFE or PSFI content does not affect the micropores, suggesting that these components only exist in the larger macropores (Gurau, et al., 2000).

Table 9 Catalyst layer material properties (Gurau, et al., 2000).

Property Value Unit

Thickness 0.01 mm

Effective porosity 0.2 -

Tortuosity 1.5 -

Vol. fraction of ionomer 0.136 -

Ionomer tortuosity in CL 1.5 -

CCL ionomer water content 14 kmol/kmol

ACL ionomer water content 10 kmol/kmol

6.5.1 Density

The density of the CLs depends on the mixture ratio of the different components in the CL. Empirical data of the CL density is not widely available and a rule of mixture solution to estimate the density is necessary. Weight fractions per square centimeter are from the work by Wang, et al. (2005) who determined the optimal Pt loading for the CL, which resulted in a power density of 0.37 𝑊/𝑐𝑚2. Density of the wet Nafion 1100 was taken from Kusoglu & Weber

(2017).

Table 10 Densities, weight fractions and total density of the CL. 20% Pt/C loading. (Wang, et al., 2005) Density (kg/m^3) Specific weight

(mg/cm^2) Weight fraction

Black carbon 2100 2 (Wang, et al., 2005) 0,526

Pt 21 450 0.4 (Wang, et al., 2005) 0,105

Nafion 1100

(λ = 14): ~1750 (Kusoglu & Weber, 2017) 1.4 (Wang, et al., 2005) 0,368 CL density

(34)

Figure 17 Density of a CL with 20% Pt/C loading as a function of water content. Based on weight fractions and properties from Figure 9, Table 5 and Table 10.

6.5.2 Specific heat capacity

The specific heat capacity of the CL is, like the density, dependent on the mixture of materials present in the CL. The rule of mixture is applied to estimate the specific heat capacity at a water content of 14 𝑚𝑜𝑙𝐻2𝑂/𝑚𝑜𝑙𝑆𝑂3−.

Table 11 Specific heat capacity of CL at λ = 14.

Specific heat capacity (J/kg,K)

Specific weight

(mg/cm^2) Weight fraction

Black carbon 710 2 0.526

Pt 130 0.4 0.105

Nafion 1100 (lambda = 14): 4136 (Wang, et al., 2005) 1.4 0.368 Total (porosity = 0.2) 1528.94 3160 3180 3200 3220 3240 3260 3280 3300 3320 0 5 10 15 20 25 De n sity (kg/m ^3)

Water content (mol H2O / mol SO3)

(35)

Figure 18 Specific heat capacity of a CL with 20% Pt/C loading as a function of water content. Based on weight fractions and properties from Figure 8, Table 5 and Table 11.

6.5.3 Thermal conductivity

Burheim et al. (2010) performed a series of experiments on a 20 % Pt/C CL to determine the thermal conductivity of the material at dry- and wet conditions at different compaction pressures. The results of the experiments were compared to a 1D thermo dynamic model, which showed an accuracy of 5 %. They found low water content (𝜆 < 25) CLs to have a thermal conductivity of 0.07 𝑊/𝐾, 𝑚 at 5 𝑏𝑎𝑟 compaction pressure and 0.11 𝑊/𝐾, 𝑚 at 15 𝑏𝑎𝑟. When more water was added to the CLs, it was observed that the water content in the CL only had a significant effect on the thermal conductivity of the CL when the water content exceeded the absorption capacity of the Nafion ionomer and floods the pores of the material. Hence, they concluded that the thermal conductivity of the CL is not significantly affected by the water content of the CL material (Burheim, et al., 2010).

Equation 11 Coefficient of thermal conductivity for a 20% Pt/C CL. (Burheim, et al., 2014) 𝑘𝐶𝐿 = 0.07 𝑊/𝑚, 𝐾

6.5.4 Mass diffusivity

The effective diffusivity of a porous media generally considers the pathway characteristics such as porosity 𝜀 and tortuosity 𝜏. However, for porous media composed of spherical particles such as the CLs, the Bruggeman correlation is commonly used to express the effective diffusivity (Gu, et al., 2010). 1524 1525 1526 1527 1528 1529 1530 1531 1532 0 5 10 15 20 25 30 Cp (J /kg, K)

Water content (mol H2O / mol SO3)

(36)

7

RESULTS

In this section are the inputs used for the simulation presented, followed by a comparison between the simulations performed by Weber & Hickner.

7.1 Inputs used in simulations

The simulations that were run with this works model had the layer inputs according to . The various velocity profiles used as inputs are illustrated in Figure 19. An electrical

Anode side Cathode side

BPP GC GDL CL MEM CL GDL GC BPP Unit Thermal conductivity, tp 350 0.206 (𝐻2) 0.22 0.07 0.27 0.07 0.22 0.0287 (𝑂2) 350 𝑊 /𝑚, 𝐾 Thermal conductivity, ip 350 0.206 (𝐻2) 2.5 0.07 0.27 0.07 2.5 0.0287 (𝑂2) 𝑊 /𝑚, 𝐾 Density 1850 𝜌𝐻2 440 3212 1771 3212 440 𝜌𝑂2 1850 𝑘𝑔 /𝑚3 Specific heat capacity 1200 𝑐𝑝𝐻2 840 1528 4189 1528 840 𝑐𝑝𝑂2 1200 𝐽 /𝑘𝑔, 𝐾 Porosity - - 0.63 0.2 - 0.2 0.63 - - − Tortuosity, tp - - 1.5 1.5 - 1.5 1.5 - - − 0 0.05 0.1 0.15 0.2 0.25 0 0.5 1 Ve locity (m /s )

Position in gas channel Velocity profile 4 Velocity profile 5 Velocity profile 6 Velocity profile 7 0 0.05 0.1 0.15 0.2 0.25 0 0.5 1 Ve loci ty (m /s )

Position in gas channel Velocity profile 1 Velocity profile 2 Velocity profile 3

Figure 19 Velocity profiles. The “position in gas channel” on the x-axis represents the node inside the gas channel. 0 is the first node and 1 is the last node. Note that these values are not at the wall boundary between the gas channel and GDL/BPP, but in the nodes making up the gas flow field only.

(37)

Table 12 Parameters used in FC simulation. Highlighted values are assumptions. *tp = through-plane, ip = in-plane.

7.2 Comparison with simulation by Weber & Hickner

A simulation of the same cell used by Hickner and Mench (2008) was performed to compare the results. The simulations were performed with as similar input parameters possible, but some of the parameters were not disclosed in the work by Weber & Hickner. The GDL type, CCL material and cell- and gas temperatures were known, but layer thicknesses and other properties were unknown. The GDL type was known to be a E-TEK LT-1400W and the catalyst layer was known to be a E-TEK 40% Pt/C coated Nafion 112 membrane. The parameter data for these layers were approximated by extracting the suggested values by Mench, in the book Fuel Cell Engines (2008). These values are shown in Table 13 and Table 14, and included in of section 7.1 Inputs Used in Simulations. The depth of the gas channels was not known and was thereby set to a standard dimension of 2mm for both anode- and cathode side GCs. The same velocity profile, Velocity profile 1 (see Figure 19), was used for all the simulations independent of current density.

Table 13 Typical thermal conductivities of the PEFC layer materials. (Mench, 2008) * Value used.

Material Measured k (W/m,K)

DuPont Nafion membrane (at 30⁰C) 0.16* ± 0.03

W.L. Gore reinforced electrolyte 0.16* ± 0.03

E-Tek ELAT carbon cloth diffusion media (LT1200-W at 33 ⁰C) 0.22* ± 0.04 Catalyst layer (0.5 mg/cm² platinum on carbon) 0.27* ± 0.05 Table 14 Typical thickness of the PEFC layers. (Mench, 2008)

* Value used.

Component Typical thickness

PEFC electrolyte 50-200* µm

PEFC bipolar plate (graphite) 2*-4 mm each PEFC gas diffusion layer 100-300* µm PEFC catalyst 5-30* µm

The simulations were run in two set ups. One set up overed the PEFC cross section where the gas channels were located (Figure 20), and the other one covered the part of the cross section in between the gas channels ribs where the BPP material was in direct contact with the GDL. To approximate the mean temperature- and concentration profiles, the average of these

(38)

Figure 20 This work's model (solid) compared to the model by Weber & Hickner (dotted). Temperature profile for the GC ribs cross section.

Figure 21 This work's model (solid) compared to the model by Weber & Hickner (dotted). Averaged temperature profile over the BPP/GC ribs.

59 60 61 62 63 64 65 66 67 0 0.2 0.4 0.6 0.8 1 Te m p era tu re (C ) Position in GDLs+MEA

Temperature profile in the gas channel section

0.2 A/cm^2 0.5 A/cm^2 0.75 A/cm^2 1 A/cm^2 1.25 A/cm^2 59 60 61 62 63 64 65 66 67 0 0.2 0.4 0.6 0.8 1 Te m p era tu re (C ) Position in GDLs+MEA

Averaged temperature profile

0.2 A/cm^2 0.5 A/cm^2 0.75 A/cm^2 1.25 A/cm^2 1 A/cm^2

(39)

Figure 22 Comparison of the water concentration profile of this work's model, expressed in RH, and the measurements of liquid water concentration by Weber & Hickner. Concentration profile for the GC ribs cross section. 100 % RH at gas inlets. Velocity profile 1.

Figure 23 Comparison of the water concentration profile of this work's model, expressed in RH, and the measurements of liquid water concentration by Weber & Hickner. Averaged concentration profile over the BPP/GC ribs. 100 RH at gas inlets. Velocity profile 1.

0.99 1 1.01 1.02 1.03 1.04 1.05 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0 0.2 0.4 0.6 0.8 1 RH (%) Liq u id w at er (m L/cm ^3) Position in GDLs + MEA

Water concentration profile in the gas channel section

Weber & Hickner This work's model

1 1.01 1.02 1.03 1.04 1.05 1.06 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0 0.2 0.4 0.6 0.8 1 RH (% ) Liq u id w at er (m L/cm ^3) Position in GDLs + MEA

Averaged water concentration profile

(40)

7.3 Temperature distribution

In this section are some illustrations over the steady state temperature distribution for a cell with no-heat-flux boundaries at various gas flow velocity profiles. The only outgoing flux is the gas channel outlets, and the ingoing fluxes are located in the CCL and gas channel inlets.

Figure 24 Temperature distribution at various velocity profiles. 1 𝐴/𝑐𝑚2, 0.6 𝑉.

Velocity profile 1: {max. ; min.} = {0.2;0.05} m/s (Top left) Velocity profile 4: {max. ; min.} = {0.1;0.05} m/s (Top right) Velocity profile 6 {max. ; min.} = {0.06;0.01} m/s (Bottom left) Velocity profile 7 {max. ; min.} = {0.01;0} (Bottom right)

(41)

7.4 Water concentration distribution

In this section is a visualization of the steady state water vapor concentrations for a cell with no-flux boundaries. The only outgoing water vapor flux is the gas channel outlets and the only ingoing water vapor fluxes are the CCL and the gas channel inlets.

Figure 25 Temperature distribution at various velocity profiles. 1 𝐴/𝑐𝑚2, 0.6 𝑉.

Velocity profile 1: {max. ; min.} = {0.2;0.05} m/s (Top left) Velocity profile 4: {max. ; min.} = {0.1;0.05} m/s (Top right) Velocity profile 6 {max. ; min.} = {0.06;0.01} m/s (Bottom left) Velocity profile 7 {max. ; min.} = {0.01;0} (Bottom right)

(42)

8

DISCUSSION

In this section, the comparison between the simulation results and the temperature- and water concentration profiles published by Weber & Hickner are discussed.

8.1 Temperature profiles

When compared to the simulations run by Weber & Hickner, this works’ model displays the same behavior with the temperature peak located in the cathode catalyst layer (CCL) and increased peak temperature with increased current density. But there are differences that raise questions about this works’ models accuracy. These differences are assessed below.

a). The cathode side features less decline in temperature through the GDL in this works model than in the model by Weber & Hickner. This is caused by the low thermal conductivity of the air in the cathode GC used in this work. But the cathode gas channel in Weber & Hickner’s work was also supplied pure oxygen, which raises some awareness whether the convection is accurately represented in this works’ model. The gas channel depth and BPP rib materials are also not known from their work, which adds to the uncertainty when comparing the two results. Evaporation of liquid water present in the cathode CL/GDL area is also a process that might have an effect in reducing the cathode side GDL temperature which is not included in this work.

b). The temperature peaks do not match well for the temperature profiles. One possible cause to this is that Weber and Hickner did not get the averaged temperature profile over the rib/GC sections in the same way as this work, or performed their simulation on only the GC section. This raises questions about how the BPP/GC rib is best accounted for in this type of 2D model. Other possible causes regard the layers included in the simulation of the cell. Firstly, the simulation by Weber & Hickner featured a MPL and this works model did not. Secondly, contact resistances were not included in this works model and the magnitude of effect these have are unclear. Thirdly, the material properties used in the two models might vary significantly.

The material properties such as thermal conductivity, specific heat capacity and density affect the temperature gradient substantially. Especially the thermal conductivity of the BPP has a potentially large impact on the temperature profile and peak temperature as the BPP-ribs are in direct contact with the GDL, allowing for a relatively heat flux through the cell compared to the gas filled ribs where the GDL is in contact with gas. The cell materials containing Nafion

(43)

8.2 Water concentration profiles

When comparing the water concentration profiles, it must be noted that it the measurements represent liquid water concentrations while the simulations by this work’s model represents relative humidity. To continue this comparison, it is assumed that the relative humidity has a strong correlation to the concentration of liquid water. It is also assumed that the abundance of liquid water due to the simulated relative humidity profile is in the same range as in the measurements.

The water simulated concentration profile aligns well with the measured data, but the concentration throughout the cathode is substantially lower in the simulation than for the measured data. The cause for this is believed to be the low diffusivity of the CCL, which causes a sharp drop in water concentration. The other differences in the slope characteristics of the concentration profile would be explained by the different material structure that is present in a real set up. Compared to the simulated set up, which features many constant values for diffusivity and hence only applies to completely uniform material structures. Especially the contact surfaces between the layers make for a disturbance in the layer

structures, and have not been accounted for in the model. There are also other mass transfer processes that affect the water distribution, for example electro-osmotic drag where water is dragged through the membrane as protons pass. This would further increase the water content at the cathode catalyst layer and likely give a better match of the water concentration profiles in the cathode side region.

(44)

9

CONCLUSIONS

The temperature profile that the model produces did not align with the temperature profile of the model by Weber & Hickner. However, there are too many uncertainties regarding the material properties used in the model by Weber & Hickner to make a valid comparison and say whether this works model is correct or incorrect. The trends and behaviors that the models display are however similar, which suggests that the characteristics given by the model are realistic. Further comparisons with simulations and measurements where the material properties are known is required to validate this model.

The water concentration profile that the model produced had similar characteristics as the measurements by Weber & Hickner. This comparison was however based on the assumption that the relative humidity profile displayed by the simulation would yield a similar liquid water concentration profile. Even though this assumption, it is not possible to say if the water concentration would be in the same range span as in the measurements. There are also characteristic differences between the measured- and simulated concentration profile, with likely causes being the structure of the contact surface between layers and the presence of other mass transfer processes such as electro-osmotic drag.

Overall, the comparison between the result of Weber & Hickner and this work’s model show that the model displays realistic characteristics and behavior, but further validation is required to verify its accuracy. If validated, the model serves as a good base for further expansions to simulate the EC PEFC concept by PSI.

(45)

10 SUGGESTIONS FOR FURTHER WORK

The model in this work is not complete in terms of mass transfer processes. Only water vapor mass transfer is modeled, through Fickian diffusion, and in the gas channels through both advection and diffusion. This diffusion considers air to be the diffusion media, which is not correct in actual fuel cell applications. The diffusion should instead consider the actual mixture of gases present at the local points of the FC to be the diffusion media. Thereby, the diffusion of other gases such as hydrogen, oxygen and nitrogen need to be included.

Although the effective thermal conductivity includes some contact resistance, the actual thermal contact resistances between each layer might also be necessary to implement.

Other significant processes to be added to the model are electro-osmotic drag and evaporation in hydrophilic lines, to simulate PSI’s EC PEFC concept.

The relative humidity is not coupled with water content of the PTFE material in the CCL and GDL layers. The effect of water content in these materials might not affect the steady state significantly but coupling the relative humidity to the water content in these layers would allow for a much more accurate representation of these materials’ properties. This is especially important to obtain the correct the density of Nafion.

Mass transfer processes such as electro-osmotic drag and diffusion due to pressure differences needs to be included for a much more accurate representation of the species concentrations. The current density should not be uniformly distributed in the CCL, as in this work’s model, but follow a realistic current density profile.

(46)

REFERENCES

Alhazmi, N. et al., 2013. The in-plane thermal conductivity and the contact resistance of the

components of the membrane electrode assembly in proton exchange membrane fuel cells,

s.l.: J Power Sources, 241 (0), pp. 136-145.

Burheim, O. S. et al., 2014. Study of thermal conductivity of PEM fuel cell catalyst layers , s.l.: International Journal of Hydrogen Energy.

Burheim, O., Vie, P., Pharoah, J. & Kjelstrup, S., 2010. Ex-situ measurements of

through-plane thermal conductivities in a polymer electrolyte fuel cell, s.l.: J Power Sources, 195, pp.

249-256.

Caulk, D. A., Brenner, A. M. & Clapham, S. M., 2012. A Steady Permeation Method for

Measuring Water Transport Properties of Fuel Cell Membranes, s.l.: Journal of the

Elechtrochemical Society.

Cochet, M. et al., 2020. Enabling High Power Density Fuel Cells by Evaporative Cooling

with Advanced Porous Media, Villigen: Journal of The Electrochemical Society.

Garcia-Salaberri, P. A., Vera, M. & Zaera, R., 2011. Nonlinear orthotropic model of the

inhomogeneous assembly compression of PEM fuel cell gas diffusion layers, s.l.:

International Journal of Hydrogen energy.

Gurau, V., Barbir, F. & Liu, H., 2000. An Analytical Solution of a Half-Cell Model for PEM

Fuel Cells, s.l.: Journal of The Electrochemical Society.

Gu, W., Baker, D., Liu, Y. & Gasteiger, H., 2010. Proton exchange membrane fuel cell

(PEMFC) down-the-channel performance model, s.l.: John Wiley & Sons.

International Energy Agency, 2018. Data and statistics. [Online] Available at: https://www.iea.org/data-and-statistics

Khandelwal, M. & Mench, M., 2006. Direct measurement of through-plane thermal

conductivity and contact resistance in fuel cell materials, s.l.: J Power Sources, 161, pp.

1106-1115.

Kusoglu, A. & Weber, A. Z., 2017. New Insights into Perfluorinated Sulfonic-Acid Ionomers, Berkeley, California: American Chemical Society.

Liu, D. & Zhong, C., 2002. Modeling of the Heat Capacity of Polymers with the Variable

Figure

Table 1 Generic PEFC layers and their function.
Figure 1 Hydrophilic lines and water supply channel in a test cell not in operation. Pure nitrogen at cathode side  and hydrogen at anode side gas channels
Figure  3  EC  PEFC  concept  sketch  (left)  and  generic  PEFC  sketch  (right)  with  gas-  and  coolant  flows
Figure 4 Computation time for equation systems of different size (number of nodes = number of equations) for  different solvers
+7

References

Related documents

ܳ ஼௢௡௩ was calculated according to eq.5 using HFS to measure the conduction heat flow ሺܳ ஼௢௡ௗ ሻ and IR-camera to measure the reflection temperature (T Reflection )

Heat generation occurs due to the electrochemical reactions at the active surfaces in the interface between the electrolyte and electrodes [55], and due to the internal

[r]

This work is limited to the study of supply energy requirements and conditions of a circular Eo5 heavy fuel oil tank described earlier in the introduction, as well as performing

Figure 3.2 shows the input and output data sets used in the system identification are 98 samples of time plot which are found from the simulation of the

Finally, in Paper E, the critical parameters for heat transfer and chemical reactions with thermosets and their influence on the peak exotherm temperature and cure time are

The methodology of this project consists mainly of laboratory work and experiments, conducted at the KTH Energy department. The study focuses on the heat

12 The adhesion force map (Fig. 4f) shows a distinct contrast between the polymer matrix and the regions where PCBM-rich domains are present, but the contrast between uncoated