DYNAMIC MODELING OF POLYMER ELECTROLYTE MEMBRANE FUEL CELL STACK WITH 1D AND 2D CFD TECHNIQUES Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classified information. Yuyao Shan Certificate of Approval: Jay M. Khodadadi Professor Mechanical Engineering Jeffrey Fergus Associate Professor Materials Engineering Song-Yul Choe, Chair Associate Professor Mechanical Engineering Stephen L. McFarland Dean Graduate School DYNAMIC MODELING OF POLYMER ELECTROLYTE MEMBRANE FUEL CELL STACK WITH 1D AND 2D CFD TECHNIQUES Yuyao Shan A Thesis Submitted to the Graduate Faculty of Auburn University in Partial Fulfillment of the Requirements for the Degree of Master of Science Auburn, Alabama August 7, 2006 iii DYNAMIC MODELING OF POLYMER ELECTROLYTE MEMBRANE FUEL CELL STACK WITH 1D AND 2D CFD TECHNIQUES Yuyao Shan Permission is granted to Auburn University to make copies of this thesis at its discretion, upon request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iv VITA Yuyao Shan, son of Xiufang Wen and Baofu Shan, was born on September 17, 1981 in Chen Ba Er Hu Qi, Nei Mongolia Province, the Peoples? Republic of China. He entered Tsinghua University in September 1999 and earned a Bachelor of Engineering degree in Thermal Science in July 2003. In August 2003, he entered Auburn University to pursue a Doctor of Philosophy degree in Department of Mechanical Engineering. In November 2005, he changed to master program. v THESIS ABSTRACT DYNAMIC MODELING OF POLYMER ELECTROLYTE MEMBRANE FUEL CELL STACK WITH 1D AND 2D CFD TECHNIQUES Yuyao Shan Master of Science, August 7, 2006 (B.S., Tsinghua University, 2003) 156 Typed pages Directed by Song-Yul Choe The PEM fuel cell is a potential candidate for use as an alternative power source in future vehicle and power conditioning applications. The electric power, as a load, continuously varies as a function of time in the system. Accordingly, the flow rates of fuels and reactants should be properly supplied to meet the dynamics of load requirement. Hence, temperature and humidity for the stack must be precisely maintained to secure an efficient, safe and durable operation of the PEM fuel cell system. As a matter of fact, the dynamic power output and efficiency profile of a PEMFC is strongly influenced by the variation of the temperature, reactant and product transfer in the fuel cell caused by a current load. In this work, a new 1-D dynamic model and a new 2-D dynamic model for the PEM fuel cell stack are proposed. Emphasis is placed on dynamic thermal effects of the PEM fuel cell stack. The 1-D model has been finished and was used to simulate the dynamic responses and performance of the stack during the start-up procedure and with vi dynamic load input. The results give interesting information on the stack temperature distribution and could be used for stack components and control design. vii ACKNOWLEDGEMENTS The author would like to thank the members of his research committee Dr. Song-Yul Choe, Dr. Jay Khodadadi, Dr. Fergus, for their constructive advices and complete support during his research. He would also like to thank Mr. Suk-Heung Song, Dr. Choe, for their help in his work. Thanks are also given to his group member, and other faculty and graduate students who were able to given their helps during this research. viii Computer software used Microsoft word 2007 pre-release, MATLAB/SIMULINK, C language ix TABLE OF CONTENTS LIST OF TABLES xi LIST OF FIGURES xii NOMENCLATURE xvi 1. INTRODUCTION .................................................................................................. 1 1.1. FUEL CELL TECHNOLOGY ............................................................................... 1 1.2. POLYMER ELECTROLYTE MEMBRANE FUEL CELL .................................. 4 1.2.1. PRINCIPLE OF OPERATION OF THE PEM FUEL CELL .................... 4 1.2.2. CHALLENGE ISSUES .............................................................................. 8 1.3. MODELING OF PEM FUEL CELLS.................................................................. 14 1.3.1. GENERAL DESCRIPTION..................................................................... 14 1.3.2. STATE-OF THE-ART OF MODELS FOR LAYERS............................. 15 1.3.3. SINGLE CELL MODEL .......................................................................... 20 1.3.4. STACK MODEL ...................................................................................... 28 1.3.5. SUMMARY.............................................................................................. 28 1.4. THESIS OBJECTIVES ........................................................................................ 29 REFERENCES ................................................................................................................. 30 2. 1-D PEM FUEL CELL STACK MODELING..................................................... 34 2.1. INTRODUCTION ................................................................................................ 34 2.2. 1-D SINGLE CELL MODELING........................................................................ 37 x 2.2.1. MODEL DESCRIPTION ......................................................................... 37 2.2.2. SIMULATION SETUP ............................................................................ 49 2.2.3. ANALYSES.............................................................................................. 49 2.3. 1-D STACK MODELING.................................................................................... 59 2.3.1. MODEL DESCRIPTION ......................................................................... 59 2.3.2. SIMULATION SETUP ............................................................................ 64 2.3.3. ANALYSIS............................................................................................... 66 2.4. SUMMARY.......................................................................................................... 88 REFERENCES ................................................................................................................. 91 3. 2-D SINGLE CELL MODELING........................................................................ 95 3.1. INTRODUCTION ................................................................................................ 95 3.2. MODEL DESCRIPTION ..................................................................................... 98 3.2.1. MODELING DOMAIN AND ASSUMPTIONS ..................................... 98 3.2.2. MODEL DESCRIPTION ....................................................................... 100 3.3. NUMERICAL SOLUTION................................................................................ 110 3.4. ANALYSES........................................................................................................ 115 3.5. CONCLUSION................................................................................................... 131 REFERENCES ............................................................................................................... 133 4. CONCLUSION AND FUTURE WORK ........................................................... 136 4.1. CONCLUSION................................................................................................... 136 4.2. FUTURE WORK................................................................................................ 137 xi LIST OF TABLES Chapter 1: Table 1-1. Comparison of three types of fuel cells???????????????3 Table 1-2. Type of models and governing equations employed??????????12 Table 1-3. Summary of parameters, variables and equations by CFD analysis????14 Chapter 2: Table 2-1. Empirical parameters for quasi 1-D models?????????????49 Table 2-2. Geometry parameters for the fuel cell???????????????53 Table 2-3. The governing equations describing phenomenon for a single cell????62 Table 2-4. Working condition for PEM fuel cell stack with low inlet RH??????65 Table 2-5. Working condition for PEM fuel cell stack with high inlet RH?????65 Chapter 3: Table 3-1. Parameters for 2-D models???????????????????114 Table 3-2. Geometry parameters for the fuel cell???????????????114 Table 3-3. Initial values????????????????????????115 xii LIST OF FIGURES Chapter 1: Fig. 1-1. Schematic simulation domain????????????????..???4 Fig. 1-2. Chemical reaction processes in a PEM fuel cell????????????6 Fig. 1-3. Polarization curve????????????????????????7 Fig. 1-4. Thin film set up for catalyst????????????????????12 Fig. 1-5. Agglomerate setup???????????????????????13 Chapter 2: Fig. 2-1 Schematic simulation domain????????????????..???37 Fig. 2-2: Voltage response by a consideration of proton concentration???????44 Fig. 2-3 Current input??????????????????????????45 Fig. 2-4 Oxygen concentration response dependent upon GDL thickness at P=1 bar and T=353?K?????????????????????????????46 Fig. 2-5: dynamic cathode pressure input??????????????????46 Fig. 2-6: Oxygen concentration response to the dynamic current input and cathode pressure at different depths??????????..?????????????47 Fig. 2-7: Oxygen concentration dynamic response for different GDL thicknesses??47 Fig. 2-8: Dynamic response of oxygen concentration for different porosities????48 Fig. 2-9: Dynamic response of oxygen concentration for different tortuosities????48 xiii Fig. 2-10: IV curve for different cell working temperature. Cell: P=1.0bar?????50 Fig. 2-11: IV curve for different cell gas pressure. Cell: T=353.15K???????50 Fig. 2-12: I-V curve for different stoichimetric number. Cell: T=353.15K, P=1.0bar?51 Fig. 2-13: IV curve for different RH at the cathode inlet. Cell: T=353.15K, P=1.0bar?52 Fig. 2-14 Simulation setup ????????????????????????52 Fig. 2-15 Simulation results a) Dynamic current, b) Dynamic temperature response of different layers of single cell, c) Dynamic voltage/power response, d) Single cell efficiency???????????????????????????????54 Fig. 2-16 Temperature gradient across the cell (Left to right: cathode coolant channel to anode coolant channel) a) 50 seconds b) 7 minutes after startup?????????56 Fig. 2-17 Analysis of transient behavior of temperature, voltage, power, and efficiency upon a current step a) Current, b) Temperature in the different layers c) Voltage and power d) Efficiency analysis???????????????????????58 Fig. 2-18: Stack schematic configuration????????????????..??61 Fig. 2-19: I-V curve for different cell working temperature???????????66 Fig. 2-20 a) Current input for the startup b) Dynamic membrane temperature c) Stack temperature distribution at the 50 second d) Stack temperature distribution at the 360 second e) Stack temperature distribution at the 450 second???????????69 Fig. 2-21 a) Dynamic voltage output b) Voltage output percentage????????72 Fig. 2-22 Dynamic behavior of membrane water uptake for a stack????????74 Fig. 2-23 a) Current input for the startup b) Dynamic membrane temperature c) Stack temperature distribution at the 50 second d) Stack temperature distribution at the 270 second e) Stack temperature distribution at the 350 second???????????76 xiv Fig. 2-24 a) Dynamic voltage output b) Voltage output percentage????????78 Fig. 2-25 Dynamic membrane water uptake?????????????????79 Fig. 2-26 a) Step current input b) Dynamic membrane temperature c) Stack temperature distribution at the 580 second d) Stack temperature distribution at the 650 second e) Stack temperature distribution at the 800 second?????????80 Fig. 2-27 a) Dynamic voltage output b) Dynamic voltage output (enlarge) c) Voltage output percentage d) Dynamic membrane water uptake ????????????82 Fig. 2-28 a) Step current input b) Dynamic membrane temperature c) Stack temperature distribution at the 400 second d) Stack temperature distribution at the 550 second e) Stack temperature distribution at the 800 second?????????84 Fig. 2-29 a) Dynamic voltage output b) Dynamic voltage output c) Voltage output percentage d) Dynamic membrane water uptake??..?????????????87 Chapter 3: Fig. 3-1 Single cell schematic configuration?????????????????99 Fig. 3-2: Stack schematic configuration???????????????????99 Fig. 3-3 Overall solution procedure????????????????????111 Fig. 3-4 Stack geometry information???????????????????112 Fig. 3-5 cell geometry information????????????????????113 Fig. 3-6 Mesh generated in the 2-cell stack?????????????????114 Fig. 3-7: Pressure drop in the cell 1????????????????????116 Fig. 3-8: Velocity profile????????????????????????117 Fig. 3-9: Dynamics of oxygen concentration at two cells???????????118 xv Fig. 3-10: Current density on the cathode side????????????????121 Fig. 3-11: Membrane potential dynamics at the cell 1?????????????123 Fig. 3-12: Stack temperature 2-D?????????????????????125 Fig. 3-13: Fluid temperature???????????????????????127 Fig. 3-14: tack temperature 1-D?????????????????????128 Fig. 3-15: Cell temperature 1-D?????????????????????130 xvi NOMENCLATURE Alphabets: Units a Catalysts area ratio m 2 /m 3 a Water activity a Partial pressure Pa b Membrane extension coefficient m 2 /m 3 c Mole concentration mole m -3 C p Thermal capacity J mol -1 K -1 D Diffusion coefficient m 2 /s F Faraday number C mol -1 I Current density A m -2 j Current density A m -3 J Mass flow rate kg m -2 k Thermal conductivity W m -1 K -1 K Gas permeability m 2 L Length m M Mole mass kg mol -1 N Mass flux kg n Electro-osmotic coefficient p Pressure (partial pressure) Pa xvii Q Heat transfer rate Q m -2 sec -1 R Universal gas constant J mol -1 K -1 R Electrical resistance Ohm S Source term S Percentage of channel saturated s Entropy J mol -1 K -1 T Temperature K Th Thickness m u Velocity m s -1 V Voltage V X Dimensionless mol concentration m -3 Greeks: ? Transfer coefficient ? Water transfer coefficient m s -1 ? Porosity ? Conductivity s m -1 ? Potential V ? Viscosity kg m -1 s -1 ? Water content ? Over-potential V ? Contact angle 90.02? ? Density kg m -3 ? Tortuosity xviii Superscripts, subscripts: a Anode c Cathode cata Catalysts ch Channel cv Control volume conv Convection cond Conduction diff Diffusion effect e Electrolyte eff Effective ele Electro osmotic effect eq Equilibrium f Fluid g Gas gdl Gas diffusion layer i, j, k Index l Liquid in Inlet condition m Mixture including solid and fluid mass Mass concentration mem Membrane mid Middle point xix OC Open circuit ref Reference value s Solid matrix sat Saturation T Thermal tota Anode total totc Cathode total v Vapor 1 CHAPTER 1: INTRODUCTION 1.1. FUEL CELL TECHNOLOGY To meet the daily life needs at advancing civilization of human beings on the earth, internal combustion engines in vehicles and combustors in industrial plants, are continuously exhausting air pollutants (for example NOx and CO 2 ) all over the world by consuming the finite fossil energy. As a result, technologies for sustainable energy sources are being researched by scientists all over the world. Intensive studies conducted by academia and industries, and government agencies have proven that the fuel cell technology is one of the potential candidates that can replace the technology based on the fossil fuels because of the highest efficiency among other alternatives and near zero emission. Additionally, the fuel cell technology can cover wide range of applications as a power source needed for portable electronics, utilities, and transportations. Among different technologies of the fuel cells, the technology based on a polymer electrolyte membrane (PEM Fuel cell) has been shown the most viable technology because of the low operating temperature and high efficiency along with high power density. However, various tests on the field have shown that the durability and reliability of the fuel cell does not meet the requirements for a power source currently used in the various applications. The most challenging issue identified is the management of water generated in the cell as a byproduct in an operating environment, where the ambient temperature 2 can drop under a freezing value. In addition, the cost for a system is at least ten times higher than current technology used. These challenging technical issues require fundamental understandings in the mechanism that associates with fluid and thermal science along with properties of the materials used for a construction of the cell. Therefore, the work presented has been approached to mathematically describe the mechanisms by applying electrochemical, fluid and thermal equations that were computationally solved. Final analyses complement understandings in the mechanisms and parameter dependency of a PEM fuel cell in operating conditions in conjunction with geometrical factors and properties of the materials. A fuel cell is a device converts the chemical energy directly into the electrical energy via chemical reaction. It consists of two electrodes, an anode to which the fuel (hydrogen for example) is fed externally and a cathode to which the oxidant is supplied, and the electrolyte which separates the two electrodes, blocks the reactants and electrons transfer, and allows the ions to flow across it. The most viable technologies among others are polymer electrolyte membrane fuel cell (PEMFC), direct methanol fuel cell (DMFC), and solid oxide fuel cell (SOFC). A comparison among these three types of fuel cell technologies are briefly summarized in table 1-1. The advantages of the PEMFC are a low operating temperature and consequently possible quick startup. However, due to the low operating temperature the catalysts need a large amount of precious material like platinum. Conversely, the highly loaded catalysts are likely to get covered and the performance of the cell drops rapidly when the fuel stream includes impurities. Particularly, the tolerance of the catalysts to carbon monoxide decreases when temperature decreases. 3 The advantage of DMFC has shown that the energy density of the liquid methanol is orders of magnitude greater than even highly compressed hydrogen. Adversely, the low efficiency associated with a high permeation of the fuel through the membrane, the sluggish dynamic behavior, and the poisonous property of the fuel has shown as drawbacks. Conversely, SOFC has revealed high efficiency, fuel and catalysts flexibility because of a high operating temperature, but slow startup process associated with a relative large amount of the thermal mass and the corresponding reliability issues on the sealing materials. Name Polymer Electrolyte membrane fuel cell Direct methanol fuel cell Solid oxide fuel cell Mobile ion Proton Proton Oxygen anion Electrolyte Solid organic polymer polyperfluorosulfonic acid Solid organic polymer polyperfluorosulfonic acid Solid zirconium oxide to which a small amount of yttria is added Working temp. 60?100 ?C, (70?C) 90?120 ?C 600?1000 ?C elec. eff. Cell: 50?70 % System: 30?50 % Cell: 20?30 % Cell: 60?65 % System: 55?60 % Power range 0.1 to 500 kW mW to 100 kW Up to 100 MW power density 2.6-5.0 kW/m 2 0.6 kW/m 2 1.5-6.5 kW/m 2 Application Automotive/small stationary Portable Stationary Startup time Seconds to minutes Seconds to minutes Hours Table 1-1. Comparison of three different types of the fuel cell technologies [43] 4 1.2. POLYMER ELECTROLYTE MEMBRANE FUEL CELL 1.2.1. PRINCIPLE OF OPERATION OF THE PEM FUEL CELL As shown in Fig. 1-1, a typical PEM fuel cell consists of several separate layers: two flow fields for hydrogen and oxygen by supplying air, two gas diffusion layers, two catalyst layers and one membrane. The hydrogen reacts with oxygen and generated electricity and water and heat as byproducts. Fig. 1-1 Schematic diagram of a fuel cell 1.2.1.1. CHEMICAL REACTION When one mole of hydrogen molecule reacts with one half mole of oxygen molecule, one mole of water molecule and heat (Q) are produced. 222 1 2 HOHOQ+?+ (1) The chemical reactions separately occur in two electrodes. When a mole of hydrogen molecule reaches to the anode catalyst, the chemical bond of the molecule is broken and the two hydrogen atoms are absorbed at the electrode surface. After the hydrogen atom 5 has lost an electron, the atom is positively charged and becomes a proton. Then, protons are transferred to the cathode through the membrane, while electrons are being transferred to the cathode via the external circuit. 2 2HPt PtH+? ? (2) 222Pt H H e +? ?? + (3) Likewise, when a half mole of oxygen molecule reaches to the cathode catalyst, the chemical bond of the oxygen is broken and then one mole of oxygen atoms is absorbed at the electrode surface. At the same time, the proton combines with the electron. Subsequently, the oxygen atom reacts with the absorbed hydrogen atom and generates one mole of water molecule. 2 2OPt PtO+? ? (4) He PtH +? +?? (5) 2 2Pt H Pt O H O?+?? (6) These processes are shown in Fig. 1-2 6 Fig. 1-2 Chemical reaction processes in a PEM fuel cell 1.2.1.2. THERMODYNAMIC ANALYSIS The ideal voltage that a PEM fuel cell can generate is calculated by using the thermodynamic analysis of the chemical reaction, where the Gibbs free energy is released. Then, the change of the Gibbs free energy is the abstraction of the change of enthalpy and entropy at a temperature, GHTS?=??? (7) If all the change of Gibbs free energy is converted to electrical energy, then we get: OC OC G GnVFV nF ?? ?=? ? = (8) where n is the mole of electrons transferred. 7 1.2.1.3. POLARIZATIONS The voltage output of a fuel cell is decreased from its ideal one because of several irreversible losses. Besides losses caused by reactants crossover through electrolyte etc, there are mainly three kinds of polarizations (Fig. 1-3): (1) activation polarization, (2) ohmic polarization, and (3) concentration polarization. The activation polarization loss is dominant at low current density. Ohmic polarization is proportional to the current density. The concentration polarization is caused by the gas transport across the gas diffusion media. Fig. 1-3 Polarization curve [42] Activation Polarization is directly related to the rates of electrochemical reactions. There is an activation barrier that must be overcome by the reacting species. Ohmic Polarization presents ohmic losses by resisting the flow of protons and electrons in the conduction regions, the proton transfer through the electrolytes and the electrons transfer through the electrodes. The former one plays a dominant role causing the losses. 8 The resistance of the membrane (Nafion) is affected by the concentration of water and the local temperature. Even though, researchers and scientists are trying to develop high temperature membranes which can conduct protons without the presence of water, Nafion is still one of the best choices. Concentration Polarization is referring to the losses of potentials because of the decrease of reactant concentrations when the reaction occurs. In fact, a gas concentration gradient is built by the flow of reactants from the gas channel to the catalyst layer, the dissolving process of reactants into the electrolyte and diffusion to the catalyst sites. This effect is dominant at high current density where the flooding occurs on the cathode side. 1.2.2. CHALLENGE ISSUES Although the fuel cells have been successfully demonstrated in both automotive and stationary power applications, there are numerous fundamental technical and logistic issues that still have to be resolved before PEM fuel cell systems are commercialized, and the quest is on for membrane and electrode materials and structures, and flow field configurations, that will enable even higher, more stable and more reliable performance. Challenges for developing PEM Fuel Cell are: (1) Cost and performance, (2) Durability and stability, (3) Balance-of-plant components including compressors, humidifiers, heat exchangers, sensors, and controls, (4) Air, Water and Thermal management and (5) Stack and system design for a better performance ? Cost is the greatest challenge that results from expensive, precious-metal catalysts and others requiring costly materials that can sustain high temperatures? The cost may be reduced by decreasing precious metal required and by developing low-cost, continuous 9 fabrication processes for fuel cell components such as membrane-electrode assemblies (MEAs) and bipolar plates. Maximal achievable power density directly translates to smaller thus less expensive fuel cells. Although maximum power density is achieved at ~0.5 V, nominal power density is typically selected at 0.6-0.7 V to get a higher efficiency. Current densities of 1.7 A/cm 2 at 0.6 V has been achieved in a single cell and 1.25 A/cm 2 at 0.6 V in a cell stack configuration that correspond to 1.0 W/cm 2 and 0.75 W/cm 2 , respectively. A technical challenge is to achieve single cell performance in multi-cell stacks, which requires cell and stack designs that ensure uniform conditions over each entire cell active area. There are many factors that affect the fuel cell performance and achievable power densities. Some of these factors include: type and thickness of the polymer membrane, electrode kinetics, i.e., electrode structure, catalyst loading and catalyst utilization, type of backing layer (structure, thickness, porosity, tortuosity, hydrophobicity), hardware resistance, main interface contact resistance, flow field configuration, and operating conditions such as temperature, pressure, flow rates, humidification of reactant gases. Another technical challenge is the need to increase fuel cell durability and dependability. Particularly, degradation of membrane performance and damage is one of the crucial factors caused by insufficient moisture and/or hotspots. Therefore, the PEM fuel cells must have effective water management systems to operate dependably and efficiently. On the other hand, there is on-going effort to develop high temperature membranes that do not require water for the transport of protons. Overall performance of a stack can be effectively analyzed by computational methods and impacts on system can be predicted. 10 Research into these areas is ongoing, and the durability of new components and designs is necessary. The complexity of balance of plant (BOP) of the fuel cell system can also be greatly reduced with the new technology developed for the fuel cell. For example, the company, UTC Fuel Cells, has come up with a new design of porous bipolar plate which has completely removed the need for gas humidifiers and best of all; the atmospheric air input is possible without any reduction in the system efficiency compared to the pressurized air fuel cell system. However, the manufacturability of the porous bipolar plates and water resided in the pores are presenting other barriers on the way to a commercialization. The fuel cell industry has traditionally focused on research and development, and proving concepts via hand made prototypes. As the technology matures, focus will shift from research to manufacturing. Efficient and automated fuel cell diagnostics should be developed to assist users in diagnosing system failures of a fuel cell system during factory acceptance testing or for fuel cell stack system control. The current methods for developing a fuel cell system and diagnosing the structures are almost purely empirical. This leads to a long, inefficient development cycle. There are needs for tools/instruments to probe the structures and understand the interactions between the components of the fuel cell and monitor the fuel cell operation on line. However, several technical barriers should be overcome before it could be used widely in a real world application. Three major areas are identified as follows: 1) thermal management for optimization of temperature control and a high efficient operation; 2) water management that includes design parameters of GDL materials and structures as 11 well as operating conditioned water flooding in the GDL and dehydration in membrane; 3) material selection, gas flow channel and coolant flow channel design for BP. Solving these problems can make the fuel cell operating under more reliable and safe mode, with high efficiency, also decrease the consumption of fuel in the subjects of transportation. Modeling and simulation are one of the main strategies to achieve these objectives. There are mainly two directions in the area of fuel cell modeling and simulation. One direction is to establish computationally efficient models including empirical modeling, equivalent circuit modeling and some mixed models for the purpose of stack performance evaluation, system control and balance of plant. The other direction is to start from physical and kinetic theories to build detailed and complete models using computational fluid dynamics (CFD), in order to understand the electrochemical and physical mechanism within the fuel cell, and in turn to improve the design and preparation of fuel cell for better performance and efficiency. Various CFD technologies have been applied to develop a fuel cell model. These efforts are led by industry, the National Laboratories and universities. The fuel cell has a cubical shape that can be mapped in a Cartesian coordinates and is mathematically described. The number of axis determines the degree of the model whether it is one, two or three-dimensional model. Basically, the model describes the major transport phenomena in components of a fuel cell that compromise of two flow fields, two gas-diffusion layers, two catalyst layers and the membrane. Others also include collector plates, cooling channels and other components in the modeling domain. Simulation can yield detailed multi-dimensional distributions of gas concentrations, current density, over-potential, temperature, and water content and ultimately polarizations curves. Major mechanistic models of PEM fuel 12 cells proposed in the literature will be briefly summarized by categorizing components and the degree of used dimensions. Table 1-2 shows types of models for each layer of a single cell and the applied governing equations. Physical Layer Membrane Catalysts GDL Bipolar plate Model type - Springer empirical model - Bernardi?s convection model - General Stefan- Maxwell Model - Percolation Theory Model - Pseudo- homogeneous - Thin film (Fig. 1-4) - Agglomerate (Fig. 1-5) - Micro-scale - Multiphas e mixture model - Separate model Governing Equations Water & Proton transfer: Nernst-Plank/ General Stefan-Maxwell; water uptake: empirical /N2 Layer Brunauer- Emmett-Teller equation/Flory-Huggins Continuity; Diffusion- reaction Continuity; Navier- Stokes; Stefan- Maxwell /Fick?s law; Darcy law Continuity; Navier- Stokes Table 1-2. Type of models and governing equations employed Fig. 1-4. Thin film set up for catalyst [6] 13 Fig. 1-5. Agglomerate setup [6] Fig. 1-4 and Fig. 1-5 show the setup for the thin film model and agglomerate model for the catalyst layer in the PEM fuel cells Table 1-3 shows a summary of the important parameters, variables and equations that are used in CFD models to describe a single cell except the gas channel, where the gas flow is described by the Navier-Stokes equation. Anode GDL Anode Catalyst Membrane Cathode Catalyst Cathode GDL Parameters D ij , ?, ?, ?, k P D i , ?, ?, ?, i o , ?, ?, k ? , k P D i , ?, ?, ?, k ? , k P D i , ?, ?, ?, i o , ?, ?, k ? , k P D ij , ?, ?, ?, k P Variables C H2 , ? s , P C H2 , ?s, ?, P C H2 , C O2 , ?, P C O2 , ?s, ?, P C O2 , x N2 , ? s , P Flux N H2 , i s , v N H2 , is, i, v N H2 , N O2 , i, v N O2 , is, i, v N O2 , N N2 , i s , v Conservatio n equations ??N H2 = 0 ? ? i s = 0 ? ? v = 0 ? ? N H2 =y H2 ? ? i s = -? ? i = j A ? ? v = 0 ? ? N H2 = 0 ? ? N O2 = 0 ? ? i = 0 ? ? v = 0 ? ? N O2 =y O2 ? ? i s = -? ? i = j C ? ? v = ? ? y w ? ? N O2 = 0 ? ? N N2 = 0 ? ? i s = 0 ? ? v ? ?. N w Charge Transport Equations Ohm?s law (is) Butler- Volmer (j A ), Ohm?s law (i s , i) Ohm?s law (i) Butler-Volmer (j C ) Ohm?s law (i s , i) Ohm?s law (i s ) Mass Transport Equations Stefan- Maxwell (H 2 , H 2 O), Henry?s law (H 2 , H 2 O) Nernst-Plank (H 2 ) Nernst- Plank (H 2 , O 2 ) Nernst-Plank (O 2 ) Stefan-Maxwell (O 2 , N 2 , H 2 O), Henry?s law (O 2 , N 2 , H 2 O) 14 Water Velocity Equations Darcy (v) Schl?gl (v) Schl?gl(v) Schl?gl(v) Darcy (v) Table 1-3. Summary of parameters, variables and equations by CFD analysis 1.3. MODELING OF PEM FUEL CELLS 1.3.1. GENERAL DESCRIPTION Numerical modeling and simulation of fuel cell systems play more and more important roles in designing and analyzing fuel cell stack and systems, whose major benefits can be summarized as follows: ? Preliminary investigations of different designs and new materials without actual construction and operation of demonstration equipment. For example, there are tremendous research efforts devoted to freezing and thawing studies. The actual investigation of the behavior of the stack takes time and resources which results in long development cycle and low efficiency. With proper modeling and simulation platform, many operating conditions can be tested and thus dramatically increase the possibility of finding durable and safe operations. ? Facilitate the development of control and diagnostic strategy for fuel cell system and fuel cell driven vehicles. It is anticipated that the fuel cell system calls for far different control algorithm and diagnostic capabilities mainly due to more difficult cooling and water demands than the traditional ICE engines. To avoid excessive deployment of expensive sensors and monitoring hardware, more efficient and sophisticated control strategies should be developed based on more realistic fuel cell system models, and 15 not just the empirical polarization map approach taken by current researchers. Distribution information can be extracted from the multilevel models of fuel cell stack system and provided to control decision components. With the continuously increasing computational power and more efficient numerical algorithms, an overall fuel cell simulation has become more feasible. The simulation can include detailed mechanistic fuel cell models and transient balance-of-plant models. Based on these models, a stack system design and optimization can be attempted. In order to fully realize the simulation potential, many crucial steps have to be taken. 1.3.2. STATE-OF-ART OF MODELS FOR LAYERS 1.3.2.1. MEMBRANE The membrane should block electrons and transfer protons from the anode to the cathode. Thus, understanding of the transport mechanism of the protons is of importance. The conductivity of the protons directly relates to the ohmic polarization. Thus, an improvement of the proton conductivity leads to an increase of the efficiency of fuel cells. Current models proposed can extensively and sufficiently describe the transportation phenomena of the protons. However, the models are using empirical parameters and data. Particular attention should be paid to experimental validation that currently presents challenging scientific issues; ? Springer et al. [1] proposed an isothermal, one-dimensional, steady state model for a partial cell including the gas diffusion layers and the membrane of Nafion 117, while the catalyst layers are treated as interfaces. The model is developed along the direction 16 perpendicular to the MEA. Even being empirical, the model can predict the water and proton transportation within the layer. They derived the water content of the membrane as a polynomial function of local vapor activity. Then the properties of the membrane including: the water diffusion coefficient and the proton conductivity of the membrane and the electro-osmotic coefficient are expressed as an empirical function of water content and temperature. Generally, this model can be utilized to investigate water transport in the membrane. However, the constraints are on the saturated membrane reflecting the flows by a single diffusion. The uniform concentration of water in the membrane and the corresponding gradient does not allow this hypothesis. Possible improvement is possible by using a two phases or hydraulic model. ? Bernardi and Vebrugge [2] proposed one-dimensional model for water transport within the membrane. The model treats the membrane as a two-phase system similar to micro porous medium, where gas and liquid channels are separate with constant porosity. By doing so, the convection of the water in the membrane caused by the pressure difference between the two electrodes is under consideration. ? Janssen [3] brought up a new model describing the transport through the membrane by using only gradients in chemical potential. Thus, this model can avoid a possible complication occurring when separating diffusive and pressure driven flow components. ? Weber [7, 11, and 19] developed a model which combined the two effects of both diffusion at the gas phase and the convection driven by capillary pressure together. The mass transportation and thermo physical properties of the membrane are derived on the basis of thermodynamic equations. As a result, these properties of membrane 17 are related to each other, which are mathematically and physically described. The calculated parameters are modified by using experimental data. 1.3.2.2. CATALYSTS Catalysts are primarily responsible for accelerate the speed of chemical reactions taking place. Complex structure of materials and chemical reaction processes requires fundamental understanding in their mechanism and phenomena. Platinum and carbon support as a basic material are currently modeled in two ways, an agglomerate or homogenous distribution. Most of the CFD modeling simply assumes the catalysts as a homogeneously distributed layer or frequently treats it as a dimensionless interface where only source terms for chemical reaction are considered. ? Springer and Gottesfeld [4] extended the cathode catalyst layer by using a pseudo- homogenous structure. They assumed that reactants are transported through a homogeneous mixture of gas pores, polymer, and carbon support. The proton transfer occurs through the polymer. And only an effective diffusion coefficient and local concentration are considered for the gas diffusion. ? Broka and Ekdunge [5] assumed the cathode catalyst layer as an agglomerate structure. Reactants are transported through gas pores in the catalyst. Upon arrival to the agglomerate, reactants are dissolved and diffuse through the surrounding polymer to reach to the reaction sites. It should be noted that the diffusive resistance of the polymer in the catalyst layer is not regarded as a limiting part for the diffusive flux of the gaseous reactants, but as the effective reaction rate. The results demonstrate the 18 capability of the agglomerate model in predicting the polarization curves until the mass transfer limited regime begins. ? Sui et al. [6] compared several types of models, a thin film structure and an agglomerate structure. The simulation with the former is conveyed under an assumption that reactant transport occurs only in the dissolved form within the polymer, i.e. the reactant gas is uniformly distributed in the pores of the catalyst layer. Although the agglomerate model is more coincidence with experimental results than the thin film one, enormous number of parameters necessary for the model adversely affects the model. ? Genevey et al. [8] presented the most comprehensive model of the cathode catalyst layer. The model is based on an agglomerate structure that includes the transport of thermal energy, gas species, liquid water, and charges. The dynamic behavior of various transport processes occurring within the catalyst layer can be predicted. The effects of the catalyst layer structural properties, porosity and reaction surface area on performance are also analyzed. ? Bultel et al. [9] proposed a one-dimensional steady state microscopic model of the catalyst layer. The dependence of geometry in the agglomerate model is analyzed. It is found that the thickness of Nafion film surrounding Carbon-Pt agglomerates is a limiting factor, particularly at medium and high current densities. However, the extension to a fuel cell model with the microscopic model seems to be very complex and impractical. A combination of macroscopic and microscopic models might be a potential solution to predict the PEMFC performance more accurately. 19 ? Wang [16] included a three-dimensional homogeneous model of the catalyst layer. The dependence of active catalyst sites are related to the liquid water concentration in the catalyst layer. By doing so the model can predict the influences of the flooding effects at the catalysts. 1.3.2.3. GAS DIFFUSION LAYER (GDL) Gas diffusion layer facilitates to transfer reactants into the catalysts and remove the products. Recently, new suggestions are proposed by adding ionomer into the carbon layer and/or embedding an additional micro-porous layer. Various models have been proposed to describe the mass transport in the layer: ? Wang, Wang and Chen [10] developed a homogeneous two phase model, where the two phases are regarded as a mixture. The Navier-Stokes equations are employed to describe the two phase flow the gas diffusion media. And the results are able to predict the liquid water concentration distribution in the layers of a fuel cell including the gas diffusion media. However, first, the interfacial dynamic mass transfer between the liquid and vapor phases is not considered, which means the interfacial mass transfer rate are regarded as infinity. Second, as no information about phase change can be derived from the calculation result, it is hard to consider the phase change into the energy conservation equation. ? Natarajan and Nguyen [12] developed a model for mass transport within the gas diffusion layer considering the two phase effect. The model describes both the gas transfer and the liquid water transfer by using the Darcy?s law. And the liquid water transfer is driven by hydraulic pressure. In their work the interfacial mass transfer 20 between vapor and liquid phases are considered. The capillary pressure is used to calculate the hydraulic pressure. And the relationship between capillary pressure and saturation is based on an interfacial drag coefficient which is assumed to be constant. ? Natarajan and Nguyen [13] proposed a pseudo three-dimensional model for water transport in the cathode GDL. First, the Stefan-Maxwell equations are used instead of the Darcy?s law, as a result their model considered the multi-component mass transfer effects. Second, the authors solve the two-dimensional model in sequence along the channel direction, so the model is able to predict some three dimensional effects. ? Meng and Wang [16] proposed a new model for water transport in the cathode GDL. A new submodel for liquid coverage at the GDL-channel interface included accounts for the water droplet emergence on the GDL surface. The interfacial model should be able to represent not only the two phase model developed to predict the cathode flooding effect on cell performance but also ultimately removes the inability of prior two phase models to correctly capture effects of gas velocity on cell performance. ? Wang and Wang [18] included thermal effects to predict water transport in the cathode GDL. However, as pointed above, the use of homogeneous two phase model cannot be used for a prediction on the phase change and consequently for inadequate description for the thermal effects. 1.3.3. SINGLE CELL MODEL The computer model for a single cell embraces gas channels, gas diffusion media and catalysts of both cathode and anode, and membrane. The complexity of the cell has been 21 attacked with helps of CFD technology to describe mass transport of protons and electrons transport. CFD technology can cover a 1-D, 2-D or 3-D domain. Challenging issues are models for components and the associated integration. Besides, exponential growth of computational time and validation of models that requires measurement and characterization of parameters and variables are another issues. In particular, the multi-phase phenomena coupled with thermal effects are yet to be solved. 1.3.3.1. 1-D MODEL ? Bernardi and Verbrugge [2] developed a one-dimensional mathematical model of a complete cell. The effect considered is perpendicular to the MEA. The model is capable of analyzing the reactant starvation and membrane dehydration. Due to the limitation of one-dimensional effects, it is hard to predict the water distribution and transportation in the GDL. Thus, this model is constrained on the range of low current density analysis because of the lacking liquid water transport. 1.3.3.2. 2-D MODEL ? Gurau, Liu, and Kakac [14] developed a two-dimensional mathematical model of a complete cell. Hence, no general form of the governing equation is used. Instead, governing equations are constructed based on layers of the fuel cell and computation is accomplished by dividing the space into three domains according to the spaces transformed. A new technique is employed to solve the equations. Firstly, the equations are separately solved and coupled to each other by an iterative technique. 22 The simulation conducted address three major issues that include the oxygen and water vapor distribution in the gas channel and gas diffusion layer. The results show the oxygen concentration drops along the channel and in the direction perpendicular to the MEA. It is noted how the temperature influences the I-V curve. A comparison of I-V characteristics at different temperature of 303K, 323K, and 353K are conveyed. However, the effect of liquid water is not under consideration. Thus, it is impossible to predict the coupled effects of the temperature and liquid water. Finally, the water velocity in the membrane layer is analyzed at three current densities, 850 A/cm 2 , 865 A/cm 2 and 900 A/cm 2 . According to the results presented, the direction of the water velocity changes when the current increases In addition, the velocity difference of the water in membrane along the gas channel mainly caused by the water accumulated in the reactant is also presented. However, the density changes due to species consumption appear to be neglected. Moreover, the connectivity of the main channel and diffusion layer involved a change of primary variables that may lead to numerical discontinuities under some operating conditions. Again, the width of the channel and portion of the gas diffusion layer hidden from the channel were neglected. ? Fuller and Newman [15] developed a two-dimensional model of a membrane electrode assembly of a fuel cell to study water and thermal management. This model includes water transport in the diffusion layer. While Wang [16] considers the mist flow, the plug flow is assumed here in the flow channel. In fact, the equilibrium sorption of water between gas phase and the polymer electrolyte depends on temperature, water, and thermal management. Particularly, the rate of heat removal was presented to be a critical parameter in the operation of the PEM fuel cells. 23 ? Nguyen and White [17] developed a two-dimensional model to predict the heat and water transport that accounted for variation in temperature and membrane hydration conditions along the flow channels. The model did not include the diffusion layer for water transport. The results show how ohmic loss in the membrane at high current densities affects a large fraction of the voltage loss in the cell. It is found that a minimization of the ohmic loss requires humidification of the anode gas. Additionally, the cathode gas must also be humidified, if pure oxygen is replaced by air. ? Musser and Wang [20] improved their model by adding heat generation by chemical reaction and others based on numerical 2-D heat transfer in a cell. The focus is placed on the analysis of thermal management of the PEM fuel cell stack, where the coolant outlet temperature influenced by the coolant flow rate is analyzed. It came to a conclusion that the outlet temperature insignificantly varies with increased coolant flow rates. In addition, the temperature of the stack is analyzed depending upon the stack geometry design, materials of the membrane and backing layers. The results show that the coolant at the inlet should have a better ability to reduce the higher temperature experienced at the anode inlet. Furthermore, Wang [21] improved the model by including the transport of gas species, momentum, and protons with the channel geometry. By treating the different terms in the governing equations as source terms, all of the equations can be uniformly solved over the entire domain, which makes developer of the codes easier to modify and parallelize on clusters of computers. Hence, the concentration discontinuity between the GDL and the catalyst layer is solved by using Henry?s law, which results from gaseous reactants dissolving in the polymer filled catalyst layer. The simulation is 24 assumed that no liquid water is present and the pores in the catalyst layers completely filled with polymer. One of the outcomes is the predictability of the overshoot of the output voltage caused by the oxygen concentration variation delay when load changes. The performance comparison between modeling results and experimental data are not considered in the mass transport limited regime. Wang [22] published a new model for a two-phase two-dimensional modeling by using a multiphase mixture model. The iteration needed for the calculation of the interface between the gas and liquid water is eliminated by mathematically manipulating models. Finally, the contour of the gas and liquid water concentration and velocity can be obtained. However, the proposed model lacks a validation that seems to be impossible with current resources. ? Siegel et al. [23] presented a two-dimensional single cell model including a catalyst model based on an agglomerate catalyst structure. Primarily, the authors focused on the fuel cell performance influenced by the catalysts properties. According to the studies, the optimization of catalyst layer porosity is required to maximize the amount of current produced at a given voltage. Larger the size of pores is, higher gets ohmic, while a smaller size leads to high mass transport losses. ? A model proposed by Singh, Lu and Djilali [24] includes the effects of liquid water transfer on the performance, but did not consider temperature influence that is improved in a following study [25]. The influences of the non-uniform temperature and pressure distributions on the liquid water and vapor fluxes in the diffusion layers are analyzed. Consequent humidification requirements are stricter than referred by isothermal and isobaric models. 25 1.3.3.3. 3-D MODEL ? Naseri-Neshat et al. [26] studied the reactants transported by both convective and diffusion transport apparatuses. Analysis shows that the geometry factors such as the width of the flow channel and gas diffusion layer are very influential on the reactants transport. In addition, varying effects in diffusion layer depending on diffusion layer thickness and channel cross-section area are also investigated. ? Shimpalee et al. [27] completed a three-dimensional straight channel model by establishing governing equations for each gas component at the inlet. It is noted that the energy equation can predict the temperature distribution in the cell and can be used to analyze thermal management of the fuel cell stack. However, only the heat source term caused by the chemical reaction among others is considered. In a revised paper [28], 20 channels with a serpentine flow path is modeled on the base of the previous model for a straight channel. The study reveals that the water evaporation and condensation produced by temperature directly change the membrane humidity, and finally influence the local current density. In a subsequent paper [29], Dutta, Shimpalee, and Van Zee used a three dimensional model to analyze the effect of a serpentine flow field on the performance. The model includes the transport of mass, momentum, gas species, and water within the polymer membrane. The catalyst layers are modeled as interfaces between the GDL and membrane, at which point various source terms are applied. The results show the velocity and density distribution within the gas channels and GDL. In addition, the variation of hydrogen consumption and water production across the surface of the MEA is 26 analyzed. The outcomes are the pressure drops caused by flow crossover from one channel to the next via the GDL and along the serpentine flow field that is less than expected for laminar flows through the gas channel alone. However, there is limitation in this model where liquid water transport is neglected. ? Um et al. [30] also developed a three-dimensional PEM fuel cell model. The model is an isothermal single phase. The focus is placed on the analysis of the reactant concentration distribution in a single cell. It is found that the force convection created by interdigitated flow field can improve mass transport of oxygen and water elimination at reaction zone, thus increasing the performance as compared to convectional flow field at a high current density. However, the calculation domain is a straight channel, where a flow turbulence occurring at a curvature of flow channels is neglected. Um et al. [31] discussed the influences of the flow patterns by conventional and interdigitated on the ohmic polarization losses in the MEA, current distribution and the oxygen concentration distribution. The results show that the performance of a cell with interdigitated gas channel is better than the conventional one, but does not allows for uniformity in the performance for different geometry. ? Ju, Meng, and Wang [32] proposed a three-dimensional non-isothermal model to rigorously account for various heat generation mechanisms, including irreversible heat because of electrochemical reactions, entropic heat, and Joule?s heat caused by electrolyte ionic resistances. The analyses show that the steady temperature distribution in the stack and an asymmetrical temperature distribution in the direction perpendicular to the MEA. 27 ? Zhou and Liu [33] coupled the governing equations for reactant mass transport and chemical reaction kinetics. The model includes the energy equation to map the temperature distribution inside a PEM fuel cell. However, the channel simulated is straight, which has 1.6 cm 2 reacting area, approximately. ? Berning, Lu and Djilali [34] demonstrated a new three-dimensional model with one gas channel, on the basis of their previous two-dimensional code. The analyses include three-dimensional effects on the reactant concentrations, current densities, temperature, and water fluxes. Significant temperature gradient existing within the cell is demonstrated by showing several degrees K within the MEA. However, no liquid water is considered. In a following work, [35], the model has been improved by including half of the coolant channel in order to analyze the temperature effects at two gas diffusion layers with the phase change. At the same time, the balance of three separate processes, temperature change, reactant gas depletion, and pressure drop inside the GDL are considered. The outcomes show that the relationship between the liquid water saturation and permeability are opposite for cathode and anode sides that is explained by the different sources of the liquid water for two electrodes. However, the catalyst layers are regarded as an interface, where the temperature influence on the water removal process and oxygen concentration in the catalyst layer is neglected. In a revised paper by Nguyen, Berning and Djilali [36], the previous model is further improved by adding a three-dimensional code for serpentine gas flow channels. A voltage to current algorithm is considered for a more realistic spatial variation of the electrochemical kinetics and a three-dimensional activity of the catalyst layer included. However, no liquid water is considered in this model. 28 1.3.4. STACK MODEL Safe and reliable operation of a fuel cell system requires a full understanding in design parameters with operating variables. Representing variables are humidity, temperature and pressure of reactants as well as the impurity. Currently, research and development of models for stacks have been rarely performed in comparison to endeavors devoted to development of models for a single cell. Some of publications are summarized as follows: ? Recently, Sundaresan [37] proposed a lumped model considering thermal effects, which facilitates to analyze influence of the endplate assembly on the startup procedure. However, no fluid dynamics are considered and all the chemical reactions are described by lumped equation, so that a realistic analysis is impossible. On the other hand, Wetton proposed a stack model that coupled the 1 D coolant channel [38] with two-dimensional MEA [39] and improved by electrical [40] and thermal [41] connection of the individual cells for a stack. However, no dynamics is considered. 1.3.5. SUMMARY Existing fuel cell models are insufficient to meet the industrial needs to understand the internal complex mechanisms of chemical reaction, water balance, reactants flows influenced by temperature fluctuation that occur in real world. According to intensive surveys conveyed, none of models are able to fulfill the needs for a model that represent the stack performance with dynamics. Thus, objectives of the research have been set to develop an integrated multilevel model capable of addressing different technological issues of materials and systems that encompass the layers such as the membrane, catalysts, flow field plates with balance-of-plant and associated controls of system. 29 Ultimately, the progress in modeling of components and system and its availability will also provide means to dramatically augment understandings in the mechanisms and physical phenomena of a cell, a stack and finally the power system that lead to advanced designs and/or diagnostic tools. 1.4. THESIS OBJECTIVES As discussed in the previous sections, most of current research focused on a cell, while the description for a stack is widely ignored. Thus, the objectives of this thesis are set to develop that considers reactants flows, chemical reactions and water balance coupled with temperature. Major benefits having developed a comprehensive PEM fuel cell stack models will be the availability of a tool that helps engineers and scientists understand the mechanisms and physical phenomena in operating environment and optimally design a fuel cell stack and system. 30 REFERENCES 1. T.E. Springer, T.A. Zawodzinski, S. Gottesfeld, (1991) Polymer Electrolyte Fuel Cell Model, J. Electrochem. Soc., Vol. 138, No. 8, pp. 2334-2341. 2. D.M. Bernardi, M.W. Verbrugge, (1992) A Mathematical Model of the Solid- Polymer-Electrolyte Fuel Cell, J. Electrochem. Soc., Vol. 139, No. 9, pp. 2477-2490. 3. G. J. M. Janssen, (2001) A Phenomenological Model of Water Transport in a Proton Exchange Membrane Fuel Cell, J. Electrochem. Soc., 148, pp. A1313-A1323. 4. T.E. Springer, S. Gottesfeld, Pseudohomogeneous Catalyst Layer Model for Polymer Electrolyte Fuel Cell, Proc. of the Symposium on Modeling of Batteries and Fuel Cells, R.E. White, M.W. Verbrugge, and J.F. Stockel, Editors, PV 91-10, pp. 197- 208, The Electrochemical Society Softbound Proceedings Series, Pennington, NJ, (b). 5. K. Broka, P. Ekdunge, (1997) Modelling the PEM fuel cell cathode, J. Applied Electrochemistry, Vol. 27, pp. 281-289. 6. P.C. Sui, L.D. Chen, J.P. Seaba, Y. Wariishi, (1999) Modeling and Optimization of a PEMFC Catalyst Layer, SAE Congress, 1999-01-0539, pp. 61-70. 7. A. Z. Weber, J. Newman, (2003) Transport in Polymer-Electrolyte Membranes I. Physical Model. Journal of the Electrochemical Society, 150, pp. A1008-A1015 8. D. Genevey, M.R. von Spakovsky, M.W. Ellis, D.J. Nelson, B. Olsommer, F. Topin, N. Montel, N.P. Siegel, (2002) Transient Model of the Heat, Mass and Charge Transfer as well as Electrochemistry in the Catalyst Layer of a PEMFC, International Mechanical Engineering Congress and Exposition ? IMECE?2002, ASME IMECE Paper No. 33322, N.Y., N.Y., November. 9. Y. Butel, P. Ozil, R. Durand, and D. Simonsson, (1995) Study of Mass Transfer within the Active Layer of P.E.M.F.C. Electrodes at the Particle Level, Electrochemical Society Proceedings, Vol. 95, pp. 34-47. 10. Z.H. Wang, C. Y. Wang and K. S. Chen, (2001) Two phase flow and transport in the air cathode of proton exchange membrane fuel cells, J. Power Sources, Vol. 94 (1), pp. 40-50. 11. A. Z. Weber, J. Newman, (2004) Transport in Polymer-Electrolyte Membranes. II. Mathematical Model. Journal of the Electrochemical Society, 151, pp. A311 12. D. Natarajan, T. Van Nguyen, (2001) A Two-Dimensional, Two-Phase, Multicomponent, Transient Model of the Cathode of a Proton Exchange Membrane Fuel Cell Using Conventional Gas Distributors, J. Electrochem. Soc., Vol. 148, No. 12, pp. A1324-A1335. 31 13. D. Natarajan, T.V. Nguyen, (2003) Three-dimensional effects of liquid water flooding in the cathode of a PEM fuel cell, J. Power Sources, Vol. 115, pp. 66-80. 14. V. Garua, H. Liu, and S. Kakac, (1998) Two-dimensional model for proton exchange membrane fuel cells, J. of AIChE, 44, pp. 2410-2422. 15. T. Fuller, and J. Newman, (1993) Water and thermal management in solid-polymer- electrolyte fuel cells, J. of Electrochemical Society, 140, pp. 1218-1225. 16. H. Meng and C.Y. Wang, (2005) New Model of Two-Phase Flow and Flooding Dynamics in Polymer Electrolyte Fuel Cells, Journal of Electrochemical Society, Vol. 152, pp. A1733-1741. 17. T. Nguyen, and R. White, (1993) A water and heat management model for proton- exchange-membrane fuel cells, J. of Electrochemical Society, 140, pp. 2178-2186. 18. Y. Wang, C. Y. Wang, (2006) A Non-isothermal, Two-phase Model for Polymer Electrolyte Fuel Cells, Journal of Electrochemical Society, Vol. 153, pp. A1193- 1200. 19. A. Z. Weber, J. Newman, (2004) Transport in Polymer-Electrolyte Membranes. III. Model Validation in a Simple Fuel-Cell Model. J. Electrochem. Soc., 151, pp. A326. 20. J. Musser, C.Y. Wang, (2000) Heat Transfer in a Fuel Cell Engine, Proceedings of NHTC'00, 34th National Heat Transfer Conference. 21. S. Um, C. Y. Wang, and K. S. Chen, (2000) Computational fluid dynamics of proton exchange membrane fuel cells, J. of Electrochemical Society, 147, pp. 4485-4493. 22. Z.H. Wang, C. Y. Wang and K. S. Chen, (2001) Two-phase flow and transport in the air cathode of proton exchange membrane fuel cells, J. Power Sources, Vol. 94 (1), pp. 40-50. 23. N. P. Siegel, M. W. Ellis, D. J. Nelson, M. R. von Spakovsky, (2003) Single domain PEMFC model based on agglomerate catalyst geometry, J. Power Sources, Vol. 115 pp. 81-89. 24. D. Singh, D.M. Lu, N. Djilali, (1999) A two-dimensional analysis of mass transport in proton exchange membrane fuel cells International Journal of Engineering Science, Vol 37 pp.431-452 25. N. Djilali, D. Lu, (2002) Influence of heat transfer on gas and water transport in fuel cells, Vol 41, pp. 29-40 32 26. H. Naseri-Neshat, S. Shimpalee, S. Dutta, W.K. Lee, and J. W. Van Zee, (1999) Predicting the effect of gas-flow channel spacing on current density in PEM fuel cells, Proc. of ASME IMECE, Nashville, TN., 4, pp. 341. 27. S. Shimpalee and S. Dutta. (2000) Numerical prediction of temperature distribution in PEM fuel cells. Numerical Heat Transfer, Part A, Vol. 38, pp. 111-128. 28. S. Shimpalee, S. Dutta, and J. W. Van Zee, (2000) Numerical prediction of local temperature and current density in a PEM fuel cell proceeding in ASME IMECE, Orlando, FL. November 5-10. 29. S. Dutta, S. Shimpalee, J.W. Van Zee, (2001) Numerical prediction of mass- exchange between cathode and anode channels in a PEM fuel cell, Int. J. Heat and Mass Transfer, Vol. 44, pp. 2029-2042. 30. S. Um, C.Y. Wang, (2000) Three-dimensional analysis of transport and reaction in proton exchange membrane fuel cells, Proc. of ASME IMECE, Orlando, FL., 1, pp. 19. 31. S. Um, C.Y. Wang, (2004) Three-dimensional analysis of transport and electrochemical reactions in polymer electrolyte fuel cells Journal of Power Sources Vol.125, pp 40-51 32. H. Ju, H. Meng, C.Y. Wang, (2005) A single-phase, non-isothermal model for PEM fuel cells, International Journal of Heat and Mass Transfer, Vol 48, pp 1301-1315 33. T. Zhou, and H. Liu, (2000), 3-D model of proton exchange membrane fuel cells, Proc. of ASME IMECE, Orlando, FL., 1, pp. 43 34. T. Berning, D. M. Lu, N. Djilali, (2002) Three-dimensional computational analysis of transport phenomena in a PEM fuel cell Journal of Power Sources, Vol. 106, Issues 1-2, pp. 284-294 35. T. Berning, N. Djilali, A 3D, Multiphase, (2003) Multi-component Model of the Cathode and Anode of a PEM Fuel Cell, Journal of The Electrochemical Society, 150 (12) A1589-A1598. 36. P.T. Nguyen, T. Berning, N. Djilali, (2004) Computational model of a PEM fuel cell with serpentine gas flow channels Journal of Power Sources, Volume 130, Issues 1-2, pp. 149-157 37. M. Sundaresan, R.M. Moore, (2005) Polymer electrolyte fuel cell stack thermal model to evaluate sub-freezing, Journal of Power Sources, Volume 145, Issue 2 pp. 534-545 33 38. B. Wetton, K. Promislow, A. ?aglar, (2004) A Simple Thermal Model of PEM Fuel Cell Stacks, Second International Conference on Fuel Cell Science, Engineering and Technology. 39. P. Berg, K. Promislow, J. St. Pierre, J. Stumper, B. Wetton, (2004) Water Management in PEM Fuel Cells J. Electrochem. Soc. Vol. 151, pp. A341. 40. G. S. Kim, J. St-Pierre, K. Promislow, B. Wetton, (2005) Electrical coupling in proton exchange membrane fuel cell, Journal of Power Sources, In Press. 41. K. Promislow, B. Wetton, A simple, (2005) mathematical model of thermal coupling in fuel cell Journal of Power Sources, In Press. 42. http://www.eere.energy.gov 43. http://www.eere.energy.gov/hydrogenandfuelcells/fuelcells/pdfs/fc_comparison_char t.pdf (July/01/2006) 34 CHAPTER 2: 1-D PEM FUEL CELL STACK MODELING 2.1. INTRODUCTION The PEM fuel cell is a potential candidate for use as an alternative power source in future vehicle and power conditioning applications. The electric power, as a load, continuously varies with a function of time in the systems. Accordingly, a corresponding chemical reaction occurs in the cell and water and heat as byproducts are produced, Transfer of reactants and water plays a significant role in the performance of a cell. Particularly, temperature variation in the individual layer of the cell strongly influences reactants transfer to the catalysts, chemical reaction and water content in the membrane, so the representing I-V polarization characteristic of a cell varies dynamically. The complexity of these interrelated physical phenomena impedes design and system engineers from optimizing design parameters and criterion that lead to high costs and engineering time. Utilization of computer models have been proven to be the most effective tool requiring a comprehensive physical model reflecting real behaviors of the fuel cell stack. Currently, models can be classified into two groups; a simple composition of passive electric circuit representing the double layer effects with a voltage source using I-V empirical curve [25, 26], and a relatively precise cell model based on CFD [27, 28, 29] whose behavior is governed by equations describing chemical reactions, fluid and 35 thermodynamic, and electric properties. These CFD models take enormous computational time necessary for a calculation of grid points and constrain the wide use of the models. Moreover, associated parameters and variables required for the models are hard to characterize and measure because of the thinness of layers during operations. Thus, these models are mostly applied to investigate parts of complex domains such as flow fields of a single cell and predict the cell performance with a constant temperature. A possible loophole has been proposed by other authors with lumped models emphasizing dynamic behavior of the stack. For example, Amphlett [2, 30] developed the first dynamic model with empirical thermal values, and Gurski [3] expanded it to consider reactant flow and coolant control. Further improvement was made by Muller [31], who calculated the temperature variation of a stack. Others proposed models calculating the temperature variation of the stack or the cell as a control volume [4, 5, 6, 7, 8, 9, 10], or two electrodes and MEA [11, 12]. However, the simulation results do not incorporate either the dynamic or transient aspects of the fuel cell system in operating environments. These drawbacks have been significantly complemented by using simplified thermodynamic models to analyze the performance of the stack, by Sundaresan [14, 32]. His model based on layers is used to analyze the start-up behavior from a sub-freezing temperature. However, the model does not fully consider several factors; 1) flow of species at the inlet must be the same as that at the outlet. Thus, no fluid dynamics are considered in the model; 2) heat source terms in both the catalysts are empirically calculated with values suggested by Wohr and Peinecke?s modeling [6, 33]. Accordingly, the anode source term is presumed as a relatively large value that should be referred to be around zero [17]. As a result, the model does not show asymmetric 36 phenomena of performance across the stack. Wetton et al. [13] proposed an explicit thermal model with the coolant channel coupled with a 1-D cell model [34], which shows an outstanding temperature gradient of the stack, but with no dynamics at all. In fact, dehydration of membranes dynamically varies with temperature, which strongly influences overall performance of cells and finally the resulting stack. Consequently, a high dynamic model for a stack based on layers of a single cell [35] has been developed by reflecting the dynamics under different operating conditions for pressure, humidity, reactants and temperature. Temperature profiles in a stack with associated effects on the stack performance show the high dependency in not only static but also dynamic behavior. On the other hand, the proton transport in the membrane and its associated ohmic losses mainly determine the characteristics of the ohmic polarization. The proton conductivity has been regarded as constant, temperature dependent [1], or temperature and water concentration dependent variables [15]. Recently, Pukrushpan [16] proposed the most comprehensive model that considers the dependence of the proton conductivity on the water concentration and temperature. However, the water concentration of the membrane obtained from the membrane relative humidity (RH) is based on an average of the anode and cathode RH. In fact, the water concentration in the anode and cathode varies rapidly for percentage, while the water concentration in the membrane does slowly because the amount of water residing in both sides is relatively less than that in the membrane [15]. In addition, the oxygen concentration in the GDL on the cathode side is continuously changing in operating environments and significantly affects the performance of the cell. Therefore, plenty of models considering multiphase multi-species have been employed to 37 investigate the transport phenomena in the GDL. However, those models do not consider the dynamics. Recently, Pukrushpan proposed a dynamic model with lumped parameters to predict the gas dynamics in a cathode electrode, which does not consider the effects in the GDL [16]. Therefore, a 1-D single phase model is added in this work to represent the dynamics present in the GDL. 2.2. 1-D SINGLE CELL MODELING 2.2.1. MODEL DESCRIPTION 2.2.1.1. MODELING DOMAIN AND ASSUMPTIONS The model has been developed on the basis of layers in a cell that consist of a MEA, two gas diffusion layers, and two gas channels sandwiched by two coolant channels, as shown in Fig. 2-1. The input variables for the model are current load, mass flow rate, the gas components fraction, temperature, pressure and relative humidity of reactants as well as the temperature and velocity of coolants at the inlets. Fig. 2-1 Schematic simulation domain 38 The main assumptions made for the new model are as follows: y Reactants are ideal gases; y There is no pressure gradient between the anode and cathode side; it means no convection but only diffusion for gas transport is considered; y There is no gas pressure drop from the inlet to the outlet of the gas channel; y The temperature gradient is linear across the layers in a fuel cell; y The thermal conductivity for the materials in a fuel cell is constant; y There is no contact resistance; y Anodic over-potential is negligible; y There is no current density gradient across the cathode catalyst layer; it means that the reactants completely reacted as soon as it reaches to the cathode catalyst layer surface. Based on these assumptions, five sub models have been developed and are described in the following sections. 2.2.1.2. MATHEMATICAL MODELING 2.2.1.2.1. ELECTROCHEMICAL MODEL Generally, the overall chemical reaction of the PEM fuel cell can be described by using the following expressions, illustrating that a chemical reaction of hydrogen and oxygen molecules produces electricity, water, and heat as a byproduct: 222 1 2 res cell HOHOQV+?++ (1) 39 The output cell voltage V cell is the difference between the open circuit voltage (OCV) E 0 and over-potentials ? and V ohm . cell OC ohm VV V?=?? (2) By neglecting the dependence of the OCV on the reactant pressure, the relationship between the OCV and the temperature can be simplified with the empirical parameter dE 0 /dT. If the reactant is ideal, its activity can be described by using the equation (3), where index i indicates H 2 and O 2 , while p i is the partial pressure of gas components, and p 0 is the overall pressure of both the anode and cathode side. Then, E 0 can be derived by modifying the Nernst equation (4). 0 i i p a p = (3) () () 22 10.5 ln 2 oc oc ref ref H O dV RT VV TT aa dT F =+ ?+ (4) The anodic over-potential is negligible; while the ? represents the over-potential of the cathode catalyst layer. Under the further assumptions that the asymmetric parameter of the reaction is 1, the Butler-Volmer equation leads to equation (5) that describes the over- potential on the cathode side. 2 2 , 0 , exp 1 Ocata eff ref cell O ref ref Hp A F jj Ap RTH ? + + ?? ?? ???? ?? =? ?? ?? ?? ?? ?? ?? (5) The ohmic over-potential V ohm is determined by the product of the current density and the proton resistance R mem according to Ohm?s law (6). ohm mem ViR= (6) 40 2.2.1.2.2. THERMAL MODEL If a cell is assembled with cubic layers whose thermo-physical properties are isotropic and constant, then according to the energy conservation equation, the total energy changes in a controlled volume equals the sum of the energy exchange at boundaries and internal energy resources. In fact, the energy exchanges at boundaries occur by three factors: a) the mass flow into each volume; b) the conduction heat transfer across the cell; c) the convection heat transfer occurring between bipolar plates with the coolant and the reactants. Thus, the thermal-dynamic behavior can be described with the following energy conservation equation: () N , cv i i mass cell cv in cell j in cv Conv cell Cond cell T i Sources Convection heat transfer Conduction heat transfer Mass flow in dT Cp c A Th m A Cp T T Q A Q A S dt =?+++ ??      (7) On the other hand, the internal energy source is composed of the entropy loss and the chemical energy required for protons to overcome the barrier of the over-potentials in both catalyst layers (Eq. 8). In addition, other heat sources are ohmic losses caused by a transport of electrons and protons in the cell. : 4 T cell cell ohm TS SiA iAR F ? ??? =?++ ?? ?? (8) In fact, the change of entropy due to the electrochemical reaction (Eq. 9 and 10) in both of the catalysts sides predominantly influences the energy sources term according to the calculation shown below. ?+ +? )()()(2 22 ptaqg eHH (9) )(2)()()(2 2221 lptaqg OHeHO ?++ ?+ (10) 41 In order to obtain the entropy change of these reactions, the zero point of semi-absolute entropy is taken as a reference according to [17]: () 0 aq sH + ?? ? ?? (11) The entropy of an electron obtained from the standard hydrogen electrode results in the following equations [17]: ( ) 0 SHE SHE SHE SHGT?=??? = (12) () () 11 2 1 65.29 2 pt g se sH Jmol K ??? ?? ? ? == ?? ? ? (13) Therefore, the entropy change of the cathode reaction is equal to the sum of that in water, oxygen and electron: () () 22 11 1 2 2 1 69.91 *205.03 2*65.29 163 2 ca ptlg SsHO sO se Jmol K ? ? ? ???? ???= ? ? ?? ???? =? ? =? (14) If the anode is assumed as a standard electrode, the anodic entropy change becomes 0. 2.2.1.2.3. PROTON CONDUCTING MODEL FOR MEMBRANE The membrane resistance is a function of the temperature and water content in a membrane layer, which is described as follows [16]: () 23 11 12 2 11 exp 303 mem mem HOSO mem Th R bbb T ? = ?? ?? ?? ?? (15) Where, the temperature T mem can be derived from the previous equation (7), while the membrane water content mem ? can be described by using the water mass concentration [15] and the mass conservation equation [16]: 42 22 23 22 , , , HOmass HO HOSO dry mem H Omass HO mem cM bc M M ? ? = ? (16) () () 2 , ,, ,, ,, ,, H O mass cell mem water ele mem a ele mem c diff mem c diff mem a dC A dm WWW W dt dt ==?++ (17) The electro-osmotic driving force created by the different electrochemical potential at the anode and cathode determines the water mass flows of W ele,mem.a and W ele,mem,c. at the boundaries of the membrane layer . In addition, the diffusion caused by the water concentration gradient at the two boundaries makes up the mass flows of W diff,mem,a and W diff,mem,c . Those relationships are described by equations (18), (19), and (20), proposed by Springer [24]. 23 23 219 0.0029 0.05 3.4*10 dHOSHOS n ?? ? =+? (18) ,, ,ele mem i water cell d i i WMAn F = (19) () ,, imid diff mem i water cell water mem CC WMAD l ? = (20) Hence, the diffusion coefficient D water and the water concentration C i is calculated from the empirical equations [24] expressed as a function of membrane water content 23 H OSO ? : () 23 11 exp 2416 303 water H O SO mem DD T ? ?? ?? =? ?? ?? ?? (21) () () 23 23 23 23 23 23 23 6 6 6 6 10 2 10 1 2 3 32 () 4.5 3 10 3 1.67 3 4.5 1.25*10 HOSO HOSO HOSO HOSO HOSO HOSO HOSO D ? ? ? ? ? ? ? ? ? ? ? ? > ? +?? ?? ? = ? >> ? ?? ? ? ? ? (22) 43 The boundary water content i ? is a function of water activity a i , which is calculated from the water vapor partial pressure: () 23 0.043 17.81 39.85 36 1 0 14 1.4 1 3 1 16.8 3 iiii ii i aaaa a ? ? +? + ?> ? =+? ? ? ? ? (23) , , vi i sat i p a p = (24) 2.2.1.2.4. PROTON CONDUCTING MODEL IN THE CATHODE CATALYSTS The dynamic behavior of a fuel cell at a load is investigated by experiments. When the output current changes abruptly, the output voltage of the fuel cell reacts with an overshoot [18]. These dynamics result from different physical phenomena of reactants and their chemical reaction in the cell, such as dynamics filling in the gas flow channel, diffusing reactants through the GDL and reacting process in the double layer at the interface of electrodes. M. Ceraolo explained the dynamic effects with a relationship between the number of mobile protons and water content [1]. As the matter of fact, when the current density increases, the hydration of the polymeric electrolyte near the cathode catalyst tends to raise as well; consequently, the proton concentration near the cathode catalyst increases rapidly. On the other hand, the proton concentration will decrease slowly at a decrease of current. Thus, the dynamics can be described by the following differential equation using the proton concentration as a variable [1]: 3 1 HHHH HH I ccc tt ? ? ?? +++ + ++ +? ???? ??+= ?? ?? ??  (25) 44 H ref cHH + ++ ????= ????  is the dimensionless proton concentration, ()? is the Heaviside function, and H ? + and H ? + are empirical parameters. Fig. 2-2 shows the calculated responses. The voltage decreases quickly when the current density increases. However, the voltage first reaches its highest value and then damps with a time constant that is associated with the proton concentration, as the current density decreases. Fig. 2-2: Voltage response by a consideration of proton concentration 2.2.1.2.5. GDL REACTANT MODEL FOR CATHODE Air contains not only oxygen but also nitrogen and water vapor. The air entering the cell diffuses through the GDL before reaching the catalyst layer. The diffusion effect is described by using the mass conservation (26) and Stefan-Maxwell equations (27): 0 g ii pN RT t y ? ?? += ?? (26) () 3 2 1 g i ik ki k cik p RT p NpN ypD ? ? = ? =? ? ? (27) 45 Hence, i, k = (1, 2, 3), where P 1 is the oxygen partial pressure, and P 2 = P sat (T) and P 3 are the water vapor and the nitrogen partial pressure, respectively. The diffusion coefficient p c D ik = D ik,eff = D ik,eff (T), and the cathode pressure of P c is the summation of the species partial pressures. The parameter ? is a constant describing the pore curvature of the GDL. The partial differential equation (PDE) systems above can further be simplified by using the following PDE [1], whereby ? is the dimensionless distance y/Th gdl : 22 2 2 2 4 OO O p pp i tF ?? ? ? ?? ? =? ?? ? (28) ()( ) 22 12 12 1 gdl sat ca sat Th p D p p D ? ? = +? (29) () g gdl ca sat RT Th p p ? ? = ? (30) In order to investigate the effects of the GDL on dynamics, simulations have been run to analyze the relationship between the GDL thickness and the dynamics of oxygen by diffusion. Fig. 2-4 shows the dynamic response of the oxygen partial pressure at the interface of the cathode GDL and cathode catalyst layer. The results show that the oxygen partial pressure drops rapidly when a step current (Fig. 2-3) is applied. Thus, the concentration over-potential increases. Accordingly, the thickness of the GDL is one of the factors influencing the dynamic response. Moreover, the steady state is reached by the thin GDL more quickly than by the thick one. Fig. 2-3 Current input 46 Fig. 2-4 Oxygen concentration responses dependent upon GDL thickness At P=1 bar and T=353K Further analysis has been undertaken to discover how the reactant partial pressure is distributed along the GDL by the given pressure (Fig. 2-5) and current step. Fig. 2-6 shows the cathode oxygen partial pressures and their responses depending upon the thickness ratio. The analysis shows that the dynamic response of the oxygen partial pressure is highly dependent upon the geometrical locations. When the cathode inlet pressure changes, the pressure at the catalysts side responds with a time delay before it has reached the steady state, which is caused by the diffusion of the reactant. Accordingly, the over-potential can not be manipulated instantly. Fig. 2-5: dynamic cathode pressure input 47 Fig. 2-6: Oxygen concentration response to the dynamic current input and cathode pressure at different depths The dynamic responses of the oxygen concentration at the catalyst layer are illustrated in Fig. 2-7. The oxygen concentration is strongly influenced by the thickness of the GDL. The thinner the layer becomes, the shorter the response time gets. On the other hand, when the inlet pressure increases (Fig. 2-5), the partial pressure at the catalysts tends to follow its increase, but the amounts of the recovered partial pressure compared to Fig. 2- 4 depend on the thickness. Therefore, the settle times to the steady state become constant regardless of the thicknesses of the GDL. Fig. 2-7: Oxygen concentration dynamic response for different GDL thicknesses In this part, the effects of material properties on the dynamics are analyzed by a given current step. Fig. 2-8 and 2-9 show the dependences on porosity and tortuosity. Generally, 48 the GDL with a high porosity allows less pressure drops than that with the lower one. However, the GDL with a high tortuosity causes higher pressure drops than the low one because of a long path for the oxygen transported. When a step current is applied to the cell, the oxygen consumption on the catalysts side will be increasing. Instantly, the high porosity enables more oxygen to transfer from the inlet to the catalysts side, and, subsequently, the oxygen partial pressure at the catalysts can quickly follows the pressure changes at the inlet. Thus, the GDL with a high porosity dynamically responds to the pressure increase. When the tortuosity increases, the dynamic response time slows from the same effects as prolonged geometry and the associated pressure drop. The settle times remain unchanged, as in the analysis for different thicknesses. Fig. 2-8: Dynamic response of oxygen concentration for different porosities Fig. 2-9: Dynamic response of oxygen concentration for different tortuosities 49 2.2.2. SIMULATION SETUP Multi-run simulations have been conducted to investigate the static and dynamic behavior of a single cell. The static behavior is analyzed by calculating the typical polarization, which takes into account variables such as, working temperature, pressure, stoichiometric number, and relative humidity (RH) at the cathode inlet. The dynamic characteristic considers two aspects, the startup and the transient response on the current as a step load. The parameters and reference data for the models chosen are shown in Table 2-1, and some of them are empirical [1, 16, 21]. Electrochemical Reaction Model Proton Conducting Model P 0 1.0 bar b 11 0.5139 [16] T ref 343.15 K b 12 0.326 [16] E ref 0.975 V [1] b 2 350 [16] dE 0 /dT 0.00027 V K -1 [1] b 0.0126 [16] A cata,eff /A cell f(i, T, P o2 )[1] n d f(C water ) [16] Gas Transport Model D W f(T, C water ) [16] D eff f(P, T) m 2 s -1 [1] Proton Concentration Model P sat f(T) bar [1] Thermal Model H ? + 5.87E-12 (m 2 A -1 ) 3 [1] H gas f(P, T) [21] C p,gas f(P, T) [21] H ? + 12.78 s [1] ? gas f(P, T) [21] Empirical Parameters Table 2-1. Empirical parameters for quasi 1-D models 2.2.3. ANALYSES 2.2.3.1. STATIC BEHAVIOR Fig. 2-10 shows the temperature dependent I-V characteristics from 333?K to 363?K with a step of 10K. As the temperature increases, the water removal will be more eased. The effects are considerably high at the range of the higher cell current, where more water is 50 produced. The increased of the voltage by increased working temperature shown in the results is approved by the CFD analysis [19]. Fig. 2-10: IV curves for different cell working temperature Cell: P=1.0bar Fig. 2-11 shows the pressure dependent polarization curve. As the pressure increases, the oxygen concentration at the surface of the catalysts layer tends to increase, too. Thus, the concentration over-potential gets lower. Otherwise, the over-potential becomes higher because of the oxygen starvation. The increased of the voltage by increased gas channel pressure shown in the results is approved by with the CFD analysis [19]. Fig. 2-11: IV curves for different cell gas pressure Cell: T=353.15K 51 Fig. 2-12 shows the stoichiometric number dependent I-V curve at a constant temperature and pressure. When the stoichiometric number is low, the removal of water at the cathode outlet flow decreases. Thus, the water concentration in the membrane layer increases. Consequently, the membrane resistance and the resulting ohmic over-potential become lower. However, the low stoichiometric number adversely affects the cathode over- potential at the high current because of the excessive water in the catalysts. Fig. 2-12: I-V curves for different stoichiometric number Cell: T=353.15K, P=1.0bar Fig. 2-13 shows the output voltage of a cell influenced by relative humidity at the cathode inlet. When the humidity increases at the cathode side, the air transported to the catalysts will be blocked, and the cathode over-potential will increase, especially at the high current range. The increased of the voltage by decreased inlet RH shown in the results is approved by the experimental data [20]. 52 Fig. 2-13: IV curves for different RH at the cathode inlet Cell: T=353.15K, P=1.0bar 2.2.3.2. DYNAMIC BEHAVIOR For a start-up of dynamic simulations, initial values are necessary for variables such as layer temperature, membrane water concentration, GDL air and oxygen concentration, and gas channel pressure. Fig 2-14 shows the setup for the simulation of a single cell from layer 1 to 11. Coolant Plate Gas Channel Gas Channel Plate Coolant Channel GDL Catalysts Membrane Catalysts GDL Channel layer 1 layer 2 layer 3 layer 4 layer 5 layer 6 layer 7 layer 8 layer 9 layer 10 layer 11 Cathode Bipolar Plate Anode Bipolar Plate Fig. 2-14 Simulation setup Geometrical and thermo-physical data for the layers are summarized in Table 2-2. Thickness m Heat conductivity W m -1 K -1 Heat capacity J kg -1 K -1 Density Kg m -3 GDL 0.0004 65 840 2000 Catalyst Layer 0.000065 0.2 770 387 53 Membrane Layer 0.000183 0.21 1100 1967 Gas Channel 0.001 52 935 1400 Plate 0.001 52 935 1400 Coolant Channel 0.001 30 935 1400 GDL Porosity 0.5 GDL Tortuosity 3.725 Bipolar Plate Contact Area Percentage 0.55 Membrane Molecular Mass 1.1 kg mol -1 Fuel Cell Area 0.0367 m 2 Fuel Cell Active Area 0.03 m 2 Table 2-2. Geometry parameters for the fuel cell 2.2.3.2.1. STARTUP The start-up temperature for the cell model is initially set to 298.15?K. The value of current density is increased continuously for the first 350 seconds in order to quickly raise the temperature to 353.15?K, which is assumed as a typical working temperature. In addition, a temperature controller is built in the simulation, as if a coolant subsystem turned on and off at this set point to extract the excessive heat produced in the cell. Fig. 2-15 shows the dynamic behavior of the temperature for different layers and voltage, as well as the efficiency during a start-up. It took 8 minutes for the cell to reach to the working temperature (Fig. 2-15a). Generally, the temperature profiles in each of the layers tend to follow the current waveform because of the associated energy losses occurring in the layers. Particularly, the temperature in the membrane and catalysts layer is the highest, which results from ohmic losses due to the membrane resistance and the heat released by the chemical reaction. The average difference of temperature between these two layers and other layers on the anode and cathode side amounts to 3?K and 2?K, respectively. Corresponding voltage and power are shown in Fig. 2-15c. When the 54 current increases, the over-potentials increase, and, subsequently, the voltage and power decrease. While the temperature is rising, the voltage fluctuates slightly during the start- up. When the current had been kept constant after the 350 seconds, the amount of water generated in the catalyst layer becomes constant. On the other hand, the continuously increased temperature leads to a high saturation pressure in the cell, which enables water residing in the catalysts to be quickly removed [23]. Otherwise, water would be flooding and blocking further influx of the oxygen into the catalysts. Therefore, the cell voltage increases rapidly. Thereafter, the water concentration in the membrane continuously decreases by the electro-osmotic force and diffusion effects, and the corresponding proton conductivity will be decreased. Thus, the cell voltage slightly drops after the temperature has reached a steady state. The overall efficiency of a cell is also strongly influenced by the variation of temperature (Fig. 2-15d). a) Dynamic current 55 b) Dynamic temperature response of different layers of single cell c) Dynamic voltage/power response d) Single cell efficiency Fig. 2-15 Simulation results 56 Fig. 2-16 shows the dynamic behavior of the temperature distribution across the fuel cell at 50 seconds and 7 minutes during a startup process. The catalyst on the cathode side shows the highest peak because of the losses associated with the over-potential being higher than on the anode side. For the first 30 seconds, the temperature rises slowly because of the slow slope of the current. The maximum difference of temperature between the GDL at the anode side and the membrane has been shown to reach 1?K. The higher the current is; the more heat is generated. Thus, the rise of temperature in the layers is accelerated. When the start-up is ended after 7 minutes, the temperature in the catalyst rises up to 353?K, approximately 2 ?K higher than the anode side catalyst. Then, the peak point of the temperature is moved from the catalyst to the membrane, which results from the dehydration of the membrane and the associated increase of losses. The dehydration is mainly caused by diffusions of water from the membrane to both sides because of higher water concentration in the membrane than in the gas channel sides. On the other hand, the increased number of protons transported takes up more water from the membrane to the cathode. Consequently, the resistance in the membrane is increased and shows the highest temperature among the layers. a) 57 b) Fig. 2-16 Temperature gradient across the cell (Left to right: cathode coolant channel to anode coolant channel) a) 50 seconds b) 7 minutes after startup 2.2.3.2.2. TRANSIENT RESPONSE In order to analyze the dynamic response on a power demand, a step current with 0.8 A cm -2 to 0.4 A cm -2 at the 600 second is applied. Fig. 2-17 shows the response of the temperature in the different layers (b), voltage (c), power output (c), and efficiency (d). The operating temperature is automatically controlled by the coolant system, the reference value for which is set to 353.15?K. The on-off control of the coolants causes a slight fluctuation of the temperature waveforms until 600 seconds. When the current suddenly decreases, the heat generated at the cathode catalyst and membrane layer decreases rapidly and leads to a temperature drop at these two layers. Then, the coolant system is turned off. The heat is transferred by the temperature gradient from the layers into the bipolar plates and stored there. Thus, the temperature of both bipolar plates tends to increase, and the temperature gradient begins to decrease. As a result, the amount of heat removed from the catalyst and membrane layers is again decreased. 58 When the temperature of the catalyst and the membrane layer reaches its lowest point, the temperature of all the layers will rise again due to the accumulated heat after the coolant is turned off. Finally, it reaches the steady state after around 10 seconds. The effects of the temperature variation on the output voltage are slightly different from the temperature profiles. When the current decreases, the voltage rapidly increases and shows a slight overshoot. Then, the voltage slowly elevates and drops back to a steady state. The first overshoot of the voltage results from the variation of the proton concentration in the cathode catalyst layer, while the first half of the voltage arc is caused by the temperature increase in all the layers and the rest by the decrease of water concentration in the membrane layer. The efficiency profile is the same as the behaviors of the output voltage, but the power remains almost the same as before the current was applied. The curve of the voltage response is comparable to an experimental result [22]. a) Current 59 b) Temperature in the different layers c) Voltage and power d) Efficiency analysis Fig. 2-17 Analysis of transient behavior of temperature, voltage, power, and efficiency upon a current step 60 2.3. 1-D STACK MODELING 2.3.1. MODEL DESCRIPTION 2.3.1.1. MODELING DOMAIN AND ASSUMPTIONS Assumptions made are as follows: y Reactants as ideal gases; y No pressure gradient between the anodic and cathodic side, which results in no convection, but only diffusion for gas transport; y Identical inlet conditions of each cell for both the cathode and anode as well as coolant channel; y No gas pressure drop along the gas channel; y Linear temperature gradient, linear across the layers in the stack; y Constant thermal conductivity of the materials in a fuel cell; y No contact resistance; y Negligible anodic over-potential; y No current density gradient across the cathode catalyst layer, which implies a complete reaction of reactants at the cathode catalyst layer surface. y The anodic stoichiometrical coefficient is 1, which indicates a complete reaction of the fuel filled in the fuel cell stack. Based on these assumptions, a layered based model for a 10-cells stack is developed. The schematic configuration of the simulation domain for the stack is shown in Fig. 2-18, where the cathode sides of the cells are located on the left hand side, while the anodic on the right hand side. 61 Fig. 2-18: Stack schematic configuration 2.3.1.2. MATHEMATICAL MODELING 2.3.1.2.1. SINGLE CELL MODEL A model for a single cell has been developed with compositions of individual layers and the performance detailed analyzed [35]. Table 2-3 summarizes the sub- models the governing the equations describing phenomenon for a single cell. Sub models Phenomenon Mathematical Equation Electrochemical model Chemical reaction Butler-Volmer Layered thermal model Temperature variation Energy conservation Membrane water balance Water transport 1-D mass conservation Catalyst layer proton model Proton concentration variation Empirical GDL O 2 diffusion model Multi-component diffusion 1-D Stefan-Maxwell Gas channel mass balance 0-D mass conservation Table 2-3. The governing the equations describing phenomenon for a single cell 62 2.3.1.2.2. STACK MODEL The model for a stack consists of models for endplates assembly, multiple single cells and coolant channels as in a typically designed stack. The stack voltage and current are obtained by series and parallel connections of outputs of individual cells, while the thermal conditions for each of the cells and coolant channels, or between coolant channels and endplate assembly, are coupled by using the energy conservation equation (Eq.31): () N , cv i i mass cell cv in cell j in cv Conv cell Cond cell T i Sources Convection heat transfer Conduction heat transfer Mass flow in dT Cp c A l m A Cp T T Q A Q A S dt =?+++ ??      (31) As a part of the source term, the change in entropy due to the electrochemical reaction in both of the catalysts sides predominantly influences the energy sources term, which can be clarified by the calculations below. The chemical reaction of the cathode catalysts can be described as, )(2)()()(2 2221 lptaqg OHeHO ?++ ?+ (32) Therefore, the entropy change of the cathode reaction is equal to the sum of that in water, oxygen and electron. The phase of water produced can be either vapor or liquid which makes a difference in the entropy change of the reaction. The value of the entropy change indirectly measured by [1] amounts to 52 J mol -1 K -1 . However, this value is quite different from other authors [27, 36, and 37]. Therefore, the energy generated in a catalyst layer for three cases are calculated and the values are compared for a better judgment. 63 The first case assumes water is generated as vapor phase. Then, the cathode entropy change can be described as follows: () () 22 11 1 2 2 1 198 *205.03 2*65.29 35.095 2 OC cptgg V nF S s H O s O s e T Jmol K ? ? ? ? ???? ? ?=? = ? ? ? ? ???? ? =? ? =? (33) The resulting energy by the entropy change is equal to: 11 35.095 c QST TJmolK ? ? =? =? (34) The second case assumes water as liquid phase. The corresponding entropy change of the reaction and the energy result in as follows: () () 22 11 1 2 2 1 69.91 *205.03 2*65.29 163 2 OC cptlg C nF S s H O s O s e T Jmol K ? ? ? ? ???? ? ?=? = ? ? ? ? ???? ? =? ? =? (35) 11 163 c QST TJmolK ?? =? =? (36) The last case assumes that the water generated is immediately evaporated at the catalysts. The corresponding energy required for the evaporation can be expressed as a product of the latent heat and the rate of mass change from liquid phase to vapor: fgfg Qhm=  (37) A comparison with the experimental value [1] indicates that most parts of water generated might be a state of vapor. In the following part of this section, simulations have been conducted for the first two cases. In contrast, the anodic side is assumed to be a standard electrode, while the membrane is assumed to provide a water-like environment. Thus, the anodic entropy change becomes almost zero. 64 ?+ +? )()()(2 22 ptaqg eHH (38) 2.3.2. SIMULATION SETUP All of the aforementioned models are coded with Matlab/Simulink and C. Simulations with different operating conditions have been conveyed to investigate the static and dynamic behavior of a single cell as well as a stack, which include the typical polarization curves at different working temperatures, a start-up behavior and transient responses of the stack on a current as a step load. The parameters and reference data for the models chosen are as shown in Table 2-1, and some of them are empirical. Currently, lack of experimental data on the water phase inside cells impedes any statement on how much of the enthalpy is generated by a water molecule, which includes decisive information on the state of water generated. Thus, the two cases of either the liquid or the vapor are selected and summarized as follows: Case 1: If the RH of the cathode inlet is low and at the same time the stoichiometric number on the cathode is high, water generated by the chemical reaction in the catalysts is immediately removed by the cathode outlet gas in the form of vapor. This state can be regarded as a vapor phase, where no liquid is formed in the stack. One of possible operating conditions is listed in Table 2-4 with other parameters necessary for a simulation: Case 2: If the RH at the anodic inlet is high and stoichiometrical number at the cathode is low, then water generated at the catalysts cannot be removed immediately in the form of vapor in the cathode outlet gas. Thus, an assumption can be made that water is mainly 65 generated and removed in the liquid phase. These cases are listed in Table 2-5 with the second condition. The simulation parameters are listed below. The start-up and transient behaviors are analyzed by using these conditions and the simulated results are discussed in the following parts. Gas inlet condition Pin (Pa)?in R.Hin fO2 or CO Cathode 141000 2.5 0.2 0.21 Anode 141000 1 0.001 or 0.5 0 Table 2-4. Working condition for PEM fuel cell stack with low inlet RH Gas inlet condition Pin (Pa)?in R.Hin fO2 or CO Cathode 141000 2.5 0.8 0.21 Anode 141000 1 0.001 or 0.5 0 Table 2-5. Working condition for PEM fuel cell stack with high inlet RH The operating conditions for a start-up and transient response are described as follows: z The start-up temperature for the cells is initially set to 298.15?K. A look-up table was created for the lowest membrane temperature among the 10 cells as an input and the current density as an output. The value of current density continuously follows the increase of the temperature in the membrane in order to quickly raise the temperature to 353.15?K that is assumed as a typical working temperature. In addition, a built-in controller regulates the temperature in the membrane to prevent dehydration. Whenever the temperature in different cells exceeds the desired working temperature, a coolant subsystem is turned on and off to extract the excessive heat produced in the stack. 66 A load current with multiple steps is applied to the stack model and the transient responses analyzed. 2.3.3. ANALYSES 2.3.3.1. STATIC BEHAVIOR Fig. 2-19 shows the temperature dependent I-V characteristics from 333?K to 363?K with a step of 10?K. As the temperature increases, the water removal is easier. The effects are considerably high at the range of the higher cell current, where more water is produced. The increased of the voltage by increased working temperature shown in the results is approved by an analysis using CFD [19]. Fig. 2-19: I-V curve for different cell working temperature 2.3.3.2. DYNAMIC BEHAVIOR 2.3.3.2.1. STARTUP 67 2.3.3.2.1.1. CASE 1: VAPOR PHASE (low cathode inlet RH & fast water removal): Fig. 2-20 shows the current load, the transient behavior of the temperature for the membrane layer of different cells, as well as the stack temperature profile at the 50, 360 and 450 seconds during the start-up. In the given start-up condition (Table 2-5); it took more than six minutes for a membrane layer with the highest temperature to reach the desired working temperature (Fig. 2-20b). Generally, the temperature of each layer rises during a start-up. Due to ohmic losses by the membrane resistances and the heat released by the chemical reaction in the cathode catalysts, the membranes and cathode catalysts show the highest value among others. In contrast, the losses on the anodic side are negligible. Therefore, the high losses of cathodic catalysts and the corresponding heat generated leads to an asymmetrical temperature distribution throughout the cell [35, 36, and 41]. In fact, the temperature of membranes rises to a reference temperature with the associated rising time strongly influenced by the location of a cell in the stack. As a result, the total heat capacity of an endplate assembly is much higher than the one of the layers in the cells. The heat generated in the end cells are quickly transferred to endplates and stored there rather than stored in the cell itself. The temperature of the membrane layers in the middle of the stack increases more rapidly than the end ones. For example, the 5th cell shows the highest temperature, while the 10th cell shows the lowest (Fig. 2-20b). The profile of the stack temperature in Fig. 2-20c, d and e shows different dynamics with an asymmetric shape that differs from the one with symmetric shape [32] or smooth curve [38]. In addition, an unsteady and abrupt drop of temperature at the end cells is observed, caused by the difference between ambient and cell temperature. The large heat 68 capacity of the bus plate and the small heat transfer coefficient of the endplate lead to a negligible heat transfer. At the 50th second after the start-up (Fig. 2-20c), the nine cells show the temperature walls at MEA except the cell number 1. In fact, the cathode catalyst layer of the cell 1 is located near the endplate and possibly conducts a large amount of heat generated toward the endplates. Compared to the 1st cell, the heat generated at the cathode catalyst in the 10th cell is being kept by two more layers of the membrane and anodic catalysts. Thus, the temperature of the cathode catalysts layer in the 1st cell becomes lower than the one of the 10th cell. At the 360th second (Fig. 2-20d), the membrane of the 5th cell shows the highest temperature at 80?C, while six cells still show a temperature wall with an asymmetric distribution. The temperature wall is generated in the 10th cell before blocking the conduction of the heat generated by the rest of the cells. The transfer of heat occurs through the left side of cells. As a result, the temperature of the left end plate assembly rises more than the one of the right end plate assembly. Subsequently, the temperature walls of cells located at the right hand side of the stack disappear. Then, the temperature of the cell rises gradually, while the 1st cell does continuously. Finally, the temperature of the cathode catalysts layer in the 10th cell becomes lower than the one of the 1st cell. At the 450th second, the numbers of temperature walls increase again. As the temperature at the end plate assembly increases, the gradient of the cells and the end plates becomes lower and thus the heat conduction from the cells decrease. As a result, the heat generated accumulates in the end cells, while the heat in the inner cells causes a temperature rise in the coolant that offsets the end cells with a temperature rise. Finally, the temperature of 69 the stack becomes more uniform. The middle cells reach a steady state, while the temperature dynamics of the end cells have not yet. a) Current input for the startup b) Dynamic membrane temperature c) Stack temperature distribution at the 50 second 70 d) Stack temperature distribution at the 360 second e) Stack temperature distribution at the 450 second Fig. 2-20 Fig. 2-21 shows a dynamic behavior of the output voltage of the individual cell and output voltage referred to the 5th cell voltage that shows the highest value. Generally, the cells located in the inner side of a stack show low losses of over-potentials because of the high temperature. The middle cells illustrate a higher output voltage than those of the end cells. The behavior of the two end cells is particularly different from the rest of the cells (Fig. 2-21a). The voltage of all the cells tends to follow the increase of temperature but those of end cells show a decrease. In fact, the lower temperature at the end cells slows down the removal of water continuously generated, which blocks the influx of the 71 reactants toward the catalysts through the GDL. Likewise, water generated in the catalyst is hard to remove, which causes a high over-potential. This curve of the output voltage is comparable to the experimental results [39]. The output voltage of the cells is shown in the Fig. 2-21b. The voltage increases until 300 seconds and remains almost constant affected mainly by dehydrated membranes and the associated low proton conductivities. On the other hand, the 10th cell shows less loss in the over-potential than the 1st cell until 300 sec, because the temperature of the cathode catalyst at the 1st cell becomes larger than the one at the 10th cell. As soon as the coolants start controlling temperature around 350 seconds, the output voltage of cells except two end cells starts decreasing until the current reaches around 0.6 (A cm -2 ). The operating condition of the anodic inlet RH is intentionally set to 0.6 from 0.001 in order to study the effects of the humidity on the membranes. As a result, the voltage of the cells starts increasing because of the increased membrane proton conductivity. The recovery speed of membrane conductivity of the individual cell becomes dependent on the cell location in the stack, whose physical reasons are described in detail in the following part. Fig. 21b shows percentage rates of the cell voltages. The voltage of the end cells amounts to half of the voltage produced by the central cells. The derivation of the end cells from the central cells varies at a rate of around 20 percent, which also worsens overall performance of the stack. The proton conductivity of membrane currently available strongly depends upon water content. Therefore, an analysis is performed to study the effect of the anodic RH and cell location on membranes. The inlet RH has been drastically changed from 0.001 to 0.6 at the 380 second. 72 a) Dynamic voltage output b) Voltage output percentage Fig. 2-21 Fig. 2-22 shows simulated results of the water uptake of membrane layers in the 10 cells during the start-up. At the low anode inlet gas humidification, water present in the membrane is only influenced by water generated in the catalysts. Thus, the water uptake continuously decreases and corresponding proton conductivity gets lower. At an application of a high humidification, the water uptake in the membranes of all the cells increases. Particularly, the cells at the right hand side show a higher rate of water uptake than its counterpart in 73 the stack. The differences in temperature between the cathodic electrodes of cells at the left-hand side of the stack and at the right-hand side leads to different saturating pressures, although the amount of water generated is the same at the load current. The lower the temperature of the cathode gas, the lower the saturation pressure. Thus, the outlet gas carries less water vapor out of cells, and consequently, the water concentration in the membrane becomes higher. A comparison of water uptake shows a large difference affected by the temperature distribution in the stack. For example, the central cells become quickly hydrated in comparison to the end cells. Under a fully consumed fuel, the anode gas gets saturated. Therefore, the electro-osmotic influence at the anodic boundary of each cell on their membranes becomes the same. However, the temperature distribution leads to much higher water concentration and diffusion coefficient in membrane layers at the central cells than the ones at the end cells. If the influence of the liquid water is neglected, the diffusion of water from the anodic electrodes to membranes in the central cells becomes higher than the one in the end cells. Conversely, the humidity at the cathodic inlet is low and the temperature is high, so the gas is not saturated. Due to the temperature difference in the cells, the saturation pressure of end cells is lower than the one of the central cells (Eq. 9) Consequently, water activity at the cathodic sides of the end cells is higher, which makes it easier to transport water out of membranes by electro-osmotic force at the cathode boundary. As a result, the recovery speed of water contents in the central cells is generally faster than the one in the end cells. However, the low water concentration in the cathode sides of the central cells 74 accelerates water diffusion from membranes to cathodes, which again slows down the recovery speed. Fig. 2-22 Dynamic behavior of membrane water uptake for a stack 2.3.3.2.1.2. CASE 2: LIQUID PHASE (high cathode inlet RH & slow water removal): Fig. 2-23 shows the current load, the transient behavior of the temperature for the membrane layers of different cells as well as the stack temperature profile at the 50, 270, and 350 second. Even with the same lookup table of current and membrane temperature for the start-up, the speed of the start-up becomes quicker than the one in the first case. It takes about 4.5 minutes for the membrane layer with the highest temperature among others to reach the desired working temperature. In fact, the water generated in the cell is removed as liquid phase, so the latent heat necessary for getting water evaporated is stored in the stack. Fig. 2-23b shows the temperature waveforms of membranes. A comparison between Fig. 2-20b and Fig. 2-23b shows two differences: 1) temperature of membranes continuously rises, even if the coolant circuit turns on; 2) the gap between the membrane temperature of the 5th cell and the 10th cell is calculated by the time the 75 coolant circuit turns on 75 seconds later. The difference in both gaps has been decreased from 8?C in the first case to 2?C in this case. Therefore, the low RH at cathode inlet in the first case enables the coolant to easily control the stack temperature at the start-up. Figs. 2-23 d, e and f show the stack temperature profile at the 50, 270, and 350 second. A comparison with the first case shows relatively higher walls of temperature with increased numbers in the stack. a) Current input for the startup b) Dynamic membrane temperature 76 c) Stack temperature distribution at the 50 second d) Stack temperature distribution at the 270 second e) Stack temperature distribution at the 350 second Fig. 2-23 77 Fig. 2-24 shows the dynamic response of the output voltage of the 10 cells on the load current and the corresponding percentage rate. A comparison with the first case shows a low output voltage at about 0.6 A cm -2 , even though the temperature of the stack and membrane conductivity are higher than before. As a matter of fact, the liquid water generated at the catalysts layers blocks the influx of oxygen and finally leads to a starvation of oxygen. This flooding effect is considered by using the empirical parameter which is function of temperature, current density, channel gas pressure, and air flow rate. a) Dynamic voltage output b) Voltage output percentage Fig. 2-24 78 Fig. 2-25 shows the water uptake in the membrane layer of the 10 cells at the start-up. The humidified cathode inlet gas slightly slows down dehydration of the membranes until the 300 second when the anode inlet gas starts to be humidified from 0.001 to 0.6. Then the water uptake of the membrane stops decreasing. A comparison with the results from the changed cathode and anode inlet gas RH shows that water transport at the boundary of membranes is dominated by the electro-osmotic drag force. Therefore, humidification on the cathode inlet gas is not sufficient to prevent dehydration, which requires additional humidification of the gas at the anodic inlet. Fig. 2-25 Dynamic membrane water uptake 2.3.3.2.2. Transient response: 2.3.3.2.2.1. Case 1: Vapor phase (low cathode inlet RH and fast water removal): Fig. 2-26 shows the load current and the simulated results for the temperature of the membrane layers and the stack at the 580, 650 and 800 seconds. 79 When the current drops with 0.15 and 0.3 A cm -2 , the temperature of all membranes drops by 1 and 1.5 degrees respectively (Fig. 2-26b). Primarily, the heat associated with the current density is responsible for the rise or drop in temperature. Therefore, the recovery behavior and final value is strongly influenced by the amount of heat dependent upon the amplitude of current density. For example, a step with a high current density (0.45 A cm -2 ) causes an immediate rise in temperature, while a low current density (0.15 A cm -2 ) takes 100 seconds for the central membrane to reach a desired working temperature of 353.15K. The simulation shows that the temperature rise of the end cells at the high current density is faster than the one at the low current. Besides, the heat depends upon the current density because the coolant controlling for a working temperature carries over the heat from the central cells and offsets the temperature of the end cells. Similarly, the transfer of the heat to the end cells gets smaller. When the amplitude of the current is small, the corresponding heat generated in the central cell becomes smaller. Those effects are illustrated in Figs. 2-26 d, e, and f. a) Step current input 80 b) Dynamic membrane temperature c) Stack temperature distribution at the 580 second d) Stack temperature distribution at the 650 second 81 e) Stack temperature distribution at the 800 second Fig. 2-26 Fig. 2-27 shows the dynamic response of the cell output voltage and membrane water uptake of the 10 cells at the same load current. The overshoot of voltage is caused by several physical factors: 1) proton concentration at the cathode catalysts [1]; 2) oxygen concentration at the catalyst [40]; 3) variation of the cell working temperature. Temperature influences the water saturation pressure (Eq. 9) and the effective areas of cathode catalysts as well as oxygen concentration. More details are discussed in the following part. Figs. 2-27 b and c show responses of water uptake on the load current. The operating condition for the RH at the anode inlet is set to 0.001 from 0.6 at about 550 second when the current density is lower than 0.6 (A cm -2 ). The dependency of the membrane water uptake on the temperature is identical with the previous analyses. 82 a) Dynamic voltage output b) c) Voltage output percentage 83 d) Dynamic membrane water uptake Fig. 2-27 2.3.3.2.2.2. Liquid phase (high cathode inlet RH & slow water removal): The step load current has been applied to study the effects of liquid water removal on the dynamics (Fig. 2-28a). Figs. 2-28 b, c, d, e, and f show the dynamic temperature behavior of the membrane layers and the stack temperature profile at the 400, 550 and 800 seconds. A comparison with the first case of operating condition shows a high temperature drop in the stack. It is caused by the large amount of heat generated in the catalysts when water is generated and removed as liquid. Consequently, the temperature cannot be controlled with the set coolant flow rate until the current density is decreased. At around 900 seconds, the temperature difference between the membrane of the central cell and two end cell amounts to about 5 and 13 respectively, le?ss than the previously calculated values of 5.5 and 15 even at the same amplitude of the current load. ? 84 Generally, the amplitude of the peak temperature at both high and low current density is higher than the first case and the numbers of temperature walls at MEA is more than before (Fig. 2-28 d, e and f). a) Step current input b) Dynamic membrane temperature 85 c) Stack temperature distribution at the 400 second d) Stack temperature distribution at the 550 second e) Stack temperature distribution at the 800 second Fig. 2-28 86 Fig. 2-29 shows the same dynamic responses of output voltage of all the cells. The amplitude of the overshoot of the output cell voltage is much larger than the first case according to Fig. 2-24a, which is strongly dependent on the temperature. An analysis shows that a sudden temperature drop with a high current causes high amplitude of the overshoot. Conversely, at the low current step, the heat generated can be more effectively extracted by the same condition of the coolant circuit, so the influence of temperature on the voltage overshoot becomes less. In addition, the waveform of voltage of the 10 cells does not have the same tendency but fluctuates with distortions at special operating conditions. For example, the output voltage of the 2nd cell at an instant can be higher than the 5th cell expected with the highest voltage. The open circuit voltage (OCV) is a function of the temperature whose derivative shows a negative value. Therefore, OCV decreases when the temperature rises. Consequently, the 2nd OCV is higher than the 5th. In addition, the high temperature causes more water to be stored as vapor at the electrodes. Particularly, the 2nd cell has lower temperature than the one in the 5th, which leads to a less water vapor. The concentration of oxygen becomes larger and the associated over-potential does smaller. On the other hand, temperature increase accelerates dehydration of the membranes (Fig. 2-29 c) that decreases the membrane proton conductivity. Therefore, corresponding proton conductivity of the 2nd cell is higher than the 5th cell. Although the liquid water in the 5th cell is less than the one in the 2nd, the overall cell performance of the 2nd cell is better because of the three reasons aforementioned. 87 a) Dynamic voltage output b) c) Voltage output percentage 88 d) Dynamic membrane water uptake Fig. 2-29 2.4. SUMMARY In this work, firstly, a new dynamic model for a PEM fuel cell is proposed, which can predict the behavior of the cell when an electric load is applied. Emphasis is placed on temperature response on the load. The model is constructed with layers that represent membrane, catalysts, gas diffusion layers, and bipolar plates. The models for the individual layers are separately and mathematically described. The results show the various effects of electric loads on the temperature, and the voltage, finally efficiency. Particularly, description of interrelated physical variables by varying temperature provides more realistic tools that can be utilized for deriving design parameters of a fuel cell stack and systems. For an example, a startup process and the associated effects are analyzed in detail. After the cell model is extended for a stack, the variables of the layers for a cell are coupled. Then, analyses have been carried out to the effects of temperature on static and 89 dynamic characteristics of cells. Two operating conditions are used for studies on the stack at a start-up and a step load. The analyses suggest a new consideration for design and operating criteria to improve the performance of cells and optimize components of the power system. In particular, proper thermal management is required to prevent dehydration of membranes and secure longevity of the cells. Likewise, optimal operations of the stack can be derived by this simulation and associated new control strategies can be developed. Major findings are summarized as follows: ? Temperature distribution through the cells is asymmetric because of the heat generated in the cathode sides are higher than in the anodic sides. This asymmetric effect coupled with the end plate finally leads to an asymmetrical temperature profile in the stack (Figs. 20b, 23b, 26b, 28b) determined by thermal conductivities of the end plates. The temperature wall in the cells blocks further conduction of heat generated by other cells contingent on their location, and strongly influences the static and dynamic characteristics; ? The heat conductivity of the membrane varies during operations because of the continuous variation of water content and the consequent swelling phenomena, which results in an increase of electrical and thermal resistance. ? The latent heat of water produced at the catalysts by the reactions can be stored in the stack and requires a proper cooling strategy at a start-up when the cathode outlet gas does not sufficiently carry all the water vapor. In this case, the time necessary for a startup is relatively short. Thus, the temperature difference between the central cells and the end cells is larger than the one with water removal during the vapor phase. 90 ? Non-uniform temperature distribution through the stack can be minimized by properly coupling the coolant for the central cells with the end cells which can offset the temperature of the end cells. ACKNOWLEDGEMENTS We are grateful to a visiting scholar, Mr. Suk-Heung Song for his support. 91 REFERENCES 1. M. Ceraolo, C. Miulli, A. Pozio (2003) Modelling static and dynamic behaviour of proton exchange membrane fuel cells on the basis of electro-chemical description. Journal of Power Sources Vol. 113, Issue 1: 131-144. 2. J.C. Amphlett, R.M. Baumert, R.F. Mann, B.A. Peppley, P.R. Roberge, T.J. Harris (1995) Performance modeling of the Ballard Mark IV solid polymer electrolyte fuel cell. 2: Empirical model development. Journal of the Electrochemical Society Vol. 142, Issue 1: 9-15. 3. S.D. Gurski, D.J. Nelson (2003) Cold start fuel economy and power limitations for a PEM fuel cell vehicle. Presented at the SAE 2003 World Congress & Exhibition, Detroit, MI. March 2003. 4. Y. Zhang, M. Ouyang, Q. Lu, J. Luo, X. Li (2004) A model predicting performance of proton exchange membrane fuel cell stack thermal systems. Applied Thermal Engineering Vol. 24, Num.4: 501-513. 5. P.R. Pathapati, X. Xue, J. Tang (2005) A new dynamic model for predicting transient phenomena in a PEM fuel cell system. Renewable Energy Vol. 30, Num. 1:1-22. 6. M. Wohr, K. Bolwin, W. Schnurnberger, M. Fischer, W. Neubrand, G. Eigenberger (1998) Dynamic modelling and simulation of a polymer membrane fuel cell including mass transport limitation. International Journal of Hydrogen Energy, Vol. 23, Num. 3: 213-218. 7. R. Eckl, W. Zehtner, C. Leu, U. Wagner (2004) Experimental analysis of water management in a self-humidifying polymer electrolyte fuel cell stack. Journal of Power Sources Vol. 138, Issues 1-2: 137-144. 8. J. Golbert, D.R. Lewin (2004) Model based control of fuel cells: (1) Regulatory control. Journal of Power Sources Vol. 135, Issue 1-2: 135-151. 9. X. Xue, J. Tang, A. Smirnova, R. England, N. Sammes (2004) System level lumped parameter dynamic modeling of PEM fuel cell. Journal of Power Sources Vol. 133, Issue 2: 188-204. 10. C. Maxoulis, Dimitrios, N. Tsnolou, G.C. Koltsakis (2004) Modeling of automotive fuel cell operation in driving cycles. Energy Conversion and Management Vol. 45, Issue 4: 559-573. 92 11. H.M. Jung, W.Y. Lee, J.S. Park, C.S. Kim (2004) Numerical analysis of a polymer electrolyte fuel cell. International Journal of Hydrogen Energy Vol. 29, Issue 9: 945- 954. 12. http://www.emmeskay.com/fuelcell.pdf. 13. B. Wetton, K. Promislow, A. Caglar (2004) A simple thermal model of a PEM fuel cell stacks. Second International Conference on Fuel Cell Science, Engineering and Technology 14. M. Sundaresan (2004) A thermal model to evaluate sub-freezing startup for direct hydrogen hybrid fuel cell vehicle polymer elctrolyte fuel cell stack and system. PhD Dissertation, University of California at Davis. 15. B. M. Eaton (2001) Cold start effects on performance and efficiency for vehicle fuel cell systems. Master Thesis, Virginia Polytechnic Institute and State University. 16. J.T. Pukrushpan, A. Stefanopoulou, H. Peng (2002) Modeling and control of fuel cell PEM stack systems. Proceedings of the 2002 American Control Conference: 3117- 3122. 17. M. Lampinen, M. Fomino (1993) Analysis of free energy and enthalpy changes for half cell reactions. Journal of Electrochemical Society. Vol. 140, No. 12: 3537-3546. 18. J.C. Amphlett, E.K. De Oliveira, R.F. Mann, P.R. Roberge, A. Rodrigues, J.P. Salvador (1997) Dynamic interaction of a proton exchange membrane fuel cell and a lead acid battery. Journal of Power Sources Vol. 65, Issues 1-2: 173-178. 19. Y.M. Ferng, Y.C. Tzang, B.S. Pei, C.C. Sun, A. Su (2004) Analytical and experimental investigations of a proton exchange membrane fuel cell. International Journal of Hydrogen Energy Vol. 29: 381-391. 20. John P. Evans (2003) Experimental evaluation of the effect of inlet gas humidification on fuel cell performance. Master Thesis, Virginia Polytechnic Institute and State University. 21. http://webbook.nist.gov. 22. P. Vie, S. Kjelstrup (2004) Thermal conductivities from temperature profiles in the polymer electrolyte fuel cell. Electrochimica Acta. Vol. 49, Issue 7: 1069-1077. 23. W. Yan, F. Chen, H. Wu, C. Soong, H. Chu (2004) Analysis of thermal and water management with temperature dependent diffusion effects in membrane of proton exchange membrane fuel cells. Journal of Power Sources Vol. 129, Issue 2: 127-137. 93 24. T.E. Springer, T.A. Zawodzinski, S. Gottesfeld (1991) Polymer electrolyte fuel cell model. Journal of the Electrochemical Society Vol. 138: 2334-2341. 25. R.M. Moore, K.H. Hauer, D. Friedman, J. Cunningham, P. Badrinarayanan, S. Ramaswamy, A. Eggert (2005) Adynamic simulation tool for hydrogen fuel cell vehicles. Journal of Power Sources, Vol. 141, Issue 2: 272-285. 26. J.T. Pukrushpan, H. Peng, A.G. Stefanopoulou (2004) Control oriented modeling and analysis for automotive fuel cell systems. Journal of Dynamic Systems, Measurement and Control, Vol. 126, Issue 1: 14-25. 27. B.R. Sivertsen, N. Djilali (2005) CFD based modeling of proton exchange membrane fuel cells. Journal of Power Sources, Vol. 141, Issue 1: 65-78. 28. U. Pasaogullari, C.Y. Wang (2005) Two-phase modeling and flooding prediction of polymer electrolyte fuel cells. Journal of the Electrochemical Society, Vol. 152, Issue 2: A380-A390. 29. D. Natarajan, T.V. Nguyen (2003) Three-dimensional effects of liquid water flooding in the cathode of a PEM fuel cell. Journal of Power Sources, Vol. 115, Issue 1: 66-80. 30. J.C. Amphlett, R.F. Mann, B.A. Peppley, P.R. Roberge, A.Rodrigues (1996) A model predicting transient responses of proton exchange membrane fuel cells. Journal of Power Sources, Vol. 61, Issue 1-2: 183-188. 31. E.A. Muller A.G. Stefanopoulou (2005) Analysis, modeling, and validation for the thermal dynamics of a polymer electrolyte membrane fuel cell system. Proceedings of FUELCELL 2005 Third International Conference on Fuel Cell Science, Engineering and Technology, ASME. 32. M. Sundaresan, R.M. Moore (2005) Polymer electrolyte fuel cell stack thermal model to evaluate sub-freezing startup. Journal of Power Sources, Vol. 145 Issue 2: 534-545. 33. V. Peinecke (1994). PhD Dissertation. University of Karlsruhe, Germany. 34. P. Berg, K. Promislow, J. St-Pierre, J. Stumper, B. Wetton (2004) Water management in PEM fuel cells. Journal of the Electrochemical Society, Vol. 151, Issue 3: A341-A353. 35. Y. Shan, S.Y. Choe (2005) A high dynamic PEM fuel cell model with temperature effects. Journal of Power Sources, Vol. 145, Issue 1: 30-39. 94 36. H. Ju, H. Meng, C.Y. Wang (2005) A single phase non isothermal model for PEM fuel cells. International Journal of Heat and Mass Transfer, Vol. 48, Issue 7: 1303- 1315. 37. S. Shimpalee, S. Dutta (2000) Numerical prediction of temperature distribution in PEM fuel cells. Numerical Heat Transfer, PartA, Vol. 38, Issue 2: 111-128. 38. J.H. Lee, T.R. Lalk (1998) Modeling fuel cell stack systems. Journal of Power Sources, Vol. 73 Issue 2: 229-241. 39. J.P. Meyers (2005) Fundamental issues in subzero PEMFC startup and operation. DOE Freeze Workshop. 40. Y. Wang, C.Y. Wang (2005) Transient analysis of polymer electolyte fuel cells. Electrochimica Acta, Vol. 50, Issue 6: 1307-1315. 41. A. Rowe, X. Li (1993) Mathematical modeling of proton exchange membrane fuel cells. Journal of Power Sources, Vol. 102, Issue 1-2: 82-96. 95 CHAPTER 3: 2-D PEM FUEL CELL STACK MODELING 3.1. INTRODUCTION Current computational models available in both the academic world and the market are either too simple or complex to particularly describe dynamic behaviors of a PEM fuel cell stack. Some authors [1, 2, and 3] simply employ empirical equations, whose coefficients are obtained by fitting a polarization curve. This approach is useful for a design of the power system, but ignores effects of temperature, water and reactants on the cell performances. Thus, it can hardly describe the complex behaviors of a stack. Conversely, two or three dimensional models using CFD techniques proposed by other authors [4, 5, 6, and 7] can capture the complexity of a single cell, but are limited to steady state analyses and unable to represent an unsteady behavior of a stack particularly taking into account the varying load. Um et al. [8] published a two-dimensional CFD model that describes the transient behavior of the bulk flow, species and electro-chemical reactions in a single cell, which was extended to an isothermal three-dimensional model [10, 11]. Likewise, Dutta et al. [9] developed a 3D model with an isothermal flow in a cell that embeds a serpentine-type gas channel. However, the simulation is conveyed by the use of commercial software package, Fluent, and they studied both the gas distribution and water generation in a single cell. 96 Um and Wang [10] developed a three-dimensional CFD model for a single cell. They intensively studied species and water removal in a straight and an inter-digitated flow channel and found enhancement of the performance at the inter-digitated shape. Furthermore, Wang [11] simulated a single cell with 36 gas serpentine channels taking a low humidity condition by using the software package of Star-CD and presented the mechanism of the species transport and the associated current density distribution. Unlike the studies above, Ju et al. [12] considered a non-isothermal condition in the three-dimensional plane and simulated a cell with a straight channel by using Star-CD. All of suggestions above, however, have been focused on description of a single cell and are still unable to describe a stack and time-varying behaviors. On the other hand, the dynamic behavior of a stack can be improved by adding a simplified thermodynamic model, which is proposed by Sundaresan [13, 14]. The model regards a cell as a composition of layers and is used to analyze the start-up behavior from a sub-freezing temperature. However, the model does not fully consider several factors; 1) Flow of species at the inlet must be the same as that at the outlet. Thus, no fluid dynamics are considered; 2) Heat source terms in both the catalysts are empirically calculated with values suggested by the Wohr and Peinecke?s model [15]. Accordingly, the anode source term is presumed as a relatively large value that in fact should be referred to be around zero [16]. As a result, the model does not show asymmetric phenomena of performance through the stack. Wetton et al. [17] proposed an explicit stack thermal model with the coolant channel coupled with a 1-D cell model [18]. It shows an outstanding temperature gradient of the stack, but with no dynamics at all. We proposed an enhanced quasi 1-D stack model [19, 20] based on the previous single cell 97 model that considers the thermal and fluid dynamics. As a result, the model proposed is capable of capturing the dynamic temperature distribution including the asymmetrical effects in the stack, but missing the water distribution that are improved by adding an empirical relationship between the flooding effect and the current density and temperature. As the matter of fact, none of current models can fully describe the stack behavior considering two phase. Therefore, the first goal has been set to develop a two dimensional stack model with a single phase. On the system aspects, the model should be a current controlled voltage source, so that the load can be easily integrated into a stack model. Moreover, the domain should be so set up that can integrate two endplates and bus plate at the anodic and cathode side with an interface plate that embeds a coolant channel as well as two bipolar plates with a basic cell unit. The mass and charge transport as well as the heat flux in the basic cell unit is described by the use of the Navier-Stokes and the potential and energy conservation equation. Likewise, the heat flux in the coolant channel is described by using both the Navier-Stokes and the energy conservation equations, while the rest of plate regions are described by the heat conservation equations. The PDEs (Partial Differential Equations) are solved by the SIMPLE [21] method, which allows for a calculation of the values at an each time step until it is converged. In the following chapters, details on the two dimensional models are summarized, which includes assumptions, simulation set-up and equations used for a description of the model proposed. Included are the procedure to solve the equations and a generation of grids. At the end, simulations are conveyed with boundary conditions used for the individual domain and simulated results are discussed. 98 3.2. MODEL DESCRIPTION 3.2.1. MODELING DOMAIN AND ASSUMPTIONS Major assumptions are made for a two-dimensional stack model as follows: 1. Reactants as ideal gases; 2. Incompressible and laminar flow; 3. Isotropic and homogeneous electrodes, catalyst layers and membrane; 4. Identical inlet conditions of each cell for both the cathode and anode as well as coolant channel; 5. Constant thermal conductivity of the materials in a fuel cell; 6. Neglect diffusions caused by multi-component; 7. No contact resistance; 8. No liquid water generated; In addition, it is assumed that a single cell has a structure of sandwiched layers shown in Fig. 3-1. The anode sides of the cells are located on the left hand side, while the cathode sides on the right hand side. The single cell domain for the model is constructed with seven different layers that are symmetrically placed at the membrane layer. A gas flow channel, a gas diffusion layer and a catalyst layer for the anode side are located at the left side of the membrane layer as well as those for cathode side with a reversed order. Thus, a stack can be easily constructed by repeating this basic unit domain and adding bipolar plates with coolant channel, interface and bus plates, and end plates shown in Fig. 3-2. Finally, a stack model is completed by coupling of the domains for the basic units and two endplate assembles and simulation can be performed; 99 Fig. 3-1 Single cell schematic configuration Fig. 3-2: Stack schematic configuration 100 3.2.2. MODEL DESCRIPTION 3.2.2.1. CHARGE TRANSPORT Protons and electrons are the positively and negatively charged ions. The proton transfer in the proton conducting regions and the electron transfer in the electronic conducting regions determines the potential distribution in a cell. The polymer as electrolyte in the membrane and the catalysts belongs to the proton conducting region, while the catalysts, GDLs and BPs including gas flow and coolant channels are regarded as electrode, which is repeating in the stack configuration. The potentials in the electrolyte are governed by the potential conservation equation according to the Ohm?s law: ()0 ee S? ? ?? ?? + = (1) jS = ? in the two catalyst layers (1a) 0S ? = in the membrane layer (1b) where the proton conductivity ? e is a function of temperature and the water content in the polymer material [24]; 23 11 100 exp(1268 ( )) (0.005139 0.00326) 303 eHOS T ??=? ? (2) According to the Butler-Volmer equation [1, 8], the current densities in the anode and cathode catalysts can be expressed by the exchange current density, reactant concentration, temperature and over-potentials according to the butler volmer equation [Eq.3 and 4]: ? ? ? ? ? ? + ? ? ? ? ? ? ? ? = ? ?? F RTX X ajj ca refH H ref aa 2/1 , ,0 2 2 (3) 101 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= ? ? RT F X X ajj c refO O ref cc exp 2/1 , ,0 2 2 (4) where the surface over-potential is defined as a difference between the electrodes and electrolyte referring to an equilibrium state. (, ) equlibrium s e OC x yVV V? =? =???? (5) ? s and ? e are the potentials of the electrons conducting solid materials and electrolyte, respectively, at electrodes and electrolyte interface. The open circuit potential at the anode is assumed to be zero, while the open circuit potential at the cathode becomes a function of a temperature as [1,9]: 0.0025 0.2329 c OC VT=+ (6) Then, the local current density for the protons can be simplified with ee I ?=? ?? (7) On the other hand, the electronic conducting region includes the entire stack except the membranes and the two endplates. In fact, the electronic conductivity in the entire electron conducting layers is at least two orders higher than the proton conductivity. Thus, the potential drop caused by the electrons transfer is negligible. Therefore, the values of the potentials in the electron conducting regions can be regarded as a single value. When the load current is applied to the model as an input, the voltage drops in the electron conducting regions vary because of the flow field being changed. This effect is reflected by using a current conservation equation at an equilibrium state of the potential field and described by the following equation. () 0 e a input ch V fjdVIL???= ? = ? in anode catalyst layers (8) 102 () 0 s c input ch V fjdVIL???= ? = ? in cathode catalyst layers (9) Accordingly, both of the electrolyte potential distribution and electrode potential (? s and ? m ) can be corrected. The resulting equations are implicit and can be solved numerically. ,( ) () 0 ee e e a input ch V fjdVIL ?+?? ?+?? = ? = ? (10) ,( ) () 0 ss s s c input ch V fjdVIL ?+?? ?+?? = ? = ? (11) Then, the new electrolyte potential is obtained by adding the correction factor calculated in the anode catalyst layer to the old potential for the domain of the catalyst and membrane, while the old potentials are obtained by solving the potential conservation equation Eq. 1. Likewise, the new electrode potentials are calculated by adding the correction factor calculated in the cathode catalyst layer to the old potential. new old ee e ?=?+?? , new old s ss ?=?+?? (12) 3.2.2.2. MASS TRANSPORT 3.2.2.2.1. ANODE SIDE The anode side includes a gas channel, a gas diffusion layer, and catalysts. The role of the layers is to provide a pathway for the transfer of the fuel and humidity from the inlet to the catalysts. Generally, the gas channels are machined into the collector plates or fabricated by stamping or injection molding. The entire network of gas channels is designed to uniformly distribute reactants across the surface of the gas diffusion layer. The flow field can be a shape of either serpentine or inter-digitated channels. Compared to the serpentine shape, the inter-digitated type can dramatically increase the mass 103 transfer in the gas diffusion layer by convection. Due to the geometry, the pressure drop along the channel is generally high and more extra power is needed to keep the same pressure at the catalysts. In this study, a straight channel is selected to represent a serpentine channel. Then, the mass transport in the channel is governed by the following equations; Mass conservation equation: () ()0=??+ ? ? u t G ?? ?? (13) Momentum conservation equation: () () ( ) u eff Supuu t u +???+??=??+ ? ? GGG G ????? ?? (14) u K S u G? ?= in the gas diffusion and catalyst layers [5] (14a) 0 u S = in the gas channel (14b) Similarly, the hydrogen concentration can be obtained by using the species conservation equation: ()( ) 2 2222 ( ) H eff H HHH X uX D X S t ? ? ? +?? =?? ? + ? G (15) 2 2 a H tota j S Fc =? in gas diffusion layers and catalyst layers [1] (15a) 2 0 H S = in gas channel (15b) where the effective diffusivity D eff H2 is, 22 1.5 eff H H D D?= , 22 OHtoa X cc= (16) 104 3.2.2.2.2. CATHODE SIDE Like the anodic side, the cathode side includes a cathode gas channel, a gas diffusion layer and a catalyst. The flow channel and the gas diffusion layer serve as a pathway for the transfer of the oxygen from the inlet to the catalysts and at the same time extract water produced from the catalysts to the outlet. The channel shape is assumed to be the same one on the anode side and a straight channel. The mass transport in the cathode side is governed by using the following equations, where the air is supplied for the oxygen: Mass conservation () ()0=??+ ? ? u t G ?? ?? (17) Momentum conservation () () ( ) u eff Supuu t u +???+??=??+ ? ? GGG G ????? ?? (18) u K S u G? ?= in the gas diffusion and catalyst layers (18a) 0 u S = in the gas channel (18b) Species conservation () ()( ) 2 2222 O eff OOOO X uX D X S t ? ? ? +?? =?? ? + ? G (19) 2 4 c O totc j S Fc = in catalyst layers (19a) 2 0 O S = in gas channel and gas diffusion layers (19b) () ()( ) 2 2222 HO eff H OHOHOHO X uX D X S t ? ? ? +?? =?? ? + ? G (20) 105 2 2 c HO totc j S Fc =? in the catalyst layers (20a) 2 0 HO S = in the gas channel and gas diffusion layers (20b) where the effective diffusivities and mole fractions of the oxygen and water are as follows; 22 1.5 eff OO D D?= , 22 1.5 eff H OHO DD?= , 22 OOtoc X cc= , 22 H OHOtoc X cc= (21) 3.2.2.2.3. WATER TRANSPORT IN MEMBRANE Water distribution in the membrane is determined by the electro-osmotic force, the diffusion and convection. At a single phase, the influence of the convection is negligible and the water concentration can be described by the following equations. 2 22 e ( )( ) HO HO HOe ed c I Dc n tF ? =?? ? ??? ? (22) The first term describes the diffusion dependent water flux and the second one does the water transport dependent upon the electro-osmotic force. And the electro-osmotic coefficient is a function of the local water content (n d ); 23 2.5 22 H OSO d n ? = (23) The water content ? H2O/SO3 is a function of the water uptake obtained from the water concentration: 2 23 2 HO e HOSO e dry H O ee c bc M ? ? = ? (24) 106 The water flux at the boundary interfaces between the catalysts and membrane is described by using the Robin type boundary condition, which presents the relationship between the water concentration in the gas phase and in the membrane [25]. 222 , () HO HO HO d bc e e eq n Jcc I F ?=?? (25) However, compared to the reactants transport phenomena in the gas diffusion layer, the water transport in the membrane is more complicated. Some experiments showed the different time scales of the hydration and dehydration processes in the Nafion [25]. It takes about tens of seconds for liquid water to get adsorbed, which is two orders of magnitude larger than for the dehydration process. This hydration-dehydration dichotomy might result from two reasons, 1) the water diffusion coefficient in the membrane is a function of local water concentration. And it varies with the changing of water concentration; 2) the Robin type water transfer coefficient ? is also affected by the water concentration at the boundaries between the membrane and catalysts. Springer [23] proposed a coefficient of water diffusion equation that is derived from experimental results: 2 11 exp 2416 303 HO e e DD T ???? ?=? ?? ?? ?? ?? (26) where 23 23 23 23 23 23 2.64227 13 1.23 7.75 11 9.5 11 6 1.23 2.5625 11 2.1625 10 14 6 HOSO HOSO HOSO HOSO HOSO HOSO e De e ee ??? ?> ? ? ?=????? ? ? ???>> ? ? (26a) 107 According to the data proposed Berg [7], the Robin type water transfer coefficient is chosen with ? l = 5 * 10 -4 when a membrane is fully immersed into liquid water, while ? g = 4.5 * 10 -6 for the membrane dried in a gaseous condition. 3.2.2.2.4. COOLANT FLOWS The coolant flow channels embedded in bipolar plates or endplates provide part of the pathway for the coolant circuits. After the heat exchange between the coolants and the heat source has occurred at the wall, the coolant releases most of the heat to the environment via the radiator. If the coolant flow channel is assumed as a straight channel, the flow of the coolant can be governed by the following equations. Mass conservation ()0u t ? ? ? +?? = ? G (27) Momentum conservation () () () u uu p u t ? ?? ? +?? =?? +?? ? ? G GG G (28) 3.2.2.3. HEAT FLUX IN A STACK The heat in the stack is produced by five different sources that include the entropy and losses caused by over-potentials at two catalysts, proton conductivity, electron conductivity, and the phase change of water. Under the condition of a single phase, there is no heat generation associated with the phase change. In addition, all the heat generated is then transferred from the source to the fluid by conduction and convection and 108 completely removed out of the stack by the coolant and gases at the channel outlet. Then, the thermal behavior of a stack can be governed by the energy conservation equation [5]: ,, ()( )() eff mpm f pf T CT uCT k TS t ??? ? +?? =?? ? + ? G GK G (29) 2 ( ) OC T e dV I Sj T dT ? ? =+ + in the catalyst layers (29a) 2 T e I S ? = in the membrane layer (29b) 0 T S = in the gas channel and gas diffusion layers (29c) where the overall density and thermal conductivity are defined as mfstack ? ??=+ , eff eff eff f stack kkk=+ (30) and the fluid mixture properties are fii i X? ?= ? , ,,fpf ipi i CC??= ? (31) 3.2.2.4. BOUNDARY CONDITIONS 3.2.2.4.1. FOR THE FLUID FLOW The fluid velocities at the inlet of the anode, the cathode, and the coolant channel are predetermined. The standard exit boundary and no-slip boundary conditions are applied to the channel exits and walls, respectively. In the species field, the inlet species concentrations are given, while the species gradients at the channel exits and walls are set to zero. All boundary conditions for the fluid flow are summarized as below. At the anode inlet 109 22222 ,,, 0, , , H in H H in H H in uvvXXTT== = = (32) At the anode outlet 222 0, 0, 0, 0, 0 HO H H XXT v u yy y y ??? ? == = = = ?? ?? (33) At the cathode inlet 222 ,, 0, , , OOOinairairn uvvXXTT== = = (34) At the cathode outlet 22 0, 0, 0, 0, 0 HO O air XX Tv u yy yy ?? ?? == = = = ?? ? (35) At the coolant inlet ,, 0, , Coolant in Coolant Coolant in uvv TT== = (36) At the coolant outlet 0, 0, 0 coolant Tv u yy ?? == = ?? (37) At the wall 0,0,0,0,0 222 = ? ? = ? ? = ? ? == y X y X y X vu OHOH (38) 3.2.2.4.2. FOR THE ELECTROLYTE POTENTIAL FIELD The boundary conditions for the gradient of the electrolyte potential field are set to zero at the left anode catalyst as well as the right cathode catalysts 0 e x ?? = ? (39) 110 0 e x ?? = ? (40) 3.2.2.4.3. FOR THE ELECTROLYTE POTENTIAL FIELD The boundary conditions for the gradient of the electrical potential field is set as zeros at the right anode catalyst boundary as well as the left cathode catalyst boundary 0 s x ?? = ? (41) 0 s x ?? = ? (42) 3.3. NUMERICAL SOLUTION First of all, the conservation equations are discretized by using the control-volume-based finite difference method. The flow solution procedure is based upon the SIMPLE algorithm [6] with a collocated grid cell centered approach. The flow chart in the Fig. 3-3 shows the procedure of how to solve the equations for a stack. When the program gets started, the first step is to generate grids for the stack. Then, the program reads an input file that includes initial values and initializes all the field values. Thereafter, calculation begins for velocity, pressure, species concentrations and phase potentials for each of the basic unit cell in the stack. Velocity and pressure fields for coolant channels are calculated in the same manner, while calculation for the energy transport is being carried out for the whole stack domain as well as in each gas and flow channel. Once the calculation of the inner loop gets converged, the calculation of the next outer loop is started. 111 Fig. 3-3 Overall solution procedure The simulation set-up for the stack includes two endplate assemblies, two coolant and flow channels, and two cells with one bipolar plate. (Fig. 3-4) 112 Fig. 3-4 Stack geometry information The simulation set-up for the cell includes anode gas channel, anode gas diffusion layers, and anode catalyst layer on the left side of the membrane, and cathode gas channel, cathode gas diffusion layers, and cathode catalyst layer on the right side of the membrane. (Fig. 3-5) 113 Fig. 3-5 cell geometry information Fig. 3-6 shows an example for meshing two-cells, where the stretched Cartesian type grids are employed. The left cell is numbered as cell number 1, and the right one as cell number 2. The input parameters used for the current simulation are summarized (see Table 3-1, 3-2 and 3-3). 114 Fig. 3-6 Mesh generated in the 2-cell stack Quantity Value Gas channel length, L 7.112 cm Oxygen diffusivity in gas 5.2197?10 -2 cm 2 s -1 Hydrogen diffusivity in gas 2.63?10 -2 cm 2 s -1 Dissolved oxygen diffusivity in active layer and membrane 2.0?10 -4 cm 2 s -1 Dissolved hydrogen diffusivity in active layer and membrane 2.59?10 -6 cm 2 s -1 Faraday constant, F 96487 C mol -1 Permeability of backing layer, K 1.76?10 -6 cm 2 Universal gas constant, R 8.314 J mol -1 K -1 Cathodic transfer coefficient 2 Anodic transfer coefficient 2 Inlet nitrogen-oxygen mole fraction 0.79/0.21 Air-side inlet pressure/fuel-side inlet pressure 5/3atm O 2 stoichiometric flow ratio 3.0 H 2 stoichiometric flow ratio 2.8 Reference exchange current density x area of anode 5.0?10 2 A cm -3 Reference exchange current density x area of cathode 1.0?10 -4 A cm -3 Total mole concentration at the anode side 66.81?10 -6 mol cm -3 Total mole concentration at the cathode side 17.81?10 -6 mol cm -3 Table 3-1. Parameters for 2-D models Thickness m Heat conductivity W m -1 K -1 Heat capacity J kg -1 K -1 Density Kg m -3 GDL 0.0004 4 840 2000 115 Catalyst Layer 0.000065 0.2 770 387 Membrane Layer 0.000183 0.21 1100 1967 Gas Channel 0.001 52 935 1400 Plate 0.001 52 935 1400 Coolant Channel 0.001 30 935 1400 GDL Porosity 0.4 Catalyst Layer Porosity 0.2 GDL Tortuosity 3.725 Bipolar Plate Contact Area Percentage 0.5 Membrane Molecular Mass 1.1 kg mol -1 Fuel Cell Area 0.0367 m 2 Fuel Cell Active Area 0.03 m 2 Table 3-2. Geometry parameters for the fuel cell Quantity Value Temperature 353 0 K Oxygen nondimensinal concentration 0.21 Cathode inlet velocity 0.334 m s -1 Anode inlet velocity 0.157 m s -1 Table 3-3. Initial values 3.4. ANALYSES Fig. 3-7 shows the simulated result of the pressure field in the cell number 1. The pressure drop on the anode flow channel is very small and negligible, simply because the viscosity of the hydrogen is smaller than that of the air on the cathode. In addition, the velocity of the anode gas is slower than that of the cathode one. It is noted that the Reynolds number for the cathode is 23.7 for the simulation and the resulting pressure drop is 10pa approximately, which is comparable to the results in reference [22]. 116 Fig. 3-7 Pressure drop in the cell 1 (Unit: Pa) Fig. 3-8 shows a velocity profile in the channels of two cells that are quite similar to each other because of the exclusion of two phases flow. All of profiles are parabolic as expected for laminar and steady flow in channel. 117 Fig. 3-8 Velocity profile Fig. 3-9 shows transient behaviors of the oxygen concentrations in the in the cathode side flow channels. When the simulation starts, the oxygen concentrations are changing, but still identical in the both channels because the temperature effects on the reaction is not so high and the reactants consumed are not so different in the two cells. However, the concentration on the left side is lower than the right side in the channel because of the oxygen being consumed by the chemical reactions. The steady state has been reached after the 1sec. 118 119 Fig. 3-9 Dynamics of oxygen concentration at two cells 120 Fig. 3-10 shows a transient behavior of the source term of the current density generated in the cathode catalysts of the cell 1 at the 0.3, 0.5, 1 and 10sec. The magnitude of the current density largely varies in the 1sec. In fact, the source term of the current density are generally influenced by three phenomena with the rise of the temperature, the distribution of the cathode overpotentials and the concentration of the reactants. The temperature in the first 10 sec is not so high and the influence is negligible. The overpotential on the left side of the cathode catalyst in the cell 1 is higher than that on the right side because the protons flow in the direction where the electrolyte potential decreases. The concentration of the reactant at the inlet side is higher than that on the outlet side. Therefore, the current density gets higher at the inlet side. However, the magnitude of the influences through and along the planes depends upon which one is dominant in an operating condition. For an example, when the concentration of the reactants decreases along channel, the overpotential gets increased because of the change of the membrane conductivity that is explained in details by Wang [12]. As a matter of fact, water generated is being accumulated at the outlet side and hydration rate in the membrane becomes higher, while the inlet side relatively gets dehydrated because of less water accumulated. Likewise, the concentration of the reactants gradually decreases through the planes from the GDL to the membrane, while the overpotential increases. As a consequence, the source term of the current density depends upon an instant which one is dominantly influenced. It is observed that the current density along the channel direction is dynamically varying within the 1sec. The current density at the outlet side decreases, while the one at the inlet 121 side increases. This effect is caused by the change of the concentration of the reactants along the channel. The results presented in this paper are not able to represent some effects associated with water concentration, overpotentials, concentration of the reactants caused by two phase flow. In addition, due to the limited calculation time of the codes developed, simulation has been run for only 10 seconds and the steady state is not reached yet. 122 Fig. 3-10 Current density on the cathode side (Unit: A/m 3 ) Fig. 3-11 shows the simulated results of a phase potential field in the electrolyte of the cell 1 representing a catalyst and membrane. The potential field has reached to a steady state at the first second, where the potential drop in the membrane along the channel direction is not high. In fact, the gradient of the membrane conductivity along the gas channel direction is negligible because of the constraints on the single phase and requires a long computational time to see the effects. Conversely, the potential in the cathode catalysts shows a high gradient, because the chemical reaction at the inlet side is high than those on the outlet side. Subsequently, protons at the inlet side are more consumed and tend to flow toward the inlet direction. When the reactant concentration is large, the source term of current density is bigger. As a result, the potential of the electrolyte at the inlet gets smaller, the over-potential gets 123 smaller. It decreases the source term of the current density, finally a equilibrium state is reached. 124 Fig. 3-11 Membrane potential dynamics at the cell 1 (Unit: Volt) Fig 3-12 illustrates a transient behavior of the temperature distribution for a stack including 2 cells at eight different instants, 0.1sec, 0.3sec, 0.5sec, 1sec, 3sec, 5sec, 8sec and 10sec. When the stack gets operated, the heat is generated and temperature rises. At the first instant, the most heat is generated in two catalyst layers when the reaction begins and at the same time the entropy change occurs. Thus, the temperature peak appears in the cathode catalyst layers of two cells instead of the membrane layer. In addition, the source term of the current density at the inlet region is higher than that of the outlet region, so the temperture along the channel direction gets decreased. At the following instants, the heat flux begins to diffuse and temperature in all the planes of the stack rises. 125 126 Fig. 3-12 Stack temperature 2-D (Unit: K) Fig. 3-13 shows the temperature distribution of the gases and fluids in the cell 1 that includes two gas flow channels. The geometry of the cell is referred to the table 3-2. The temperature profiles in both of channels are not symmetrical and the temperature on the cathode side is higher than the anode side because of the larger amount of heat generated in the cathode than in the anode. Due to the large heat transfer area of the GDL by the porosity, the temperature gradient between the catalyst and the GDL is relatively low, while the temperature drop at the interface between the GDL and the flow channel is high. 127 Fig. 3-13 Fluid temperature Fig. 3-14 shows the dynamics of the temperature distribution. Firstly, the temperature dynamically changes at different instants when the times go by from 0.1sec, 0.3sec, 0.5sec, 1sec, 3sec, 5sec, 8sec, and 10sec. It is observed that the temperature in the cells at the very beginning rapidly rises from an initial value of 353K because of the heat generated in the catalysts by chemical reaction and in the membranes by ohmic losses. The shapes of the rising temperature in both of cells are identical. It is shown that the peak value in the catalyst layer of the cell 1 gradually becomes higher than the one in the cell 2. In fact, the endplate assembly used for the stack should be comparably thick because of the mechanical requirement for robustness and represents two large heat sinks. However, the distance between the layer with the heat source and the end assemblies in 128 the stack is different in a typical stack construction and consequently heat transfer properties are not identical to each other. Finally, the temperature profile through the plane gets asymmetrical and the one on the left side cell becomes higher. Details are explained in the chapter II. 129 Fig. 3-14 tack temperature 1-D Fig. 3-15 shows the temperature distribution in the cell 1 and 2, which allows for a direct comparison of how the temperature profiles are dynamically changing at different instants. It is shown that the shapes of the temperature profile at the 0.1sec are the same. The endplate effect is to recognize at the 0.2sec, where the cell 1 shows a low temperature at the anode side. At the 0.5sec and 1sec, the gap between two cells gets larger. The temperature of the cell 1 on the cathode side is lower than that of the cell 2, while the temperature of the cell 1 on the anode side is higher than the one of the cell 2. 130 At the 1sec and 3sec, the gaps on the both sides get larger and the peak values of temperature for both cells begin to deviate. The temperature of the cell 1 becomes higher than the cell 2. At the 8sec and 10sec, the behavior of the temperature profiles remains as the same as before. As a consequence, the performance of the cell 1 becomes better than the cell 2. However, there are some constraints on the analyses conveyed. Firstly, the simulation time has not been sufficiently long enough. Secondly, effects of two phase flow have not considered in this study. Therefore, it is necessary to improve the model and the simulation that do take account those factors. 131 Fig. 3-15 Cell temperature 1-D 3.5. CONCLUSION In this work, a new dynamic model for the PEM fuel cell stack considering the fluid and thermal characteristics is proposed. Emphases are placed on numerically solving the equations and effect of the temperature distribution on the stack performance, which varies dynamically during operations. The model developed assumes a structure sandwiched by layers that include membrane, catalysts, gas diffusion layers, bipolar plates and endplate assembly. The domains specified for the modeling are determined by the way of physically working principles rather than the layers used commonly. The 132 separate setup of modeling domains provided an easiest way to understand the physics involved and consequently reduce the development time for the codes. The model provides distribution of pressure and reactants, current density, temperature and potentials in a stack that vary dynamically. It is particularly interesting for designing a stack and BOPs, where the interacting influences from the neighboring cells are considered. For example, calculations of temperature distribution across the stack can contribute to develop a management strategy for the coolant circuit. Temperature distributions through the stack is asymmetrical and varies dynamically, which partially are confirmed by the simulation using the 1-D model described in previous chapter. The analyses and modeling proposed laid a ground work that can be utilized to design a PEM fuel cell stack and BOPs, which secures a safe operation and increases performance. Future work will include an expansion of the model for two phase flow in a three dimensional geometry. ACKNOWLEDGEMENTS I am very grateful to Dr. Hyung-il Choi for his advice to this work. 133 REFERENCES 1. M. Ceraolo, C. Miulli, A. Pozio (2003) Modelling static and dynamic behaviour of proton exchange membrane fuel cells on the basis of electro-chemical description. Journal of Power Sources Vol. 113, Issue 1: 131-144. 2. J.T. Pukrushpan, A. Stefanopoulou, H. Peng (2002) Modeling and control of fuel cell PEM stack systems. Proceedings of the 2002 American Control Conference: 3117- 3122. 3. S.D. Gurski, D.J. Nelson (2003) Cold start fuel economy and power limitations for a PEM fuel cell vehicle. Presented at the SAE 2003 World Congress & Exhibition, Detroit, MI. March 2003. 4. T. Nguyen, and R. White, (1993) A water and heat management model for proton- exchange-membrane fuel cells, J. of Electrochemical Society, 140, pp. 2178-2186. 5. D. Natarajan, T.V. Nguyen, (2003) Three-dimensional effects of liquid water flooding in the cathode of a PEM fuel cell, J. Power Sources, Vol. 115, pp. 66-80. 6. D. Singh, D.M. Lu, N. Djilali, (1999) A two-dimensional analysis of mass transport in proton exchange membrane fuel cells International Journal of Engineering Science, Vol 37 pp.431-452 7. V. Garua, H. Liu, and S. Kakac, (1998) Two-dimensional model for proton exchange membrane fuel cells, J. of AIChE, 44, pp. 2410-2422. 8. S. Um, C. Y. Wang, and K. S. Chen, (2000) Computational fluid dynamics of proton exchange membrane fuel cells, J. of Electrochemical Society, 147, pp. 4485-4493. 9. S. Dutta, S. Shimpalee, J.W. Van Zee, (2001) Numerical prediction of mass- exchange between cathode and anode channels in a PEM fuel cell, Int. J. Heat and Mass Transfer, Vol. 44, pp. 2029-2042. 10. S. Um, C.Y. Wang, (2004) Three-dimensional analysis of transport and electrochemical reactions in polymer electrolyte fuel cells Journal of Power Sources Vol.125 pp 40-51 11. Y. Wang and C.Y. Wang, (2005) Simulation of Flow and Transport Phenomena in a Polymer Electrolyte Fuel Cell under Low-Humidity Operation, Journal of Power Sources, Vol. 147, pp. 148-161. 134 12. H. Ju, H. Meng, C.Y. Wang, A single-phase, non-isothermal model for PEM fuel cells, International Journal of Heat and Mass Transfer, Vol 48, pp 1301-1315 13. M. Sundaresan, R.M. Moore, Polymer electrolyte fuel cell stack thermal model to evaluate sub-freezing, Journal of Power Sources, Volume 145, Issue 2 pp. 534-545 14. M. Sundaresan (2004) A thermal model to evaluate sub-freezing startup for direct hydrogen hybrid fuel cell vehicle polymer elctrolyte fuel cell stack and system. PhD Dissertation, University of California at Davis. 15. M. Wohr, K. Bolwin, W. Schnurnberger, M. Fischer, W. Neubrand, G. Eigenberger (1998) Dynamic modelling and simulation of a polymer membrane fuel cell including mass transport limitation. International Journal of Hydrogen Energy, Vol. 23, Num. 3: 213-218. 16. M. Lampinen, M. Fomino (1993) Analysis of free energy and enthalpy changes for half cell reactions. Journal of Electrochemical Society. Vol. 140, No. 12: 3537-3546. 17. B. Wetton, K. Promislow, A. Caglar (2004) A simple thermal model of a PEM fuel cell stacks. Second International Conference on Fuel Cell Science, Engineering and Technology 18. P. Berg, K. Promislow, J. St-Pierre, J. Stumper, B. Wetton, (2004) Water management in PEM fuel cells, Journal of the Electrochemical Society, 151 A341- A353 19. Y. Shan, S.Y. Choe, (2005) A high dynamic PEM fuel cell model with temperature effects. Journal of Power Sources, Vol. 145, Issue 1: 30-39. 20. Y. Shan, S.Y. Choe, (2005) Modeling and simulation of a PEM fuel cell stack considering temperature effects. Journal of Power Sources, in press. 21. S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980. 22. S.M. Senn, D. Poulikakos, (2004) Laminar mixing, heat transfer and pressure drop in tree-like microchennel nets and their application for thermal management in polymer electrolyte fuel cells, Journal of Power Sources Vol. 130, pp 178-191 23. T.E. Springer, T.A. Zawodzinski, S. Gottesfeld (1991) Polymer electrolyte fuel cell model. Journal of the Electrochemical Society Vol. 138: 2334-2341. 24. Brandon M. Eaton (2001) One Dimensional, transient model of heat, mass, and charge transfer in a Proton Exchange Membrane. the Virginia Polytechnic and State University. Master Thesis. 135 25. P. Berg, K. Promislow, (2003) Modeling water uptake of proton exchange membrane, Nanotech Vol. 3: 493-496 136 CHAPTER 4: CONCLUSION AND FUTURE WORK 4.1. CONCLUSION Research has been focused on development of a dynamic temperature model and its analyses. In this work, a new 1-D dynamic and a new 2-D dynamic model for the PEM fuel cell stack are proposed. The 1-D model has been applied to analyze dynamic behavior of a fuel cell and stack. For exsample, a start-up and a dynamically varying load are taken to demonstrate the validity of the model. Further improvement has been made by developing a 2D stack model with a single phase. The 2D model includes: 1. Proton transport in the electrolyte; 2. Water diffusion in the electrolyte; 3. Dynamic water transport at the interface between the electrolyte and the electrode 4. Thermal model for the stack components; 5. Thermal model for the fluid in the channels; 6. Fluid dynamics in the gas and coolant channel; 7. Fluid dynamics in the gas diffusion layer and the catalysts; 8. Water transport in membrane governed by hydraulic pressure; The main discovering in this thesis is the dynamic and asymmetrical temperature distribution in the stack, the temperature peaks in the MEA, and the endplate effects on 137 the temperature distribution in the stack. And based on this result, new design of control of coolant circuit will be designed in the future. 4.2. FUTURE WORK However, current model does not consider two phase flow, so that the major issues like water balance in the cell cannot be fully described. Therefore, following improvements are necessary to fully describe the dynamic behaviors of the stack; ? Two phase transport in the gas channel and gas diffusion layer; ? Phase change in the thermal model. This effect is important to predict the temperature distribution in the stack at high current density. And it is one of the critical issues for the water management. ? Water transport in electrolyte governed by hydraulic pressure. This effect is required to predict the potential distribution in the membrane. Consequently, it is important to precisely predict the source of the current density