Numerical investigation of transport component design
effect on a proton exchange membrane fuel cell
Mu-Sheng Chiang
a,b, Hsin-Sen Chu
a,∗aDepartment of Mechanical Engineering, National Chiao Tung University, Hsin-Chu 30010, Taiwan, ROC
bDepartment of Mechanical Engineering, Nan Kai Institute of Technology, 568 Chung Cheng Road, Tsao-Tun , Nantou 54243, Taiwan, ROC Received 12 December 2005; received in revised form 16 January 2006; accepted 25 January 2006
Available online 2 March 2006
Abstract
A numerical investigation of the transport phenomena and performance of a proton exchange membrane fuel cell (PEMFC) with various design parameters of the transport component is presented. A three-dimensional fuel cell model, incorporating conservations of species, momentum, as well as current transport, is used. The Bulter–Volmer equation that describes the electrochemical reaction in the catalyst layer is introduced; the activation overpotential connects the solid phase potential field to that of the electrolyte phase. Through cell performance simulation with various channel aspect ratios and gas diffusion layer (GDL) thicknesses, a slender channel is found suitable for cells operating at moderate reaction rate, and a flat channel produces more current at low cell voltage. Plots of transverse oxygen concentration and phase potential variation indicate that these oppositely affect the local current density pattern. The relative strengths of these two factors depend on the transport component position and geometry, as well as on the cell operating conditions. Consequently, the curves of cell output current density demonstrate that the optimal GDL thickness increases as the cell voltage decreases. However, at the lowest considered cell voltage of 0.14 V, optimal thickness decreases as that of a thick GDL. The oxygen deficiency caused by long traveling length and clogging effect of liquid water reverses this relationship.
© 2006 Elsevier B.V. All rights reserved.
Keywords: Proton exchange membrane fuel cell; Channel aspect ratio; Gas diffusion layer; Overpotential; Cell performance; Reactant concentration
1. Introduction
Fuel cell technology greatly impacts industry development and has attracted recent attention because the shortage of fos-sil oil and environmental protection have become increasingly important issues. The proton exchange membrane fuel cell (PEMFC) particularly has excellent high efficiency and zero emission characteristics and can be operated at low temperature. It is a commonly proposed substitute, therefore, to the internal combustion engine[1].
The PEMFC uses pressurized reactants that circulate through grooved bipolar plate channels to provide cell energy. Oxidation occurs on the anode side, releasing energetic electrons and pro-tons. Meanwhile, the cathode reduces oxygen with electrons and protons to form water. Electric energy is released during elec-trons transport through external loading. Structural simplicity,
∗Corresponding author. Tel.: +886 3 571 2121x55115; fax: +886 49 2566725. E-mail addresses: [email protected], [email protected] (H.-S. Chu).
compact operating principles, and product neatness are the most attractive features.
PEMFC performance is dictated by several factors such as transport component geometry and morphology, and operating conditions. Giorgi et al.[2]investigated diffusion layer porosity influence on effective catalyst activity for the cathode reduction reaction. Thickness and diffusion layer structure effects on low Pt-loading electrode performance for PEMFC were also con-ducted experimentally. An optimal thickness of diffusion layer was found due to its lower electrical resistance. Decreased per-formance with the thickest diffusion layer was attributed to a long reactant transport passage and to the flooding problem. Jordan et al.[3]reported the carbon support and electrode load-ing effects on PEMFC performance. Various parameter impacts, such as diffusion layer thickness, PTFE content and morphology, as well as operating temperature and humidification condition on cell performance were addressed. Lee et al.[4]asserted from experimental results that an optimal bolt torque was obtained for a soft commercial diffusion layer because of porosity and electrical contact resistance changes. Their results also revealed
0378-7753/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jpowsour.2006.01.086
Nomenclature
c molar concentration (mol m−3)
C form drag constant
D diffusivity (m2s−1)
F Faraday constant, 96,500 C mol−1
g gravitational acceleration (m s−2) H height (m) i current density (A m−2) k relative permeability K permeability (m2) L length (m) m catalyst loading (mg m−2)
M molecular weight (kg mol−1)
n number of species
P pressure (Pa)
q number of electron transferred
r rate constant
R universal constant, 8.314 J mol−1K−1
s saturation level
S source term
T temperature (K)
U mixture velocity (m s−1)
v pore volume of porous medium (m3)
V Voltage (V)
w mass fraction of species
W width of domain (m)
x molar fraction
y switching function
Greek symbols
β transfer coefficient
γ water content in membrane
δ tortuosity ε dry porosity η overpotential (V) θ contact angle µ dynamic viscosity (N s m−2) ν kinematic viscosity (m2s−1)
ξ stoichiometry of hydrogen or oxygen
σ electric conductivity (S m−1) ς surface tension (Nm−1) φ potential (V) ϕ concentration dependence Superscripts c capillary eff effective value
Subscripts a anodic c cathodic cat catalyst cel cell cha channel con condensation eva evaporation F Forchheimer term g gaseous phase
j volumetric transfer current density
k anode or cathode l liquid phase m momentum
mas per unit mass of carbon supported catalyst particle mem membrane phase
p index for phases plt platinum s species sat saturation sol solid phase v per unit volume w water
α index for species
that a cell with high bore torque corresponding to a thin diffusion layer generates more current at a moderate reaction rate, and a moderate bore torque performs better at a high reaction rate for this kind of diffusion layer.
Several elaborate mathematical models were developed to provide qualitative insights into the transport phenomena in a PEMFC, besides the experimental works. He et al. [5] pre-sented a two-dimensional, two-phase, multi-component trans-port model investigating gas and liquid water hydrodynamic effects on PEMFC performance employing an inter-digitated flow field. Parametric study on electrode thickness effect indi-cated that within an electrode finite thickness, average cur-rent density increases monotonically with electrode thickness increase, no matter what the operating voltage. Furthermore, beyond this finite thickness, cell performance decreases. This phenomenon was explained by reactant transport and liquid water removal viewpoints. The effects of channel and shoulder width ratio (C/S ratio) revealed that cell performance enhances with a greater C/S ratio. Natarajan and Nguyen [6] proposed a half-cell, transient model for the cathode using the conven-tional flow field to evaluate design parameter influences on fuel cell transport processes. Contrary to He et al.[5], they found, as diffusion layer thickness decreases at lower cathode overpo-tential, the reactant transport is not affected by the liquid water and at some diffusion layer thickness, current density loss at the shoulder region outweighs channel region gain, leading to an optimum diffusion layer thickness. However, a thinner diffusion layer thickness results in better performance and an optimum thickness does not exist at a higher reaction rate.
In the two-dimensional half-cell model performed by West and Fuller[7], effects of rib sizing and geometry of GDL on cur-rent and water distributions within the cell were reported. Special attention was given to the relation between rib width and water content in the membrane. Kulikovsky et al.[8]conducted the parametric studies with high and low values of solid phase elec-trode conductivity elucidating the possibility of locally reducing
catalyst loading in the catalyst layer. Their study concluded that at low solid phase conductivity, cell reaction concentrates mainly at the rib region catalyst layer and the catalyst can be removed from the channel region catalyst layer with little impact on total cell performance. Simplified descriptions of the kinetic and mass transport model developed from Jeng et al.[9]concluded that GDL effectiveness decreases with reaction rate and increases with flow channel width. Also, when GDL porosity is low, cell performance decreases with an increase in GDL thickness. Sun et al. [10,11]conducted a series of investigations on catalyst layer structural parameter influences with an improved cross-the-channel model. Various design and operating factors such as
C/S ratio, orthotropic GDL conductivity, and electrode
compres-sion effects were investigated in the model while charge transfer and oxygen diffusion were simultaneously taken into account. Gas and electron transport plays an essential role in cathode performance in that both electric conductivity and GDL thick-ness could be vital parameters for transport component optimal design, as conclusions indicated.
Most groups researching cell performance simulation have recently begun applying commercial or self-developed compu-tational fluid dynamic (CFD) algorithms in their study. Hence, a more rigorous three-dimensional simulation has been made pos-sible for investigating the complex transport phenomena. The species and properties in the physical domain can be described in greater detail, and numerical simulation results are more reli-able than one- or two-dimensional models[12–17].
The objective of this study is to present the interactions between transport component design and PEMFC performance. A three-dimensional, multi-species and two-phase model is used to investigate channel aspect ratio and GDL thickness effects on cell performance. Local resolution of current generation domi-nant mechanisms is also emphasized, revealing electrochemical crucial reaction factors, such as reactant concentration and acti-vation overpotential.
2. Mathematical formulations
Fuel cell configuration is symmetric with a typical straight flow field, so only half of the domain is considered here. The model considers nine components, as Fig. 1 illustrates. The anode and cathode catalyst layers (CLs) and membrane are con-sidered separate entities, despite the ultra-thinness of the CLs. Two GDLs are compressed by current collector ribs 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 follow-ing transport phenomena:
• three-dimensional convection and diffusion in flow channels and porous media;
• electrochemical reactions in the catalyst layers of the anode and the cathode;
• multi-species mass transport;
• formation and transport of liquid water;
Fig. 1. Physical domain of interest, the channel height is in Z-direction and the channel width is in Y-direction.
• electronic current transport through bipolar plate and elec-trode;
• positive ionic current transport through the membrane and catalyst layer.
2.1. Basic assumption
Basic assumptions are made to simplify actual cell condi-tions in the theoretical model and thus facilitate the modeling approach of transport component influence on transport phe-nomena and cell performance. The following are the most impor-tant:
• the Reynolds number of the fluid is below 100 due to low mixture velocity, laminar flow is considered;
• the GDLs, CLs and membrane are modeled as porous media; • each component has isotropic transport properties;
• the gaseous phase of the working fluid behaves as an ideal gas;
• the electric potential on the outer surface of each bipolar plate is constant;
• no charge accumulates in the electrodes and the domain is electrically neutral;
• the system operates in a steady and isothermal state.
2.2. Governing equations
A unified domain approach is employed to avoid tedious boundary condition appointments on the component interface in
the following formulations. Instead, proper transport properties such as porosity and permeability are designated on each dis-tinct domain. This method can significantly simplify the model system equations and numerical processes.
The continuity equation is used to describe the mass conser-vation of mixture throughout the domain:
∇(εeffρU) = 0
(1) where the mixture density ρ is the volume-weighted average of the phase mass concentration for the consideration of the two-phase flow[18]. It can be expressed as
ρ = ρls + ρg(1− s) (2)
where s stands for the saturation level in porous medium, repre-senting the volume fraction of the pore occupied by the liquid phase and can be given as
s = vl vl+ vg
(3) The generalized Navier–Stokes equation is introduced to rep-resent momentum conservation of the mixture; a source term is included to consider the additional drag forces in the porous medium:
∇(ρεeffUU) = −εeff∇P + ∇(εeffµeff∇U) + εeffρkg + Sm
(4) where εeff represents the effective porosity given by
εeff= ε(1− s), µeff is the effective viscosity of the mixture, ρ k
the kinetic density, Smis the sum of Darcy and Forchheimer drag
forces associated with the morphology (porosity and permeabil-ity) of the porous medium:
µeff= ρls + ρg(1− s) (krl/νl)+ (krg/νg) (5) ρk = ρlkrlν νl + ρ gkrg ν νg (6) Sm= −(εeff) 2U √ K µeff √ K + ε effCFρ|U| (7) 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[18]:
krl= s3
(8)
krg= (1 − s)3 (9)
The mass fraction of each species in the gas mixture can be given as follows: ∇ εeffρw αUkrgνgν = ∇[ρDα(εeff)δ∇wα]+ Ss (10)
Notably, the diffusivity is modified from that of the free stream value by applying the Bruggemann correction. The parameter δ is the tortuosity of the porous medium. The source
term Ss defines the creation or consumption of species at the
electrode catalyst sites and is given by
SO2 = Sj,cMO2 4000F (11) SH2O= −Sj,c MH2O 2000F + Sl (12) SH2 = Sj,aMH2 2000F (13)
At the electrode catalyst sites, reactants undergo an electro-chemical 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 Eqs.(11)–(13):
Sj,k = Aviref,k cα cα,ref ϕα × exp β k,aFηk RT − exp −β k,cFηk RT (14) where Av 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=
Amasmplt
Hcat (15)
where Amasis the catalyst area per unit mass of catalyst
parti-cle, mpltrepresents the catalyst loading and Hcatstands for the
catalyst layer thickness.
The term iref,k in Eq. (14)is the exchange current density,
which characterizes the catalyst layers. Due to sluggishness of oxygen reduction reaction at the cathode, its value is sev-eral orders of magnitude smaller than that at the anode. Hence, the cathode exhibits significant activation overpotential during cell operation. The experimental results of Parthasarathy et al. yielded the following relation[19]:
iref,O2 = 104exp 3.507−4001 T (16) Two charged species, electrons and protons, are transported in the fuel cell and subjected to individual driving forces deter-mined by the potential gradient. The negatively charged elec-trons 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 expres-sions for electrical current conservation in the catalyst layer: ∇(−σsol∇φsol)= −Sj,k (17)
∇(−σmem∇φmem)= Sj,k (18)
Outside the catalyst layer, no current sink or source is pre-sented and the right-hand sides of Eqs.(17) and (18)are equal to zero, suggesting that no species is created or consumed; elec-trical 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[20]:
σmem= (0.005139γ − 0.00326) exp 1268 1 303− 1 T (19) An empirical relationship between the water content in the membrane and the partial pressure of the water is
γ =
0.043+ 17.81a − 39.85a2+ 36a3, for 0 < a ≤ 1 12.6+ 1.4a, for 1 < a≤ 3
(20) where a is the water activity and is given by
a = xwppsat (21)
In the above equation, the saturation pressure varies with temperature and can be determined from the thermodynamic table or using the following empirical expression:
log10psat= −2.1794 + 0.02953T − 9.1837 × 10−5T2
+ 1.4454 × 10−7T3 (22)
where T is in unit of◦C and psatis in unit of bars. This formulation
can be utilized with a relative error less than 1.5% in the range of 20–100◦C.
During fuel cell operation, water partial pressure in the elec-trode may exceed its saturation pressure if the local water con-centration 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 overpo-tential because the diffusion species are blocked. Additionally, extremely small pores in the fuel cell porous media cause the cap-illary 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 formu-lated because the real liquid–gas interface configuration is not tractable. The generalized Richards equation, originally devel-oped for the oil recovery field and adapted by Wang and Becker-mann[21]to study two-phase flow transport in capillary porous medium, is applied: ∇ εeffρUk rlνlν + ∇Nl = ∇ εeffDc∇s − Kkrlkrg(ρl− ρg)g krlνg+ krgνl + Sl (23)
where Nlrepresents the liquid water flux due to electro-osmosis
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 gravi-tational separation between the liquid and gas phases of water. The variable Dcis the capillary diffusion coefficient and can be expressed as a function of saturation level, surface tension, and permeability[22,23]: Dc=− Kkrlkrg ς cos θc εeff K 1/2 (−3.789s2+3.338s−0.966) krlνg+ krgνl (24) where krl, krgare the relative permeabilities associated with pore
space reduction by the co-existence of multiphase fluids. The last term in Eq.(23)is introduced to account for liquid water condensation or evaporation, and can be expressed as[5]:
Sl= Mlrcon εeffxw
RT (xwp − psat)y
+ revaεeffsρl(xwp − psat)(1− y) (25)
where y is a constant with unity or zero value depending on water species condensation or evaporation scenarios in the porous medium.
2.3. Boundary conditions
Model boundaries fall into three categories as indicated by
Fig. 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 gra-dients. The Neumann conditions are also assigned on IBs, 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 perfor-mance. At the inlet ports of the channels, velocities, thermo-dynamic states, and mass fractions of mixture are specified according to the desired operating scenarios and mixture of inter-est. Based on the stoichiometric flow ratio, the inlet velocity of the reactant is
vk= ξIWcelLchaRTk,in qFWchaHcha 1−Psat,w,kPk xα,inPk (26)
The mass fractions of the species on the CBs are specified according to operating pressure and fully saturated humidifi-cation 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
wi= xiMi
jxjMJ (27)
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,
Table 1
Cell geometries and operating conditions of base model
Geometry and condition Value Unit
Domain length 50 mm
Domain width (half) 0.8 mm
Gas channel width (half) 0.4 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 (N2/O2) 79/21
Reactant relative humidity 100%
Reactant stoichiometry 3
Cell back pressure 2 atm
Operating temperature 353 K
and a cell total overpotential is set on the cathode BP. Mean-while, 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 conduct-ing 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:
Vcel= Voc+ ηtot (28)
where ηtotis the total cell overpotential and Vocis the open circuit
voltage given by[17]:
Voc= 0.2329 + 0.0025T (29)
2.4. Method of solution
The nonlinearity and coupling in the model formulation eliminate an analytical solution possibility. The forgoing trans-port equations are converted into a conservation form of convection–diffusion conservation equations and numerically solved using the commercial computational fluid dynamic code CFD-ACE+, based on the control volume formulation and the SIMPLE algorithm[24–28].
Three structure mesh systems—60× 30 × 12, 80 × 40 × 16, and 90× 48 × 18 are constructed to explore numerical result dependence on computational cell numbers. Base model geome-tries and operating conditions are listed in Table 1. Model
Table 2
Electrochemical parameters and transport properties
Parameters and properties Value Unit Sources Porosity of the diffusion and catalyst
layer
0.4 [29]
Porosity of the membrane 0.28 [29]
Permeability of the diffusion and catalyst layer
2.3E-11 m2 [26]
Permeability of the membrane 1.0E-18 m2 [26]
Tortuosity of the diffusion and catalyst layer
1.5 [30]
Tortuosity of the membrane 3 [30]
Condensation rate constant 100 s−1 [5]
Evaporation rate constant 100 atm−1s−1 [5]
Concentration dependence of H2 0.5 [31]
Concentration dependence of O2 1 [31]
Transfer coefficients at anode (anodic and cathodic)
0.5 [30,32]
Transfer coefficients at cathode (anodic and cathodic)
1.5 [30,32]
Electrical conductivity of electrode 114 S m−1 [33]
Catalyst loading 0.4 mg cm−2
Catalyst surface area per unit mass 100 m2g−1 [34]
Contact angel 0 ◦ [18]
Exchange current density for anode reaction
1.4E5 A cm−3 [35]
Reference concentration of oxygen 3.39E-6 mol cm−3 [35]
Reference concentration of hydrogen 5.64E-5 mol cm−3 [35]
component electrochemical parameters and transport proper-ties, obtained mainly from open literature[29–35], are listed in
Table 2. Numerical calculations are performed on a Pentium IV 2.4 GHz PC with a 1G RAM. The mesh system (80× 40 × 16) is adopted because current density values are satisfactory with an error range of 2% when using base model geometries and param-eters of these two tables, see Table 3. Based on the reported operating conditions and component geometries of Um and Wang[17], data inTable 4show that the present model to be con-sistent with previous studies. However, as pointed by Sivertsen and Djilali[36], this is not sufficient to validate computational model, comparison with experimental data providing detailed information of key quantities is essential in the future.
3. Results and discussion
Transport component design effects on physical quantity dis-tributions and PEMFC performance are analyzed by numerical
Table 3
Results of three computation grid systems based on the parameters inTables 1 and 2
Cell voltage (V) Grid
60× 30 × 12 80× 40 × 16 90× 48 × 18
Current density (A cm−2)
CPU time Relative error (100%)
Current density (A cm−2)
CPU time Relative error (100%) Current density (A cm−2) CPU time 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 0.67 0.433 6113 3.81 0.444 31270 1.44 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
Table 4
Comparison of present study with results from literature based on the parameters in Ref.[17]
Cell voltage (V) Method
Numerical[17] Present study Experiment[17]
Current density (A cm−2) Relative error (100%) Current density (A cm−2) Relative error (100%) Current density (A cm−2)
0.47 0.998 8.71 0.889 3.16 0.918
0.57 0.856 23.70 0.659 4.77 0.692
0.67 0.462 4.74 0.428 11.75 0.485
0.77 0.179 26.03 0.219 9.50 0.242
0.87 0.035 12.90 0.032 3.23 0.031
method. The influences of channel aspect ratio and GDL thick-ness are described in detail in the following discussion.
3.1. Effects on cell performance
Flow channel geometry is important for optimizing cell per-formance with delivery of sufficient reactant to the reaction sites. Various configurations with different channel aspect ratios are investigated to explore the effects of channel geometry design. The aspect ratio (AR) is defined as the height (in Z-direction) of the channel divided by the width (in Y-direction) of the channel. Channel height and width are modified according to the AR, such that the channel cross-sectional area is fixed. Hence, the reactant mass flow rate at the inlet and the channel hydraulic diameter are the same in all case studies. To maintain a fixed cell width, the shoulder width is also changed. Detailed chan-nel geometries for each considered AR are listed in Table 5.
Fig. 2(a) shows the relationships between output current den-sities and AR at moderate and low cell voltages. Clearly, cell performance dependence on the channel aspect ratio varies with operating conditions. The cell output current density is larger at moderate cell voltage and performance declines rapidly in the high reaction rate region, when the AR is large as in the case of 1.5. However, despite poor performance at cell potentials of 0.42–0.62 V when the channel is flat, such as at AR = 0.5, current density exceeds that obtained with other channel geometries at 0.14 V.
Fig. 2(b) displays cell current density curves with five GDL thicknesses. A thinner diffusion layer is generally believed to reduce the reactant vertical diffusion path such that more oxy-gen arrives at the catalyst layer to gain more current. In contrast, cell performance varies with increased thickness of the diffu-sion layer because the reactant is transported through a longer distance and the flooding problem occurs. However, the data in
Fig. 2(b) reveals that optimal thickness increases from 152m at
Table 5
Channel and shoulder geometries for various ARs (unit: mm)
AR Channel height Half channel width Half shoulder width
0.5 0.57 0.57 0.23
0.75 0.69 0.46 0.34
1 0.080 0.40 0.40
1.25 0.089 0.36 0.45
1.5 0.98 0.33 0.47
0.62 V to 254m at 0.22 V, beyond which potential, it decreases. The tendency is slightly different from the work of Jeng et al.[9]
with a monotonously decreased performance when GDL thick-ness increases for a low-porosity GDL. By using a high-porosity GDL, however, the trends coincide.
Fuel cell output is controlled by three degrading mechanisms—loss associated with reactant activation energy at the electrodes; loss associated with transport of the current in the electrically conducting media, and loss associated with the transport of reactant gases in the cell. These are called activation overpotential, ohmic overpotential, and concentration overpo-tential, respectively. The AR and GDL thickness influences relative magnitudes of these three irreversibilities according to the operating voltage, as is obvious from the data inFig. 2. Local variations of these factors may also play an important
Fig. 2. Current density distribution at various cell voltages as function of (a) AR; (b) GDL thickness.
role on the trend transition shown in the figure. A fixed cell overpotential is notably designated on the cathode BP surface in the present study. The current conservation equation solution provides detailed information of potential distribution through-out the domain. Therefore, ohmic overpotential and activation overpotential can be resolved and compared with the oxygen concentration for different operating voltage. Without general-ity loss, representative locations in the cell domain are assigned to investigate the effects of these two parameters and various transport properties on cell performance variation in the follow-ing sections.
3.2. Effects on membrane conductivity and phase potential
A series of demonstrations involving the most important transport quantities, such as phase potentials, reactant concen-trations at selected locations, and cell voltages are performed to elucidate the cause of the forgoing results. All of the data are mirrored to yield results for a complete channel and a pair of half shoulders adjacent to the channel. Notably, the negatively charged electrons move from low potential to high potential and the positively charged protons move oppositely. The main con-cern in the potential variation discussion is that the passages with lowest electron potential increase from the cathode boundary, and the ionic potential decrease from the anode catalyst layer. Therefore, ohmic overpotential decreases, and the absolute value of the electrochemical reaction driving force, the activation over-potential, in Eq.(14)increases.Fig. 3plots and compares the
Fig. 3. Distributions of potential and conductivity of membrane in the transverse direction of the middle X–Y plane for three values of AR at cell voltages of (a) 0.62 V; (b) 0.14 V. (() AR = 0.5, () AR = 1 and () AR = 1.5). The short bars are the interfaces between channel and shoulder.
membrane potential and conductivity in the transverse direc-tion of the middle X–Y plane of the membrane for AR = 0.5, 1, and 1.5, and cell voltages of 0.62 and 0.14 V at Z-coordinate of 2.152 mm. The short bars on the plot indicate the interfaces between channel and shoulder for each AR. At a cell voltage of 0.62 V, the membrane conductivity patterns in the channel region and the shoulder differ, according to the AR value. At the mem-brane location of the channel region, conductivity is greater at a lower AR. Membrane conductivity is a function of water activ-ity, so increasing the channel transverse dimension at AR = 0.5 facilitates water transport at the anode through the channel to the membrane, causing water activity in the central region to exceed that in the shoulder region. Trends are reversed at the other AR values because of local current density effect. Membrane phase potential variations are strongly related to the conductivities; a lower conductivity is responsible for a larger membrane ohmic loss at a moderate reaction rate at which the local change in cur-rent density is expected to be small. However, the data inFig. 3(b) indicate that at a high reaction rate, higher membrane conduc-tivity locations exhibit a larger membrane ohmic loss because the local current density varies markedly at this cell voltage. The ohmic law and the fact that a large current density variation outweighs the trivial local conductivity fluctuation demonstrate that ohmic loss in the membrane phase is consistent with local current density.
Fig. 4plots the membrane phase conductivity and potential at GDL thicknesses of 152, 203 and 356m and cell voltages of 0.62 and 0.14 V, to examine GDL thickness and cell
volt-Fig. 4. Distributions of potential and conductivity of membrane in the transverse direction of the middle X–Y plane for three values of GDL thickness at cell voltage of (a) 0.62 V; (b) 0.14 V. ((), () and () are denoted as 152, 203 and 356m).
age effects on the various transport properties. The data were obtained from the same position as inFig. 3. A thick GDL sta-bilizes transverse variation of local conductivity. Reducing this parameter leads to a highly non-homogeneous local conductiv-ity distribution at a cell voltage of 0.14 V. The data inFig. 4
indicates that conductivity is related to membrane potential in a manner similar to that inFig. 3. At a low reaction rate, these two variables vary oppositely, but at a high reaction rate, they vary similarly.
3.3. Effects on solid phase potential and activation overpotential
Solid phase potential and activation overpotential in the trans-verse direction for various AR at the interface between the cathode catalyst layer and the GDL are described inFig. 5. The data inFig. 5(a) reveal that the shoulder area exhibits a small ohmic overpotential from the outer surface of the cathode BP. The longer electron passage in the catalyst layer beneath the channel region corresponds to greater potential variation on the way to this location. Consequently, the activation overpotential absolute value in the shoulder region exceeds that in the channel region, indicating that the electrochemical reaction driving force is stronger there. This finding of non-uniform activation over-potential is consistent with that of Kulikovsky et al.[8]and Sun et al.[11]. Comparing the potential variation effects in the three
Fig. 5. Transverse distributions of solid phase potential and activation over-potential at the interface between the cathode catalyst layer and the GDL for three values of AR at cell voltage of (a) 0.62 V; (b) 0.14 V. (() AR = 0.5, () AR = 1 and () AR = 1.5). The short bars are the interfaces between channel and shoulder.
selected designs of the channel AR, show that a slender chan-nel design exhibits a smaller ohmic overpotential, and a larger absolute activation overpotential at all positions of interest. This result is explained by the wide electron transport passage along the shoulder height and the reduced average distance between the shoulders and the channel center, leading to a reduced current resistance. Variation of these two potentials at 0.14 V appears ini-tially to be similar to that inFig. 5(a). The channel with AR = 1.5 has a stronger driving force of electrochemical reaction in the shoulder region as it has the smallest membrane and solid phase ohmic overpotentials compared to the other two channel designs.
Fig. 6compares the transverse variation in the solid phase potential and the activation overpotential at the location between the cathode GDL and the catalyst layer at three GDL thick-nesses and two operating potentials. The GDL thickness effects on ohmic overpotential is clearly demonstrated in the data in
Fig. 6(a). Despite the fact that the catalyst layer shoulder area exhibits minor ohmic overpotential in the design with a thickness of 152m, the potential increases abruptly toward the channel central region, because the height and cross-sectional area of the electron transverse transport passage are small in this thinnest design. For the design with the thicker GDL, the larger ohmic overpotential in the catalyst layer shoulder region is problematic, but the moderate potential variation toward the central chan-nel is advantageous. Accordingly, the activation overpotential variation exhibits the same tendency as that of the solid phase potential. These situations are similar at the two operating
volt-Fig. 6. Transverse distributions of solid phase potential and activation overpo-tential at the interface between the cathode catalyst layer and the GDL for three GDL thicknesses at cell voltage of (a) 0.62 V; (b) 0.14 V. ((), () and () are denoted as 152, 203 and 356m).
age, but the absolute values and variations are greater in the 0.14 V case.
3.4. Effects on water saturation level
Channel geometry effects on the transverse saturation level of the cathode GDL along the channel direction is presented in
Fig. 7. The flat channel advantage is evidenced by liquid water accumulation. The relatively short distance through which water is transported at AR = 0.5 leads to the low saturation level in the channel direction even at a high reaction rate. That means the fast diffusion associated with the flat channel design causes the local partial pressure of water vapor to be low. In contrast, sat-uration level at the rear section of the channel with AR = 1.5 is high, as the water vapor cannot easily escape from GDL under the shoulder. In this scenario, effective pore space in the porous medium is reduced considerably, and more mass transport over-potential is activated, so cell performance is drastically degraded at low cell voltage.
Fig. 8plots the saturation level in the transverse direction of the cathode GDL at three positions x = 0.05, 0.025 and 0.045 m, to elucidate GDL thickness effects on liquid water distribution. The plot demonstrates that the thinner GDL design has the lowest saturation level since water vapor produced in the cell reaction easily transports to the channel and outside the cell. On the con-trary, the longer path and water vapor lower diffusion rate in the design with the 356m-thick GDL results in greater water
Fig. 7. Effect of channel aspect ratio on transverse water saturation of the cath-ode GDL along channel direction at cell voltage of (a) 0.62 V; (b) 0.14 V. (() AR = 0.5, () AR = 1 and () AR = 1.5). Right half portion is the enlargement of the crowded lines.
Fig. 8. Effect of GDL thickness on transverse water saturation of the cathode GDL along channel direction at cell voltage of (a) 0.62 V; (b) 0.14 V. ((), () and () are denoted as 152, 203 and 356 m). Right half portion is the enlargement of the crowded lines.
vapor concentration with a high probability of over-saturation and liquid water formation, contributing to concentration over-potential and performance reduction, especially at x = 0.45 m at a cell potential of 0.14 V.
Note that temperature variation is not considered in this study, model result may over predict saturation level at high reaction rate because the saturation pressure would be greater accom-panied by larger local temperature. However, the membrane water content adjacent to the high temperature region is ele-vated due to isothermal assumption. This could enhance ionic current transport and compensate the loss from over-predicted saturation level. Nevertheless, this highlights the necessity of conducting a non-isothermal analysis to interpret these mutual-coupled effects.
3.5. Effects on oxygen concentration and local current density
Fig. 9depicts oxygen mass fraction distribution and local cur-rent density under the same condition and at the same location as inFig. 5. Intuitively, the flat channel supplies more oxygen to the catalyst layer and a high current density is expected. However, calculations reveal entirely different trends between channel and shoulder region catalyst layers. The design with AR = 0.5 gen-erates the lowest local current density in the channel region, despite its having the largest oxygen mass fraction, and the expected relationship between reactant concentration and the
Fig. 9. Transverse distributions of oxygen mass fraction and local current density at the interface between the cathode catalyst layer and the GDL for three values of AR at cell voltage of (a) 0.62 V; (b) 0.14 V. (() AR = 0.5, () AR = 1 and () AR = 1.5).
reaction rate is not observed except at the shoulder center. On average, the slender design with AR = 1.5 generates more cur-rent at moderate cell voltage. This phenomenon is explained by the fact that in this scenario, the reaction is relatively slow and the high reactant concentration is not as important as low ohmic overpotential and high activation overpotential provided by slender channel geometry such as AR = 1.5. This design has a wider rib zone than those of other designs, providing a small increase and variation in the solid phase potential as well as large absolute activation overpotential, causing high current density and favorable cell performance. At a high reaction rate as shown inFig. 9(b), the dominant mechanism of local current density transits from electric potential at the channel region, to oxygen concentration at the shoulder region at AR = 1.5. This transition is also found when AR = 1. The design with AR = 0.5 exhibits the same local current density variation trend as in the case of 0.62 V. At AR = 1.5, a high reactant concentration need outweighs a high activation overpotential need, so the expected relationship between the concentration and the local reaction rate appears earlier. The large shoulder area width hinders reactant transport, resulting in a relatively low level of local oxygen mass fraction and a sharp decline in the local current density under the shoul-der area. In contrast, the sufficient oxygen provided by the flat channel such as AR = 0.5 causes most of the region to exhibit a potential controlled state and on average, produces a greater current than that generated by other designs.
GDL thickness effects on transverse distribution of oxygen concentration and current density is plotted inFig. 10. Exhibiting
Fig. 10. Transverse distributions of oxygen mass fraction and current density at the interface between the cathode catalyst layer and the GDL for three GDL thicknesses at cell voltage of (a) 0.62 V; (b) 0.14 V. ((), () and () are denoted as 152, 203 and 356m).
a trend opposite to that of the potential, a thinner GDL provides more oxygen to the channel region catalyst layer. The thickest design is associated with a greater concentration in the shoul-der region. These variations are related to the vertical depth and transverse cross-section of reactant delivery. In the case in which GDL thickness is 152m, the vertical path is short and oxygen concentration is high at the catalyst layer beneath the channel. Nevertheless, the transverse transport cross-section is reduced and oxygen concentration falls substantially at the shoulder region. With reference to local current density dis-tribution, mechanisms of activation overpotential mechanism and oxygen concentration have different effects in different regions. The activation overpotential dominates the reaction in the catalyst layer under the channel region because oxygen con-centration fulfills the requirement for electrochemical reaction. Accordingly, the current density trend is consistent with that of the activation overpotential throughout the entire region in the 356m case and in most of the region in the 203 and 152 m cases at 0.62 V. A peak point in the concave pattern of the local current density appears only near the central shoulder region in the 152m cases, suggesting that from this point to the central shoulder, oxygen deficiency forces the current density to drop according to the oxygen mass fraction, and local performance is dominated by oxygen concentration. Nevertheless, most current density gain in the shoulder region catalyst layer arises from higher activation overpotential for the design with the 152 m-thick GDL, compensating the loss at the central region and generating more current than other designs at 0.62 V.
Similar result can be found in the work of Sun et al.[10]which investigated the GDL thickness effects at total call overpotentials of−0.4, −0.5 and −0.6 V. The thinner GDL thickness design exhibits the slowest reaction rate at the catalyst site under the channel region. This phenomenon was interpreted by a longer electron traveling length to the channel region than to the shoul-der region, such that fewer electrons would participate in the electrode reaction in the region under the channel.
At a high reaction rate with a cell voltage of 0.14 V as shown inFig. 10(b), oxygen concentration is lower than inFig. 10(a) because the reaction rate is higher, but transport depth and transverse cross-section effects on oxygen concentration are the same. However, local current density distributions exhibit quite different patterns. The plot depicts clear peak points in the con-cave pattern for all three GDL thicknesses. As stated previously, the 356m design generates more current than the other two designs at the central channel region catalyst layer because the activation overpotential dominates the reaction. This situation changes rapidly toward the shoulder region since oxygen con-centration is extremely low. The figure also shows that the peak current density points move near the shoulder in the remain-ing two cases. A thinner GDL corresponds to a latter peak point from the channel central and a larger peak current den-sity. Beyond these points, current density variation is nearly the same as oxygen concentration. At the peak points, the design with the 152m-thick GDL clearly has the largest local cur-rent density; unfortunately, oxygen concentration falls sharply toward the shoulder center and current density is the lowest of the three designs. Consequently, the 152m-thick case does not exhibit the best cell performance. On the other hand, despite the fact that the 203m-thick design has a moderate local current density output at the channel center, peak point, and shoulder center, its local current density curve envelopes the largest area with the X–Y plane and becomes the optimal GDL thickness at low cell voltage. This finding is quite different from the result of Sun et al.[10]that no significant variation of average current density was found when GDL thickness reduced from 0.35 to 0.15 mm. The discrepancy may be attributed to different oper-ating conditions or cell geometries between these two models.
4. Conclusions
A multi-dimensional, multi-component, computational fluid dynamic model was employed to elucidate transport phenom-ena and electrochemical reaction in PEMFCs. Two important transport component design parameters—channel aspect ratio and GDL thickness, are investigated in detail. The simulation and discussion support the following conclusions:
1. The designs of channel aspect ratio and GDL thickness affect the physical property distributions. At a high reaction rate, these two parameters have a strong influence on cell perfor-mance.
2. Three mechanisms of cell irreversibility are resolved locally from the model and are considered to determine variation in the macroscopic cell current density and performance.
3. Transverse current density distribution is governed by both electron conduction and activation overpotential in channel region or by reactant concentration in shoulder region of cat-alyst layer. The relative strengths of these two factors depend on the operating voltage and the transport geometry. 4. At a moderate reaction rate, the transverse direction current
density in most regions is a convex function of position and is influenced by solid phase potential and activation overpo-tential, making the peak point of the concave pattern current density either close to the shoulder region or non-existent. Thus, a geometry design that facilitates electron transport such as a large channel aspect ratio or a thin GDL thickness causes a relatively larger current in the catalyst layer. 5. At a low cell voltage, the largest reaction rate location is close
to the channel center, therefore, the electrochemical reaction in the majority regions is dominated by the reactant transport and a smaller channel aspect ratio design is preferred. 6. Large GDL thickness positively increases oxygen transport
under the shoulder region; therefore, output current density elevates according to cell voltage decrease. However, at the lowest considered cell voltage of 0.14 V, oxygen deficiency caused by long traveling length and liquid water clogging effect reverses this relationship.
Acknowledgements
The authors would like to acknowledge the National Science Council, Taiwan, the Republic of China, for financially support-ing this research under Contract No. NSC 93-2212-E-009-001. Furthermore, the authors would like to thank Dr. Sheng-Rui Jian, “Institute and Department of Electrophysics, National Chiao Tung University”, for his helpful discussion.
References
[1] H. Tsuchiya, O. Kobayashi, Int. J. Hydro. Energy 29 (2004) 985–990. [2] L. Giorgi, E. Antolini, A. Pozio, E. Passalacqua, Electrochim. Acta 43
(1998) 3675–3680.
[3] L.R. Jordan, A.K. Shukla, T. Behrsing, N.R. Avery, B.C. Muddle, M. Forsyth, J. Power Sources 86 (2000) 250–254.
[4] W.K. Lee, C.H. Ho, J.W. Van Zee, M. Murthy, J. Power Sources 84 (1999) 45–51.
[5] W. He, J.S. Yi, T.V. Nguyen, AIChE J. 46 (2000) 2053–2064.
[6] D. Natarajan, T.V. Nguyen, J. Electrochem. Soc. 148 (2001) A1335–A1342. [7] A.C. West, T.F. Fuller, J. Appl. Electrochim. 26 (1996) 557–565. [8] A.A. Kulikovsky, J. Divisek, A.A. Kornyshev, J. Electrochem. Soc. 146
(1999) 3981–3991.
[9] K.T. Jeng, S.F. Lee, G.F. Tsai, C.H. Wang, J. Power Sources 138 (2004) 41–50.
[10] W. Sun, B.A. Peppley, K. Karan, J. Power Sources 144 (2005) 42–53. [11] W. Sun, B.A. Peppley, K. Karan, Electrochim. Acta 50 (2005) 3359–3374. [12] S. Dutta, S. Shimpalee, J.V. Van Zee, J. Appl. Electrochem. 30 (2000)
135–146.
[13] S. Dutta, S. Shimpalee, J.V. Van Zee, Int. J. Heat Mass Trans. 44 (2001) 2029–2042.
[14] T. Berning, D.M. Lu, N. Djilali, J. Power Sources 106 (2002) 284–294. [15] T. Berning, N. Djilali, J. Power Sources 124 (2003) 440–452.
[16] L. Wang, A. Husar, T. Zhou, H. Liu, J. Hydro. Energy 28 (2003) 1263–1272. [17] S. Um, C.Y. Wang, J. Power Sources 125 (2004) 40–51.
[18] Z.H. Wang, C.Y. Wang, K.S. Chen, J. Power Sources 94 (2001) 40–50. [19] A. Parthasarathy, S. Srinivasan, A.J. Appleby, J. Electrochem. Soc. 139
[20] T.E. Springer, T.A. Zawodzinski, S. Gottesfeld, J. Electrochem. Soc. 138 (1991) 2334–2342.
[21] C.Y. Wang, C. Beckermann, Int. J. Heat Mass Transfer 36 (1993) 2747–2758.
[22] C.Y. Wang, P. Cheng, Advance in Heat Transfer, Acadmic Press, San Diego, 1997.
[23] L. You, H.T. Liu, Int. J. Heat Mass Transfer 45 (2002) 2277–2287. [24] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New
York, 1980.
[25] J.J. Hwang, C.K. Chen, R.F. Savinell, C.C. Liu, J. Wainright, J. Appl. Electrochem. 34 (2004) 217–224.
[26] CFD-ACE(U)TM User Manual, CFD Research Corp., Huntsville, AL, 2004.
[27] H.C. Liu, W.M. Yan, C.Y. Soong, F. Chen, J. Power Sources 142 (2005) 125–133.
[28] C.Y. Soong, W.M. Yan, C.Y. Tseng, H.C. Liu, F. Chen, H.S. Chu, J. Power Source 143 (2005) 36–47.
[29] S. Um, C.Y. Wang, K.S. Chen, J. Electrochem. Soc. 147 (2000) 4485–4493. [30] S. Mazumder, J.S. Cole, J. Electrochem. Soc. 150 (2003) A1503–A1509. [31] C. Marr, X. Li, J. Power Sources 77 (1999) 17–27.
[32] F.B. Weng, A. Su, G.B. Jung, Y.C. Chiu, S.H. Chan, J. Power Sources 145 (2005) 546–554.
[33] B. Hum, X. Li, J. Appl. Electrochem. 34 (2004) 205–215. [34] http://www.etek-inc.com/home.php.
[35] P.T. Nguyen, T. Berning, N. Djilali, J. Power Sources 130 (2004) 149–157. [36] B.R. Sivertsen, N. Djilali, J. Power Sources 141 (2005) 65–78.
