Chapter 1 Introduction
1.3 Scope of Present Study
According to the results of reference [16], low current density is observed at the bending area of flow channels, where low water content and high temperature are presented. The non-uniform distribution of water and temperature will seriously cause mechanical damage to the membrane and shorten its lifetime. Therefore, the flow pattern design in this research will focus on the improvement of cell performance and the increase of distribution
uniformity of current density, temperature and water content.
This research applies a commercial package software, CFD-ACE+, for the full-cell simulation. The investigation mainly concerns on the improvement of flow-field design by changing the bend angles and widths in flow channels.
Simulations are performed in a steady, two-phase and three-dimensional model with different flow channel patterns at different operating temperatures.
This numerical investigation consists of two parts. The first one is to study the effects of bend angle on PEMFC performance with three different flow patterns at two different operating temperatures, 323K and 353K.
Comparisons will be given afterward to discuss how the performance and variables distributions change with bend angle. Another part of the investigation is to study the influences of bend width under the same operation conditions as the first part. The performance and variable distributions will be compared and analyzed as well. Finally, the advantages of those configurations will be discussed.
Fig. 1.1 Cross section of a PEMFC illustrating steps of electrochemical reaction and major structure.
O2
O2
Anode Flow Channel H2
Membrane
Cathode Flow Channel H2→2H++2e
-1/2O2+2e-+2H+→H2O GDL
O2
H2
H2
GDL Catalyst
layer Catalyst
layer
Load 2e
-2e -2H+
+H2O +H2O
+H2O +H2O
Fig. 1.2 Schematic of fuel cell I-V curve with three losses that influence the curve.
Current density (A/cm2) Activation
region
Ohmic region
Concentration region
Cell voltage (V)
Theoretical reversible voltage
Fig. 1.3 Pictorial illustration of the effect of a leakage current loss on overall fuel cell performance.
No parasitic loss
FC with parasitic loss Theoretical reversible voltage
Current density (A/cm2)
Cell voltage (V)
Chapter 2
Simulation Model
The convenience of computational fluid dynamic (CFD) provides a low-cost way to study PEMFC and enables researchers to precisely understand the flow and temperature fields inside PEMFC. Besides, the distribution of chemical parameters can be studied and the effect of each parameter can also be further compared. Based on the governing equations that describe the correlations of various properties can be developed and solved by numerical methods, providing the necessary information for PEMFC design.
2.1 Model Descriptions
Figure 2.1 shows the three different serpentine flow patterns under investigation of bend angle (a) 45–135, (b) 60–120 and (c) 90–90, respectively.
The first number in each pattern represents the first angle of the channels, and the second number represents the second angle. These two are combined to generate a complete u-turn flow channel. In addition, the patterns for investigation on the effect of bend width are shown in figure 2.2. Patterns (a) 1007, (b) 90–90 and (c) 1020, are considered. The first two numbers, 10’s, represent the channel width of 1 mm, and the last two numbers, 07 and 20, represent the bend width of 0.7 mm and 2.0 mm, respectively. 90–90 is the most conventional flow pattern applied in most researches and is considered as a reference for comparison in this study. According to the results from Scholta et al. [14], the channel and rib dimension ratio of this model is 0.5 which would lead to superior performance. Based on the results from Chiang and Chu [15],
the AR is chosen to be 1 which provides the most uniform membrane conductivity. Table 2.1 shows the detailed geometries of current model. The electrochemical parameters and transport properties from Springer et al. [5] and Mazumder and Cole [11, 12] are listed in Table 2.2.
Table 2.1 Geometries of Model
Geometry Vlaues
Reaction area 25 (cm2), 5×5
Channel Width 1 (mm)
Channel Depth 1 (mm)
Rib Width 2 (mm)
GDL Thickness 0.3 (mm)
Catalyst Layer Thickness 10 (μm) Membrane (Nafion117) Thickness 163 (μm)
Table 2.2 Electrochemical parameters and transport properties
Parameter Values Electrochemical parameters
Transfer coefficients at anode 0.5 Transfer coefficients at cathode 1.5
Concentration dependence of H2 0.5 Concentration dependence of O2 1 GDL and catalyst layer
GDL and catalyst layer porosity 0.4
GDL and catalyst layer permeability 1.76×10-11
Tortuosity 1.5 Electrode conductivity of catalyst
layer 53 (1/ohm m)
Effective surface to volume ration of
catalyst layers 1000
Membrane
Membrane porosity 0.28
Membrane permeability 1.8×10-18
Tortuosity 5 Electrode conductivity of membrane 1×10-20 (1/ohm m)
Bipolar Plates
Electrical conductivity Temperature dependence Reference Temperature 293 (K)
Electrical conductivity at reference
Temperature 3.5×10-18 (ohm m)
Temperature coefficient at reference
Temperature 5×10-4
2.2 Model Assumptions
The following assumptions are utilized to simplify simulation conditions.
1. All fluids are treated as continuums.
2. Gases are ideal gases throughout the entire model.
3. The velocity and temperature of reactant species are distributed uniformly over the channel inlet.
4. The flows are incompressible and laminar due to small Reynolds number.
5. The gas diffusion layer (GDL), catalyst layers and membrane are all isotropic porous media with pores uniformly distributed.
6. Gravity is neglected.
7. No fuel crossover.
8. Evaporation and condensation of water is neglected.
9. Dissolution of gases in the liquid is neglected.
10. The phenomenon of double charge layer is not considered.
2.3 Governing Equations
The steady, two-phase and three-dimensional model is basically modified from the one proposed by Mazumder and Cole [11, 12]. The flow pattern has been developed from a single channel to a full-size serpentine channel with 25cm2 reaction area. The membrane domain is built for Nafion 117, and the data is available from Springer et al. [5].
The PEMFC model consists of five volume-averaged conservation equations, which are mass, momentum, energy, species and electric current conservation equations. All the equations are applied in every domain of the model. In addition, fluxes across any two different domains are automatically continuous.
2.3.1 Continuity and Momentum Equations
The continuity equation and momentum conservation equations are as follows
( ) (
U)
0t εeffρ +∇⋅ εeffρ =
∂
∂ (2-1)
( ) ( ) ( )
where ρ is the fluid density, U is the fluid velocity, p is the pressure, τ is the shear stress tensor, B is the body force, μ is the dynamic viscosity of the fluid, εeff is the effective porosity of the medium, and represents the ratio of the volume occupied by the pores to the total volume of porous medium while the permeability κ, is the square of the effective volume to surface area ratio of porous medium. The last term in Eq. (2-2) represents Darcy’s drag force imposed by pore walls on the fluid, and leads to a significant pressure drop across the porous medium. In a purely fluid region, εeff→1 and κ→∞, and the standard Navier-Stokes equation is recovered.
For high accuracy, mix kinetic theory is chosen to estimate the viscosity value of mixture in this research, and is expressed as
ij viscosity of species i and Φij is the dimensionless quantity. The dynamic viscosity data input form [29] is linearly interpolated for each species. And the dimensionless quantity is given by
2
where Mi is the molecular weight of species i, T is the temperature in Kelvin, σi is the characteristic diameter of the molecular in Angstroms and Ωμ is the collision integral and is given by
)
where T* is the dimensionless temperature and is defined as
= Ε
∗ kT
T (2-6)
where k is the Boltzmann’s constant and E is the characteristic energy. The details are given by Incropera [30], Wilke [31] and Hirschfelder et al. [32].
2.3.2 Energy Equation
The temperature field of each domain can be obtained from solving the energy conservation equation. The energy conservation equation can be written as
where h and hs are the gas mixture enthalpy and solid phase enthalpy, respectively. ρs is the solid phase density, jT is the transfer current density which will be discussed further in the section below, (S/V)eff is a direct representation of catalyst loading, η is the overpotential between solid phase and electrolyte phase of the electrode, and S‧
h is the enthalpy source due to phase change. The heat flux, q, is comprised of contributions due to thermal conduction and species diffusion without considering radiation effect, and is written as
i enthalpies defined as sum of the enthalpy of formation and sensitive enthalpy,
and Ji is the diffusion flux of species i. The thermal conductivity, keff, of the porous medium is an effective thermal conductivity of the pores and solid combined region and is expressed as
S
where kS and kF are the thermal conductivities of the solid and fluid regions, respectively. The second term of on the right hand side (RHS) of Eq. (2-7) is the effect of viscous dissipation, the third term on the RHS is the work done by pressure, the last three terms in Eq. (2-7) represent electrical work, Joule heating and energy interactions due to phase change. The irreversible losses due to reaction (conversion of chemical energy to heat) manifest through the last term on the right hand side of Eq. (2-8) because the definition of enthalpy includes both the enthalpy of formation as well as sensitive enthalpy. Mix kinetic theory is also chosen for calculating the thermal conductivity across the whole domain. The data for evaluating the thermal conductivity is available from [33]. Mix Jannaf method is applied for calculating the specific heat and estimated as form:
4
where R is the value that the universal constant divided by each species’
molecular weight, and a0, a1, a2, a3 and a4 are constant vary with species kind.
2.3.3 Species Equation
The species conservation equation is estimated as form:
(
eff i)
i iwhere Yi are the mass fraction of i-th species, Ji are the mass diffusion flux of species i and ω‧
i are the mass production/destruction rates of i-th species in the gas phase. The species mass diffusion rate is expressed as
j
the first term represents Fickian diffusion due to concentration gradients, and the other three terms are correction terms necessary to satisfy the Stefan-Maxwell equations for multicomponent species study. In this case, Di are the effective diffusion coefficient of species i within porous medium, and depend on the porosity and tortuosity, τ, of the medium
Deff=εeffτD (2-13) where the value of tortuosity is usually 1.5 in Eq. (2-13), resulting in the so-called Bruggeman model [11].
The volumetric production/destruction rate of a species i, ω‧
i, due to heterogeneous surface reaction is treated by performing a diffusion flux at the interface of the GDL and catalysts, and is obtained from the surface flux in discrete form expressed as
eff
where Yi denotes the mass fraction at the interface of GDL and catalyst layer, while YPi denotes the species mass fraction in the catalyst layer and δ is the diffusion length scale. Usually, the diffusion length scale cannot be obtained from grid in the case of porous media, because the region under consideration is already much smaller than the size of a control volume.
Typically, the diffusion length scale varies locally from primary pores to secondary pores. With the consideration of continuum mechanism, Knudsen
effect is neglected in tiny pore region, and the diffusion length scale is assumed to be equal to the average pore size in this case.
However, in the case of electrochemical reactions, the volumetric production/destruction rate of a given species is expressed as
F reactants, respectively. F is the Faraday constant. Therefore, reaction-diffusion balance is necessarily required to be written as
F
Eq. (2-16) clearly states that the reactants/products flowing through the interface where chemical reaction takes place by diffusion is equal to the volumetric production/destruction rate of a specific species from chemical reaction.
Maxwell-Stefan Model is applied for multi-species diffusion when over and above three species are involved in a diffusion process. At low density, multi-component gas diffusion can be approximated by the Maxwell-Stefan equation:
2.3.4 Electrochemical Reaction
At the electrode catalyst layers, current is supplied from electrochemical reactions that reactant species participate at anode and cathode. Hydrogen is oxidized at the anode and oxygen is reduced at the cathode. These two reactions are driven by the potential difference between the solid phase and electrolyte phase, called the activation overpotential. The Bulter-Volmer
equation describes this phenomenon and is expressed as
where αa,j and αc,j are kinetic constants determined from experimentally generated Tafel plots, αk,j are the concentration exponents of the k-th species for the j-th step. [Λk] represents the average interfacial molar concentration of the k-th species, and the subscript ref refers to concentration values at the reference state at which the reference current density is prescribed. The electrode potential, η, is the potential difference between the solid phase (ΦS) and the ionic phase (ΦF) of the electrode.
The term j0,j in Eq. (2-18) is the reference current density with units of A/m3. As a matter of fact, the reference current density at the cathode is several orders of magnitude smaller than the value at the anode. Both the reference current densities at anode and cathode are expressed as [34]
Anode: j0=1×109 (2-19)
Cathode: j0=3×105 exp[0.014189(T-353)] (2-20) 2.3.5 Current Equation
The continuity equation of current in any material leads to:
0 i =
⋅
∇ (2-21)
where i is the current density vector. When the material is a porous electrode, the current density is split into two parts: one flowing though the polymer electrolyte (ionic phase) denoted as iF and the other flowing through the solid parts (electronic phase) of the porous matrix denoted as iS. Fig. 2.3 simply shows the transfer of current within porous media. Eq. (2-18) can be written as
0 i iF +∇⋅ S =
⋅
∇ (2-22)
During the electrochemical reactions, electrons are either transferred from ionic phase to solid phase or vice versa shown in figure 2.3. Thus, Eq. (2-19) can be rewritten as
It is generally assumed that each phase has its own electrical potential. In the ionic phase, the electrical current is denoted as electrolytic term while the current is denoted as electronic term in solid phase. Chemical reaction can only occur at the interface of electrolyte and solid when there’s interaction between these two current components. Under this assumption, application of ohm’s law to Eq. (2-23) yields
where ΦS and ΦF are the electrical potentials of the ionic and solid phases, respectively.
2.3.6 Liquid Water Transport
The foundation of this two-phase model is an additional governing equation for the transport and formation of liquid water, given by
l
where s is the saturation of liquid water, and is defined as the ratio of liquid volume to the volume of pores. ρl is the density of liquid water and ζ is the relative mobility of liquid water. εd is the dry porosity which represents the porosity of the medium in the dry state. In this model where liquid water effect
is considered, pores in porous media are occupied by liquid, and the value of effective porosity decreases as the amount of liquid water in the medium increases. Therefore, εeff has the form
)
The first term in Eq. (2-25) represents the storage of liquid water in a certain domain. And the second term represents transport of liquid water by pressure-driven advection.
Membrane water content, λ, is given as the ratio of the number of water molecules to the number of charge (SO−3H+) sites. It shows how intense the electrochemical reaction is, and dictates the electrical conductivity of membrane.
The value of membrane water content is determined by the activity of water vapor at the catalyst/membrane interface. The activity in the water phase, a, is xwP/Psat. xw is the mole fraction of water vapor and the saturation pressure of water is denoted by Psat in units of atm, which can be computed from the curve-fitted expressions provided by Springer et al. [5]
3
Membrane water content is function of water activity for the Nafion 117 membrane by weighing membranes equilibrated above aqueous solutions of various lithium chloride concentrations, and the experimental relationship of water content versus water activity, a, used in this model is
⎪⎩
Eq. (2-28) is employed to calculate the value of λ based on solution of water vapor outside the membrane domain.
When water profile in the membrane has been determined, the membrane electrical resistance and the potential drop across it can easily be obtained. The electrical conductivity of Nafion 117 was proposed by Springer et al. [5], and is given by
Once the value of water content drops less than one, the membrane conductivity is assumed to be constant.
In a PEMFC, electro-osmotic water drag moves liquid water from anode to cathode. As this water built up at cathode, back diffusion occurs, resulting in the transport of water from cathode back to anode. This phenomenon is strongly dominated by capillary diffusivity. Capillary diffusivity depends on factors like the contact angle and saturation, and can be written in terms of the capillary pressure head as
ds
where pc is the capillary pressure head, κl is the permeability of the liquid and μl is the dynamic viscosity of the liquid. Based on the correlation that Springer et al. [5] proposed, capillary diffusivity is function of water content and temperature, and is expressed as [5]
⎪⎩
T ]
In Eq. (2-32), s is considered as a constant with value of 0.0126.
2.3.7 Concentration Loss
During the electrochemical reaction, performance drops when reactant species concentrations are deficient on the reaction surfaces, especially when operating a PEMFC at low operating voltage due to high formation of water that causes clogging. The total concentration loss is expressed as
R reactant molar concentration in catalyst layers and α is the mass transfer coefficient expressing how variations in electrical potential across reaction interfaces changes the reaction rate. The value of α depends on the reaction and electrode material.
One of concentration loss is due to reactant species depletion in the catalyst layers called Nerst potential changes, and it has the form
R
The second way that concentration contributes to concentration loss is due to reaction kinetics, and it has the form
R
Different bend angles and widths of a channel lead to different minor losses;
therefore, fluid velocity through bend decreases. Figure 2.4 illustrates a schematic of the PEMFC mass transport model, indicating that the reactant
concentration on the reaction surface is dominated by fluid velocity in the flow channel. The correlation between gas velocity in the flow channel and gas density on reaction surfaces is expressed as [3]
C the Sherwood number. According to the Eq. (2-37), the geometry of flow patterns is one of the parameters that dictate the performance of PEMFC.
Hence, high performance can be achieved from flow pattern design.
2.4 Boundary Conditions
The governing equations for the current PEMFC model are elliptic, partial differential equations. Hence, boundary conditions are required for all boundaries in the computational domain.
The boundary conditions are presented as follows:
Anodic inlet of the channel:
u=uin, T=Toperating , Y= YH2+ YH2O=1, j=0 (2-38) Anodic outlet of the channel:
P=Patm, j=0 (2-39)
Cathodic inlet of the channel:
u=uin, T=Toperating , Y= YO2+YH2O=1, j=0 (2-40) Cathodic outlet of the channel:
P=Patm, j=0 (2-41)
Outer surface of bipolar plate at Anode:
T=Toperating, V=0 (2-42)
Outer surface of bipolar plate at Cathode:
T=Toperating, V= ηtot=Vre+Vcell (2-43) where ηtot is the total overpotential and Vre is the reversible voltage given by reference [33]:
Vre=0.2329+0.0025T (2-44)
Walls:
T=Toperating, j=0 (2-45)
(a)
(b) (c) Fig. 2.1 Flow field pattern in PEMFC (a) 45-135; (b) 60-120; (c) 90-90,
respectively.
Z X
Y
(a)
(b) (c) Fig. 2.2 Flow field pattern in PEMFC (a) 1007; (b) 1010; (c) 1020, respectively.
Z X
Y
Fig. 2.3 Schematic of electrochemical reaction and the transfer current within a porous catalyst-containing conductor [33].
Fig. 2. 4 Schematic of the PEMFC mass transport model [3].
uin
conv
JC diff
JG
A B
Y GDL
Flow channel
HC
HG
Z
Catalyst layer
Chapter 3
Numerical Methods
3.1 Introduction to CFD-ACE+ software
CFD-ACE+ is a set of computer applications for multi-physics computational analyses. The programs provide integrated geometry and grid generation software, a graphical user interface for preparing the model, a computational solver for performing the simulation, and an interactive visualization software for examining and analyzing the simulation results.
Based on the finite-volume approach [35], governing equations mentioned in chapter 2 are solved by discretization methods which will be described in the following sections.
3.2 Numerical Method for CFD-ACE+
The CFD-ACE+ employs the finite-volume method to discretize the
The CFD-ACE+ employs the finite-volume method to discretize the