- Email: [email protected]
1D and 3D numerical simulations in PEM fuel cells
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
1D and 3D numerical simulations in PEM fuel cells D.S. Falca˜o a, P.J. Gomes b, V.B. Oliveira a, C. Pinho a, A.M.F.R. Pinto a,* a Centro de Estudos de Feno´menos de Transporte, Departamento de Engenharia Quı´mica, Faculdade de Engenharia da Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal b Laboratory of Separation and Reaction Engineering, Departamento de Engenharia Quı´mica, Faculdade de Engenharia da Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal
article info
abstract
Article history:
The potential of fuel cells for clean and efficient energy conversion is generally recognized.
Received 30 April 2011
The proton-exchange membrane (PEM) fuel cells are one of the most promising types of fuel
Received in revised form
cells. Models play an important role in fuel cell development since they enable the under-
21 June 2011
standing of the influence of different parameters on the cell performance allowing a systematic
Accepted 26 June 2011
simulation, design and optimization of fuel cells systems. In the present work, one-
Available online 31 July 2011
dimensional and three-dimensional numerical simulations were performed and compared with experimental data obtained in a PEM fuel cell. The 1D model, coupling heat and mass
Keywords:
transfer effects, was previously developed and validated by the same authors [1,2]. The 3D
PEM fuel cell
numerical simulations were obtained using the commercial code FLUENT e PEMFC module.
Water transport
The results show that 1D and 3D model simulations considering just one phase for the
Modeling studies
water flow are similar, with a slightly better accordance for the 1D model exhibiting
CFD
a substantially lower CPU time. However both numerical results over predict the fuel cell performance while the 3D simulations reproduce very well the experimental data. The effect of the relative humidity of gases and operation temperature on fuel cell performance was also studied both through the comparison of the polarization curves for the 1D and 3D simulations and experimental data and through the analysis of relevant physical parameters such as the water membrane content and the proton conductivity. A polarization curve with the 1D model is obtained with a CPU time around 5 min, while the 3D computing time is around 24 h. The results show that the 1D model can be used to predict optimal operating conditions in PEMFCs and the general trends of the impact on fuel cell performance of several important physical parameters (such as those related to the water management). The use of the 3D numerical simulations is indicated if more detailed predictions are needed namely the spatial distribution and visualization of various relevant parameters. An important conclusion of this work is the demonstration that a simpler model using low CPU has potential to be used in real-time PEMFC simulations. Copyright ª 2011, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
1.
Introduction
Fuel cells are extensively being studied today because of their potential as an alternate energy source for a wide range
of applications. Fuel cells have unique technological attributes: efficiency, absence of moving parts and very low emissions. In particular, the proton-exchange membrane (PEM) fuel cells are today in the focus of interest as one of the
* Corresponding author. Tel.: þ351225081675; fax: þ351225081449. E-mail address: [email protected] (A.M.F.R. Pinto). 0360-3199/$ e see front matter Copyright ª 2011, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2011.06.133
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
12487
Fig. 1 e Schematic representation of the experimental set-up.
most promising technology in development in power generation. The necessity of formulation of adequate mathematical models is based on the need to reach a deeper understanding of internal processes (reactions, mass transport, heat transport) occurring in a fuel cell. These processes cannot be observed directly in experiments and mathematical model could be applied to obtain hints for the development of optimal control and operating strategies. One of the most important operational issues of PEMFC is the water management in the cell [3,4]. Water content of the membrane is determined by the balance between water production and three water transport processes: electro-osmotic drag of water (EOD), associated with proton migration through the membrane; back diffusion from the cathode to the anode and diffusion of water to/from the oxidant/fuel gas streams. Understanding the water transport in the PEM is a key issue to study. In
order to achieve optimal fuel cell performance, an adequate water balance should be ensured because membrane should remain hydrated for sufficient proton conductivity, while cathode flooding and anode dehydration should be avoided [5e7]. Different models were developed in the last decade to describe several water transport mechanisms through the membrane such as Springer et al. [8] using a diffusion model, Bernardi and Verbrugge [9] considering a hydraulic permeation model and Kulikovsky [10] developing a semi analytical model 1D þ 1D taking into account oxygen and water transport across the cell and deriving an expression for the limiting current density. Scott and Argyropoulos [11] developed a onedimensional model of the current and potential distribution in a porous electrode. De Francesco et al. [12] used an analytical approach to calculate the membrane resistance. More recently, there was a trend to apply CFD methods to fuel cells modeling [13,14]. Due to cell configuration, transport
12488
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
phenomena are three-dimensional, so, some works taking in account three dimensions arose [15e22]. In this study, 1D a 3D numerical simulations are presented and compared, in order to verify advantages and disadvantages of each type of modeling. The one-dimensional model was already developed in a previous work [1]. This model considers the effects of coupled heat and mass transfer, along with the electrochemical reactions occurring in PEMFC. The model can be used to predict the hydrogen, oxygen and water concentration profiles in the anode, cathode and membrane, as well as the temperature profile across the cell. The model was validated with published experimental data. This model is also used and validated in another reference [2]. The 3D numerical simulations were obtained in a commercial software FLUENT with the PEMFC module. In the present work, the influence of two-phase flow effects, as well as the influence of reactants humidity and operating temperatures is studied recurring to the models and to experimental data. Besides polarization curves, the presented models simulate several important parameters with impact on fuel cell performance, such as membrane water content and protonic conductivity.
2.
Experimental set-up
2.1.
Apparatus
A schematic drawing of the experimental apparatus used in this work is shown in Fig. 1. Pure hydrogen (humidified or dry) as fuel and air (humidified or dry) as oxidant are used. The pressure of the gases is controlled by pressure regulators (Air e Norgreen 11400, H2 e Europneumaq mod. 44-2262-241) and flow rates controlled by flow meters (KDG e Mobrey). The reactants humidity and temperatures are monitored by adequate humidity and temperature probes (Air e Testo, H2 e Vaisala). The humidification of air and hydrogen gases is conducted in Erlenmeyer flasks by a simple bubbling process. To control the humidification temperature each Erlenmeyer flask is thermally isolated and surrounded with an electrical resistance (50 W/m) activated by an Osaka OK 31 digital temperature controller. The same procedure is applied along the connecting pipes from the humidification point up to the entrance of the fuel cells to guarantee the temperature stabilization of each reacting gas flow, as well as to control the operating temperature of the fuel cell. For the measurement and control of the cell electrical output, an electric load reference LD300 300 W DC Electronic Load from TTI is used. This load was connected to a data acquisition system composed by Measurement Computing boards installed in a desktop computer. The used data acquisition software was DASYLab.
2.2.
Fuel cell design
In the present work, all the components of the PEMFC were “in house” designed, except the MEA. A Nafion 111 from Dupont
Fig. 2 e Representation of the serpentine flow channels.
membrane with 25 cm2 active surface area is used. The channels configuration used for the anode and cathode flow channels is represented in Fig. 2. The channels depth is 1 mm for hydrogen side and 1.5 mm for air side. Table 1 summarizes the cell design specifications.
2.3.
Experimental conditions
The studied operating conditions were the reactants humidity and the operating temperature. All the experiments were conducted at atmospheric pressure. Conditions are presented in Tables 2 and 3.
3.
Numerical simulations
3.1.
One-dimensional modelling
The 1D model presented here for comparison with the 3D numerical simulations was already developed and validated. The detailed description of the model can be found in references [1] and [2]. The main assumptions are recalled here: - Steady-state and one-dimensional heat and mass transport; - Ideal gas mixture assumed;
Table 1 e Fuel cell design parameters. Gas difusion layer Catalytic layer
Channels configuration Membrane
Anode Cathode Anode
Thickness: 0.0195 cm Thickness: 0.0195 cm Catalyst Loading: 0.5 mgPt/cm2 Thickness: 0.006 cm Cathode Catalyst Loading: 0.5 mgPt/cm2 Thickness: 0.006 cm Anode Serpentine Cathode Serpentine Nafion 111, Thickness: 0.00254 cm
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
12489
Table 2 e Operating conditions used to analyse reactants humidity effect. Operating temperature (K) 333
H2 flow rate
Air flow rate
Calculated for Calculated for z ¼ 1 for 1 A/cm2 z ¼ 1 for 1 A/cm2
H2 RH (%)
Air RH (%)
95 0 95
95 95 0
- Heat and mass transport through the gas diffusion and catalyst layers assumed to be a diffusion-predominated process (negligible convection effects); - Effective Fick models for the mass transport in the diffusion layers and membrane are considered; - The thermal energy model is based on the differential thermal energy conservation equation (Fourier’s law); - The thermal conductivity for all the materials is assumed to be constant; - The anodic and cathodic overpotential is constant through the catalyst layers; - Heat generation or consumption is considered in the catalyst layers; - The heat generated by the current flow (Joule effects) through each cell layer is ignored; - Water transport through the membrane assumed to be a combined effect of diffusion and electro-osmotic drag; - Membrane proton conductivity is a function of l, the number of water molecules per ionic group; - Local equilibrium at interfaces is represented by partition functions; - Reactions in the catalyst layers are considered as homogeneous; - Kinetics of the cathode is described by a Tafel expression; - Kinetics of the anode is described by a Tafel expression; - Anode and cathode flow channels are treated as a continuous stirred tank reactor (CSTR), so, the composition and temperature inside the channels are uniform; - Anode and cathode streams act as heat transfer fluids removing heat from the cell at the exit temperatures.
3.2.
Fig. 3 e Representation of the flow channels Mesh.
3.2.1.
Geometry and mesh generation
The geometry of the fuel cell includes the MEA and the channel plates. Terminal plates and collector plates are not considered, due to simplification (shorter simulation times). The used mesh is not uniform, with elements type Hex/Wedge Cooper and around 1,600,000 elements. Mesh in the channels is more refined, because the channel dimensions are very small. Due to this dimensional difference between channels and the other fuel cell components, the mesh generation is complicated. A representation of the mesh in the channels is shown in Fig. 3. Concerning mesh division, in terms of depth, the anode channel is divided in 2 (depth of 0.7 mm) and the cathode channel is divided in 5 (depth of 1.5 mm). In the MEA, the mesh is also built with high refinement, because the layers are very thin. The scheme used in mesh division is presented in Fig. 4. The mesh resolution can affect the solution precision and therefore a grid independency test was done. In order to do that, a more refined mesh was created, with 2,400,000
Three-dimensional model formulation
The 3D numerical simulations were obtained with commercial code FLUENT, in particular using the module PEFC. Details of the implementation and resolution of equations can be found in [23]. Geometry and mesh of the fuel cell had to be generated.
Table 3 e Operating conditions used to analyse operating temperature effect. Operating temperature (K) 298 313 333
H2 flow rate
Air flow rate
Calculated for Calculated for z ¼ 1 for 1 A/cm2 z ¼ 1 for 1 A/cm2
H2 RH (%)
Air RH (%)
95
95
Fig. 4 e General indication of the grid of the MEA mesh.
12490
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
elements and simulations with this new mesh were done. The results show variations around 0.03% relatively to those obtained with regular mesh. So, it could be concluded that results were independent of grid refinement.
3.2.2.
Boundary conditions
Boundary conditions need to be defined for PEM fuel cell simulations, based on problem specification. In this model, the following boundary conditions were defined as: - Inlet (mass flow inlet type) and outlet (pressure outlet type) zones for the anode and cathode gas channels - Surfaces (wall type) representing anode and cathode channels The following continuum zones were also required: -
Anode and cathode flow channels (fluid zone) Anode and cathode current collectors (solid zone) Anode and cathode diffusion layers (fluid zone) Anode and cathode catalyst layers (fluid zone) Membrane (fluid zone)
boundary, there is a zero flux boundary condition for the membrane phase potential on all outside boundaries. For the solid phase potential, there are external boundaries on the anode and the cathode side that are in contact with the external electric circuit. Through these boundaries passes the electrical current generated in the fuel cell. On all other external boundaries there is a zero flux boundary condition solid phase potential. Anode potential was set to zero, so the positive value prescribed on the cathode side is the cell voltage. The transfer currents, or the source terms in Equations (1) and (2), are non-zero only inside the catalyst layers and are computed as: - For the solid phase, RS ¼ Ra (<0) on the anode side and RS ¼ þRc (>0) on the cathode side. - For the membrane phase, RM ¼ þRa (>0) on the anode side and RM ¼ Rc (<0) on the cathode side. Source terms on Equations (1) and (2) are also called exchange current density and have the generic Tafel formulation [23]:
Boundary conditions definition was done in Gambit (the mesh generator software used). After geometry and mesh generation and boundary and continuum zones established, the mesh was found to be suitable to be exported to FLUENT.
Ra ¼
jref a
Rc ¼
jref c
CACL H2
!ga
ref
CH2
CCCL O2
!gc
aa Fh =RT e a
ref
CO2
eac Fhc =RT
(3)
(4) ref
3.2.3.
Model equations
In the PEM fuel cells module of Code FLUENT, two electrical potential fields are solved: one in the membrane and catalytic layers and the other in catalytic layers, diffusion layers and current collectors. Surface reactions occurring on the porous catalyst region are solved and the reaction diffusion balance is applied to compute the rates. Based on a cell voltage prescribed, the current density value is computed.
3.2.3.1. Electrochemistry model. The most important electrochemistry aspects are hydrogen oxidation and oxygen reduction reaction rates. In this model, these electrochemical processes are treated as heterogeneous reactions taking place on the catalyst surfaces. The driving force behind these reactions is the surface overpotential: the difference between the phase potential of the solid (ES) and the phase potential of the electrolyte/membrane (EM). Therefore, two potential equations are solved, one for electrons transport and the other for hydrogen protons transport: VðsS VES Þ þ RS ¼ 0
(1)
VðsM VEM Þ þ RM ¼ 0
(2)
where R is the volumetric transfer current [A/m3], s is the electrical conductivity [S/m] and the indices S and M refers to, respectively, solid phase and membrane phase. For the resolution of these equations resolution, several boundary conditions need to be defined. There are two types of external boundaries: those through which an electrical current passes and those through which no current passes. As no protonic current leaves the fuel cell through any external
where CACL H2 is the hydrogen concentration in catalytic layer, CH2 is the reference hydrogen concentration, CCCL O2 is the oxygen ref concentration in catalytic layer, CO2 is the reference oxygen ref ref concentration, ja and jc are, respectively, the anode and cathode volumetric reference exchange current density, aa and ac are, respectively, anode and cathode transfer coefficient, ga and gc are, respectively, the anode and cathode concentration dependence, and ha and hc are, respectively, the anode and cathode overpotentials. These overpotentials are the difference between the solid and membrane potentials. The gain in electric potential from crossing from the anode to the cathode side can then be taken into account by subtracting the open circuit voltage, E0, on the cathode side: ha ¼ ES EM
(5)
hc ¼ ES EM E0
(6)
Adequate combination of Equations (1)e(6) allows the resolution of the two potential fields.
3.2.3.2. Current and mass conservation. The volumetric source terms (S) for the species equations [kg/m3s] and energy equation [W/m3] are given in Equations (7)e(9) [23]: MH SH2 ¼ 2 Ra 2F
(7)
MO SO2 ¼ 2 Rc 2F
(8)
MH2 O Rc 2F
(9)
SH2 O ¼
where M is the molecular weight [kg/mol].
12491
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Additional volumetric sources (Senergy) to the energy equation implemented in the FLUENT - PEMFC model include ohmic heating (Rohm), heat of formation of water (Hreac), electric work and latent heat of water (Hphase). Senergy ¼ I2 Rohm þ Hreac þ hRa;c þ Hphase
(10)
The electrochemical reactions are treated as surface reactions in the two catalyst layers, and it is assumed that the diffusive flux of any reacting species is balanced by its rate of production: rDj Mj yj;sur yj;cen r ¼ Ra;c d nF
(11)
where r is the specific reacting surface area of the catalyst layer [m1], yj, sur is the mass fraction of species j at the reactant surface, yj, cen is the mass fraction of species j at the cell center and d is the average distance between the reaction surfaces and the cell center [m]. The left hand side of Equation (11) represents the diffusive flux at the reacting surface and the right hand side represents the rate of mass generation. The average distance from the cell center to the reacting surface is estimated as d ¼ 1/r. Equation (11) is used to obtain the surface values of H2 and O2 concentrations, applying a Newtonian solution procedure. These surface values are then used to compute the rates in Equations (3) and (4).
3.2.3.3. Liquid water formation, transport, and its effects. Since PEM fuel cells operate under relatively low temperature (<100 C), the water vapor may condense to liquid water, especially at high current densities. While the existence of the liquid water keeps the membrane hydrated, it also blocks the gas diffusion passage, reduces the diffusion rate and the effective reacting surface area and hence the cell performance. To model the formation and transport of liquid water, FLUENT uses a saturation model. In this approach, the liquid water formation and transport is governed by the following conservation equation for the volume fraction of liquid water, s, or the water saturation: vðerl sÞ þ V$ðrl vl sÞ ¼ rw vt
(12)
where the subscript, l, represents liquid water, vl represents liquid velocity, which is equivalent to the gas velocity inside the gas channel and rw is the condensation rate, which is defined as: PH O Psat MH2 O ; ½srl rw ¼ cr max ð1 sÞ 2 RT
scosqc pc ¼ 0;5 1;417ð1 sÞ 2;12ð1 sÞ2 þ1;263ð1 sÞ3 ; for qc < 90+ K e (15) scosqc pc ¼ 0;5 1; 417s 2; 12s2 þ 1; 263s3 ; for qc > 90+ K e
(16)
where qc is the contact angle and s is the surface tension [N/m2]. Equation (12) is used to model various physical processes such as condensation, vaporization, capillary diffusion, and surface tension. The clogging of the porous media and the flooding of the reaction surface are modeled by multiplying the porosity and the active surface area by (1-s).
3.2.4.
Solution procedure
Solver default definitions from FLUENT are not appropriate to reach solution convergence. Convergence is here considered when the difference between two consecutive iterations is less than 0.001 for the residuals of continuity, velocity and concentrations equations and for UDFs and less than 0.0001 for the residuals of the energy equation. Multigrid cycle was changed for F-Cycle for all equations. Stabilization method was changed for BCGSTAB for species concentrations, water saturation, electric and protonic potential. Terminal restrictions were also changed for 0.001 for species concentrations, water saturation and for 0.0001 for electric and protonic protential. Maximum number of cycles was increased for 50. Under-relaxation factors were changed to the following values: 0.45 for pressure, 1 for density, 0.9 for body forces, 0.3 for momentum, 0.9 for species concentration, 0.2 for energy, 1 for electric potential, 0.95 for protonic potential and 0.5 for water saturation and water content. The discretization method was changed to second order (leading to more precise solutions) for all equations.
3.2.5.
Model parameters
The model parameters are already defined in the PEM Fuel Cells module but could be set by the user. All default parameters were not changed, except those related with the kinetics, presented in Table 4. Reference concentrations were set to the concentrations used in each experiment.
(13)
where rw is added to the water vapor equation, as well as the pressure correction. This term is not applied inside the membrane. The condensation rate constant is defined as cr ¼ 100/s. Inside the highly-resistant porous zones, the use of the capillary diffusion term allows us to replace the convective term in Equation (12): vð˛rl sÞ Ks3 dpc Vs ¼ rw þ V$ rl ml ds vt
where K is the absolute permeability and pc the capilarity pressure. Depending on the wetting phase, the capillary pressure is computed as a function of s (the Leverett function) [23]:
(14)
Table 4 e Model kinetic parameters. Parameter E0 ref ja ref jc aa ac ga gc
Value
Reference
0:0025 T þ 0:2329 V 1 109 A/m3 2 104 exp [0.014189 (T-353)] A/m3 0.5 0.2 0.5 1.0
[18] [19] [19] [20] [21] [22] [22]
12492
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
1.2 Exp
Mod 3D Two-Phase
Mod 3D One-Phase
Mod 1D
Voltage (V)
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Current Density (A/cm2 )
Fig. 5 e Simulated and experimental polarization curves for 333 K, both gases RH of 95%, H2/Air flow rates with stoichiometry flow ratios za [ 1 and zc [ 2 for 1 A/cm2.
4.
Results and discussion
As already referred, the main objective of the present work was the comparison between 1D and 3D modeling, bearing in mind the demonstration that a simpler model using low CPU has potential to be used in real-time PEMFC simulations. Particular attention was devoted to the impact of considering two-phase flow effects in the 3D simulation results. The effect of the relative humidity of gases and operation temperature in fuel cell performance was also studied. Preliminary tests showed that the results using Tafel or Buttler-Volmer kinetic equations were very similar, consequently, the Tafel equation leading to lower computing times was used in all the simulations here presented.
4.1.
Two-phase flow effects
In Fig. 5, the 3D model simulations considering one and two phases for water flow and the 1D model simulations are
compared with experimental data, to evaluate the impact of two-phase flow effects on simulation results. The data plotted refer to an operation temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios of za ¼ 1 and zc ¼ 2 for 1 A/cm2. As can be seen from Fig. 5, the 3D simulation results considering two-phase flow effects compare very well with the experimental data. For this operating condition, the formation of liquid water is likely to occur, partially filling the pores of the electrodes structure and leading to a lower fuel cell performance. The consideration of only water vapor, in the 1D and 3D models (without two-phase flow effects) lead to higher predicted fuel cell performances, as expected, but reproduce however worse the experimental data. It is worth mentioning that these 1D and 3D model simulations considering just one phase are similar with a slightly better accordance for the 1D model exhibiting a substantially lower CPU time. In PEMFC literature, validation is commonly performed against polarization curves and using a significant number of fitting parameters. Furthermore, there are still many technical difficulties in obtaining detailed and reliable in situ experimental data. Consequently, verification of model performance based on polarization curves will persist but should however be complemented with an analysis of the physical results obtained with the model. Given the complexity of the transport phenomena involved and the uncertainty in determining some of the physico-chemical parameters, it is any case unrealistic to expect accurate quantitative agreement between experiments and computational model simulations [24]. So, it is essential to analyze the results in terms of the physical parameters predicted by the models. Three-dimensional numerical simulations enable the prediction of the distribution and visualization of various relevant parameters influencing the fuel cell behavior, namely, the mass or mole fractions of the different species in the channels and the membrane water content. These parameters change with voltage variations and therefore, the
Fig. 6 e Contours of hydrogen mole fraction in an XY plane at half of the anode channels depth for a) 0.8 V, b) 0.6 V. Operating conditions: operating temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios of za [ 1 and zc [ 2 for 1 A/cm2.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
12493
Fig. 7 e Contours of the oxygen mole fraction in an XY plane at half of the cathode channels depth for a) 0.8 V, b) 0.6 V. Operating conditions: operating temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and zc [ 2 for 1 A/cm2.
contours for two values of voltage (0.8 and 0.6 V) are presented in Fig. 6 for an operation condition revealing good agreement with experimental data (3D simulation considering two-phase flow effects). This figure shows the contours of hydrogen mole fraction along the XY plane of anode channels, for a value of Z corresponding to half the channels depth, for the two selected voltage values. It is evident from Fig. 6, that hydrogen mole fraction changes along the channels, mainly for low voltage values. More hydrogen is consumed for lower voltage values (corresponding to increasing current densities). This result is clearly, in discordance with the assumption made in the
development of the 1D model i.e. that the composition inside the channels is uniform. This assumption is acceptable for high voltage values because the mole fraction, effectively, does not change significantly (variations between 0.74 and 0.85), but leads to higher errors for low voltage values. It should be noted that, at the contour plan, hydrogen enters in the cell at the left side, so, higher values are expected and obtained at the cell inlet (left) and lower at the outlet (right). The contours of the oxygen mole fraction are represented in Fig. 7, once again along the XY plane of anode channels, for a value of Z corresponding to half the channels depth and for
Fig. 8 e Contours of water mole fraction in an XY plane at half of the cathode channels for depth a) 0.8 V, b) 0.6 V. Operating conditions: operating temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios of za [ 1 and with zc [ 2 for 1 A/cm2.
12494
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Fig. 9 e Contours of membrane water content in an XY plane at half of the membrane thickness for a) 0.8 V, b) 0.6 V. Operating conditions: operating temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and with zc [ 2 for 1 A/cm2.
The analysis of the membrane water content is essential in the evaluation of the water management in the cell. Fig. 9 shows the water membrane content (in terms of parameter l, the number of water molecules per ionic group) along the membrane centered XY plane, for the two selected voltage values. It is evident from the contours represented in Fig. 9 that the membrane water content distribution is not uniform for both voltage values being this non-homogeneity more significant for the higher voltage (lower current density) showing l values between 6 and 15 while for the lower voltage the l values fall between 10 and 15. The developed one-dimensional model predicts for the same operating conditions, average values for l around 13 and 14 for the voltages of 0.8 V and 0.6 V. These values are close to the higher values of l obtained in the 3D simulations. The possibility of obtaining the water content
1.2
1.0
Voltage (V)
the two selected voltage values. The plotted contours show variations along the cathode channel. Air enters in counter flux with hydrogen being at the cathode side the cell inlet at the right side of the contour plan and the cell outlet at the left side. For the lower voltage value corresponding to a higher current density and consequently higher kinetic rates, a decrease in the oxygen mole fraction along the channels is evident while for the higher value of the voltage, the contour distribution is more uniform (variations between 0.168 and 0.174) as for the case of the contour of hydrogen mole fractions. Fig. 8 shows the contours of the water mole fraction along the selected XY plane of the cathode channels for the same voltage values. This representation is relevant because at this side of the cell water is formed by the oxygen reduction and must be efficiently removed from the cell. For the voltage values selected, the water mole fraction does not follow a specific trend perhaps, due to the complexity of the water transport phenomena involved. For the higher value of voltage (situation a) the water mole fraction decreases along the cathode channel accompanying the slight increase of the oxygen mole fraction). For the voltage of 0.6 V (higher current density), the water mole fraction first increases and then decreases near the channel outlet. For this voltage, the back-diffusion mechanism could become dominant, in the final zone of the channel, due to a progressive increase of the water concentration. In fact, some authors [25] measured the water transport coefficient along the channel and reported the existence of positive and negative values along the same channels. It is worth mentioning that the negative water transport values indicate that the water flux occurs from the cathode to the anode. It should also be noted, when comparing the two situations in Fig. 8 that for the higher current density (lower voltage) the values of the water mole fraction are in general higher due to higher water production rates and to the occurrence of an increased water electro-osmotic drag.
H2 Hum, Air Hum Exp
H2 Hum, Air Hum Mod
H2 Dry, Air Hum Exp
H2 Dry, Air Hum Mod
H2 Hum, Air Dry Mod
0.8
0.6
0.4
0.2
0.0 0.0
0.1
0.2
0.3
0.4
0.5
0.6
Current Density (A/cm2 )
Fig. 10 e Comparison between experimental and 3D numerical predictions of the polarization curves for 333 K, when humidified both gases with an RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios of za [ 1 and zc [ 2 for 1 A/cm2.
12495
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Fig. 11 e Contours of membrane water content in an XY plane at half of the membrane thickness for 0.6 V, for three humidification situations: a) both gases fully humidified, b) air fully humidified and dry hydrogen and c) dry air and hydrogen fully humidified. Operating conditions: operating temperature of 333 K, H2/Air flow rates calculated with stoichiometry flow ratios of za [ 1 and zc [ 2 for 1 A/cm2.
visualization all over the membrane is an asset of the 3D numerical simulations while the 1D model gives an average value for l the plane XY, for a given value of Z. The ideal distribution of water content should be uniform all over membrane, avoiding dry-out or flooding zones which affect negatively the cell behavior. For these two conditions, a dryer zone tends to appear at the cathode outlet.
4.2.
Reactants relative humidity effect
The effect of reactants relative humidity on cell performance was studied at an operating temperature of 353 K, for three humidification conditions: both gases fully humidified, air fully humidified and dry hydrogen and dry air and hydrogen fully humidified. The 3D numerical predictions of the fuel cell performance are compared to the experimental curves in Fig. 10. Experimentally, the experiment with dry air could not
be performed due to difficulties in the stabilization of the fuel cell. Simulations reveal a better fuel cell performance for the situation when both gases enter in the cell fully humidified, mainly for high current densities. Curiously, for lower current density values, the predictions for fuel cell performance are slightly better when one of the reactants enters dry in the cell. These results are somewhat unexpected, because, for low current density values, the water production rate is lower, so, the fuel cell performance should be enhanced by a better humidification. However, it seems that the humidification of just one stream is sufficient to maintain the membrane hydrated. Probably, there is an initial membrane humidification provided by back-diffusion mechanism in the case of humidified air and by electro-osmotic drag in the case of humidified hydrogen. For high current density values, the water electro-osmotic drag is very intense, leading eventually
1.2 298 K Exp 313 K Exp 333 K Exp
1.0
298 K Model 313 K Model 333 K Model
Voltage (V)
0.8
0.6
0.4
0.2
0.0 0.0
Fig. 12 e Membrane water content along the membrane predicted by the 1D model for 0.6 V and for the three studied humidification conditions. Operating conditions: operating temperature of 333 K, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and zc [ 2 for 1 A/cm2.
0.1
0.2
0.3
0.4
0.5
0.6
Current Density (A/cm2)
Fig. 13 e Comparison between experimental data and 3D numerical simulations of polarization curves for the three temperatures studied a) 298 K b) 313 K and c) 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and and zc [ 2 for 1 A/cm2.
12496
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Fig. 14 e Protonic conductivity along the membrane predicted by the 1D model for 0.6 V and for the three studied temperatures. Operating conditions: both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and zc [ 2 for 1 A/cm2.
to membrane dehydration, and, consequently, to poorer performances when only one humidified stream is used. This decrease in cell performance is more noticeable when dry air is used, probably because the produced water is dragged by the dry air stream, leading to a more pronounced dehydration in membrane. The 3D simulations are close to experimental data when both gases are fully humidified. For the dry air case, the simulations deviate from experimental data in particular for lower current density values. The contours of water membrane content (in terms of parameter l) along the membrane centered XY plane, for the three studied conditions and 0.6 V, are presented in Fig. 11. The distribution of the membrane water content is not uniform all over the membrane and is substantially higher when both gases enter in the cell humidified. The lowest values are obtained when dry air is used. These results are in accordance with the 3D simulated polarization curves for this operating condition (0.6 V), corresponding to better performances when both gases are humidified and lower for the
case of dry air. For this latter case, the water content values are extremely low revealing dehydration of the membrane. This fact is in accordance with the experimental difficulties in reproducing this operating condition. The average values of the water content along the membrane for the same value of voltage predicted by the 1D model are displayed in Fig. 12 for the different humidification conditions. As can be seen the trend is similar: higher average water content values are obtained when both gases are humidified and lower values for the case of dry air. When comparing the average water content at the center of the membrane (z ¼ 0.013 cm) with the information given by the water contours in Fig. 11, it is evident that the predictions provided by the simpler model 1D are higher that the 3D simulations for the cases of one humidified stream. This discrepancy is somewhat expected, since the 1D model predicts higher performances (recall Fig. 5) corresponding to higher values of the current densities for the same voltage value and consequently enhanced predicted average water contents.
4.3.
Operating temperature effect
The operating temperature effect was studied for 298 K, 313 K and 333 K. The comparison between 3D numerical simulations and experimental data are presented in Fig. 13. As can be seen from Fig. 13, the 3D predicted curves are close to the experimental data. The curves are close to each other showing a slight increase of fuel cell performance with increasing temperature for the entire range of current densities covered. For this type of membranes (very thin), the improvement of performance with increasing temperature is not perceptible in experimental data. A well-known impact of an increase in the cell temperature is the enhancement of the membrane conductivity and, consequently, the decrease of the membrane resistance to the transport of hydrogen protons. Following this idea, the predictions of the average values of the proton conductivities along the membrane and the contours of protonic conductivity for the membrane
Fig. 15 e Contours of protonic conductivity in an XY plane at half of the membrane thickness given by the 3D numerical simulations, for 0.6 V and for the three studied temperatures a) 298 K b) 313 K and c) 333 K. Operating conditions: both gases RH of 95%, H2/Air flow rates calculated with za [ 1 and with zc [ 2 for 1 A/cm2.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
middle plane are presented in Figs. 14 and 15 for the three studied conditions. The simulation results from the two Figures show, effectively, an increase in the proton conductivity with increasing temperatures. The curves in Fig. 14 also reveal an increase of the proton conductivity along the membrane accompanying the increase in water content. The contours in Fig. 15 put in evidence the non-homogeneity of proton conductivity in the studied plane. The ideal situation should correspond to high values of protonic conductivity with a uniform distribution all over the membrane in order to facilitate the protons flux. The comparison between the average value at the center of the membrane (z ¼ 0.013 cm) taken from Fig. 14 with the information given by the contours in Fig. 15, puts in evidence the accordance between 1D and 3D predictions. For 298 K, the contour exhibits proton conductivity values between 5.6 U1 m1 and 5.9 U1 m1, while the 1D model predicts an average value for the proton conductivity in the membrane middle plane of 6.3 U1 m1. For 313 K, the contour gives values between 6.2 U1 m1 and 7.3 U1 m1, while the predicted average value in the centered plane is 7.3 U1 m1. Finally, for 333 K, the 3D predicted values are between 7 U1 m1 and 10.8 U1 m1, while the 1D model predicts a proton conductivity of 8.6 U1 m1. One can conclude that the developed one-dimensional model gives adequate predictions for the proton conductivity expending a much lower CPU computing time.
5.
Conclusions
In the present work a comparison was performed between 1D and 3D numerical simulations in PEM fuel cells. The onedimensional model used was previously developed and validated [1]. The 3D numerical simulations were obtained using the commercial code FLUENT-PEMFC module. Particular attention was devoted to the impact of considering two-phase flow effects in the 3D simulation results. The 3D simulations considering one and two phases for water flow and the 1D model simulations were compared with experimental data. The results showed that 1D and 3D model simulations considering just one phase are similar, with a slightly better accordance for the 1D model exhibiting a substantially lower CPU time. However both numerical results over predict the fuel cell performance. Actually, the 3D simulations reproduce very well the experimental data when two-phase flow effects are considered. The formation of liquid water (considered in this case) tends to block the pores, leading to a decrease in fuel cell performance. Three-dimensional numerical simulations of oxygen, hydrogen and water mole fractions along the channels were analyzed. The distribution of mole fractions for the three species is not uniform along the channels in particular for higher current densities (corresponding to higher electrochemical reaction rates). The water mole fraction does not follow a defined trend, showing for certain operating conditions an increase and afterwards a decrease near the channel outlet. This finding is in accordance with published results
12497
showing the change of the water dominant mechanism from the cell inlet to the cell outlet [24]. The effect of the relative humidity of gases and operation temperature on fuel cell performance was also studied both through the comparison of the polarization curves for the 1D and 3D simulations and experimental data and through the analysis of relevant physical parameters such as the water membrane content and the proton conductivity. The results showed that, as expected, the fuel cell performance is improved when both gas streams are humidified. The 1D and 3D simulations corroborate this finding, revealing higher water content values for this condition. Fuel cell performance slightly increases with increasing temperature, being this trend perceptible in the simulated polarization curves and confirmed by the protonic conductivity predictions. The protonic conductivity values are higher for higher temperature values. The water content predictions by the simpler 1D model are in agreement with the 3D numerical simulations, except for the cases when one stream enters dry in the cell. In these cases, the simpler model predicts higher water content values. Concerning the predictions of the proton conductivities, the 1D and 3D results are in agreement. The presented results showed that the developed 1D model (coupling heat and mass transfer effects) can be used to predict optimal operating conditions in PEMFCs and the general trends of the impact of several important physical parameters in fuel cell performance such as membrane water content or the membrane protonic conductivities. The use of the 3D numerical simulations is indicated if more detailed predictions are needed. An asset of this type of simulations is the prediction of the spatial distribution and visualization of various relevant parameters. However simulation times are still too high. A polarization curve with the 1D model is obtained with a CPU time around 5 min, while the 3D computing time is around 24 h. Besides this disadvantage, before the generation of the 3D numerical simulations, it is necessary to create the mesh for the whole cell, (in this case, GAMBIT was the mesh generator used), a very difficult task due to the discrepancy between the channels dimensions and the rest of the cell geometry. An important conclusion of this work is the demonstration that a simpler analytical model using low CPU has potential to be used in real-time PEMFC simulations. Empirical models can always be used for quick predictions, however, they are «validated» against polarization curves using a significant number of fitting parameters. The smart use of «old fashioned» analytical or semi analytical models using a rationally chosen set of parameters to obtain useful reduced models is put in evidence in the present work.
Acknowledgements The partial support of “Fundac¸a˜o para a Cieˆncia e Tecnologia Portugal” through scholarship SFRH/BD/28166/2006 is gratefully acknowledged. POCI (FEDER) also supported this work via CEFT.
12498
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Nomenclature
CACL H2 ref
CH2 CCCL O2 ref
CO2 cr E0 EM ES ref ja ref jc K r rw Ra Rc RM
RS s S vl Yj, cen Yj, sur
hydrogen concentration at anode catalytic layer, mol/m3 hydrogen reference concentration, mol/m3 oxygen concentration at anode catalytic layer, mol/ m3 oxygen reference concentration, mol/m3 condensation rate constant, s1 open circuit voltage, V phase potential of the membrane, V phase potential of the solid, V anode volumetric transfer current, A/cm3 cathode volumetric transfer current, A/cm3 absolute permeability, m/s specific reacting surface area, m1 condensation Rate, m/s anode volumetric transfer current, A/m3 cathode volumetric transfer current, A/m3 membrane volumetric transfer current, A/m3 solid volumetric transfer current, A/m3 water saturation source term liquid velocity, m/s mass fraction of species j at the cell center mass fraction of species j at the reactant surface
Greek letters anode transfer coefficient aa cathode transfer coefficient ac d average distance between reaction surfaces and cell center, m anode concentration dependence ga cathode concentration dependence gc anode overpotential, V ha cathode overpotential, V hc l membrane water content q contact angle, s surface tension, N/m2 sM membrane electric conductivity, S/m solid electric conductivity, S/m sS anode stoichiometric ratio za cathode stoichiometric ratio zc
references
[1] Falca˜o DS, Oliveira VB, Rangel CM, Pinho C, Pinto AMFR. Water transport through a PEM fuel cell: a one-dimensional model with heat transfer effects. Chemical Engineering Science 2009;64(9):2216e25. [2] Falcao DS, Rangel CM, Pinho C, Pinto AMFR. Water transport through a proton-exchange membrane (PEM) fuel cell operating near ambient conditions: experimental and modeling studies. Energy & Fuels 2009;23(1):397e402.
[3] Eikerling M, Kornyshev AA. Modelling the performance of the cathode catalyst layer of polymer electrolyte fuel cells. Journal of Electrochemical Society 1998;453(1e2):89e106. [4] Eikerling M, Kornyshev AA, Kucernak AR. Water in polymer electrolyte fuel cells: friend or foe? Physics Today 2006;59(10): 38e44. [5] Baschuk JJ, Li X. Modelling of polymer electrolyte membrane fuel cells with variable degrees of water flooding. Journal of Power Sources 2000;86(1e2):181e96. [6] BiyIkoglu A. Review of proton exchange membrane fuel cell models. International Journal of Hydrogen Energy 2005; 30(11):1181e212. [7] Ji M, Wei Z. A review of water management in polymer electrolyte membrane fuel cells. Energies 2009;2(4):1057e106. [8] Springer TE, Zawodzinski TA, Gottesfeld S. Polymer electrolyte fuel cell model. Journal of Electrochemical Society 1991;138(8):2334e42. [9] Bernardi DM, Verbrugge MW. A mathematical model of the solid-polymer-electrolyte fuel cell. Journal of Electrochemical Society 1992;139(9):2477e91. [10] Kulikovsky AA. Semi-analytical 1D þ 1D model of a polymer electrolyte fuel cell. Electrochemistry Communications 2004; 6(10):969e77. [11] Scott K, Argyropoulos P. A current distribution model of a porous fuel cell electrode. Journal of Electroanalytical Chemistry 2004;567(1):103e9. [12] De Francesco M, Arato E, Costa P. Transport phenomena in membranes for PEMFC applications: an analytical approach to the calculation of membrane resistance. Journal of Power Sources 2004;132(1e2):127e34. [13] Wang ZH, Wang CY, Chen KS. Two-phase flow and transport in the air cathode of proton exchange membrane fuel cells. Journal of Power Sources 2001;94(1):40e50. [14] Um S, Wang CY, Chen KS. Computational fluid dynamics modeling of proton exchange membrane fuel cells. Journal of Electrochemical Society 2000;147(12):4485e93. [15] Ju H, Meng H, Wang C-Y. A single-phase, non-isothermal model for PEM fuel cells. International Journal of Heat and Mass Transfer 2005;48(7):1303e15. [16] Berning T, Lu DM, Djilali N. Three-dimensional computational analysis of transport phenomena in a PEM fuel cell. Journal of Power Sources 2002;106(1e2):284e94. [17] Le AD, Zhou B. A general model of proton exchange membrane fuel cell. Journal of Power Sources 2008;182(1):197e222. [18] Le AD, Zhou B. A generalized numerical model for liquid water in a proton exchange membrane fuel cell with interdigitated design. Journal of Power Sources 2009;193(2):665e83. [19] Meng H. A three-dimensional PEM fuel cell model with consistent treatment of water transport in MEA. Journal of Power Sources 2006;162(1):426e35. [20] Wu H, Li X, Berg P. On the modeling of water transport in polymer electrolyte membrane fuel cells. Electrochimica Acta 2009;54(27):6913e27. [21] Dawes JE, Hanspal NS, Family OA, Turan A. Threedimensional CFD modelling of PEM fuel cells: an investigation into the effects of water flooding. Chemical Engineering Science 2009;64(12):2781e94. [22] Jang J-H, Yan W-M, Li H-Y, Tsai W-C. Three-dimensional numerical study on cell performance and transport phenomena of PEM fuel cells with conventional flow fields. International Journal of Hydrogen Energy 2008;33(1):156e64. [23] Fluent. Fluent 6.3-fuel cells module manual; 2006. [24] Sivertsen BR, Djilali N. CFD-based modelling of proton exchange membrane fuel cells. Journal of Power Sources 2005;141(1):65e78. [25] Liu F, Lu G, Wang C-Y. Water transport coefficient distribution through the membrane in a polymer electrolyte fuel cell. Journal of Membrane Science 2007;287(1):126e31.
Available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/he
1D and 3D numerical simulations in PEM fuel cells D.S. Falca˜o a, P.J. Gomes b, V.B. Oliveira a, C. Pinho a, A.M.F.R. Pinto a,* a Centro de Estudos de Feno´menos de Transporte, Departamento de Engenharia Quı´mica, Faculdade de Engenharia da Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal b Laboratory of Separation and Reaction Engineering, Departamento de Engenharia Quı´mica, Faculdade de Engenharia da Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal
article info
abstract
Article history:
The potential of fuel cells for clean and efficient energy conversion is generally recognized.
Received 30 April 2011
The proton-exchange membrane (PEM) fuel cells are one of the most promising types of fuel
Received in revised form
cells. Models play an important role in fuel cell development since they enable the under-
21 June 2011
standing of the influence of different parameters on the cell performance allowing a systematic
Accepted 26 June 2011
simulation, design and optimization of fuel cells systems. In the present work, one-
Available online 31 July 2011
dimensional and three-dimensional numerical simulations were performed and compared with experimental data obtained in a PEM fuel cell. The 1D model, coupling heat and mass
Keywords:
transfer effects, was previously developed and validated by the same authors [1,2]. The 3D
PEM fuel cell
numerical simulations were obtained using the commercial code FLUENT e PEMFC module.
Water transport
The results show that 1D and 3D model simulations considering just one phase for the
Modeling studies
water flow are similar, with a slightly better accordance for the 1D model exhibiting
CFD
a substantially lower CPU time. However both numerical results over predict the fuel cell performance while the 3D simulations reproduce very well the experimental data. The effect of the relative humidity of gases and operation temperature on fuel cell performance was also studied both through the comparison of the polarization curves for the 1D and 3D simulations and experimental data and through the analysis of relevant physical parameters such as the water membrane content and the proton conductivity. A polarization curve with the 1D model is obtained with a CPU time around 5 min, while the 3D computing time is around 24 h. The results show that the 1D model can be used to predict optimal operating conditions in PEMFCs and the general trends of the impact on fuel cell performance of several important physical parameters (such as those related to the water management). The use of the 3D numerical simulations is indicated if more detailed predictions are needed namely the spatial distribution and visualization of various relevant parameters. An important conclusion of this work is the demonstration that a simpler model using low CPU has potential to be used in real-time PEMFC simulations. Copyright ª 2011, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
1.
Introduction
Fuel cells are extensively being studied today because of their potential as an alternate energy source for a wide range
of applications. Fuel cells have unique technological attributes: efficiency, absence of moving parts and very low emissions. In particular, the proton-exchange membrane (PEM) fuel cells are today in the focus of interest as one of the
* Corresponding author. Tel.: þ351225081675; fax: þ351225081449. E-mail address: [email protected] (A.M.F.R. Pinto). 0360-3199/$ e see front matter Copyright ª 2011, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2011.06.133
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
12487
Fig. 1 e Schematic representation of the experimental set-up.
most promising technology in development in power generation. The necessity of formulation of adequate mathematical models is based on the need to reach a deeper understanding of internal processes (reactions, mass transport, heat transport) occurring in a fuel cell. These processes cannot be observed directly in experiments and mathematical model could be applied to obtain hints for the development of optimal control and operating strategies. One of the most important operational issues of PEMFC is the water management in the cell [3,4]. Water content of the membrane is determined by the balance between water production and three water transport processes: electro-osmotic drag of water (EOD), associated with proton migration through the membrane; back diffusion from the cathode to the anode and diffusion of water to/from the oxidant/fuel gas streams. Understanding the water transport in the PEM is a key issue to study. In
order to achieve optimal fuel cell performance, an adequate water balance should be ensured because membrane should remain hydrated for sufficient proton conductivity, while cathode flooding and anode dehydration should be avoided [5e7]. Different models were developed in the last decade to describe several water transport mechanisms through the membrane such as Springer et al. [8] using a diffusion model, Bernardi and Verbrugge [9] considering a hydraulic permeation model and Kulikovsky [10] developing a semi analytical model 1D þ 1D taking into account oxygen and water transport across the cell and deriving an expression for the limiting current density. Scott and Argyropoulos [11] developed a onedimensional model of the current and potential distribution in a porous electrode. De Francesco et al. [12] used an analytical approach to calculate the membrane resistance. More recently, there was a trend to apply CFD methods to fuel cells modeling [13,14]. Due to cell configuration, transport
12488
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
phenomena are three-dimensional, so, some works taking in account three dimensions arose [15e22]. In this study, 1D a 3D numerical simulations are presented and compared, in order to verify advantages and disadvantages of each type of modeling. The one-dimensional model was already developed in a previous work [1]. This model considers the effects of coupled heat and mass transfer, along with the electrochemical reactions occurring in PEMFC. The model can be used to predict the hydrogen, oxygen and water concentration profiles in the anode, cathode and membrane, as well as the temperature profile across the cell. The model was validated with published experimental data. This model is also used and validated in another reference [2]. The 3D numerical simulations were obtained in a commercial software FLUENT with the PEMFC module. In the present work, the influence of two-phase flow effects, as well as the influence of reactants humidity and operating temperatures is studied recurring to the models and to experimental data. Besides polarization curves, the presented models simulate several important parameters with impact on fuel cell performance, such as membrane water content and protonic conductivity.
2.
Experimental set-up
2.1.
Apparatus
A schematic drawing of the experimental apparatus used in this work is shown in Fig. 1. Pure hydrogen (humidified or dry) as fuel and air (humidified or dry) as oxidant are used. The pressure of the gases is controlled by pressure regulators (Air e Norgreen 11400, H2 e Europneumaq mod. 44-2262-241) and flow rates controlled by flow meters (KDG e Mobrey). The reactants humidity and temperatures are monitored by adequate humidity and temperature probes (Air e Testo, H2 e Vaisala). The humidification of air and hydrogen gases is conducted in Erlenmeyer flasks by a simple bubbling process. To control the humidification temperature each Erlenmeyer flask is thermally isolated and surrounded with an electrical resistance (50 W/m) activated by an Osaka OK 31 digital temperature controller. The same procedure is applied along the connecting pipes from the humidification point up to the entrance of the fuel cells to guarantee the temperature stabilization of each reacting gas flow, as well as to control the operating temperature of the fuel cell. For the measurement and control of the cell electrical output, an electric load reference LD300 300 W DC Electronic Load from TTI is used. This load was connected to a data acquisition system composed by Measurement Computing boards installed in a desktop computer. The used data acquisition software was DASYLab.
2.2.
Fuel cell design
In the present work, all the components of the PEMFC were “in house” designed, except the MEA. A Nafion 111 from Dupont
Fig. 2 e Representation of the serpentine flow channels.
membrane with 25 cm2 active surface area is used. The channels configuration used for the anode and cathode flow channels is represented in Fig. 2. The channels depth is 1 mm for hydrogen side and 1.5 mm for air side. Table 1 summarizes the cell design specifications.
2.3.
Experimental conditions
The studied operating conditions were the reactants humidity and the operating temperature. All the experiments were conducted at atmospheric pressure. Conditions are presented in Tables 2 and 3.
3.
Numerical simulations
3.1.
One-dimensional modelling
The 1D model presented here for comparison with the 3D numerical simulations was already developed and validated. The detailed description of the model can be found in references [1] and [2]. The main assumptions are recalled here: - Steady-state and one-dimensional heat and mass transport; - Ideal gas mixture assumed;
Table 1 e Fuel cell design parameters. Gas difusion layer Catalytic layer
Channels configuration Membrane
Anode Cathode Anode
Thickness: 0.0195 cm Thickness: 0.0195 cm Catalyst Loading: 0.5 mgPt/cm2 Thickness: 0.006 cm Cathode Catalyst Loading: 0.5 mgPt/cm2 Thickness: 0.006 cm Anode Serpentine Cathode Serpentine Nafion 111, Thickness: 0.00254 cm
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
12489
Table 2 e Operating conditions used to analyse reactants humidity effect. Operating temperature (K) 333
H2 flow rate
Air flow rate
Calculated for Calculated for z ¼ 1 for 1 A/cm2 z ¼ 1 for 1 A/cm2
H2 RH (%)
Air RH (%)
95 0 95
95 95 0
- Heat and mass transport through the gas diffusion and catalyst layers assumed to be a diffusion-predominated process (negligible convection effects); - Effective Fick models for the mass transport in the diffusion layers and membrane are considered; - The thermal energy model is based on the differential thermal energy conservation equation (Fourier’s law); - The thermal conductivity for all the materials is assumed to be constant; - The anodic and cathodic overpotential is constant through the catalyst layers; - Heat generation or consumption is considered in the catalyst layers; - The heat generated by the current flow (Joule effects) through each cell layer is ignored; - Water transport through the membrane assumed to be a combined effect of diffusion and electro-osmotic drag; - Membrane proton conductivity is a function of l, the number of water molecules per ionic group; - Local equilibrium at interfaces is represented by partition functions; - Reactions in the catalyst layers are considered as homogeneous; - Kinetics of the cathode is described by a Tafel expression; - Kinetics of the anode is described by a Tafel expression; - Anode and cathode flow channels are treated as a continuous stirred tank reactor (CSTR), so, the composition and temperature inside the channels are uniform; - Anode and cathode streams act as heat transfer fluids removing heat from the cell at the exit temperatures.
3.2.
Fig. 3 e Representation of the flow channels Mesh.
3.2.1.
Geometry and mesh generation
The geometry of the fuel cell includes the MEA and the channel plates. Terminal plates and collector plates are not considered, due to simplification (shorter simulation times). The used mesh is not uniform, with elements type Hex/Wedge Cooper and around 1,600,000 elements. Mesh in the channels is more refined, because the channel dimensions are very small. Due to this dimensional difference between channels and the other fuel cell components, the mesh generation is complicated. A representation of the mesh in the channels is shown in Fig. 3. Concerning mesh division, in terms of depth, the anode channel is divided in 2 (depth of 0.7 mm) and the cathode channel is divided in 5 (depth of 1.5 mm). In the MEA, the mesh is also built with high refinement, because the layers are very thin. The scheme used in mesh division is presented in Fig. 4. The mesh resolution can affect the solution precision and therefore a grid independency test was done. In order to do that, a more refined mesh was created, with 2,400,000
Three-dimensional model formulation
The 3D numerical simulations were obtained with commercial code FLUENT, in particular using the module PEFC. Details of the implementation and resolution of equations can be found in [23]. Geometry and mesh of the fuel cell had to be generated.
Table 3 e Operating conditions used to analyse operating temperature effect. Operating temperature (K) 298 313 333
H2 flow rate
Air flow rate
Calculated for Calculated for z ¼ 1 for 1 A/cm2 z ¼ 1 for 1 A/cm2
H2 RH (%)
Air RH (%)
95
95
Fig. 4 e General indication of the grid of the MEA mesh.
12490
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
elements and simulations with this new mesh were done. The results show variations around 0.03% relatively to those obtained with regular mesh. So, it could be concluded that results were independent of grid refinement.
3.2.2.
Boundary conditions
Boundary conditions need to be defined for PEM fuel cell simulations, based on problem specification. In this model, the following boundary conditions were defined as: - Inlet (mass flow inlet type) and outlet (pressure outlet type) zones for the anode and cathode gas channels - Surfaces (wall type) representing anode and cathode channels The following continuum zones were also required: -
Anode and cathode flow channels (fluid zone) Anode and cathode current collectors (solid zone) Anode and cathode diffusion layers (fluid zone) Anode and cathode catalyst layers (fluid zone) Membrane (fluid zone)
boundary, there is a zero flux boundary condition for the membrane phase potential on all outside boundaries. For the solid phase potential, there are external boundaries on the anode and the cathode side that are in contact with the external electric circuit. Through these boundaries passes the electrical current generated in the fuel cell. On all other external boundaries there is a zero flux boundary condition solid phase potential. Anode potential was set to zero, so the positive value prescribed on the cathode side is the cell voltage. The transfer currents, or the source terms in Equations (1) and (2), are non-zero only inside the catalyst layers and are computed as: - For the solid phase, RS ¼ Ra (<0) on the anode side and RS ¼ þRc (>0) on the cathode side. - For the membrane phase, RM ¼ þRa (>0) on the anode side and RM ¼ Rc (<0) on the cathode side. Source terms on Equations (1) and (2) are also called exchange current density and have the generic Tafel formulation [23]:
Boundary conditions definition was done in Gambit (the mesh generator software used). After geometry and mesh generation and boundary and continuum zones established, the mesh was found to be suitable to be exported to FLUENT.
Ra ¼
jref a
Rc ¼
jref c
CACL H2
!ga
ref
CH2
CCCL O2
!gc
aa Fh =RT e a
ref
CO2
eac Fhc =RT
(3)
(4) ref
3.2.3.
Model equations
In the PEM fuel cells module of Code FLUENT, two electrical potential fields are solved: one in the membrane and catalytic layers and the other in catalytic layers, diffusion layers and current collectors. Surface reactions occurring on the porous catalyst region are solved and the reaction diffusion balance is applied to compute the rates. Based on a cell voltage prescribed, the current density value is computed.
3.2.3.1. Electrochemistry model. The most important electrochemistry aspects are hydrogen oxidation and oxygen reduction reaction rates. In this model, these electrochemical processes are treated as heterogeneous reactions taking place on the catalyst surfaces. The driving force behind these reactions is the surface overpotential: the difference between the phase potential of the solid (ES) and the phase potential of the electrolyte/membrane (EM). Therefore, two potential equations are solved, one for electrons transport and the other for hydrogen protons transport: VðsS VES Þ þ RS ¼ 0
(1)
VðsM VEM Þ þ RM ¼ 0
(2)
where R is the volumetric transfer current [A/m3], s is the electrical conductivity [S/m] and the indices S and M refers to, respectively, solid phase and membrane phase. For the resolution of these equations resolution, several boundary conditions need to be defined. There are two types of external boundaries: those through which an electrical current passes and those through which no current passes. As no protonic current leaves the fuel cell through any external
where CACL H2 is the hydrogen concentration in catalytic layer, CH2 is the reference hydrogen concentration, CCCL O2 is the oxygen ref concentration in catalytic layer, CO2 is the reference oxygen ref ref concentration, ja and jc are, respectively, the anode and cathode volumetric reference exchange current density, aa and ac are, respectively, anode and cathode transfer coefficient, ga and gc are, respectively, the anode and cathode concentration dependence, and ha and hc are, respectively, the anode and cathode overpotentials. These overpotentials are the difference between the solid and membrane potentials. The gain in electric potential from crossing from the anode to the cathode side can then be taken into account by subtracting the open circuit voltage, E0, on the cathode side: ha ¼ ES EM
(5)
hc ¼ ES EM E0
(6)
Adequate combination of Equations (1)e(6) allows the resolution of the two potential fields.
3.2.3.2. Current and mass conservation. The volumetric source terms (S) for the species equations [kg/m3s] and energy equation [W/m3] are given in Equations (7)e(9) [23]: MH SH2 ¼ 2 Ra 2F
(7)
MO SO2 ¼ 2 Rc 2F
(8)
MH2 O Rc 2F
(9)
SH2 O ¼
where M is the molecular weight [kg/mol].
12491
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Additional volumetric sources (Senergy) to the energy equation implemented in the FLUENT - PEMFC model include ohmic heating (Rohm), heat of formation of water (Hreac), electric work and latent heat of water (Hphase). Senergy ¼ I2 Rohm þ Hreac þ hRa;c þ Hphase
(10)
The electrochemical reactions are treated as surface reactions in the two catalyst layers, and it is assumed that the diffusive flux of any reacting species is balanced by its rate of production: rDj Mj yj;sur yj;cen r ¼ Ra;c d nF
(11)
where r is the specific reacting surface area of the catalyst layer [m1], yj, sur is the mass fraction of species j at the reactant surface, yj, cen is the mass fraction of species j at the cell center and d is the average distance between the reaction surfaces and the cell center [m]. The left hand side of Equation (11) represents the diffusive flux at the reacting surface and the right hand side represents the rate of mass generation. The average distance from the cell center to the reacting surface is estimated as d ¼ 1/r. Equation (11) is used to obtain the surface values of H2 and O2 concentrations, applying a Newtonian solution procedure. These surface values are then used to compute the rates in Equations (3) and (4).
3.2.3.3. Liquid water formation, transport, and its effects. Since PEM fuel cells operate under relatively low temperature (<100 C), the water vapor may condense to liquid water, especially at high current densities. While the existence of the liquid water keeps the membrane hydrated, it also blocks the gas diffusion passage, reduces the diffusion rate and the effective reacting surface area and hence the cell performance. To model the formation and transport of liquid water, FLUENT uses a saturation model. In this approach, the liquid water formation and transport is governed by the following conservation equation for the volume fraction of liquid water, s, or the water saturation: vðerl sÞ þ V$ðrl vl sÞ ¼ rw vt
(12)
where the subscript, l, represents liquid water, vl represents liquid velocity, which is equivalent to the gas velocity inside the gas channel and rw is the condensation rate, which is defined as: PH O Psat MH2 O ; ½srl rw ¼ cr max ð1 sÞ 2 RT
scosqc pc ¼ 0;5 1;417ð1 sÞ 2;12ð1 sÞ2 þ1;263ð1 sÞ3 ; for qc < 90+ K e (15) scosqc pc ¼ 0;5 1; 417s 2; 12s2 þ 1; 263s3 ; for qc > 90+ K e
(16)
where qc is the contact angle and s is the surface tension [N/m2]. Equation (12) is used to model various physical processes such as condensation, vaporization, capillary diffusion, and surface tension. The clogging of the porous media and the flooding of the reaction surface are modeled by multiplying the porosity and the active surface area by (1-s).
3.2.4.
Solution procedure
Solver default definitions from FLUENT are not appropriate to reach solution convergence. Convergence is here considered when the difference between two consecutive iterations is less than 0.001 for the residuals of continuity, velocity and concentrations equations and for UDFs and less than 0.0001 for the residuals of the energy equation. Multigrid cycle was changed for F-Cycle for all equations. Stabilization method was changed for BCGSTAB for species concentrations, water saturation, electric and protonic potential. Terminal restrictions were also changed for 0.001 for species concentrations, water saturation and for 0.0001 for electric and protonic protential. Maximum number of cycles was increased for 50. Under-relaxation factors were changed to the following values: 0.45 for pressure, 1 for density, 0.9 for body forces, 0.3 for momentum, 0.9 for species concentration, 0.2 for energy, 1 for electric potential, 0.95 for protonic potential and 0.5 for water saturation and water content. The discretization method was changed to second order (leading to more precise solutions) for all equations.
3.2.5.
Model parameters
The model parameters are already defined in the PEM Fuel Cells module but could be set by the user. All default parameters were not changed, except those related with the kinetics, presented in Table 4. Reference concentrations were set to the concentrations used in each experiment.
(13)
where rw is added to the water vapor equation, as well as the pressure correction. This term is not applied inside the membrane. The condensation rate constant is defined as cr ¼ 100/s. Inside the highly-resistant porous zones, the use of the capillary diffusion term allows us to replace the convective term in Equation (12): vð˛rl sÞ Ks3 dpc Vs ¼ rw þ V$ rl ml ds vt
where K is the absolute permeability and pc the capilarity pressure. Depending on the wetting phase, the capillary pressure is computed as a function of s (the Leverett function) [23]:
(14)
Table 4 e Model kinetic parameters. Parameter E0 ref ja ref jc aa ac ga gc
Value
Reference
0:0025 T þ 0:2329 V 1 109 A/m3 2 104 exp [0.014189 (T-353)] A/m3 0.5 0.2 0.5 1.0
[18] [19] [19] [20] [21] [22] [22]
12492
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
1.2 Exp
Mod 3D Two-Phase
Mod 3D One-Phase
Mod 1D
Voltage (V)
1.0
0.8
0.6
0.4
0.2
0.0 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Current Density (A/cm2 )
Fig. 5 e Simulated and experimental polarization curves for 333 K, both gases RH of 95%, H2/Air flow rates with stoichiometry flow ratios za [ 1 and zc [ 2 for 1 A/cm2.
4.
Results and discussion
As already referred, the main objective of the present work was the comparison between 1D and 3D modeling, bearing in mind the demonstration that a simpler model using low CPU has potential to be used in real-time PEMFC simulations. Particular attention was devoted to the impact of considering two-phase flow effects in the 3D simulation results. The effect of the relative humidity of gases and operation temperature in fuel cell performance was also studied. Preliminary tests showed that the results using Tafel or Buttler-Volmer kinetic equations were very similar, consequently, the Tafel equation leading to lower computing times was used in all the simulations here presented.
4.1.
Two-phase flow effects
In Fig. 5, the 3D model simulations considering one and two phases for water flow and the 1D model simulations are
compared with experimental data, to evaluate the impact of two-phase flow effects on simulation results. The data plotted refer to an operation temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios of za ¼ 1 and zc ¼ 2 for 1 A/cm2. As can be seen from Fig. 5, the 3D simulation results considering two-phase flow effects compare very well with the experimental data. For this operating condition, the formation of liquid water is likely to occur, partially filling the pores of the electrodes structure and leading to a lower fuel cell performance. The consideration of only water vapor, in the 1D and 3D models (without two-phase flow effects) lead to higher predicted fuel cell performances, as expected, but reproduce however worse the experimental data. It is worth mentioning that these 1D and 3D model simulations considering just one phase are similar with a slightly better accordance for the 1D model exhibiting a substantially lower CPU time. In PEMFC literature, validation is commonly performed against polarization curves and using a significant number of fitting parameters. Furthermore, there are still many technical difficulties in obtaining detailed and reliable in situ experimental data. Consequently, verification of model performance based on polarization curves will persist but should however be complemented with an analysis of the physical results obtained with the model. Given the complexity of the transport phenomena involved and the uncertainty in determining some of the physico-chemical parameters, it is any case unrealistic to expect accurate quantitative agreement between experiments and computational model simulations [24]. So, it is essential to analyze the results in terms of the physical parameters predicted by the models. Three-dimensional numerical simulations enable the prediction of the distribution and visualization of various relevant parameters influencing the fuel cell behavior, namely, the mass or mole fractions of the different species in the channels and the membrane water content. These parameters change with voltage variations and therefore, the
Fig. 6 e Contours of hydrogen mole fraction in an XY plane at half of the anode channels depth for a) 0.8 V, b) 0.6 V. Operating conditions: operating temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios of za [ 1 and zc [ 2 for 1 A/cm2.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
12493
Fig. 7 e Contours of the oxygen mole fraction in an XY plane at half of the cathode channels depth for a) 0.8 V, b) 0.6 V. Operating conditions: operating temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and zc [ 2 for 1 A/cm2.
contours for two values of voltage (0.8 and 0.6 V) are presented in Fig. 6 for an operation condition revealing good agreement with experimental data (3D simulation considering two-phase flow effects). This figure shows the contours of hydrogen mole fraction along the XY plane of anode channels, for a value of Z corresponding to half the channels depth, for the two selected voltage values. It is evident from Fig. 6, that hydrogen mole fraction changes along the channels, mainly for low voltage values. More hydrogen is consumed for lower voltage values (corresponding to increasing current densities). This result is clearly, in discordance with the assumption made in the
development of the 1D model i.e. that the composition inside the channels is uniform. This assumption is acceptable for high voltage values because the mole fraction, effectively, does not change significantly (variations between 0.74 and 0.85), but leads to higher errors for low voltage values. It should be noted that, at the contour plan, hydrogen enters in the cell at the left side, so, higher values are expected and obtained at the cell inlet (left) and lower at the outlet (right). The contours of the oxygen mole fraction are represented in Fig. 7, once again along the XY plane of anode channels, for a value of Z corresponding to half the channels depth and for
Fig. 8 e Contours of water mole fraction in an XY plane at half of the cathode channels for depth a) 0.8 V, b) 0.6 V. Operating conditions: operating temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios of za [ 1 and with zc [ 2 for 1 A/cm2.
12494
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Fig. 9 e Contours of membrane water content in an XY plane at half of the membrane thickness for a) 0.8 V, b) 0.6 V. Operating conditions: operating temperature of 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and with zc [ 2 for 1 A/cm2.
The analysis of the membrane water content is essential in the evaluation of the water management in the cell. Fig. 9 shows the water membrane content (in terms of parameter l, the number of water molecules per ionic group) along the membrane centered XY plane, for the two selected voltage values. It is evident from the contours represented in Fig. 9 that the membrane water content distribution is not uniform for both voltage values being this non-homogeneity more significant for the higher voltage (lower current density) showing l values between 6 and 15 while for the lower voltage the l values fall between 10 and 15. The developed one-dimensional model predicts for the same operating conditions, average values for l around 13 and 14 for the voltages of 0.8 V and 0.6 V. These values are close to the higher values of l obtained in the 3D simulations. The possibility of obtaining the water content
1.2
1.0
Voltage (V)
the two selected voltage values. The plotted contours show variations along the cathode channel. Air enters in counter flux with hydrogen being at the cathode side the cell inlet at the right side of the contour plan and the cell outlet at the left side. For the lower voltage value corresponding to a higher current density and consequently higher kinetic rates, a decrease in the oxygen mole fraction along the channels is evident while for the higher value of the voltage, the contour distribution is more uniform (variations between 0.168 and 0.174) as for the case of the contour of hydrogen mole fractions. Fig. 8 shows the contours of the water mole fraction along the selected XY plane of the cathode channels for the same voltage values. This representation is relevant because at this side of the cell water is formed by the oxygen reduction and must be efficiently removed from the cell. For the voltage values selected, the water mole fraction does not follow a specific trend perhaps, due to the complexity of the water transport phenomena involved. For the higher value of voltage (situation a) the water mole fraction decreases along the cathode channel accompanying the slight increase of the oxygen mole fraction). For the voltage of 0.6 V (higher current density), the water mole fraction first increases and then decreases near the channel outlet. For this voltage, the back-diffusion mechanism could become dominant, in the final zone of the channel, due to a progressive increase of the water concentration. In fact, some authors [25] measured the water transport coefficient along the channel and reported the existence of positive and negative values along the same channels. It is worth mentioning that the negative water transport values indicate that the water flux occurs from the cathode to the anode. It should also be noted, when comparing the two situations in Fig. 8 that for the higher current density (lower voltage) the values of the water mole fraction are in general higher due to higher water production rates and to the occurrence of an increased water electro-osmotic drag.
H2 Hum, Air Hum Exp
H2 Hum, Air Hum Mod
H2 Dry, Air Hum Exp
H2 Dry, Air Hum Mod
H2 Hum, Air Dry Mod
0.8
0.6
0.4
0.2
0.0 0.0
0.1
0.2
0.3
0.4
0.5
0.6
Current Density (A/cm2 )
Fig. 10 e Comparison between experimental and 3D numerical predictions of the polarization curves for 333 K, when humidified both gases with an RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios of za [ 1 and zc [ 2 for 1 A/cm2.
12495
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Fig. 11 e Contours of membrane water content in an XY plane at half of the membrane thickness for 0.6 V, for three humidification situations: a) both gases fully humidified, b) air fully humidified and dry hydrogen and c) dry air and hydrogen fully humidified. Operating conditions: operating temperature of 333 K, H2/Air flow rates calculated with stoichiometry flow ratios of za [ 1 and zc [ 2 for 1 A/cm2.
visualization all over the membrane is an asset of the 3D numerical simulations while the 1D model gives an average value for l the plane XY, for a given value of Z. The ideal distribution of water content should be uniform all over membrane, avoiding dry-out or flooding zones which affect negatively the cell behavior. For these two conditions, a dryer zone tends to appear at the cathode outlet.
4.2.
Reactants relative humidity effect
The effect of reactants relative humidity on cell performance was studied at an operating temperature of 353 K, for three humidification conditions: both gases fully humidified, air fully humidified and dry hydrogen and dry air and hydrogen fully humidified. The 3D numerical predictions of the fuel cell performance are compared to the experimental curves in Fig. 10. Experimentally, the experiment with dry air could not
be performed due to difficulties in the stabilization of the fuel cell. Simulations reveal a better fuel cell performance for the situation when both gases enter in the cell fully humidified, mainly for high current densities. Curiously, for lower current density values, the predictions for fuel cell performance are slightly better when one of the reactants enters dry in the cell. These results are somewhat unexpected, because, for low current density values, the water production rate is lower, so, the fuel cell performance should be enhanced by a better humidification. However, it seems that the humidification of just one stream is sufficient to maintain the membrane hydrated. Probably, there is an initial membrane humidification provided by back-diffusion mechanism in the case of humidified air and by electro-osmotic drag in the case of humidified hydrogen. For high current density values, the water electro-osmotic drag is very intense, leading eventually
1.2 298 K Exp 313 K Exp 333 K Exp
1.0
298 K Model 313 K Model 333 K Model
Voltage (V)
0.8
0.6
0.4
0.2
0.0 0.0
Fig. 12 e Membrane water content along the membrane predicted by the 1D model for 0.6 V and for the three studied humidification conditions. Operating conditions: operating temperature of 333 K, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and zc [ 2 for 1 A/cm2.
0.1
0.2
0.3
0.4
0.5
0.6
Current Density (A/cm2)
Fig. 13 e Comparison between experimental data and 3D numerical simulations of polarization curves for the three temperatures studied a) 298 K b) 313 K and c) 333 K, both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and and zc [ 2 for 1 A/cm2.
12496
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Fig. 14 e Protonic conductivity along the membrane predicted by the 1D model for 0.6 V and for the three studied temperatures. Operating conditions: both gases RH of 95%, H2/Air flow rates calculated with stoichiometry flow ratios za [ 1 and zc [ 2 for 1 A/cm2.
to membrane dehydration, and, consequently, to poorer performances when only one humidified stream is used. This decrease in cell performance is more noticeable when dry air is used, probably because the produced water is dragged by the dry air stream, leading to a more pronounced dehydration in membrane. The 3D simulations are close to experimental data when both gases are fully humidified. For the dry air case, the simulations deviate from experimental data in particular for lower current density values. The contours of water membrane content (in terms of parameter l) along the membrane centered XY plane, for the three studied conditions and 0.6 V, are presented in Fig. 11. The distribution of the membrane water content is not uniform all over the membrane and is substantially higher when both gases enter in the cell humidified. The lowest values are obtained when dry air is used. These results are in accordance with the 3D simulated polarization curves for this operating condition (0.6 V), corresponding to better performances when both gases are humidified and lower for the
case of dry air. For this latter case, the water content values are extremely low revealing dehydration of the membrane. This fact is in accordance with the experimental difficulties in reproducing this operating condition. The average values of the water content along the membrane for the same value of voltage predicted by the 1D model are displayed in Fig. 12 for the different humidification conditions. As can be seen the trend is similar: higher average water content values are obtained when both gases are humidified and lower values for the case of dry air. When comparing the average water content at the center of the membrane (z ¼ 0.013 cm) with the information given by the water contours in Fig. 11, it is evident that the predictions provided by the simpler model 1D are higher that the 3D simulations for the cases of one humidified stream. This discrepancy is somewhat expected, since the 1D model predicts higher performances (recall Fig. 5) corresponding to higher values of the current densities for the same voltage value and consequently enhanced predicted average water contents.
4.3.
Operating temperature effect
The operating temperature effect was studied for 298 K, 313 K and 333 K. The comparison between 3D numerical simulations and experimental data are presented in Fig. 13. As can be seen from Fig. 13, the 3D predicted curves are close to the experimental data. The curves are close to each other showing a slight increase of fuel cell performance with increasing temperature for the entire range of current densities covered. For this type of membranes (very thin), the improvement of performance with increasing temperature is not perceptible in experimental data. A well-known impact of an increase in the cell temperature is the enhancement of the membrane conductivity and, consequently, the decrease of the membrane resistance to the transport of hydrogen protons. Following this idea, the predictions of the average values of the proton conductivities along the membrane and the contours of protonic conductivity for the membrane
Fig. 15 e Contours of protonic conductivity in an XY plane at half of the membrane thickness given by the 3D numerical simulations, for 0.6 V and for the three studied temperatures a) 298 K b) 313 K and c) 333 K. Operating conditions: both gases RH of 95%, H2/Air flow rates calculated with za [ 1 and with zc [ 2 for 1 A/cm2.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
middle plane are presented in Figs. 14 and 15 for the three studied conditions. The simulation results from the two Figures show, effectively, an increase in the proton conductivity with increasing temperatures. The curves in Fig. 14 also reveal an increase of the proton conductivity along the membrane accompanying the increase in water content. The contours in Fig. 15 put in evidence the non-homogeneity of proton conductivity in the studied plane. The ideal situation should correspond to high values of protonic conductivity with a uniform distribution all over the membrane in order to facilitate the protons flux. The comparison between the average value at the center of the membrane (z ¼ 0.013 cm) taken from Fig. 14 with the information given by the contours in Fig. 15, puts in evidence the accordance between 1D and 3D predictions. For 298 K, the contour exhibits proton conductivity values between 5.6 U1 m1 and 5.9 U1 m1, while the 1D model predicts an average value for the proton conductivity in the membrane middle plane of 6.3 U1 m1. For 313 K, the contour gives values between 6.2 U1 m1 and 7.3 U1 m1, while the predicted average value in the centered plane is 7.3 U1 m1. Finally, for 333 K, the 3D predicted values are between 7 U1 m1 and 10.8 U1 m1, while the 1D model predicts a proton conductivity of 8.6 U1 m1. One can conclude that the developed one-dimensional model gives adequate predictions for the proton conductivity expending a much lower CPU computing time.
5.
Conclusions
In the present work a comparison was performed between 1D and 3D numerical simulations in PEM fuel cells. The onedimensional model used was previously developed and validated [1]. The 3D numerical simulations were obtained using the commercial code FLUENT-PEMFC module. Particular attention was devoted to the impact of considering two-phase flow effects in the 3D simulation results. The 3D simulations considering one and two phases for water flow and the 1D model simulations were compared with experimental data. The results showed that 1D and 3D model simulations considering just one phase are similar, with a slightly better accordance for the 1D model exhibiting a substantially lower CPU time. However both numerical results over predict the fuel cell performance. Actually, the 3D simulations reproduce very well the experimental data when two-phase flow effects are considered. The formation of liquid water (considered in this case) tends to block the pores, leading to a decrease in fuel cell performance. Three-dimensional numerical simulations of oxygen, hydrogen and water mole fractions along the channels were analyzed. The distribution of mole fractions for the three species is not uniform along the channels in particular for higher current densities (corresponding to higher electrochemical reaction rates). The water mole fraction does not follow a defined trend, showing for certain operating conditions an increase and afterwards a decrease near the channel outlet. This finding is in accordance with published results
12497
showing the change of the water dominant mechanism from the cell inlet to the cell outlet [24]. The effect of the relative humidity of gases and operation temperature on fuel cell performance was also studied both through the comparison of the polarization curves for the 1D and 3D simulations and experimental data and through the analysis of relevant physical parameters such as the water membrane content and the proton conductivity. The results showed that, as expected, the fuel cell performance is improved when both gas streams are humidified. The 1D and 3D simulations corroborate this finding, revealing higher water content values for this condition. Fuel cell performance slightly increases with increasing temperature, being this trend perceptible in the simulated polarization curves and confirmed by the protonic conductivity predictions. The protonic conductivity values are higher for higher temperature values. The water content predictions by the simpler 1D model are in agreement with the 3D numerical simulations, except for the cases when one stream enters dry in the cell. In these cases, the simpler model predicts higher water content values. Concerning the predictions of the proton conductivities, the 1D and 3D results are in agreement. The presented results showed that the developed 1D model (coupling heat and mass transfer effects) can be used to predict optimal operating conditions in PEMFCs and the general trends of the impact of several important physical parameters in fuel cell performance such as membrane water content or the membrane protonic conductivities. The use of the 3D numerical simulations is indicated if more detailed predictions are needed. An asset of this type of simulations is the prediction of the spatial distribution and visualization of various relevant parameters. However simulation times are still too high. A polarization curve with the 1D model is obtained with a CPU time around 5 min, while the 3D computing time is around 24 h. Besides this disadvantage, before the generation of the 3D numerical simulations, it is necessary to create the mesh for the whole cell, (in this case, GAMBIT was the mesh generator used), a very difficult task due to the discrepancy between the channels dimensions and the rest of the cell geometry. An important conclusion of this work is the demonstration that a simpler analytical model using low CPU has potential to be used in real-time PEMFC simulations. Empirical models can always be used for quick predictions, however, they are «validated» against polarization curves using a significant number of fitting parameters. The smart use of «old fashioned» analytical or semi analytical models using a rationally chosen set of parameters to obtain useful reduced models is put in evidence in the present work.
Acknowledgements The partial support of “Fundac¸a˜o para a Cieˆncia e Tecnologia Portugal” through scholarship SFRH/BD/28166/2006 is gratefully acknowledged. POCI (FEDER) also supported this work via CEFT.
12498
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 1 2 4 8 6 e1 2 4 9 8
Nomenclature
CACL H2 ref
CH2 CCCL O2 ref
CO2 cr E0 EM ES ref ja ref jc K r rw Ra Rc RM
RS s S vl Yj, cen Yj, sur
hydrogen concentration at anode catalytic layer, mol/m3 hydrogen reference concentration, mol/m3 oxygen concentration at anode catalytic layer, mol/ m3 oxygen reference concentration, mol/m3 condensation rate constant, s1 open circuit voltage, V phase potential of the membrane, V phase potential of the solid, V anode volumetric transfer current, A/cm3 cathode volumetric transfer current, A/cm3 absolute permeability, m/s specific reacting surface area, m1 condensation Rate, m/s anode volumetric transfer current, A/m3 cathode volumetric transfer current, A/m3 membrane volumetric transfer current, A/m3 solid volumetric transfer current, A/m3 water saturation source term liquid velocity, m/s mass fraction of species j at the cell center mass fraction of species j at the reactant surface
Greek letters anode transfer coefficient aa cathode transfer coefficient ac d average distance between reaction surfaces and cell center, m anode concentration dependence ga cathode concentration dependence gc anode overpotential, V ha cathode overpotential, V hc l membrane water content q contact angle, s surface tension, N/m2 sM membrane electric conductivity, S/m solid electric conductivity, S/m sS anode stoichiometric ratio za cathode stoichiometric ratio zc
references
[1] Falca˜o DS, Oliveira VB, Rangel CM, Pinho C, Pinto AMFR. Water transport through a PEM fuel cell: a one-dimensional model with heat transfer effects. Chemical Engineering Science 2009;64(9):2216e25. [2] Falcao DS, Rangel CM, Pinho C, Pinto AMFR. Water transport through a proton-exchange membrane (PEM) fuel cell operating near ambient conditions: experimental and modeling studies. Energy & Fuels 2009;23(1):397e402.
[3] Eikerling M, Kornyshev AA. Modelling the performance of the cathode catalyst layer of polymer electrolyte fuel cells. Journal of Electrochemical Society 1998;453(1e2):89e106. [4] Eikerling M, Kornyshev AA, Kucernak AR. Water in polymer electrolyte fuel cells: friend or foe? Physics Today 2006;59(10): 38e44. [5] Baschuk JJ, Li X. Modelling of polymer electrolyte membrane fuel cells with variable degrees of water flooding. Journal of Power Sources 2000;86(1e2):181e96. [6] BiyIkoglu A. Review of proton exchange membrane fuel cell models. International Journal of Hydrogen Energy 2005; 30(11):1181e212. [7] Ji M, Wei Z. A review of water management in polymer electrolyte membrane fuel cells. Energies 2009;2(4):1057e106. [8] Springer TE, Zawodzinski TA, Gottesfeld S. Polymer electrolyte fuel cell model. Journal of Electrochemical Society 1991;138(8):2334e42. [9] Bernardi DM, Verbrugge MW. A mathematical model of the solid-polymer-electrolyte fuel cell. Journal of Electrochemical Society 1992;139(9):2477e91. [10] Kulikovsky AA. Semi-analytical 1D þ 1D model of a polymer electrolyte fuel cell. Electrochemistry Communications 2004; 6(10):969e77. [11] Scott K, Argyropoulos P. A current distribution model of a porous fuel cell electrode. Journal of Electroanalytical Chemistry 2004;567(1):103e9. [12] De Francesco M, Arato E, Costa P. Transport phenomena in membranes for PEMFC applications: an analytical approach to the calculation of membrane resistance. Journal of Power Sources 2004;132(1e2):127e34. [13] Wang ZH, Wang CY, Chen KS. Two-phase flow and transport in the air cathode of proton exchange membrane fuel cells. Journal of Power Sources 2001;94(1):40e50. [14] Um S, Wang CY, Chen KS. Computational fluid dynamics modeling of proton exchange membrane fuel cells. Journal of Electrochemical Society 2000;147(12):4485e93. [15] Ju H, Meng H, Wang C-Y. A single-phase, non-isothermal model for PEM fuel cells. International Journal of Heat and Mass Transfer 2005;48(7):1303e15. [16] Berning T, Lu DM, Djilali N. Three-dimensional computational analysis of transport phenomena in a PEM fuel cell. Journal of Power Sources 2002;106(1e2):284e94. [17] Le AD, Zhou B. A general model of proton exchange membrane fuel cell. Journal of Power Sources 2008;182(1):197e222. [18] Le AD, Zhou B. A generalized numerical model for liquid water in a proton exchange membrane fuel cell with interdigitated design. Journal of Power Sources 2009;193(2):665e83. [19] Meng H. A three-dimensional PEM fuel cell model with consistent treatment of water transport in MEA. Journal of Power Sources 2006;162(1):426e35. [20] Wu H, Li X, Berg P. On the modeling of water transport in polymer electrolyte membrane fuel cells. Electrochimica Acta 2009;54(27):6913e27. [21] Dawes JE, Hanspal NS, Family OA, Turan A. Threedimensional CFD modelling of PEM fuel cells: an investigation into the effects of water flooding. Chemical Engineering Science 2009;64(12):2781e94. [22] Jang J-H, Yan W-M, Li H-Y, Tsai W-C. Three-dimensional numerical study on cell performance and transport phenomena of PEM fuel cells with conventional flow fields. International Journal of Hydrogen Energy 2008;33(1):156e64. [23] Fluent. Fluent 6.3-fuel cells module manual; 2006. [24] Sivertsen BR, Djilali N. CFD-based modelling of proton exchange membrane fuel cells. Journal of Power Sources 2005;141(1):65e78. [25] Liu F, Lu G, Wang C-Y. Water transport coefficient distribution through the membrane in a polymer electrolyte fuel cell. Journal of Membrane Science 2007;287(1):126e31.