CHAPTER 1 INTRODUCTION
1.4 Outlines
This dissertation is organized as follows. Chapter 2 introduces the model development and numerical method for solving the transport phenomenon and electrochemical reaction of the PEMFCS. Basic assumptions employed in the modeling process are also presented. Chapter 3 addresses the effects of existing temperature and humidification level gradients on cell boundaries are explored by the developed model. Polarization curves are presented for various scenarios of these water and thermal management schemes. The mechanisms that cause the variation of the performance curves are discussed. Also, the interrelations among the physical properties and cell performance are elucidated. The influences of transport component
design on the local physical properties as well as cell performance are described in Chapter 4. Two important design parameters, channel aspect ratio and gas diffusion layer thickness are examined in great detail. Transverse distributions of various model variables as well as local reaction rate are plotted and discussed. Chapter 5 introduces the investigation of the electrochemical reaction and performance of PEMFCs with a novel cathode flow channel shape. Through the assignment of various values of shoulder/channel (S/C) ratio at cathode outlet port, the channel configurations can be classified as convergent, straight and divergent shapes. Variations of various model transport variables as well as cell performance are proposed and discussed.
Mechanism that dominates local cell reaction is examined through the comparison of the variation trends among local current density, oxygen concentration and potential fields. In Chapter 6, the conclusions of this investigation are drawn and suggestions for future study are proposed.
e
-e
-e
-e
-e
-Cathode Reaction
Anode Reaction
H
+H
+H
+H
2H
2H
2O
2O
2O
2e
-e
-e
-e
-e
-e
-H
2O H
2O H
2O Cathode Flow Field
Anode Flow Field PEM
e
-e
-e
-Fig. 1.1 Operating principle of a single PEMFC
End plate
Currentcollector
Bipolar plate Diffusion media
Catalyst layer
Membrane Catalyst layer
Diffusion media
Bipolar plate
Current
collector
End plate
Fig. 1.2 Basic components of a single PEMFC
Electrode Kinetics
Dominated Ohmic Polarization Dominated
Cell V o ltage
Mass Transport Dominated
Current Density
Fig. 1.3 The dominant mechanisms of cell performance
Anode PEM Cathode
Total Overpotential
Cell Potential Open Circuit Voltage
Solid Phase Potential
Membrane Phase Potential
Activation Overpotential
Fig. 1.4 Various overpotentials and available potential in fuel cells
CHAPTER 2
MATHEMATICAL MODELING
The computational fluid dynamic (CFD) method plays an important role on the resolution of the flow and temperature fields of the thermal-fluid problems. As the physical domain in a unit PEMFC is extremely thin, usually within the mm scale for the main components, it is extremely difficult to measure the local values of various important physical properties such as species concentration, velocity, temperature and current density. Based on the conservation law, several governing equations describing the inter-correlated relations of various properties can be developed and solved by numerical algorithm, providing essential and fundamental information to the engineer working on the technique development of fuel cell.
In this chapter, the general modeling framework employed in the subsequent study is presented. Equations depicting the electrochemical reactions as well as current transport are incorporated with the conservation equations of the CFD.
Relations featuring the reaction kinetics, membrane conductivity and the transport parameters of the reactant and the product are also provided along with proper boundary conditions to complete the system of the investigation model.
2.1 Model Description
The general structure is periodical and symmetric for a fuel cell with straight flow channel. Therefore, as illustrated in Fig. 2.1, the physical region containing nine major components with one single channel and two half shoulders are considered in modeling. The anode, and cathode catalyst layers (CLs) and membrane are regarded as separate entities, despite the ultra-thinness of the CLs. Two GDLs are compressed by current collector ribs(shoulders) on each cell side to deliver reactants to the reaction sites under the ribs. The spaces in the bipolar plate (BP) grooves, commonly called the anode and cathode channels, are passages in which gases are transported throughout the cell. This work describes the following transport phenomena and reaction in a PEMFC:
z Three-dimensional convection and diffusion in flow channels and porous media;
z Electrochemical reactions in the catalyst layers of the anode and the cathode;
z Multi-species mass transport;
z Formation and transport of liquid water;
z Electronic current transport through bipolar plates and electrodes;
z Positive ionic current transport through the membrane and catalyst layers.
2.2 Basic Assumptions
Basic assumptions are made to simplify actual cell conditions in the theoretical model and thus facilitate the modeling approach of transport component and operating condition influences on transport phenomena and cell performance. The following are the most important:
z The Reynolds number of the fluid is below 100 due to low mixture velocity, laminar flow is considered;
z The GDLs, CLs and membrane are modeled as porous media;
z Each component has isotropic transport properties;
z The gaseous phase of the working fluid behaves as an ideal gas;
z The electric potential on the outer surface of each bipolar plate is constant;
z No charge accumulates in the electrodes and the domain is electrically neutral;
z The system operates in a steady state.
2.3 Governing Equations
A unified domain approach is employed to avoid tedious boundary condition appointments on the component interface in the following formulations. Instead, proper physical and transport properties such as porosity and permeability are designated on each distinct domain. This method can significantly simplify the model system equations and numerical processes. Also the boundary conditions at the
interfaces between different layers are not required (Gurau et al., 1998; Um et al., 2000; 2004; You et al., 2002; Wang et al., 2003).
2.3.1 Continuity Equation
The continuity equation is used to describe the mass conservation of mixture throughout the domain:
(
εeffρ)
=0⋅
∇
U
(2-1)where the mixture density ρ is the volume-weighted average of the phase mass concentration for the consideration of the two-phase flow (Wang et al., 2001). It can be expressed as:
(2-2)
where s stands for the saturation level in porous medium, representing the volume fraction of the pore occupied by the liquid phase and can be given as:
g
2.3.2 Momentum Equation
The generalized Navier-Stokes equation is introduced to represent momentum conservation of the mixture; a source term is included to consider the additional drag forces in the porous medium (Mazumder et al., 2003; Wang et al., 1993):
(2-4)
(
ρεeff)
=−εeff∇P+∇⋅⎜⎝⎛εeffμeff∇ ⎟⎠⎞+εeffρk +Sm⋅
∇
UU U g
where represents the effective porosity given by
;
is the effective viscosity of the mixture;) Darcy and Forchheimer drag forces associated with the morphology (porosity and permeability) of the porous medium:
Sm
The momentum equation automatically becomes the Darcy equation in the porous area where the magnitudes of convection and diffusion term are extremely small. In Eq. (5), the relative permeabilities of the liquid and gas phase can be expressed as a function of saturation (Wang et al., 2001):
(2-8)
2.3.3 Mass Transport
The mass fraction of each species in the gas mixture can be given as follows:
α
Notably, the Bruggemann correction is applied in this equation to consider the mass transport in porous medium. The parameter δ is the tortuosity of the porous
medium. The multi-component diffusion coefficient is a function of species concentration and binary diffusion coefficient, the latter can be written as (Cussler, 1997; Wangard et al., 2001):
Dαβ
The source term defines the creation or consumption of species at the electrode catalyst sites and is given by:
α
2.3.4 Electrochemical Reaction
At the electrode catalyst sites, reactants undergo an electrochemical reaction.
Hydrogen is oxidized and oxygen reduced at the anode and the cathode, respectively.
These two reactions are driven by the potential difference between the solid phase and the electrolyte phase, called the activation overpotential η . The activation overpotential in the anode tends to energize the electrons lost by hydrogen; at the cathode, hydrogen ions and electrons react with oxygen to form water. The
Bulter-Volmer equation describes this important phenomenon and is related to the source terms in Eq. (2-12)-(2-14):
⎥⎦
where is the effective reaction surface area of the catalyst particle and is closely related to the morphology and platinum loading of the catalyst layer. An expression for the relationship among these parameters is written as:
AV
where is the catalyst area per unit mass of catalyst particle, represents the catalyst loading and stands for the catalyst layer thickness.
Amas mplt
Hcat k ,
iref
The term in Eq. (2-15) is the exchange current density, which characterizes the catalyst layers. Due to sluggishness of oxygen reduction reaction at the cathode, its value is several orders of magnitude smaller than that at the anode.
Hence, the cathode exhibits significant activation overpotential during cell operation.
The experimental results of Parthasarathy yielded the following relation (Parthasarathy et al., 1992):
⎟⎠
2.3.5 Potential Fields
Two charged species, electrons and protons, are transported in the fuel cell and subjected to individual driving forces determined by the potential gradient. The negatively charged electrons flow through the catalyst layer, the gas diffusion layer, and the bipolar plate; the positively charged protons flow goes through the catalyst layer and the membrane. The catalyst layer is a source of charges at the anode and a sink at the cathode. The assumption of electro-neutrality yields the following expressions for electrical current conservation in the catalyst layer:
(2-18)
Outside the catalyst layer, no current sink or source is presented and the right hand sides of Eqs. (2-18) and (2-19) are equal to zero, suggesting that no species is created or consumed; electrical current conservation is therefore easily understood.
In the membrane, the ionic conductivity σmem is strongly related to the water contentγ , defined as the ratio of the number of water molecules to the number of the charge sites (Springer et al., 1991):
⎭⎬
An empirical relationship between the water content in the membrane and the partial pressure of the water is:
where a is the water activity and is given by:
In the above equation, the saturation pressure varies with temperature and can be determined from the thermodynamic table or using the following empirical expression: with a relative error less than 1.5 percent in the range of 20 to 100 .
oC
Psat
oC oC
2.3.6 Liquid Water Transport
During fuel cell operation, water partial pressure in the electrode may exceed its saturation pressure if the local water concentration is high. Therefore, liquid water is possibly formed and occupies the electrode pore space. Operating the cell at a high reaction rate may cause severe mass transport overpotential because the diffusion species are blocked. Additionally, extremely small pores in the fuel cell porous media cause the capillary force to dominate liquid water transport in the hydrophilic surfaces.
Smaller pores correspond to stronger capillary forces. However, the actual expression of this force cannot be formulated because the real liquid-gas interface configuration is not tractable. The generalized Richards equation, originally developed for the oil
recovery field and adapted by Wang et al. (1992) to study two-phase flow transport in capillary porous medium, is applied :
)
where represents the liquid water flux due to electro-osmosis drag in the membrane, it is associated with proton transport and is a function of local current density and the electro-osmotic drag coefficient. The first term on the right-hand side describes the capillary pressure effect, whereas the second term is the gravitational separation between the liquid and gas phases of water. The variable D is the capillary diffusion coefficient and can be expressed as a function of saturation level, surface tension, and permeability (Wang and Beckermann, 1997; You and Liu, 2002):
Nl
where , are the relative permeabilities associated with pore space reduction by the co-existence of multiphase fluids.
krl krg
The last term in Eq. (2-24) is introduced to account for liquid water
condensation or evaporation, and can be expressed as (He et al., 2000):
)
where y is a constant with unity or zero value depending on water species condensation or evaporation scenarios in the porous medium.
2.3.7 Energy Equation
The temperature field of the cell domain can be obtained from solving energy equation which describes the thermal energy conservation. Except for inlet and outlet boundaries, enthalpy generation is considered from three sources; joule heating due to current transport in the medium, electrochemical reaction irreversibility due to over-potential at electrodes, and latent heat of water due to phase change. Many physical properties are closely related to temperature in the cell domain. First, mixture velocities at channel inlet ports for a fixed reactant stoichiometry and a given cell geometry are primarily decided by their thermodynamic states subsequently dictated by mixture temperatures. Secondly, electrochemical reaction kinetic in the catalyst layer and membrane conductivity are also influenced by local temperatures. Thirdly, the temperature field is also responsible for water formation or evaporation.
Consequently, it has certain effects on reactant concentration and membrane water content. The generalized steady state energy equation can be expressed as:
h
where the first two terms on the right hand side represent the conduction energy and reactant enthalpy flux, the third term is the irreversible viscous dissipation. The fourth
and the fifth terms describe the electrical related thermal effects where Sj is the volumetric current density in catalyst layer that can be expressed by the Bulter-Volmer equation. The last term results from the consideration of phase change when water evaporation or condensation occurs in the cell domain and is presented in the next section following the phase change formulation. Detail expansion of this equation in scalar form is expressed in App. A.
2.4 Boundary Conditions
2.4.1 Flow Field Boundary Conditions
Model boundaries fall into three categories as indicated by Fig. 2.1 - symmetric boundaries (SBs), impermeable boundaries (IBs), and channel boundaries (CBs). The SB conditions are imposed on both transverse sides of the model and physical quantities such as mass flux or momentum flux have zero gradients. The Neumann conditions are also assigned on IB, the top (or bottom) and side surfaces of the channels adjacent to the bipolar plates. Additionally, non-slip conditions are applied for the IB velocity fields.
Unlike these simple boundary conditions, the CB conditions importantly determine operating conditions and cell performance. At the inlet ports of the channels,
velocities, thermodynamic states, and mass fractions of mixture are specified according to the desired operating scenarios and mixture of interest. Based on the stoichiometric flow ratio, the inlet velocity of the reactant is:
k
The mass fractions of the species on the CB are specified according to operating pressure and fully saturated humidification condition. Dalton’s law and the ideal gas law yield the molar fraction of individual species. The mass fraction, used in
the preceding formulation, is given by:
=
∑
2.4.2 Potential Field Boundary Conditions
Specific solid phase potential is assigned to each bipolar plate outer surface in addition to these transport variables. This value is generally set to zero on the anode BP outer surface, and a cell total overpotential is set on the cathode BP. Meanwhile, a membrane phase potential zero gradient is applied at the interface between GDL and CL, representing a protonic current absence through this interface because of the lack of a conducting medium. This method yields the local current density in the cell from ohmic law; it also accurately yields the catalyst layer activation overpotential. The cell potential can be obtained from the following expression:
tot l oc
ce V
V = +η (2-30)
where is the total cell overpotential and is the open circuit voltage given by (He et al, 2000):
2.5.1 Calculation Procedure
The nonlinearity and coupling in the model formulation eliminate an analytical solution possibility. The forgoing transport equations are numerically solved using the commercial computational fluid dynamic code CFD-ACE+, based on the control volume formulation and the SIMPLE algorithm (Patankar, 1980; CFD-ACE TM User Manual, 2004; Liu et al., 2005; Soong et al., 2005). The numerical flow chart of present investigation is shown in Fig. 2.2. The calculation is regarded to be converged if the normalized residual of each variable is less than 10e-4.
2.5.2 Model Validation
Three structure mesh systems- 60x30x12, 80x40x16, and 90x48x18 are constructed to explore numerical result dependence on computational cell number.
Base model geometries and operating conditions are listed in Table 2.1. Model component electrochemical parameters and transport properties, obtained mainly from open literature (Um et al., 2000; Mazumder and Cole, 20003, Marr and Li, 1999;
Weng et al., 2005; Hum and Li, 2004; Nguyen et al., 2004;
http://www.etek-inc.com/home.php.), are listed in Table 2.2. Numerical calculations are performed on a Pentium IV 2.4GHz PC with a 1G RAM. The mesh system (80x40x16) is adopted because current density values are satisfactory with an error range of 2% when using base model geometries and parameters of these two Tables, please see Table 2.3. Based on the reported operating conditions and component geometries of Wang et al.’s report (Wang, Husar, Zhou and Liu, 2004), data in Figure 2.3 show that the present model to be consistent with previous studies.
Table 2.1 Cell geometries and operating conditions of base model
Geometry and condition Value Unit
Domain length 50 mm
Domain width 1.6 mm
Gas channel width 0.8 mm
Domain height 4.303 mm
Gas channel height 0.8 mm
Diffusion layer height 0.254 mm
Catalyst layer height 0.01 mm
Membrane height 0.175 mm
Molar ratio of cathode side dry air (N /O ) 2 2 79/21
Reactant relative humidity 100%
Reactant stoichiometry 3
Cell back pressure 2 atm
Operating temperature 353 K
Table 2.2 Electrochemical parameters and transport properties
Parameters and Properties Value Unit Sources
(Um et al., 2000) Porosity of the diffusion and catalyst layer 0.4
(Um et al., 2000)
Porosity of the membrane 0.28
(CFD-ACE TM User Manual, 2004) Permeability of the diffusion and catalyst layer 2.3E-11 m2
(CFD-ACE TM User Manual ,2004)
Permeability of the membrane 1.0E-18 m2
(Mazumder and Cole, 2003) Tortuosity of the diffusion and catalyst layer 1.5
(Mazumder and Cole, 2003)
Tortuosity of the membrane 3
(He et al., 2000)
Condensation rate constant 100 sec-1
(He et al., 2000)
-1 sec-1
Evaporation rate constant 100 atm
(Marr and Li, 1999) 0.5
Concentration dependence of H2
(Marr and Li, 1999)
Concentration dependence of O2 1
(Mazumder and Cole, 2003) Transfer coefficients at anode (anodic and cathodic) 0.5
(Mazumder and Cole, 2003) Transfer coefficients at cathode(anodic and cathodic) 1.5
(Hum and Li, 2004) Electrical conductivity of electrode 114 S m-1
Catalyst loading 0.4 mg cm-2
(Nguyen et al., 2004) Exchange current density for anode reaction 1.4E5 A cm-3
(Nguyen et al., 2004) Reference concentration of oxygen 3.39E-6 mol cm-3
(Nguyen et al., 2004) Reference concentration of hydrogen 5.64E-5 mol cm-3
Table 2.3 Results of three computation grid systems based on the parameters in Table 2.1 and Table 2.2
90x48x18
Grid 60x30x12 80x40x16
Cell
Voltage current density (A/cm
cpu
time relative current density (A/cm
cpu
time relative current cpu time
error error density
(A/cm
(V) 2) (100%) 2) (100%) 2)
0.47 0.866 17324 3.12 0.885 43989 0.99 0.894 105114
0.57 0.645 14657 3.84 0.661 37643 1.43 0.670 89178
6113 3.81 31270 1.44
0.67 0.433 0.444 0.450 79827
0.77 0.239 4688 3.51 0.245 24386 1.33 0.248 56722
0.87 0.075 3467 2.67 0.077 20043 1.00 0.077 47062
1 Computational domain 2 Cathode bipolar plate 3 Cathode flow channel 4 Cathode gas diffusion layer 5 Cathode catalyst layer 6 Membrane
7 Anode catalyst layer 8 Anode gas diffusion layer 9 Anode bipolar plate 10 Anode flow channel
Fig. 2.1. Physical and computational domains considered in this study
START
Set up grid system and initial values
Impose BCs
Refresh model variables Calculate transport properties and
correlation functions
Calculate velocity and pressure fields
Calculate other scalar variables
Fig. 2.2. Numerical flow chart of current study Post processing
Converged? no
yes
STOP
0 0.4 0.8 1.2 1.6 Current Density (A/cm2)
0.4 0.6 0.8 1
Cell Potential (V)
Exp. result of Wang et al. at 323K Model result of Wang et al. at 323K Current model result at 323K Exp. result of Wang et al. at 343K Model result of Wang et al. at 343K Current model result at 343K
Figure 2.3 Comparison of current model results with Wang et al. at 323K and 343K
CHAPTER 3
EFFECTS OF TEMPERATURE AND HUMIDIFICATION LEVELS ON THE PERFORMANCE OF PROTON EXCHANGE MEMBRANE FUEL CELLS
3.1 Introduction
The water and thermal management plays an important role on the transport phenomena and electrochemical reaction of the PEMFCs. When fuel cell operates with a proper control of these two factors, the membrane is able to conduct irons effeciently because of the higher water content and electric conductivity. Moreover, elevated cell temperature is beneficial to the enhancement of ionic phase conductivity as well as the prevention of liquid water formation if the counter effect on membrane water content can be ignored. Unfortunately, these three factors are coupled and mutual interlocked during the operation of PEMFCs, increasing the difficulty of reaching an optimal operation of the water and thermal management. Generally, the operating temperature and humidification level are two important factors of implementing the cell water and thermal management. In the following sections, cell polarization curves obtained from scenarios of various temperature gradient conditions as well as humidification schemes are presented. Prediction results
The water and thermal management plays an important role on the transport phenomena and electrochemical reaction of the PEMFCs. When fuel cell operates with a proper control of these two factors, the membrane is able to conduct irons effeciently because of the higher water content and electric conductivity. Moreover, elevated cell temperature is beneficial to the enhancement of ionic phase conductivity as well as the prevention of liquid water formation if the counter effect on membrane water content can be ignored. Unfortunately, these three factors are coupled and mutual interlocked during the operation of PEMFCs, increasing the difficulty of reaching an optimal operation of the water and thermal management. Generally, the operating temperature and humidification level are two important factors of implementing the cell water and thermal management. In the following sections, cell polarization curves obtained from scenarios of various temperature gradient conditions as well as humidification schemes are presented. Prediction results