Control of a Fuel Delivery System for Polymer Electrolyte Membrane Fuel Cells based on a Two-Phase Transient Model by Jinglin He A dissertation submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama August 6, 2011 Key words: PEM fuel cell, Fuel delivery system, Two-phase, Transient model, Ejector, Controls Copyright 2011 by Jinglin He Approved by Song-Yul Choe, Chair, Associate Professor of Mechanical Engineering Sushil H. Bhavnani, Alumni Professor of Mechanical Engineering Jeyhoon M. Khodadadi, Alumni Professor of Mechanical Engineering Andrew B. Shelton, Assistant Professor of Aerospace Engineering ii Abstract Polymer electrolyte membrane fuel cells (PEMFCs) are a type of fuel cell that converts the chemical energy released by the reaction of hydrogen (fuel) and oxygen into electrical energy and generates water and heat. The fuel delivery system (FDS) is designed to supply hydrogen from a storage tank to the fuel cell stack and, in some designs, reuses the exhausted fuel. In this research work, a hybrid FDS used in fuel cell vehicles is proposed, which uses an ejector and a blower dependent upon loads that circulate unconsumed hydrogen to increase efficiency of fuel usage. In addition, stoichiometric ratio (SR) of the hydrogen, defined as the ratio of the supplied hydrogen flow rate to the consumed by the reaction in cells, should be maintained a constant to prevent fuel starvation at abrupt load changes. Moreover, the hydrogen pressure imposed to the stack should follow any change of the cathode pressure to prevent large pressure difference across thin membranes. Furthermore, liquid water, impurities and contaminant species in the anode gas flow channels should be purged out in time to prevent flooding and catalysts poisoning in cells. Design of model-based controls for the FDS is a challenging issue. A transient two-phase model of a single cell was developed for the design of controls, which considered the phase changes of water and two phase flows in cells. The model was experimentally validated by a segmented single cell that allows for measurements of current distributions and visualization of liquid water in gas flow channels. The experiment results of I-V curves shown that the air iii humidity in gas flow channels had larger influence on the cell performance than the air flow rates did. The images of liquid distribution in the channels indicated that most liquid water was accumulated near the outlet of gas flow channels and the amount of liquid water in the channels was affected by the air humidity and flow rates. The I-V curves and liquid water amount variation in channels have the similar trend with the simulation results of the transient model of the single cell. An anode model of a stack including two-phase phenomena was developed based on the transient model of the single cell, which was integrated to a set of control oriented models of FDS components. The integrated model was analyzed and linearized to develop a state feedback controller with integral and observer (SFB), which was compared with other two classic controls such as the proportional and integral (PI) and static feed-forward (SFF) controllers. It was found that the FDS could not be stabilized because of the liquid water accumulation in the system and cells without purging process. A dynamic purging process based on the time integral of stack current was designed and implemented to control the liquid water amount in the system. The simulation results of SFB, PI and SFF controllers with FDS model shown that the SFB controller had the best tracking and rejection performance on the control of the supplied hydrogen pressure and stoichiometric ratio under the disturbance of step change of stack current and purging process. In the simulation results, the liquid water was found in the anode side of fuel cells and manifolds in FDS. The amount of liquid was also effectively limited in a small range to prevent flooding in FDS and cells. iv Acknowledgments I would like to express my sincere gratitude to my advisor Dr. Song-Yul Choe, for his guidance, support and encouragement in all stages of my graduate program. I would also like to thank Dr. Sushil H. Bhavnani, Dr. Jeyhoon M. Khodadadi and Dr. Andrew B. Shelton for serving on my committee, and Dr. Steve R. Duke for his assistance as my university reader. Special thanks to Jiahua Li for the assistance in my experimental work. Thanks to my other colleagues in the Advanced Propulsion Research Lab for their help and knowledge, especially Jongwoo Ahn, Russell Green, Taylor Wingo, Matthew Thames and Kyle Malinowski. Special thanks also to Dr. Lewis N. Payton, Mr. John Marcell and Mr. Shannon Price for their assistance in machining parts, electrical work and operation of high performance computer cluster of engineering school. Thanks also to Zhi-Yuan (Jerry) Chen of Hephas Energy for his help on fabrication of fuel cell components. Finally, I would like to thank my parents and brother for their love, encouragement and support to my work and life. Thanks my parents-in-law for their time, energy and funding to help my family. I also thank my dear wife, Xinying Tian for her love, patience, encouragement, support and understanding throughout my whole graduate school experience. Special thanks to my sweet daughter, Caroline Ho, for her encouragement by smiling, lovely hugs and kisses every day. You are the sense of my life. v Table of Contents Abstract ........................................................................................................................................... ii Acknowledgments.......................................................................................................................... iv List of Tables ................................................................................................................................. ix List of Figures ................................................................................................................................. x Nomenclature ............................................................................................................................... xiv Chapter 1 Introduction .................................................................................................................... 1 1.1 Background ........................................................................................................................... 1 1.1.1 Working principle of PEM fuel cells ............................................................................. 3 1.1.2 Fuel delivery system ...................................................................................................... 7 1.2 Literature Review.................................................................................................................. 9 1.2.1 Review of the fuel cell modeling ................................................................................. 11 1.2.2 Review of FDS modeling and control ......................................................................... 13 1.3 Objectives and Dissertation Structure ................................................................................. 13 Chapter 2 A Two-Phase Transient Model of a PEM Fuel Cell .................................................... 15 2.1 Electrochemical Equations of PEM Fuel Cells ................................................................... 18 2.1.1 Open circuit voltage ..................................................................................................... 20 vi 2.1.2 Activation losses .......................................................................................................... 22 2.1.3 Ohmic losses ................................................................................................................ 24 2.2 One Dimensional Model of Current Densities.................................................................... 25 2.3 Coolant Channel.................................................................................................................. 28 2.4 Channel Plate ...................................................................................................................... 29 2.5 Gas Flow Channel ............................................................................................................... 32 2.6 Gas Diffusion Layer ............................................................................................................ 40 2.7 Catalyst Layer ..................................................................................................................... 46 2.8 Membrane ........................................................................................................................... 56 Chapter 3 Experiment Setup for the PEM Fuel Cell Model Validation ....................................... 62 3.1 Fuel Cell Assembly ............................................................................................................. 64 3.1.1 Cathode flow field for the liquid water observation .................................................... 65 3.1.2 Segmentation of anode flow field ................................................................................ 68 3.2 Air and Hydrogen Supply Systems ..................................................................................... 70 3.3 Instrumentation ................................................................................................................... 72 Chapter 4 Model Integration, Validation and Comparison ........................................................... 78 4.1 Development of the Transient Model ................................................................................. 78 4.2 Boundary Conditions of the Transient Model .................................................................... 81 4.3 Comparison between Simulation and Experimental Results .............................................. 83 4.3.1 Measurement of I-V curves by current sweeping ........................................................ 84 vii 4.3.2 Validation of the model and analysis ........................................................................... 97 Chapter 5 Modeling of Fuel Delivery System ............................................................................ 103 5.1 Manifolds .......................................................................................................................... 105 5.2 Blower ............................................................................................................................... 108 5.3 Ejector ............................................................................................................................... 112 5.3.1 Mathematical modeling ............................................................................................. 114 5.3.2 CFD validation ........................................................................................................... 120 5.4 Pressure Regulator ............................................................................................................ 127 5.5 Flow Control Valve........................................................................................................... 128 5.6 Purge Valve ....................................................................................................................... 129 Chapter 6 Fuel Delivery System Control .................................................................................... 131 6.1 Anode Model for Controller Design ................................................................................. 131 6.1.1 Gas flow channels ...................................................................................................... 133 6.1.2 Gas diffusion layers ................................................................................................... 136 6.1.3 Membrane .................................................................................................................. 140 6.2 Analysis of Integrated System of FDS and Stack ............................................................. 142 6.3 Design of State Feed-back Control with an Observer ...................................................... 146 Chapter 7 Analysis of Simulation Results of FDS with Feed-back Controllers ......................... 154 7.1 Response of Step Change Current .................................................................................... 156 7.2 Response of Multistep Change Current ............................................................................ 158 viii Chapter 8 Conclusion .................................................................................................................. 167 References ................................................................................................................................... 170 Appendix A Parameters of the transient model of the fuel cell ................................................. 177 ix List of Tables Table 3.1 Geometry of the designed single cell. .......................................................................... 65 Table 3.2 Errors of the sensors and data for validation. .............................................................. 72 Table 4.1 Operating conditions for RH experiments. .................................................................. 84 Table 4.2 Operating conditions for dynamic experiments. .......................................................... 99 Table 5.1 Blower map function parameters. .............................................................................. 110 Table 5.2 Comparison of ejector MATLAB model and CFD model ........................................ 125 Table 5.3 Constants for the LPR model. .................................................................................... 128 Table 6.1 Model parameters. ..................................................................................................... 143 x List of Figures Figure 1.1 Basic structure of a single fuel cell............................................................................... 3 Figure 1.2 Microstructure of the CL coated on the membrane. ..................................................... 5 Figure 1.3 Exploded view of a fuel cell stack (abridged to 4 cells for clarity). ............................. 6 Figure 1.4 Schematic diagram of a fuel cell system. ..................................................................... 7 Figure 1.5 Simple system configuration of a FDS without recirculation. ..................................... 8 Figure 1.6 FDS with a blower recirculation................................................................................... 9 Figure 1.7 A hybrid FDS. ............................................................................................................ 10 Figure 2.1 Structure diagram of a single fuel cell. ....................................................................... 16 Figure 2.2 Layers and state variables in half cell and membrane. ............................................... 18 Figure 2.3 A polarization curve to show the typical voltage losses of a PEM fuel cell. ............. 19 Figure 2.4 I-V curves of segments and a cell at time t. ................................................................ 27 Figure 2.5 The heat transfer between the CP with its neighbor channels and layer. ................... 29 Figure 2.6 The heat transfer between the GFC and CP. .............................................................. 30 Figure 2.7 Heat and mass transfer in GFC in x-y plane. .............................................................. 33 Figure 2.8 Mass and heat transfer in GDL along y-direction. ..................................................... 41 Figure 2.9 Heat and mass transfer of CL with GDL and membrane. .......................................... 47 Figure 2.10 Water phase change in CLs of PEM fuel cells. ........................................................ 48 Figure 2.11 Water content function for different temperatures [42, 49]. .................................... 51 xi Figure 2.12 Dissolved water transport and heat transfer between CL and membrane. ............... 56 Figure 2.13 Calculation flow chart of the cell voltage and current densities. ............................. 58 Figure 2.14 Flow chart of the mass and heat transfer for the transient model of the half cell. ... 59 Figure 2.15 Flow chart of the transient model of the PEM. ........................................................ 60 Figure 3.1 An assembled single fuel cell with segmented anode and transparent cathode. ........ 63 Figure 3.2 An exploded view of the single fuel cell. ................................................................... 64 Figure 3.3 An exploded view of the cathode part. ....................................................................... 66 Figure 3.4 The graphite plate engraved with GFCs. .................................................................... 67 Figure 3.5 A graphite block engraved with 5 GFCs. ................................................................... 68 Figure 3.6 An exploded view of anode part. ................................................................................ 69 Figure 3.7 Basic schematic of air supply system. ........................................................................ 70 Figure 3.8 Basic schematic of hydrogen supply system. ............................................................. 71 Figure 3.9 Diagram of cell instrumentation. ................................................................................ 73 Figure 3.10 Schematic diagram of the DAQ system. .................................................................. 74 Figure 3.11 The fuel cell and experiment setup pictures. ............................................................ 77 Figure 4.1 Subsystem of the two-phase transient model. ............................................................ 79 Figure 4.2 Function structures of FC_Update_Model. ................................................................ 80 Figure 4.3 Final transient model of the fuel cell in MATLAB/Simulink. ................................... 81 Figure 4.4 I-V curves of the fuel cell for different RH of air at cathode GFC inlet. .................... 85 Figure 4.5 Current distribution along GFCs for different air RH. ............................................... 88 Figure 4.6 Liquid water distribution along GFCs for different inlet RH of air. .......................... 89 Figure 4.7 I-V curves for the cases of different flow rates of air and hydrogen. ......................... 90 Figure 4.8 I-V curves of 10 segments for the cases of different flow rates. ................................ 92 xii Figure 4.9 Liquid water distribution in GFCs for different gas flow rates. ................................. 94 Figure 4.10 Effects of RH on pressure drop along GFCs at cathode side. .................................. 95 Figure 4.11 Effects of air flow rates on pressure drop along GFCs. ........................................... 96 Figure 4.12 I-V curves of model predication and experiments for different air RH. ................... 98 Figure 4.13 Liquid water saturation in GFCs at cathode side of model predication. .................. 99 Figure 4.14 The cell current of the fuel cell used in the experiments and simulation. .............. 100 Figure 4.15 The cell voltage responses of experiment and simulation result. ........................... 101 Figure 5.1 Schematic diagram of FDS. ...................................................................................... 104 Figure 5.2 Blower map of the mass flow rate and pressure ratio. ............................................. 111 Figure 5.3 Basic structure of an ejector. .................................................................................... 113 Figure 5.4 Working mode of ejector. ......................................................................................... 117 Figure 5.5 Contours of static pressure in an ejector. .................................................................. 122 Figure 5.6 Contours of Mach number in an ejector. .................................................................. 122 Figure 5.7 Pressure variation at the centerline of the ejector. .................................................... 123 Figure 5.8 Mach number variation along the centerline of the ejector. ..................................... 124 Figure 5.9 Primary flow rate comparison between MATLAB and CFD models. ..................... 126 Figure 5.10 Secondary flow rate comparison between the MATLAB and CFD models. ......... 126 Figure 5.11 Pressure and flow rate variation curve. ................................................................... 127 Figure 5.12 Diagram of integrated model of FDS. .................................................................... 130 Figure 6.1 Schematic diagram for mass transport in half a cell of anodic side. ........................ 132 Figure 6.2 The diagram of integrated model of fuel cell stack. ................................................. 141 Figure 6.3 Control and state variables at steady state for different current densities. ............... 145 Figure 6.4 Measurable output at steady state for different current density. .............................. 148 xiii Figure 6.5 State feedback control with integral and observer. .................................................. 149 Figure 6.6 The state of purge valve. .......................................................................................... 152 Figure 7.1 SFF control of FDS. ................................................................................................. 154 Figure 7.2 PI control of FDS...................................................................................................... 155 Figure 7.3 SFB control of FDS. ................................................................................................. 156 Figure 7.4 Step change response of medium mode controller. .................................................. 157 Figure 7.5 Step change response of high mode controller. ........................................................ 158 Figure 7.6 Multi-step changes of stack current, cathode pressure and purging state. ............... 159 Figure 7.7 The supply manifold pressure response for the multi-step current change. ............. 161 Figure 7.8 The SRH2 response of multi-step current change. ..................................................... 162 Figure 7.9 The water activity responses in manifolds of multi-step change of stack current. ... 164 Figure 7.10 The water activity responses in fuel cells of multi-step change of stack current. .. 165 Figure 7.11 Diagram of the control system for a real FDS and fuel cell stack. ......................... 166 xiv Nomenclature a Species activity aw Water activity A Area (m2) c Specific heat capacity (J kg-1 K-1) C Gas concentration (kg m-3) dh Hydraulic diameter (m) D Diameter (m) or diffusion coefficient (m2 s-1) F Faraday constant (96,487 C mol-1) Gibbs free energy of formation per mole (J mol-1) Enthalpy of formation per mole (J mol-1) hl Liquid water transfer coefficient (kg m-2 s-1 Pa-1) hlv Latent heat of water phase change (J mol-1) hm Mass transfer coefficient (m s-1) H Height (m) i Local current density (A m-2) I Current density (A m-2) j Gas flux (mol m-2 s-1) or liquid flux (kg m-2 s-1) J Rotational inertia (kg m2) xv m Mass (kg) M Mach number or molecular weight (kg mol-1) ndrag Electro-osmotic drag coefficient Ncell Cell number p Pressure (Pa, or bar) q Heat flux (J m-2 s-1) Q Volume flux (SLPM) R Gas constant (J kg-1 K-1) Ri Volumetric current density (A m-3) rm Mass source (mol m-3 s-1) rT Heat source (J m-3 s-1) s Liquid saturation Entropy of formation per mole (J mol-1) S Reduced liquid saturation T Temperature (K) t Time (s) U Velocity (m s-1) V Volume (m3) or voltage (V) W Width (m) or Mass flow rate (kg s-1) Greek Symbols x Mole fraction y Mass fraction xvi ? Specific heat ratio ? Thickness (m) ? Porosity ? Overpotential (V) ? Water content ? Viscosity (kg s-1 m-1) ? Mass density (kg m-3) ? Electrical or ionic conductivity (S?m?1) ? Angular velocity (rad s-1) ? Torque (N m) Subscripts an Anode bl Blower bm Blower motor ca Cathode cl Catalyst layer con Condensation crit Critical ec Evaporation and condensation eff Effective ej Ejector em Ejector manifold xvii evp Evaporation f Fluid cv Flow control valve g Gas gfc Gas flow channel gdl Gas diffusion layer i Index of gas species im Immobilized in Inflow or inlet ion Ionomer l Liquid water lpr Low-pressure regulator max Maximum value oc Open circuit out Outflow or outlet p Primary pem Membrane purge Purge valve ref Reference rea Reactant rib Channel rib rm Return Manifold s Solid or Secondary xviii sat Saturation sm Supply manifold tp Two-phase v Vapor w Water Abbreviations BOP Balance of Plant CC Coolant Channel CCM Catalyst Coated Membrane CCP Coolant Channel Plate CFD Computational Fluid Dynamics CL Catalyst Layer CP Channel Plate DMFC Direct Methanol Fuel Cell FCS Fuel Cell System FCV Fuel Cell Vehicle FDS Fuel Delivery System GCP Gas Channel Plate GDL Gas Diffusion Layer GFC Gas Flow Channel HOR Hydrogen Oxidation Reaction xix ICE Internal Combustion Engine LQI Linear-Quadratic Integral LQR Linear-Quadratic Regulator LQG Linear-Quadratic Gaussian MCFC Molten Carbonate Fuel Cell MEA Membrane Electrode Assembly MIMO Multi-Input Multi-Output OCV Open Circuit Voltage ORR Oxygen Reduction Reaction PEM Polymer Electrolyte Membrane PEMFC Polymer Electrolyte Membrane Fuel Cell PI Proportional?Integral PID Proportional?Integral?Derivative RH Relative humidity SFB State Feed-back SFF Static Feed-forward SISO Single Input Single Output SOFC Solid Oxide Fuel Cell SR Stoichiometric Ratio 1 Chapter 1 Introduction 1.1 Background A fuel cell is an electro-chemical energy device that converts chemical energy of source fuel into electricity. Fuel cells are recognized as one of the most promising technologies that can potentially meet requirements of future energy resources, where emissions can be substantially reduced while conversion efficiency can be increased. Fuel cells can be applied to various power applications that include portable electronics, transportation, and stationary power. The power range for the applications is typically 50-100 W for consumer electronics, 1-5 kW for residential power and backup power generators, 50-125 kW for vehicles, and 1-200 MW or more for centralized stationary power generators [1]. The fuel cells can be classified into different technologies based on the types of electrolytes, polymer electrolyte membrane fuel cells (PEMFCs), solid oxide fuel cells (SOFCs), molten carbonate fuel cells (MCFCs), direct methanol fuel cells (DMFCs), and others. Among these various technologies, the PEMFC is widely considered as one of the promising candidates to replace the internal combustion engine (ICE) for future vehicles because of its relatively high efficiency, high power density, near to zero emissions, low working temperature, fast start-up and quiet operation. 2 However, the primary challenges for the commercialization of the fuel cell vehicles are cost, durability and hydrogen storage [1]. The PEM fuel cells should have the durability of 5000 hours at a cost of $30/kW that are comparative to ICEs [2]. Storage for hydrogen should be able to store the fuel that covers a 300-mile driving range at least across all vehicle platforms [1]. In PEM fuel cells, chemical reactions take place at particles deposited in catalyst layers (CLs). Oxidation reactions in anode sides separate hydrogen molecules and produce electrons and protons, while reduction reactions in cathode sides separate oxygen and combine with the protons migrated across membrane from anode to cathode and the electrons through outer circuit. Heat and water are produced as byproducts. Cells are connected in series or parallel and become a stack that delivers a higher power. In order to operate a fuel cell stack continuously and reliably, fuels should be supplied according to the demand of loads and produced heat and water should be removed and maintained properly. The corresponding systems are called balance of plant (BOP) that mainly include three subsystems, fuel delivery system (FDS) for supplying hydrogen, air supply system for supplying oxygen and thermal management system for rejection of produced heat [3]. The author?s major research has been focused on the FDS. Design of the FDS should consider different aspects, not only supply of hydrogen from a tank that usually contains a highly pressurized hydrogen, but also supply of the fuel to the stack with a flow rate following demand of load and a reduced pressure that should not exceed a specified limit of the fuel cell components. In addition, the exiting hydrogen from the stack should be reused to increases overall efficiency of fuel usage [4]. Purging cells is also required to remove water, impurities and contaminants present at anode side 3 1.1.1 Working principle of PEM fuel cells A schematic diagram of a basic structure of a single PEM fuel cell is shown in Figure 1.1, which consists of gas channel plates (GCPs), GDLs at cathode and anode sides, and the polymer electrolyte membrane (PEM) and CLs. When oxygen and hydrogen gases are supplied to a cell, the gases flow through gas flow channels (GFCs) and then some of them diffuse through the GDLs and reach the catalysts. Hydrogen on the anode side breaks into protons and electrons. This reaction is called hydrogen oxidation reaction (HOR). The electrons released from the HOR move via the anode GDL, GCP, and loads that connect the anode and cathode, and finally reach the CL at the cathode side. Conversely, protons mitigate across the membrane, and then combine with electrons and separated oxygen on the cathode side and produce water. This reaction at the CL at cathode side is called oxygen reduction reaction (ORR). Figure 1.1 Basic structure of a single fuel cell. Lo a d A n o d e g as chann el p late A n o d e g as dif f u sion la y er A n o d e ca t aly st lay er P o ly mer ele ctroly t e memb r ane Cat h o d e ca t aly st lay er Cat h o d e g as dif f u sion la y er Cat h o d e g as chann el p late H 2 O 2 H 2 O H 2 O H + H 2 O 2 A n o d e S ide Cat h o d e S ide 4 Since the GCPs serve to provide pathways for reactants, pattern of the channels should be designed carefully to distribute the reactant gas uniformly over the GDLs. In addition, the plate should be electrically, thermally high conductive and corrosion resistive since the plate contacts fuels, oxidants, impurities and water. Preferred materials are graphite, stainless steel, coated metal, polymer-carbon and polymer-metal plates [5]. GDLs are very thin layers made of porous materials either the carbon fiber paper or carbon cloth mixed with hydrophobic materials like Teflon. The thickness is on the order of 100?m. The GDL provides paths for transport of reactant gases from GFCs to catalysts, removal of water from the CL to GFCs, conduction of electrons and transfer of the heat from the membrane and catalysts to GCPs [6]. The polymer electrolyte membrane coated with catalysts at both surfaces is named as membrane electrode assembly (MEA). The thickness of the MEA is on the order of 10?m. The membrane is generally made of ionomers that allow for transport of protons while blocking electrons flow and being impermeable to the reactant gases. The most commonly used membrane materials for the PEM fuel cells are Nafion that has been developed by DuPont. The Nafion should be fully humidified to keep high proton conductivity. Catalysts in CL are fabricated by a mixture of the ionomers and platinum or ruthenium supported by carbons. The carbon-supported platinum (Pt/C) is widely accepted as catalysts, but the cost is still too high. Currently, minimum Pt loading reaches 0.15mg cm-2 [3] without losing any performance. The ionomers bind Pt/C particles, which forms agglomerates, as shown in the Figure 1.2. 5 Figure 1.2 Microstructure of the CL coated on the membrane. Chemical reactions in the cell take place at the surface of particles where the three-phase boundary of Pt/C, ionomers and reactant gases are present [7] as shown in Figure 1.2. Although there are some internal reactions and side reactions, the basic electro-chemical reactions include the HOR, ORR and overall reaction written into following three equations: Anode side HOR: -2 2e2HH ?? ? Cathode side ORR: O2H4e4HO 2-2 ??? ? Overall reaction: O2HO2H 222 ?? The electrons released by the half reaction of HOR enter into the carbon and then are conducted to GDLs, while the protons go into the ionomers in CL and migrate across the membrane to the cathode side. At the CL of cathode side, the electrons enter into the carbon from I o n o mer P t/C H 2 e H + w a t er fil m O 2 e H + H 2 O w a t er fil m GDL CL PEM CL GDL Pt Carb o n Pt Carb o n A n o d e Ca th o d e 6 the cathode GDLs and then combine with the protons and oxygen to form water at the surface of platinum particles where the three-phase boundaries are present. Generally, maximum power of a single cell is not enough to meet demand of loads. Therefore, multiple cells are connected in series, which is shown in the Figure 1.3. Figure 1.3 Exploded view of a fuel cell stack (abridged to 4 cells for clarity). The unit cell is constructed from components including a gasket, GDLs, a MEA, and bipolar plates. The bipolar plates also allow for connection of two adjacent cells, whereby electrons flow from the anode to the cathode side of neighboring cells. Coolant channels (CCs) and GFCs are 7 engraved on bipolar plates. The highly conductive current collector plates connect to terminals of the loads, and the end plates are designed to connect reactant gases and coolant lines. 1.1.2 Fuel delivery system A fuel cell system on vehicles consists of a stack and BOP, as shown in Figure 1.4, which involves the air supply system (red lines), FDS (green lines) and cooling system (blue lines). Figure 1.4 Schematic diagram of a fuel cell system. A typical air supply system consists of a blower or compressor and a humidifier, which is also a heat exchanger. The blower or compressor is used to pressurize the ambient air supplied to the stack. The humidifier (for example, a membrane type of humidifier) is used to humidify and heat the air by using the exhaust wet air from the outlet of the cathode of the stack. S H y drogen T a nk P re s s ure R e gul a t or Ej e c t or P urge V a l v e Fuel C e l l S t a c k P neum a t i c C ontr ol V a l v e M o t o r A m bi e nt A i r C om pre s s or E x haus t Ga s t o T a i l P i pe T her m ost a t C ool a nt B y pas s C ool a nt R e s e rv oi r R a di a t or & Fan C ool a nt P um p A node E x haus t C a t hode E x haus t M e m bra ne H um i di f i e r 8 A cooling system in Figure 1.4 consists of a coolant pump, a thermostat valve, a radiator and a reservoir. The thermostat valve can be advantageously used to bypass the coolant if the temperature is not very high, so that the stack temperature can be kept in an elevated level. There are different configurations of FDS [4]. A basic FDS is composed of a tank connected with a pressure regulator and a pneumatic control valve that serve to stabilize the hydrogen pressure and flow rate at the outlet of the hydrogen tank as shown in Figure 1.5. A purge value is opened and closed periodically to remove excessive liquid water and impurities at anode side. Figure 1.5 Simple system configuration of a FDS without recirculation. The system above shows an obvious drawback that the effective use of the fuel cannot be accomplished because of no pathway for the reuse of the unconsumed fuel. To improve the efficiency of fuel usage, the excessive hydrogen is delivered into the stack, and the preferred stoichiometric ratio (SR), defined as the ratio of hydrogen flow rate supplied to stack to that consumed by the reaction in cells, should be larger than 1. Therefore, a recirculation loop is added to recapture the exhaust hydrogen from the anode channel to supply line as shown in Figure 1.6, where the blower is used to control the mass flow rate of recirculation. S H y drogen T a nk P re s s ure R e gul a t or P urge V a l v e Fuel C e l l S t a c k P neum a t i c C ontr ol V a l v e A node E x haus t 9 Figure 1.6 FDS with a blower recirculation. Another FDS as shown in the Figure 1.4 uses an ejector as a pump for recirculation line that uses the relatively high pressure fuel of the supply line. In addition, the ejector offers a distinguishing advantage that the entrained ratio, indicating the ratio of recirculation flow rate to the supply flow rate, can be kept constant under specific operating conditions. 1.2 Literature Review As discussed in previous sections, requirements for FDS are as follows: [1] Supply sufficient hydrogen to the anode of the fuel cell stack to prevent any shortage of hydrogen when fuels are consumed at a suddenly load change and supplied with a delay. [2] It is also desirable for unconsumed hydrogen to be circulated to a supply line to increase the efficiency of fuel usage. [3] Excessive hydrogen supply can reduce the response time at an increased power demand and eases water removal. [4] Since the large pressure difference between anode and cathode can potentially lead to a damage of the thin membranes, the pressure of hydrogen at anodic side should follow the air pressure at the cathode side not to exceed allowed pressure limit. S H y drogen T a nk P re s s ure R e gul a t or P urge V a l v e Fuel C e l l S t a c k P neum a t i c C ontr ol V a l v e A node E x haus t 10 [5] The liquid at anode side may cause flooding and block channels. In addition, the impurities lead to contamination of catalysts. Both these excessive liquid water and impurities should be removed to prevent decrease of fuel cell performance. Thus, controls of the fuel flow rate, pressure and purging under a rapidly varying current requests are three important objectives that should be considered in the FDS. In this work, the author has proposed new FDS used on fuel cell vehicles as shown in Figure 1.7. This FDS is a hybrid system that is made of two supply and two recirculation lines. The supply lines are operated based on the load demand. At a relatively low load demand, the supply line with a low pressure regulator accounts for the supply of hydrogen. The other line with a duty-control solenoid valve is used to supply additional flow at a relatively high load demand. An ejector and a blower both serve to circulate the exhausted hydrogen at the high load demand. The ejector is a passive device, while the blower is actively used to control the recirculation flow rate. Figure 1.7 A hybrid FDS. S C om pre s s e d H y drogen T a nk H i gh P re s s ure R e gul a t or P urge V a l v e Fuel C e l l S t a c k A node E x haus t Low P re s s ure R e gul a t or S B l ow e r E j e c t orD ut y - C ontr ol S ol e noi d V a l v e H y drogen S t ora ge S y s t e m 11 In the hybrid FDS, operations of each controllable actuator, duty-control solenoid valve or blower, can affect the hydrogen pressure and flow rate at the inlet of the fuel cell stack simultaneously. In fact, the control targets of hydrogen flow rate, pressure and purging have coupling effects when manipulating these two actuators. 1.2.1 Review of the fuel cell modeling Analysis of the FDS in Figure 1.7 and optimization of controls need development of dynamic models of the fuel cell stack. Modeling is also a scientific tool to understand mechanisms and advance development of fuel cells. Reviews about the fuel cell modeling have been published by Weber and Newman [7], Yao et al. [8], Faghri and Guo [9], Djilali [10], Haraldsson [11], Wang [12], Cheddie and Munroe [13], Atilla B?y?ko?lu [14], Ma and Ingham [15], and Siegel [16]. Different model dimensionalities are discussed in these reviews including 0-D, 1-D across the membrane, 2-D with one of other two directions, and 3-D models. In most published articles, the computational fluid dynamics (CFD) methods are used to solve the models by utilizing the computational power of computers. Two-phase and transient models are discussed in these reviews as the challenging parts of fuel cell modeling. Transport of species including protons, electrons, water, heat and reactants in the cells are solved using CFD techniques that allows for calculation of the species in a high resolution [17- 29]. Models proposed by Chen [18], Dutta [19] and Maharudrayya [23] assumed the flow as single-phase and simplified the transport modeling, while two-phase phenomenon of water was considered by Berning [17], Sivertsen [26], Wang [27], Le and Zhou [22]. Most models using CFD techniques are based on the average properties of the two-phase flow because the model itself was derived from pseudo-one-phase equations in the same view of Weber and Newman [7]. 12 Some other authors proposed simple lumped or one-dimensional models. Pei et al. [30] proposed a model that is in an integration form and reflects a single-phase pressure gradient at the anode side. Chang et al. [31] and Baschuk et al. [32] calculated flow distribution and pressure drop in manifolds and flow channels by a hydraulic network based on continuity and momentum equations of single-phase flow. In comparison to the CFD approaches, the models allow for quick calculation of pressure without detailed calculations, but their applications were constrained because of neglected liquid water effects. Recently, Yi?s group [33] presented a one- dimensional two-phase model that considered flow distribution and pressure gradients. The model was used to investigate water management along the flow channel. However, the pressure drop along the channel was not considered. Rodatz et al. [34] introduced an integration form of two-phase pressure gradient model that was initially proposed by Argyropoulos et al. [35] for analysis of a direct methanol fuel cells (DMFC). This model used a homogenous two-phase theory that assumed the same velocity of the gas and liquid water. Most CFD models are embedded in commercial CFD software, which is very computationally intensive, while zero dimensional lumped models do not represent details of the fuel cells. Both of them cannot be used for design of controls and analysis. It is also true that most control-oriented lumped models do not take account into the liquid water effect in GFC and GDL. Particularly, the two-phase dynamic characteristics, especially in the GFCs, are neglected in most of the control oriented models. Therefore, a two-phase transient model of a fuel cell model is presented and coded in MATLAB/Simulink by M-code which is easy to be converted into a stack model for the control design and simulation of BOP and FDS. 13 1.2.2 Review of FDS modeling and control Review of recent publications on FDS shows that there are a few articles that discuss the FDS and the associated controls. Most researchers modeled the FDS as a part of the whole fuel cell system. A model proposed by Pukrushpan et al. [36] was for a simple configuration that considered a flow control valve and purging valve without a recirculation line and was used to design controls. Bao et al. [37] used an ejector as a recirculation pump, but the recirculation flow rate was not actively controlled. Karnik and Sun [38] built a FDS model that considered one supply line and one ejector recirculation line to improve water management in the fuel cell, where a flow control valve was used to supply hydrogen flow. However, three control objectives aforementioned could not be simultaneously met by using only one actuator. He et al. [39] proposed a controller for a configuration of the FDS that includes two recirculation pumps and were able to control the two pumps and other actuators, but did not consider liquid water effect and purging cells. 1.3 Objectives and Dissertation Structure The objectives of this research are to complete a dynamic model and design associate control strategies for the hybrid FDS. The basic structure of the dissertation is as follows: 1) Introduction This chapter includes the research background, motivation and objectives. 2) A two-phase transient model of a PEM fuel cell This chapter presents a quasi two-dimensional model that considers two-phase, non- isothermal and transient behavior. 3) Experiment setup for the PEM fuel cell model validation In this chapter, a single cell and an associated test station are designed for validation experiments for the two-phase transient model. 14 4) Model integration, validation and comparison In this chapter, static and dynamic experiment data of the designed single cell are collected and compared with simulation results. The characteristics are I-V curves, liquid water distribution and dynamic response of cell voltage. 5) Modeling of fuel delivery system This chapter includes static and dynamic models of FDS components, as shown in Figure 1.7 such as manifolds, valves, a blower and an ejector. A complete FDS model was developed by connecting these components 6) Fuel delivery system control In this chapter, formation of the control objectives, design of linear state feed-back controller controls along with an observer as well as gain optimization are described. 7) Analysis of simulation results of FDS with feed-back controllers In this chapter, simulation results of the FDS including state feedback control are compared with those of the classic controls such as the proportional and integral (PI) and static feed-forward (SFF) controllers and its tracking performance and disturbance rejection are analyzed. 8) Conclusion 15 Chapter 2 A Two-Phase Transient Model of a PEM Fuel Cell In this chapter, a one dimensional two-phase transient model of a single PEM fuel cell was developed, which was used to study steady state performance and transient responses of step changes of load current. The single cell model considers electrochemical reactions, species mass and heat transport as well as the two-phase phenomenon in the GFCs, GDL and CLs of at cathode and anode sides of cells. The schematic structure of the modeled single fuel cell is shown in Figure 2.1, which includes the cathode, anode and membrane. A half cell shown in the figure is composed of cooling channel plates (CCPs), a gas channel plate (GCP), a gas diffusion layer (GDL), a catalyst layer (CL), and the membrane that is located between the cathode and anode side. The coolant plate embeds several coolant channels (CCs). The reactant gases flow in the gas flow channel (GFC) that is engraved on the surface of the gas channel plate (GCP). Due to the symmetric structure of the fuel cell, a mathematical model for a half cell was constructed and used for both of cathode and anode side. 16 Figure 2.1 Structure diagram of a single fuel cell. The cell is divided into several layers in the x-y plane as shown in Figure 2.2. In each layer, a one-dimensional transient model along channel direction (z-direction) is built, where the state variables such as gas concentrations and temperature are updated at every calculation time step. The mass and heat transport considering two-phase phenomenon in each section and z-directions are considered in the transient model. Although separated CCPs between cells are used in some designs, it is assumed that the temperatures of CCP and the GFP are the same because of the relatively high heat conductivity of the solid materials. Thus, these two plates are regarded as a channel plate (CP) without temperature gradient. In the CCs, it is assumed that the coolant fluid is incompressible and its temperature is the only state variable, where the convective heat transfer between the coolant fluid and channel pate is the driven force to change the coolant temperature dynamically when the operating condition of the coolant at the channel inlet is constant. Cath od e CC P Cath od e G DL PEM A no de CL A no de G DL A no de G CP A no de CCP y z x Cath od e G CP Cath od e CL 17 In the GFC, concentrations of gas species, amount of liquid water and temperature are the state variables that were updated. The phase change of water affects the dynamic change of the amount of liquid water, reactant gas transport and heat transfer. The GDL for a half cell is divided into three sub-layers (GDL1, GDL2 and GDL3) evenly along y-direction to study the gradient of state variables along this direction. In each of these sub-layers, the gas species concentrations, temperature and amount of liquid water are the state variables. The mass transport of gas species is assumed to be mainly driven by the gradient of species concentration gradient (diffusion), while the liquid water transport is driven by the capillary pressure gradient. The heat conduction is considered as the main heat transport mechanism in these layers. Due to the condensation and evaporation of water, the mass and heat sources are included in the model of GDLs. In the CL, the gas species concentrations, liquid water amount, temperature and dissolved water amount in the ionomers are state variables. Phase changes among the water vapor, liquid water and dissolved water in the ionomers are considered in the model of CL. In the membrane, the dissolved water amount and temperature are the state variables, where transport of dissolved water is driven by back diffusion and electro-osmotic, and transport of heat is only by conduction. Half a single cell in the x-y plane is shown in Figure 2.2. The number of the layers for half a cell is 7 in addition to one layer for membrane, so that a total number of the layers for a single cell becomes 15. In z-direction, the cell is segmented to present effects in this direction. Since each of the layers in x-y plane can be regarded as a one dimension transient model along the z- direction, each of the state variables in the x-y plane becomes a n element state vector, where n 18 presents the total number of segments along z-direction. In following sections, governing equations for each section of both cathode and anode side of the cell are derived. Figure 2.2 Layers and state variables in half cell and membrane. 2.1 Electrochemical Equations of PEM Fuel Cells Performance of a fuel cell is measured by a polarization curve (I-V curve), as shown in Figure 2.3. The x-axis presents the current density in the active area of MEA, and y-axis does the cell voltage measured from the terminals. When the current density of fuel cells increases, the cell voltage drops because of different types of losses, as shown in Figure 2.3. The losses are caused by fuel crossover and internal T cf T CP C i , G FC s G FC T G FC C i , G D L 1 s G D L 1 T G D L 1 C i , G D L 2 s G D L 2 T G D L 2 C i , G D L 3 s G D L 3 T G D L 3 C i , CL s CL T CL ? CL ? P EM T P EM T cf y x 19 current, activation, resistances and mass transport. Consequently, the resulting terminal voltage is lower than the theoretical reversible voltage defined as ideal open circuit voltage. Figure 2.3 A polarization curve to show the typical voltage losses of a PEM fuel cell. The real open circuit voltage (OCV) is lower than the ideal one caused by of fuel crossover and internal currents. The crossover of fuels from anode toward the cathode takes place because of pores present in the membrane. In addition, internal currents can flow through membranes because the membrane cannot perfectly block electrons flow toward to cathode. T h e o r e t i ca l Reve r si b l e V o l ta g e ( No L o ss) Rap i d d r o p d u e a cti va ti o n l o sse s A l m o st l i n e a r d r o p d u e t o o h m i c l o sse s Rap i d d r o p a t h i g h cu r r e n t d e n si ty d u e t o m a ss t r a n sp o r t l o sse s Cur r e n t Den si ty Cel l V o l ta g e Op e n Ci r cu i t V o l ta g e ( OCV ) In i ti a l o p e n c i r cu i t d r o p d u e t o fu e l cr o sso ve r /i n te r n a l cu r r e n t l o sse s 20 The activation losses occur because activation energy is needed to start the initial reactions in fuel cells toward the formation of water and electricity. The activation losses increase apparently at low current density while it almost keeps constant for high current density if other conditions do not change. Ohmic losses are an important factor for performance of a cell at medium and high current densities since it increases almost linearly with the current density. The losses are mainly caused by transport resistances of electrons in electrodes and protons in the membrane. When the current becomes high, losses caused by mass transport is dominant because the reactant gases at reaction sites cannot be supplied sufficiently due to the faster consumption rate of reactant gases and much liquid water condensed in the channels and pores at high current densities and increased mass transport resistance in cells. 2.1.1 Open circuit voltage In the current model, the fuel crossover and internal current effects are neglected, and open circuit voltage (OCV) of fuel cell is assumed to be equal to the theoretical reversible voltage. The change of Gibbs free energy is used to calculate the theoretical reversible voltage of PEM fuel cells. For the reaction of PEM fuel cell: OHO21H 222 ?? The change of Gibbs free energy of this chemical process is determined by, ? ? ? ? ? ? 222 2 1 OfHfOHff gggg ???? (2-1) 21 where is the change of Gibbs free energy of this process per mole product, and is the Gibbs free energy of formation per mole. If the process is reversible, which means all change of Gibbs free energy is converted into electric energy that is a product of charge and terminal voltage. Because of the number of electrons produced in the HOR, the change of the Gibbs free energy becomes as follows: revf FVg 2??? (2-2) Thus, the theoretical reversible voltage is obtained by: FgV f rev 2??? (2-3) However the value of is affected by temperature and species pressures or states. The effects can be reformulated using the definition of activity as: 0P Pa? (2-4) where P is the partial pressure of a reactant, and P0 is the standard pressure (1 bar). Then the change of Gibbs free energy of a cell reaction can be calculated by: ?? ? ? ? ?? ? ? ? ?? OH OH ff a aaRTgg 2 22 2 1 0 ln?? (2-5) where is the change of Gibbs free energy at standard pressure. If the pressure of gas species of H2 and O2 are given in bar, and the H2O is assumed in liquid water, which means its activity is unit, the Eq. (2-3) becomes: 22 ?? ? ??? ? ????? ? ??? ? ???? 210210 2222 ln2ln22 OHr e vOH fr e v PPFRTVPPFRTFgV ? (2-6) where is the theoretical reversible voltage at standard pressure, which is a function of temperature. Based on the definition of the change of Gibbs free energy, the is given by: F sThV ffrev 2 000 ?? ??? (2-7) where and are the changes of enthalpy and entropy at standard pressure, which are assumed to be constant for a working temperature range. Thus the change of has a linear relationship with respect to the change of temperature. By substituting the values of , , and = 1.229V at standard states (1bar and 25?C) into Eq. (2-6), the theoretical reversible voltage or open circuit voltage (OCV) of PEM fuel cells is given by [40]: ])(l n [4)15.298(1085.0229.1 23 22 HOOC PPF RTTV ????? ? (2-8) where T, PO2 and PH2 are the temperature (K) and partial pressures (bar) of oxygen and hydrogen at reaction sites. 2.1.2 Activation losses Activation losses are the main reason of initial rapid voltage losses in PEM fuel cells. These losses represent the activation energy to overcome an energy barrier and force the half reactions to proceed. These losses are termed as activation overpotentials. The activation overpotentials of the half reactions at cathode and anode have the relationship of local current densities approximated by the Butler-Volmer (B-V) equation given as [28]: 23 ? ????? ?????? ??????????? ??? ? ??? RT FRTFCCisi cacacaan r e fO CLOcaCLcaca ???? e x pe x p)1( , ,,05.2, 2 2 (2-9) ? ????? ?????? ??????????? ??? ? ??? RT FRTFCCisi ancaanan r e fH CLHanCLanan ???? e x pe x p)1( , ,,05.2, 2 2 (2-10) where the ica and ian are the current density at cathode side and anode side, i0 is the exchange current density, ? is the activation overpotentials, CO2, CH2 and T are the oxygen, hydrogen concentrations and temperature in the CL, ?ca is the cathode activation overpotential, ?ca and ?an are the asymmetry parameters, and s is the liquid water saturation in the CL, defined as the liquid water volume fraction to the volume of void space in the porous layer. The term of liquid saturation accounts for the effect of liquid water on the performance of fuel cell in CL. The temperature effect on the exchange current density is taken account by [41]: ??? ? ??? ? ?? ? ??? ? ? ?? TTREii r e fr e f 11e x p,00 ? (2-11) where i0,ref is the exchange current density at reference temperature Tref, and ?E is the activation energy at reference state. The asymmetry parameters are typically assumed to be 0.5, and the Eqs. (2-9) and (2-10) can be rearranged into a general form as: ? ???????? ??? ? ??? RTFCCisi r e fr e a CLr e aCL ??s i n h)1(2 ,,05.2 (2-12) In order to obtain the activation losses, the Eq. (2-12) needs to be rearranged to give ?: 24 ?? ? ??? ? ? ?? ?? 0, ,5.21 2)1(s i n h iiCCsFRT CLr e a r e fr e aCL?? (2-13) Therefore, the reactant gas concentrations, liquid water amount in CL and the current density will affect the activation overpotentials of half cell reactions. Based on the properties of hyperbolic function, the activation losses due to activation overpotentials increase more rapidly at low current density. At high current density, the reactant gas concentrations in CL decrease because of the increase resistance of mass transport, and liquid water amount in CL increases for the condensation of saturation water vapor especially at the cathode side. Both factors will increase the activation overpotentials based on Eq. (2-12) by which the mass transport losses are taken account in. 2.1.3 Ohmic losses In PEM fuel cells, there are also voltage losses caused by the transport of charges. The two types of charge particles are electrons and ions, both of whom lead to the voltage loss because of the electrical or ionic resistance of all fuel cell components including the electrolyte, the CL, the GDL, bipolar plates, interface contacts, and terminal connections. The reduction of voltage due to the resistance for charge transport is called ohmic losses. The ohmic losses include the ionic loss (in membrane and CLs) and electronic losses in other components at cathode and anodes sides given by: ?? ? ??? ? ? ???? e l e canCLi o n anCLcaCLi o n caCLPEMi o nPEMo h m i c ri ,, ,,, ,, ? ?? ?? ?? (2-14) 25 where ?ohmic is the ohmic losses, i is the current density (A m-2), ? is the thickness of layers and ? is the conductivity of materials. The first term in the parentheses at right hand side is the ionic specific area resistance (? m-2) of proton migration in membrane, the second and third terms is the ionic specific area resistance in CLs at cathode and anode side, and the final term is the total electronic specific resistance in other components which are assumed to be constant. In Eq. (2-14), ?ion is the ionic conductivity in the membrane and CLs, which is determined by the ionomer volume fraction (?ion), water content (?) that is used to describe the dissolved water amount in the ionomers, and temperature as [42]: ? ? ? ????? ?????? ??? Ti o n 130311 2 6 8e x p326.05 1 3 9.05.1 ??? (2-15) Hence, the cell voltage of a fuel cell can be given by its open circuit voltage, cathode and anode activation overpotentials, and the ohmic losses as: oh m i cancaOCc e l l VV ??? ???? (2-16) where VOC is the open circuit voltage, ?ca and ?an are the over-potentials at cathode and anode sides, and ?ohmic is the total ohmic losses in the fuel cell. The activation overpotential at anode side can be neglected comparing that at cathode side because of the fast reaction rate, low activation energy and high exchange current density for the half reaction at anode side. 2.2 One Dimensional Model of Current Densities For analysis of fuel cell systems, stack current is the input variable, while terminal cell voltage, power and current density distributions are the interested outputs. As explained in the previous section, the terminal voltage primarily is a function of the open circuit voltage, cathode 26 activation overpotential and ohmic losses that are determined by local current densities and species states such as reactant gas concentrations, temperatures and water amount. Since the local current densities are unknown and the cell voltage cannot be directly calculated by the given total current directly, a numerical method based on interpolation and iteration is applied to obtain the local current density distributions along the z-direction from the given total cell current. From Eq. (2-13), when the all state variables in the CL of cathode side are known at time t, the activation overpotentials of the CL at the cathode side is only a function of the local current density. Thus, the local cell voltage is a function of the local current density given by Eq. (2-16) when all the state variables of the cell are known, and the local I-V curves of all segments along z-direction of the cell can be drawn as depicted in Figure 2.4. If the cell is divided into n segments along z-direction, there are n local I-V curves obtained from the local cell voltage as shown in Figure 2.4. Due to high electronic conductivities of electrodes at cathodes and anodes, the local voltages between both GCPs for all segments are all equal to the cell voltage. The local current density of each segment for a given cell voltage is obtained by interpolating the local I-V curve as shown in Figure 2.4. The total cell current can is obtained by integrating the local current density along z-direction as: ?? ?? n kCLCLc e l l idzWi d zWI (2-17) where the Icell is the total current of the fuel cell, i is the local area current density, and WCL is the width of the CL. By changing the total cell voltage and integrating the current densities along z- direction, a new I-V curve of total cell voltage to total cell current can be drawn as shown in 27 Figure 2.4. Then the cell voltage and local current density distribution can be obtained by interpolating the total I-V curve and the segment I-V curves for a given cell current. Figure 2.4 I-V curves of segments and a cell at time t. In the following sections, the models of the each layers of the segmented fuel cell are proposed to solve all the state variables based on the mass and energy conservation principles. I V I V I V In te g rate th e cu rr e n t d e n sities o f a ll se g m e n ts i 1 i 2 i n I ce l l V c e l l z 28 2.3 Coolant Channel All the CCs are assumed to have the same temperature distribution and inlet operating conditions. The coolant fluid is incompressible, and the outlet flow rate is equal to the inlet one. There is only one state variable, coolant fluid temperature that is governed by: CCh CCCPcfcfcfcfcfcf dqzTcjtTc , ,4???????? (2-18) where ?cf, ccf, jcf,in and Tcf,in are the density, specific heat capacity, mass flux and temperature of coolant fluid, respectively. The dh,CC and qCP,CC are the hydraulic diameter of the CCs and the surface heat flux from the wall of the GCP to the fluid in channels. The hydraulic diameter of the channel for a rectangular section is calculated by its width (W) and height (H) of cross section as [43]: HWWHd h ?? 2 (2-19) The heat flux from CP to CC is obtained by the convective heat transfer equation of laminar internal flows as: )( , ,,, CCCP CCh CCfCCfCCCP TTdkNuq ?? (2-20) where Nuf,CC and kf,cc are the Nusselt number and heat conductivity of the CC respectively. For the channels with rectangular section, the Nusselt number can be interpolated by the aspect ratio of the channels based on the experiment data as the table shown in [44]. 29 2.4 Channel Plate Convective heat transfers between the CP and the channels (GFCs and CCs) and conductive heat transfer at the contacting interface between GDL and CP are considered in the channel plate model as shown in Figure 2.5. Natural convection is neglected in the current model. Figure 2.5 The heat transfer between the CP with its neighbor channels and layer. Since there is no mass transfer in channel plate, the temperature of CP is only the state variable considered: CP CCCCCCCPCC CP G F CG F CCPG F Ctp CP r i bCPG D LCPCPCPCPCP A HWqnA HWqAWqzTkztTc )(2)2( ,,,,1 ??????????? ???????? (2-21) The first term at right hand side is the heat flux gradient along z-direction in CP, the second term is the heat transfer from GDL to CPs, the third term is that from GFC to CP, and the final T CP q CP , CC CC CC CP GDL q CP , CC q G FC , CP q GDL , CP q GDL , CP G FC 30 term is that taken away by the coolant fluid, where the ncc is the number of the CCs for the fuel cell. The Wrib is the rib width of the CP contacting the GDL, and the ACP is the area of the cross section of the channel plate in x-y plane. The heat flux from the GDL to the CP is also governed by the heat conduction equation as: )(2 111 1,1 CPGDLGDLCPCPGDL CPGDLCPGDL TTkk kkq ??? ?? (2-22) Due to the presence of water vapor in GFC, condensation and evaporation may occur on the wall of the channel. It is assumed that the latent heat released by the condensation or absorbed by the evaporation only conducts to or from CP because of the low heat conductivity of the gas in the channels. Hence, the heat flux from GFC to CP are caused by the gas convective heat transfer from channel to plate and that of phase change on the channel wall as follows: CPG F CecCPG F CgCPG F Ctp qqq ,,,,,, ?? , where )(,,, CPG F CG F CgCPG F Cg TThq ?? (2-23) Figure 2.6 The heat transfer between the GFC and CP. T CP CP h G F C ( T G F C - T CP ) G FC q e c , G F C , C P 31 where qsp,GFC,CP is the convective heat flux from GFC to CP because the gas flow, qec,GFC,CP is the heat flux caused by condensation and evaporation of water in GFC, and hg,GFC is the convective heat transfer coefficient of the gas flow in the channel given by: G F ChnG F C G F CgG F CgG F Cg ds kNuh s , ,,, 1)1( ?? (2-24) where kg,GFC is the thermal conductivity of the gas flow in the channel, sGFC is the liquid water saturation defined as the volume fraction of liquid water in a control volume. The Nusselt number, Nug,GFC, is a constant related to the aspect ratio of the GFC. The liquid water amount in gas flow channel also can affect the heat transfer between gas flow channels and channel walls. It is assumed that the heat transfer coefficient is decreased with the increase of water amount by a power function of void space volume (1-sGFC). The heat flux caused by water phase change on the GFC wall is calculated by: CPG F CecwlvCPG F Cec jhq ,,,,, ? (2-25) where the jw,ec,GFC,GFP is the molar flux of the water condensed on the channel wall of CP, and hlv is the latent heat caused by the phase change that is expressed as a function of the CP temperature as [45], 41036 23 )15.273(1098.8)15.273(1054.2 )15.273(1044.3)15.273(9.4145070 ?????? ?????? ?? ? TT TTh lv (2-26) The molar flux of water phase change due to condensation or evaporation on the wall surface of the CP is calculated by: 32 ???? ? ?? ???? s a tvw a l ls a tG FCve v pw s a tvw a l ls a tG FCvvc o nww a l lecw CCifCCsk CCifCCsxkj ),( ),)(1( ,,, ,,,,, (2-27) where kw,con and kw,evp are the condensation rate and evaporation rate at wall, and Cv,GFC and Csat,wall are the vapor concentration in the channel and the saturated concentration near the channel wall, respectively, calculated as: GFC GFCvGFCv RTpC ,, ? , w al l w al lsatw al lsat RTTpC ][, ? (2-28) From the above equation, the saturation density near the channel wall is determined by the wall temperature, Twall, that is equal to the temperature of the CP for the condensation or evaporation on the GFC wall of the CP. The psat is the saturation pressure near the channel wall and is a function of the wall temperature [42]: ? ?? ? ? ? 372510 15.273104454.115.273101837.9 15.2730 2 9 5 3.08206.2])[(l o g ?????? ??? ?? TT TTp s a t (2-29) when Eq. (2-27) is used to calculate the condensation flux on the wall of GFC on the CP, the temperature of CP, TCP, can be considered as the wall temperature. 2.5 Gas Flow Channel Water vapor transported from the MEA or supplied with reactant gases may be condensed in GFCs by variation of concentration or temperature. In addition, liquid water may be transported from the GDL to the GFC by capillary force. When water droplets appear and coalesce in the 33 channels, they finally form a film, slug or annular flow in the GFC dependent on the gas flow velocities. Thus, two-phase flow should be considered in the modeling of GFC. Figure 2.7 Heat and mass transfer in GFC in x-y plane. As shown in the Figure 2.7, the gas species diffuse to or from the GDL, and water vapor concentration is affected by phase change process on the channel walls as well as surface of GDL. For the gas species in the GFC, gas concentrations are governed by: G F CimG F CGDLG F CiG F CziG F CiG F C rHjzjtCs ,,1,,,,,)1( ????????? (2-30) where Ci,GFC is the molar concentration of gas species i, sGFC is the liquid water saturation in GFCs. The ri,GFC is the source term of species in the channels. The ji,GFC is the superficial molar fluxes defined and calculated as: CC CC CP GDL GFC T CP C i ,G F C , s G F C , T G F C h GFC ( T GFC - T CP ) h GFC ( T G D L - T GFC ) He a t f l ux Gas sp e ci e s f l u x L i q u i d f l u x 34 G F CgG F CiG F CG F CziG F Czi UCANj ,,,,,, ?? (2-31) where Ni,GFC and Ug,GFC are the molar flow rate of gas species i and gas flow superficial velocity in the GFC along z-direction, and AGFC is the area of cross section of GFC. The ji,GFC,GD1 is the superficial molar flux of gas species i from GFC to GDL1. For the fuel cell in the validation experiments, the Reynolds number at the inlet of GFC is estimated in the range of 50~400 at cathode side, and no more than 40 at anode side for SR in the range of 1~10. Hence, the gas flow in GFC is assumed to be laminar flow, and the superficial velocity of the gas flow, Ug,GFC in Eq. (2-31) is assumed to be proportional to the pressure gradient. In fact, the calculation of pressure gradient of two-phase flow in small or mini-channels is still unknown because of change of flow pattern in the channels and the associated complex flow. Thus, the superficial species velocity of two-phase flow in the GFC is simplified as: zpskU G F C G F CgG F CG F CG F Cg ? ??? ,, 1][ ? (2-32) where kGFC[sGFC] is a flow coefficient that is a function of the liquid saturation in the channel. When the pressure gradient is constant, an increasing liquid water amount in the GFC decreases the superficial flow velocity. When the superficial flow rate is constant, the pressure gradient (absolute value) or pressure drop along the channels tends to follow the increase of liquid water amount. In this work, a power function is assumed for the flow coefficient function as: 2)1()( , snG F CG F CgG F CG F C sksk ?? (2-33) 35 where ksp,GFC is the flow coefficient for the single phase flow (gas flow) when there is no liquid water in the channels that is a function of the channel geometric shape and wall condition such as roughness. The gas viscosity in Eq. (2-32) is determined based on species viscosities and molar fractions by the Wilke mixture viscosity equation [46]and results in: ? ? ? ?? ?? ?? ? 5.0 225.05.0 3,2,1 3,2,1 , /12/4 //1, ji ijji iji j jij iicag MM MMw h e r exx ???? ? ? ? ? ??? ? ?? , (2-34) where 1,2 and 3 refer to oxygen, nitrogen and water vapor in cathode GFC. The viscosities of the species are a function of temperature using the power law [47], 69.05 2 7 3109 1 9.12 ???????? ? TO? , 67.05 2 7 3106 6 3.12 ???????? ? TN? , 15.15 3501012.1 ???????? ? Tv? (2-35) Likewise, the gas mixture viscosity for the flows in the anode GFC was also calculated. The total pressure of the gas species in GFC in Eq. (2-32) is calculated by the state equation: ?? i G F CuG F CiG F C TRCp , (2-36) In Eq. (2-30), the gas species superficial molar flux ji,GFC,GDL1 is determined by the concentration difference between the GFC and GDL (GDL1). Since the liquid water may appear on the surface of the GDL in GFC by liquid transport or condensation, part of the surface is covered by the liquid water that blocks the gas species transport path in GDL. Thus, the liquid water effect for the gas transport on the surface of GDL is considered in the calculation of ji,GFC,GDL1 as: 36 ))(1( 1,,c o v,1,,1,, GDLiG F CilGDLG F CmGDLG F Ci CCrhj ??? (2-37) where hm,GFC,GDL1 is the mass transfer coefficient between the GFC and GDL (GDL1), and rl,cov is the area ratio of surface covered by liquid water to total surface of GDL in the GFC. The mass transfer coefficient in Eq. (2-37) is determined by the coefficients in GFC and GDL1, respectively, as: 1,, 1,,1,, GDLmG F Cm GDLmG F CmGDLG F Cm hh hhh ?? (2-38) It is assumed the gas species is diffused from the GFC to the surface of GDL. In GFC, the mass transfer coefficient is obtained by: G F Ch G F CeffiG F CG F Cm d DShh , ,,, ? (2-39) where ShGFC is the Sherwood number of GFC, and Di,eff,GFC is the effective diffusion coefficient of gas species i in GFC. Another mass transfer coefficient for the gas species is that in GDL that is derived by the Fick?s law [47] in porous media as: 1 1,,1, 2 GDL GDLeffiGDLm Dh ?? (2-40) In Eq. (2-39) and (2-40), the effective diffusion coefficients of gas species i in the GFC and GDL1 are calculated by [12]: 37 ? ? ?? ij e f fji jie f fi D xxD ,, , 1)1( where ? ? ?? ? ??? ? ? ??? ? 1 ,, 5.25.1,, 111 ikjie f fji DD sD ? (2-41) where xi is the molar fraction of gas species in gas phase, and Di,j,eff is the effective binary diffusion coefficient of gas species i that is a function of porosity (?) of porous materials, liquid water saturation (s) in void space, binary diffusion coefficient Di,j and Knudsen diffusion coefficient Dk,i [12]. These two diffusion coefficients can be calculated by [12]: 5.1 0 0,,0, ????????????????? TTppDD jiji , 2 1 , 83 ?? ? ??? ? ?? i p o r eik MRTdD ? (2-42) where D0,i,j is the reference binary diffusion coefficient at the reference pressure p0 and temperature T0. The Knudsen diffusion effect is only considered in small size of pores. Because of no porous structure in the GFC, the effect is neglected in the GFC. The porosity and the liquid saturation, s, in Eq. (2-41) are set to 1 and 0 in order to calculate the effective diffusion coefficient of gas species in GFC. The area ratio of GDL surface is covered by liquid water in GFC and affected by the water contact angle on the surface, flow pattern, liquid water amount, and others. This phenomena is simplified using an area ratio, rl,cov that is a power function of liquid saturation in GFC as: ),1m in ( 3co v,co v, snG F Cll sKr ? (2-43) where Kl,cov is a coefficient, and ns3 is assumed to be a constant in current model. 38 In GFC, the mass source term is the liquid water that results from phase change of water vapor on the walls of CP and surface of GDL: G F C G D LG F Cecw G F CG F C G F CG F CCPG F CecwG F CecwG F Cv HjWH WHjrr ,,,,,,,,, )2( ?????? (2-44) where the first term of the right hand side is the part from condensation molar flux on the CP wall, and the second term is that on GDL wall. jw,ec,GFC,CP and jw,ec,GFC,GDL are condensation flux in channels that are calculated by Eq. (2-27), where the wall temperature, Twall, is substituted by TCP and TGDL1, respectively. Liquid water present in the GFC comes from condensed water vapor in the GFC or carried by air and transport from the GDL by capillary force. Change of the amount of liquid water is expressed using the liquid saturation s: G F CecwOHG F CG F CG D LlG F CzlG F Cl rMHjzjts ,,,1,,, 2????????? (2-45) where jl,z,GFC is the liquid water superficial mass flux along GFC that is assumed to be entrained by the viscous force of gas flow and calculated by: G F CglG F CgG F CG F CltpG F Czl Us sCj ,, 2 ,, 1 ? ?? ?? ? ??? ? ? ?? (2-46) The Ctp is the correction factor that is as a constant in current model. In the second term of right hand side of Eq. (2-45), jl,GDL1,GFC is the liquid water mass flux from GDL1 to GFC driven by capillary force. When the capillary force in the GDL1 becomes large, the liquid water leaves the GDL1 and enters the GFC. Because of high hydrophobicity of 39 the GDL materials, it is assumed that the liquid water cannot enter the GDL from GFC. Based on the Darcy?s law in porous materials, the liquid water flux is calculated by: ? ? G D Lc r i tcG D LclG D L G D LrG D LlG G FG D Ll pp KKj ,,,1 ,,1, 2 ??? ??? (2-47) where KGDL is the absolute permeability of the GDL, Kr,GDL and pc,gdl are the relative permeability and capillary pressure of GDL that are a function of liquid saturation in the GDL as [20]: 3)( imr ssK ?? (2-48) ? ? ? ? ? ??? ? ??? ? ????????? ???????????? ? ? ? 90263.112.2417.1cos 90)1(263.1)1(12.2)1(417.1cos 32 5.0 32 5.0 ccw ccw c SSSK SSSK Sp ???? ???? (2-49) where the s is the liquid saturation, sim is the immobilized liquid saturation, and S is the reduced liquid saturation in GDL defined as [20]: ?? ??? ?? ????? im imim ss sssssS 0,0 1,1 (2-50) The last variable in GFC is the gas temperature, Tg,GFC. Based on the assumption that the phase change heat flux only is transferred from or to the CP, there is no phase change heat term in the energy equation of gas mixture in GFC as follows: 40 G F C G F CG D L G F CG F C G F CG F CCPG F CgG F CzG F Cg i ivi H qHW HWqzqtTsCC , ,,,,, 2)1( ???? ??????? (2-51) where Cv,i is the specific heat capacity at constant volume of gas species i in the GFC, and Tg,GFC is the gas temperature in GFC. At the right hand side of Eq. (2-51), the first term is the gradient of heat flux along z-direction given by: ?? i G F CgipG F CziG F Cz TCjq ,,,,, (2-52) where ji,z,GFC and Cp,i are the molar flux and specific heat capacity at constant pressure, respectively, of the gas species in the GFC. The heat flux from GFC to CP in the second term at right hand side of Eq. (2-51) is calculated by Eq. (2-23), and the last term is the volumetric enthalpy change rate caused from the heat transfer from GDL to GFC is dominated by conduction by neglecting the heat transfer of gas diffusion and liquid water transport. Thus, the heat flux from GDL to GFC is: )(,, G F CGDLG F CgG F CGDL TThq ?? (2-53) where hg,GFC is the convective heat transfer coefficient calculated by Eq. (2-24). 2.6 Gas Diffusion Layer In the GDL, the pressure difference and the convective mass transfer are negligible. The gas species transfer and the liquid water are driven by the gradient of concentration and capillary force, respectively. As shown in the Figure 2.2, the GDL is divided into three layers (GDL1, GDL2 and GDL3). The gas concentration for each layer is described as follows: 41 imyo u tiyinizii r jjzjtCs ,,,,,,)1( ? ????????? ?? (2-54) where ? is the porosity of the GDL, s is the liquid water saturation, ji,z is the molar flux along z- direction. The ji,in,y and ji,out,y is the inflow and outflow molar flux of gas species i at the interface for each layers (GDL1, GDL2 and GDL3) as shown in Figure 2.8. The last term at the right hand side is the mass source term that is only applied for vapor concentration. Figure 2.8 Mass and heat transfer in GDL along y-direction. In each of the three layers of GDL, the z-direction gas species molar diffusion flux ji,z in above equation is driven by the gradient of concentration gradient as: zCDj i effizi ???? ,, (2-55) C i , G D L 1 s G D L 1 T G D L 1 C i , G D L 2 s G D L 2 T G D L 2 C i , G D L 3 s G D L 3 T G D L 3 He a t f l u x Gas sp e ci e s f l u x L i q u i d f l u x GDL1 GDL2 GDL3 G FC CP 42 where the effective diffusion coefficient geffiD, is calculated by the Eq. (2-41) including the Knudsen diffusion effect. The interface gas species of GDL1, 2 and 3 are shown in Figure 2.8. For the GDL1, the inflow gas species flux at its boundary is that from the GFC to GDL1. Since the ribs between GFCs on CP block a part area of GDL that transfers the gas species, the ji,in,y of GDL1 in Eq. (2-54) has a relationship with ji,GFC,GDL1 by Eq. (2-37) as: G D L G FCG D LG FCiG D Lyini W Wjj 1,,1,,, ? (2-56) where WGFC is the width of GFC, and WGDL1 is the width of the GDL that is the same as the total width of the fuel cell. The outflow gas species flux is that from GDL1 to GDL2, ji,GDL1, GDL2, and is also the inflow gas species flux of GDL2. Similarly, the outflow gas species ji,GDL2,GDL3 is also the inflow gas species flux of GDL3. The molar flux of gas species between the GDL layers is calculated by the Fick?s law. For example, the molar flux of species i between the GD1 and GDL2 is given by: ? ?2,1,12,,12, iiimi CChj ?? , where 12,,21,, 2,,1,,12,, 2 ?? e f fie f fi e f fie f fiim DD DDh ?? (2-57) where Di,eff is the effective diffusion coefficients if corresponding GDL layer, hm,i,12 is the mass transfer coefficient, and ? is the thickness of the layer. Since the CL has porosity, the equation above is used to calculate the outflow gas species flux of GDL3 ji,GDL3,CL as: ? ? CLiG D LiG D LCLe f fiCLG D Le f fi CLe f fiG D Le f fiCLG D Li CCDD DDj ,3,3,,3,, ,,3,,,3, 2 ??? ?? (2-58) 43 In GDL, the source term caused by condensation and evaporation in the pores is determined by [48]: ? ? ? ??? ? ?? ??? ? s a tvs a tv OH l e v p s a tvs a tv u v c o n ecw ppifppMsk ppifppTR xsk r ),( ),()1( 2 , ? ? ? (2-59) where rw,ec is the water condensation source term, the vapor source term rm,v = -rw,ec in GDL. The condensation rate kcon and evaporation rate kevp are assumed to be constant. The liquid water change in GDL considering the capillary force and water phase change is as follows: ecwyo u tlyinlzll rjjzjts ,,,,,, ????????? ??? (2-60) where jl,z is the z-direction liquid water mass flux calculated by: zSSDj capzl ??? )(, , where ? ????? ???? SpKKSD cl r e la b slc a p ??)( (2-61) where S is the reduced liquid saturation in GDL, Dcap is the capillary liquid water diffusion coefficient that is a function of reduced liquid water saturation. The liquid water fluxes across the boundaries of GDL1, 2 and 3 are also shown in the Figure 2.8. For the GDL1, the jl,out,y is the flux from GDL1 to GFC, which is calculated using the jl,GDL,GFC in Eq. (2-47) as: 44 GDL G F CG F CGDLlGDLyo u tl W Wjj ,,1,,, ? (2-62) If the liquid water flux is also assumed to be driven by the capillary force between the three layers of GDL, the flux becomes as follows: )( 2112,12, SShj ll ?? , where 12,21, 2,1,12,, 2 ?? c a pc a p c a pc a pil DD DDh ?? (2-63) where the Dcap is a function of the reduced liquid saturation. At the interface of GDL3 and CL, it is assumed that the capillary pressure is continuous, while the liquid saturation is discontinuous because of the different liquid water contact angle, porosities and permeability in GDL and CL from the capillary pressure function of Eq. (2-49). Hence, the liquid water flux from CL to GDL is determined by: )( 3,,,,3,, GDLcCLcGDLCLlGDLCLl pphj ?? , where CLG D LpcG D LCLpc G D LpcCLpcil DD DDh ?? 3,3, 3,,12,, 2 ?? (2-64) where the hl,CL,GDL is the liquid water transfer coefficient from CL to GDL3 based on capillary pressure difference, and Dpc,CL is the capillary-pressure-based diffusion coefficient defined as: l relabslpc KKD ???? (2-65) The last term in Eq. (2-60) is the liquid water mass source term caused by condensation according to Eq. (2-59). Since gas, liquid and solid are present in the GDL, it is assumed that the three phase mass have the same temperature that is governed by: 45 ? ? Tyo u tyinzssllipi rqqzqtTccsCCs ????????????? ? ?????? ,,, )1()1( (2-66) where qz is the z-direction heat flux dominated by conduction as: zTkq z ???? (2-67) where the k is the average heat conductivity of three phases as: ))(()1( 2/1 glgs kkskkk ????? ?? , where ?? iig xkk (2-68) where ks, kg and kl is the heat conductivity of solid, gas and liquid phases, and kg is the average value of gas species conductivity based on molar fraction. In the right hand side of Eq. (2-66), the qin,y and qout,y are the input heat flux and output heat flux of each layer as shown in Figure 2.8. For GDL1, the input heat flux is that transferred between the sub layers of GDL, which is calculated by: )(2 121221 211,2 TTkk kkq ??? ?? (2-69) The output heat flux qout,y of GDL1 includes the heat from GDL1 to the GFC and GDL1 to CP, which is given by: G D L r i bCPG D LG F CG F CG D LecG F CG D LCPG F C W WqWqqq ,,,,,1 )( ???? (2-70) where qGDL,GFC and qGDL,CP are calculated by Eq. (2-53) and (2-22). 46 Since it is assumed that the latent heat due to the phase change of the liquid water in GFC on the surface of GDL comes from or goes to GDL1, the qec,GDL,GFC term is added in above equation, which is calculated by: GDLG F CecwlvG F CGDLec jhq ,,,,, ?? (2-71) where hlv is the latent heat of water, and jw,ec,GFC,GDL is the condensation flux on the surface of the GDL in GFC calculated by Eq. (2-27). Likewise, Eq. (2-69) is used to calculate the heat fluxes between GDL2, GDL3 and CL. The last term at right hand side of Eq. (2-66) is the heat source term caused by the phase change of water in GDL, which is obtained by the mass source term of water as: GDLecwlvGDLT rhr ,,, ? (2-72) where rw,ec,GDL is the water condensation mass source term calculated by Eq. (2-59). 2.7 Catalyst Layer The reactant gases (hydrogen at anode side and oxygen at cathode side) are consumed in CL by chemical reactions that produce water and heat as byproducts. In the modeling of CL, assumptions are made that no gas species can cross membranes except water that can be transported through the membrane in the form of dissolved water in the ionomers. A schematic diagram indicating the mass and heat transfer between CL and its neighbor layers is shown in Figure 2.9. 47 Figure 2.9 Heat and mass transfer of CL with GDL and membrane. Mechanism of water phase change and transport is shown in Figure 2.10. Water existing in CL can be in form of three phases, liquid water, water vapor and dissolved water in the ionomers. Once water is generated, transport to the ionomers in CL takes place and then to voids in form of vapor by desorption and adsorption processes and in form of liquid water by uptake and release. The amount of vapor and liquid water in CL can be described by vapor concentration (Cv) and liquid saturation (s), respectively. The amount of dissolved water in the ionomers is described by water content (?) that is defined as the ratio of the number of water molecules to the number of charge ( ) sites in the ionomers [42]. The dissolved water concentration in the ionomers is given by: EWC ionion ionw ????, (2-73) where ?ion is the volume fraction of the ionomers in CL, ?ion is the ionomer density (kg m-3) and EW is the equivalent weight (kg mol-1), which is the ionomer weight containing one mole in membrane. He a t f l u x Gas sp e ci e s f l u x L i q u i d w a t e r f l u x GD L3 C i , CL , s CL , T CL , ? CL Di sso l v e d w a t e r f l u x CL PEM 48 Hence, there are other source terms for the mass and heat conservation because of reactions and water phase changes. The concentration of gas species in the CL is given by: CLimCL CLGDLiCLziCLiCLCL rjzjtCs ,,,,,,,)1( ????????? ?? (2-74) where ji,z,CL is the molar flux of gas species along z-direction in the CL that is also calculated using Eq. (2-55) by considering porosity of the CL. The ji,GDL,CL is the mass flux from GDL3 to CL calculated by Eq. (2-58). Figure 2.10 Water phase change in CLs of PEM fuel cells. Since half a reaction takes place in cathode and anode CLs, the reactant species (hydrogen and oxygen) are consumed in the CLs, and two extra source terms are considered for the reactant gas concentration as: D is s olv ed w at er in ionom er W at er v apor Liquid w at er C ondensatio n Ev aporation 49 ? ? ? ??? ? ??? ??? s i d ea n o d eFiFRr s i d ec a t h o d eFiFRr CL anani CLHm CL cacai CLOm ? ? 22 ,44 , ,, , ,, 2 2 (2-75) where Ri (A m-3) is the volumetric current density in corresponding CL. As shown in Figure 2.10, the water vapor source term involves two parts as: CLadwCLecwCLv rrr ,,,,, ??? (2-76) where rw,ec,CL is the source term caused by the evaporation and, and rw,ad,CL is that caused by the adsorption and desorption between water vapor and dissolved water in the ionomers. The mass source term rw,ec,CL is calculated using Eq. (2-59) for water condensation and evaporation in porous. The source term due to the adsorption and desorption is given by [20]: ?? ??? ????? ??? ? 1)1)()(1( )1)(( ,, ,, , RHandifRHs EWk ifsCk r eqveqvM E Mde s eqveqvvad s adw ????? ???? (2-77) where RH is relative humidity (RH) of vapor, kads and kdes are the adsorption and desorption constant respectively. The equilibrium water content ?v,eq is determined by the water activity aw in the CL at different temperatures obtained from the experimental data as [42, 49]: ? ? ?? ??? ? ???? ???? ? 122 311814 1,0.3685.3981.1704.0 32 ]30[ w ww wwww C a aa aaaa ?? (2-78) 50 ? ? ?? ??? ? ???? ???? ? 18.16 311792.3216.9 1,116.14168.103.0 32 ]80[ w ww wwww C a aa aaaa ?? where aw is defined as the ratio total water amount of vapor and liquid water to that of saturated vapor in the void space. It is equal to the RH in CL for unsaturated vapor and larger than 1 when liquid water appears in pores of CL. Thus, the water activity is calculated by the vapor concentration, liquid saturation and temperature as [42]: ? ? )()1( 22 Tp TRssCMa s a t OHlvOHw ???? (2-79) Thus, the water content for other values of temperature is obtained by interpolating the function values above at given specific water activity as [50]: ? ? ][50 15.303][][],[ ]30[]30[]80[ wCwCwCw aTaaTa ??? ???? ???? (2-80) where T is the temperature (K) of the CL. The water content functions are shown in (6-2), where the curve of water content at 60 ?C is obtained by an interpolation from the curves of water content at 30 ?C and 80 ?C. From Eqs (2-78) and (2-80), the water content curve is divided into three phases as shown in Figure 2.11. The transition phase of the curve is used to remove the discontinuity of the water content curve at the saturation point (aw=1). The ?vmax and ?lmax are the water content values at the water activity aw=1and aw?3, respectively. 51 Figure 2.11 Water content function for different temperatures [42, 49]. For the liquid water phase in CL, the liquid water saturation equation is given by: CLurwCLecwCL CLGDLlCLzlCLlCL rrjzjts ,,,,,,,, ???????? ??? (2-81) where jl,z,CL is the liquid water flux along z-direction in CL calculated by Eq. (2-61), the rw,ur is the water source term caused by the liquid water uptake and release by the ionomers, which is 0 1 2 3 4 0 5 10 15 20 25 Water activity a w ? Water content at 30 ? C Water content at 80 ? C Water content at 60 ? C ? l m a x ? v m a x Single phase Two-phase Transition phase 52 considered as the mechanism of phase change between liquid water in pores of CL and the dissolved water in the ionomers of CL. Calculation of the water mass source term by liquid water uptake and release is performed based on a cluster-network model proposed by Weber and Newman[51]. When the liquid water condensed on the ionomer surface in CL, the equilibrium water content for the ionomers should be in the range from maxv? to maxl? as shown in Figure 2.11depending on fraction of expanded hydrophobic channels in the ionomers. The equilibrium water content under occurrence of liquid water is related to liquid water pressure in CL [51] as: ?? ? ? ? ? ? ?? ? ? ? ? ? ?? ? ? ? ? ? ?? ? ? ? ? ? ?? ??? ? ??? ? ? ???? ? 23.0 )1025.1l n (c o s2ln 1)(21 9 m a xm a xm a x , l i o nw vlveql pe r f ?? ???? (2-82) where erf is the error function, ?ion is the water contact angle on the ionomers and pl is the liquid water pressure in the hydrophobic channel that is given by: cgl ppp ?? (2-83) where pg is the gas pressure, and pc is the capillary pressure calculated by the Eq. (2-49). Thus, uptake of liquid water by the ionomers occurs if the liquid equilibrium water content in void space of CL is higher than that dissolved in the ionomers of CL, and the liquid water in the ionomers releases from the ionomers to the void space in the inverse case. The source term of uptake and release is given by [20]: 53 ????? ?? ??? eqleqlr eqleqluurw ifk ifkr ,, ,,, )( )( ???? ???? (2-84) where ku and kr is the water uptake and release rates, respectively. The analysis above shows that the amount change of water dissolved in the ionomers of CL will affected by the water absorption and water uptake process. Thus, change of the water content can be described using the dissolved water change in the ionomers of the CL: urwadwr e a c twCL PEMCLCLzCLi o ni o n rrrjzjtEW ,,,,,,, ?????????? ???? ?? (2-85) where the EW, ?ion and ?ion are the equivalent weight, the volume fraction and the density of the ionomers, respectively. The z-direction of the dissolved water flux in the ionomers is only forced by the gradient of dissolved water concentration in terms of the water content as: zDj ionw ???? ?? , , where io nio nio nio nw DEWD ,5.1, ???? (2-86) where D?,ion is given by [52]: ????? ?? ????? ??? ?? o t h e r w i s eee ifeeD T T i o n /2 4 3 68 /2 4 3 628.07 , )1161(1017.4 30)1(101.3 ? ? ? ? ?? (2-87) The second term in Eq. (2-85) presents the dissolved molar flux of the water between the ionomers in CL that has the same property as that of membrane. The water flux driven by diffusion and electro-osmotic force is: 54 Finhj d r a gPEMCLPEMCLPEMCL ?)(,,,, ???? ?? , (2-88) where h? is the water content transfer. In the second term at right hand side of above equation, the negative symbol is for cathode side while positive symbol is for anode side, and the ndrag is the electro-osmotic water drag coefficient, and i is the local area current density (A m-2). The water content transfer coefficient between CL and membrane is given by: P E MCL P E MCLP E MCL hh hhh ,, ,,,, ?? ??? ??? , where ][ ][,][, 2 ?? PEMwDh ? (2-89) where index [] represents the CL or membrane, PEMwD][, is the dissolved water diffusion coefficient calculated from Eq. (2-86), and ?[ ] is the thickness of the layer along y-direction. The volume fraction of the ionomers ?ion in membrane used in Eq. (2-86) is equal to 1. The electro-osmotic water drag coefficient is a function of water content at the interface between the CL and membrane layer as [51]: ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? o t h e r w i s e ?? ?? . ??if ?if? n v,l, v, v,d r a g 151 11 1 m a xm a x m a x m a x (2-90) The water content at the interface of CL and membrane is: CLM E M CLM E MM E MCLer ?? ????? ???i nt (2-91) 55 The last term of the Eq. (2-85) is the mass source terms of dissolved water produced by reaction, adsorption and desorption, and uptake and release. Since no water is generated in anode CL, the reaction term of dissolved water is zero for anode side while the term for cathode side is: CL cacair ea ctw FiFRr ?22 ,, ?? (2-92) where Ri,ca (A m-3) is the volumetric current density at the cathode side. The equation for temperature should consider the additional source term as: ? ? o h m i cTr e a c tTecT GDLCLCLPEMCLz ssllivi rrrqqzq t TccsCCs ,,, ,,, , )1()1( ????????? ? ????? ? ? ????? (2-93) where qz,CL is the z-direction heat flux calculated by Eq. (2-67), qPEM,CL is the heat flux from membrane to CL calculated by Eq. (2-69). The qCL,GDL is the input heat flux of GDL3. The equation above shows three heat source terms in the CL that are water phase change, reaction and Joule heating. The latent heat at a phase change in CL is calculated by Eq. (2-72). The rT,react is the heat source caused by reactions, and rT,ohmic is the Joule heating that are given by: ? ????? ?? ??NFSTRr ir eac tT , (2-94) 2, 2 , CLiohmicT ir ?? (2-95) 56 where Ri is the volumetric current density (A m-3), ?S is the entropy change of half reaction at cathode or anode, T is the temperature in the CL reaction, ? is the activation overpotential of the CL, and N is the number of electrons per one mole reaction, 4 for cathode side and 2 for anode side. i is the local area current density (A m-2). 2.8 Membrane A simple schematic diagram for water transport and heat transfer is depicted in Figure 2.12. Figure 2.12 Dissolved water transport and heat transfer between CL and membrane. Since no gas and liquid phase of species are present in membranes, the dissolved water may transfer from CL to the membrane. Governing equation for water content can be expressed by: PEM caPEMCLcaPEMCLPEMzPEMPEM jjzjtEW ??? ??? ,,,,,,,, ???????? (2-96) He a t f l u x Di sso l v e w a t e r f l u x C L (c at hode) PEM ? P EM T P EM C L (anode) 57 where j?,z,PEM is the dissolved water flux along z-direction calculated by the water content gradient as Eq. (2-86) and (2-87). The j?,CL,PEM is the dissolved water flux from the CL to membrane at cathode or anode side obtained by Eq. (2-88). The temperature for membranes considering the Joule heating caused by proton conductivity in the membrane is: PEMo h m i cTPEM caCLPEMcaCLPEMPEMzPEMi rqqzqtTc ,,,,,,, ?????????? ?? (2-97) The qz,PEM is the conduction heat flux along z-direction which is calculated by the temperature gradient, while qPEM,CL is the heat flux from membrane to CL. The total heat capacity of the membrane is calculated by: P E MP E MlpOHP E MP E Mp ccEW Mc ???? ??? ,2 (2-98) The effective heat conductivity of the membrane is also a function of water content as: lOHPEMPEM llOHPEMPEMPEMPEMe f f MEW kMkEWk ??? ??? // )/()/( 2 2, ? ?? (2-99) The Joule heating is given by: PEMPEMh ir ? 2, ? (2-100) where i is the local area current density (A m-2) across the membrane, and ?PEM is the ionic conductivity in membrane. 58 Based on the section 2.1 and 2.1, the calculation processes of the cell voltage and local current densities are shown in Figure 2.13. If all state variables of the fuel cell are known, the I-V curves of all segments along z-direction and the I-V curve of the single cell can be obtained. Then the cell voltage and local current densities are calculated by interpolating these I-V curves. Figure 2.13 Calculation flow chart of the cell voltage and current densities. The modeling process of the mass and heat transfer in a half cell is shown in Figure 2.14. At the very beginning of the simulation, the initial state variables are inputted as the state variables of the fuel cell transient model. Then time derivatives of state variables of the half cell are obtained from local current densities and the state variables. The new state variables are updated by the integral of the derivatives of theses state variables. The fluxes between the CL and PEM are also outputted to the PEM model. S t a t e v a r iab le o f c a t a ly s t lay e r a t c a t h o d e s ide a n d P E M L o c a l I - V c u r v e s f o r n s e g m e n t s C e ll I - V c u r v e b y int e g r a t ing t h e loc a l c u r r e n t d e n s it ies I c e l l V c e l l L o c a l c u r r e n t d e n s it y i 59 Figure 2.14 Flow chart of the mass and heat transfer for the transient model of the half cell. A ll s t a t e v a r iab les ( x ) o f h a lf c e ll q C P ,C C ( E q . 2 - 20) ( E q . 2 - 18) t T cf ? ? q G D L 1,,C P ( E q . 2 - 22) q t p,G F C ,C P ( E q . 2 - 23) ( E q . 2 - 21) t T CP ? ? ( E q . 2 - 30 ) , ( E q . 2 - 4 5 ) , ( E q . 2 - 51) t C G FCi ? ? , t s G F C ? ? t T G F Cg ? ? , For GD L 1 , GD L 2 a n d GD L 3 ( E q . 2 - 5 4 ) , ( E q . 2 - 6 0 ) , ( E q . 2 - 66) t C GDLi ? ? , t s GDL ? ? t T G D L ? ? ( E q . 2 - 7 4 ) , ( E q . 2 - 8 1 ) , ( E q . 2 - 8 5 ) , ( E q . 2 - 93) t C CLi ? ? , t s CL ? ? t T CL ? ? t CL ? ?? j i ,z,GF C ( E q . 2 - 31) j i ,G F C ,G D L 1 ( E q . 2 - 37) r i ,G F C ( E q . 2 - 44) j l ,z,GF C ( E q . 2 - 46) j l , G D L 1,G F C , ( E q . 2 - 47) r w ,e c ,G F C ( E q . 2 - 45) q z,GF C ( E q . 2 - 52) q g,G F C ,C P ( E q . 2 - 23 ) q G D L ,G F C ( E q . 2 - 53) j i ,z,GD L ( E q . 2 - 55) j i ,i n,y , j i ,out ,y ( E q . 2 - 57 ) r m ,i ,G D L ( E q . 2 - 59) j l ,z,GD L ( E q . 2 - 61) j l ,i n,y , j l ,out ,y ( E q . 2 - 64) r w ,e c ,G D L ( E q . 2 - 59) q z,GD L ( E q . 2 - 67) q i n,y , q out ,y ( E q . 2 - 69) r T ,G D L ( E q . 2 - 72) j i ,z, CL ( E q . 2 - 55) j i ,G D L , CL ( E q . 2 - 57 ) r m ,i ,C L ( E q . 2 - 5 9 , 7 5 , 7 7 ) j l ,z,C L ( E q . 2 - 61) j l ,G D L ,C L ( E q . 2 - 64) r w ,e c ,C L ( E q . 2 - 59) r w, ur ,C L ( E q . 2 - 84) j ?,z ,C L ( E q . 2 - 86) j ? ,C L ,P EM ( E q . 2 - 88) r w ,re ac t ( E q . 2 - 92) r w ,ad ( E q . 2 - 77) r w ,ur ( E q . 2 - 84) q z,C L ( E q . 2 - 67) q P EM ,C L , q C L ,G D L ( E q . 2 - 69) r T ,e c ( E q . 2 - 72) r T ,re ac t ( E q . 2 - 94) r T ,ohm i c ( E q . 2 - 95) D e r iv a t iv e s o f s t a t e v a r iab les ( ) j ? ,C L ,P EM , q P EM ,C L x? L o c a l c u r r e n t d e n s it y i I n it ial s t a t e v a r iab les ( x 0 ) o f h a lf c e ll I n t e g r a l ( ) ? x? Update state variables for next time step. 60 The calculation process of the membrane is shown in Figure 2.15, where the fluxes from anode side and cathode side to the membrane, the state variables of PEM and local current densities are the inputs of the membrane model. The derivatives of the state variables are calculated based on the mass and energy conservation principles, and then are integrated with time to obtain the new state variables. Figure 2.15 Flow chart of the transient model of the PEM. In above sections, the one dimensional transient models in all layers across membrane were proposed which were based on the mass and energy conservation principles. There are total 15 layers across membrane and 53 state variables in the single fuel cell model. The two-phase (gas H a lf c e ll m o d e l a t a n o d e s ide H a lf c e ll m o d e l a t c a t h o d e s ide L o c a l c u r r e n t d e n s it y i j ? ,C L ,P EM , q P EM ,C L a t A n o d e S ide j ? ,C L ,P EM , q P EM ,C L a t C a t h o d e S ide j ? , zP EM ( E q . 2 - 86) q z ,P EM ( E q . 2 - 67) r T ,ohm i c ,P EM ( E q . 2 - 100) S t a t e v a r iab les ( x ) o f P E M ( E q . 2 - 9 6 ) , ( E q . 2 - 97 t T PEM ? ? t PEM ? ?? D e r iv a t iv e s o f s t a t e v a r iab les ( ) o f P E M x? 61 and liquid phases) effect were involved in the modeling processes of the channel plate, gas flow channel, gas diffusion layer and catalyst layer. In the catalyst layer, the dissolved water in the ionomers was considered in the phase transition with the water vapor and liquid water. The calculation processes of the transient model involve three parts, cell voltage model, half cell model and membrane model shown in Figure 2.12, Figure 2.13 and Figure 2.14, respectively. The transient model of the single fuel cell will was implemented in MATLAB/Simulink, and then was validated by experiments based on a single cell that was specially designed. A control- oriented fuel cell stack model was obtained by simplifying the transient model of the single fuel cell that can be used for integration of FDS components and associated design of controllers. 62 Chapter 3 Experiment Setup for the PEM Fuel Cell Model Validation In order to validate the model for a cell, a single fuel cell and a test station are designed and constructed. Experimental data for static and dynamic behaviors are collected and compared with simulation results. The cell at anode side is segmented to measurement of currents, and the design of cathode GFP enables visualization of liquid water along channels. Operations for the cell are controlled by flow rates of air and hydrogen, respectively. Segmentation is a design that can be used for examining the steady and transient states of current and water distribution in the fuel cells. Mench et al. [53] used cell segmentation to study the current density distribution in a fuel cell with a serpentine channel flow pattern. Yoon et al. [54] investigated the liquid water distribution and flooding phenomenon using a segmentation design of a cell. Strickland et al. [55] used a three by three segmentation of anode CP to study the current distribution in a fuel cell with active water management. A detail review about the fuel cell segmentation is well summarized by P?rez et al. [56]. Visualization allows for study of two-phase phenomena, especially by visually observing the liquid water formation, departure and finally removal. There are several techniques proposed by some authors. Use of transparent materials [57-60], neutron radiography [61-65] and X-ray [66- 69] are the three main methods utilized to observe and quantify liquid water distribution along the GFCs and other layers. 63 Our design for a single cell is shown in Figure 3.1. The GCP at the anode side is divided into 10 segments and electrically insulated between each other to measure the current density distribution along the channels separately. On the other hand, the GFCs at the cathode side are so designed as to be transparent, which enables the one to observe flow of liquid water along flow pattern of the channels. Five straight GFCs at both sides are engraved on the surface of GCPs. Figure 3.1 An assembled single fuel cell with segmented anode and transparent cathode. I-V curve of the cell is obtained by simultaneously measuring the terminal voltage at different currents produced by an electronic load connected at the terminal and controlled by LabVIEW program. Due to the measurement of current of each segment, the I-V curve of each segment can be also measured. Images of liquid water distribution in flow channels are captured by a camera that periodically estimates the liquid water amount. The dynamic response of the cell on current steps was also measured in the transient experiments. 64 3.1 Fuel Cell Assembly The design cell has a size of 1cm?10cm that results in an active area 10cm2 catalyst coated membrane electrolyte assemble (MEA). 5 straight GFCs (1mm width, 0.6mm depth and 100mm length) were machined on the surfaces of GCPs. The cell consists of a cathode part, an anode part, two GDLs, a catalyst coated membrane (CCM) and gaskets as shown in Figure 3.2. Figure 3.2 An exploded view of the single fuel cell. Air Air H y d r o g en H y d r o g en 65 The cathode and anode parts mainly include the GFCs, current collectors and end plates. The end plate gasket is used to seal the gap between the end plates. The membrane gaskets are used to seal the catalyst coated membrane at cathode and anode sides. The MEA is a commercial catalyst coated membrane (NAFION catalyst coated membrane NR-212). The membrane size is 160 mm ? 40 mm, and the 100 mm ? 10 mm catalyst areas are coated on the surfaces of the membrane at both sides with a loading of 0.3 mg Pt cm-2 on each side. The summary of the fuel cell system components is shown in Table 3.1. Table 3.1 Geometry of the designed single cell. Name Number Scale Active area 1 cm ? 10 cm Membrane area 1.6 cm ? 16 cm Straight channel 5 1 mm ? 0.6 mm ?100 mm Segment 10 along channels 3.1.1 Cathode flow field for the liquid water observation Cathode GFCs are designed using a transparent endplate that enables to observe flow of liquid water. Five channels are engraved on the CP through the whole plate for visualization across the transparent end plate. 66 Figure 3.3 An exploded view of the cathode part. As shown in Figure 3.3, the cathode part consists of a graphite GCP, a current conductor and an end plate. Five straight channels with 1mm rib between them are engraved through the whole graphite plate except ends of the channels for connections. The graphite plate is assembled on the transparent end plate made of acrylic, and five GFCs (1mm width, 0.6mm depth) are formed by the fitting of graphite plate and the ribs on the end plate as shown in the section drawing of Figure 3.4. 67 Figure 3.4 The graphite plate engraved with GFCs. At the end of the cathode graphite plate along the channels, the channels are not machined through the whole plate and only 0.6 mm depth channels are engraved to connect all channel ribs. The copper foil (3M 1182 Tape) that is electrically conductive and adhesive is used on both sides to decrease contact resistance between the graphite plate and copper current collector. The current collector made of copper is fixed by the eight sets of aluminum screws and nuts on the end plate that serve as wire connecting terminals (positive terminals) with the load circuit connected to the cell. The inlet and outlet manifolds are machined on the surface of acrylic end plate. Two Swagelok male tube connectors are screwed on the end plate as the air inlets and outlet. Formation, departure and removal of liquid water produced by reactions can be observed through the acrylic plate directly by human eyes or image record devices. Periodic image capture of the liquid water can help understand the liquid water distribution along flow pattern transition in the channels, which is used to validate the simulation results of the two-phase. 68 3.1.2 Segmentation of anode flow field GCP on the anode side is segmented by ten graphite blocks to measure current distribution along channel, as shown in Figure 3.5. Figure 3.5 A graphite block engraved with 5 GFCs. The anode part consists of ten graphite blocks, current conductor and end plate as shown in Figure 3.6. The ten graphite blocks are attached together by the air dry vanish (Dolph Synthite AC-43) on the side surfaces that insulates the electron transport between neighboring graphite blocks. Five straight channels (1mm?0.6mm) are engraved on the surface of each graphite block, and the GFCs are formed by connecting the ten blocks. Like the cathode part, 10 copper adhesive foils are used to decrease the contact resistance between the graphite blocks and the10 current collectors. 20 sets of aluminum screws and nuts are used to fix the current collectors on the end plate, and at the same time to conduct the currents of 10 segments to the load as terminals. The end plate at anode side is also made of acrylic. The tube connectors are screwed on the end plate to connect the hydrogen supply system. Electrically insulated surface Electrically insulated surface 69 Figure 3.6 An exploded view of anode part. While operating the cell, the electrons generated in the CL at anode side and are transported to each graphite blocks through the GDL. Ten wires connected to the aluminum screws conduct the currents of the 10 segment current collectors separately, which are measured by 10 current sensors. In this way, I-V curves of ten blocks along the channels can be obtained. 70 3.2 Air and Hydrogen Supply Systems The air and hydrogen supply systems are designed to feed the reactant gases from gas cylinders to the fuel cell. The air flow and humidity at cathode side can be changed by the air supply system, while the hydrogen flow is controlled by a mass flow controller. The operating parameters of the reactant gases at the inlets and outlets of cathode and anode side are monitored by a LabVIEW program residing in a PC. Schematic diagram for the air supply system is shown in Figure 3.7. Figure 3.7 Basic schematic of air supply system. Fue l C ell Un it R H , T P H u m id ity an d T e m p e ra tu re Se n s o r Pres s u re Se n s o r Dif f e re n tia l Pres s u re Se n s o r C h e ck Valv e Bu b b le Hu m idif ie r C h e ck Valv e Chec k Valv e Air Cylin d e r dP MFC MFC Pres s u re G au ge Man u al Valv e R H , T H u m id ity an d T e m p e ra tu re Se n s o r Ma s s Flow Con tro lle r Ma s s Flow Con tro lle r 71 Two mass flow controllers are used to control the total mass flow rate and relatively humidity of supplied air. The mass flow controller line is connected to a bubble humidifier supplies wet air, while the second mass flow controller line is used to supply dry air. The air entering the cell is the mixed one from these two lines and humidity of this air is determined by a ratio of the wet air flow rate to the dry air. The operating condition of the air entering the cell is measured by a humidity/temperature transmitter and a pressure transmitter. At the outlet, temperature and humidity of exhausted air are measured by a second humidity/temperature transmitter. In addition, pressure drop between the inlet and outlet at the cathode side is also measured by a differential pressure transmitter. Likewise, hydrogen supply system is shown in Figure 3.8. The mass flow rate of supplying hydrogen is controlled by a mass flow controller. The hydrogen operating conditions at the inlet are measured by a thermocouple and a pressure transmitter. At the outlet, temperature, humidity and pressure drop are measured using the same sensors as the cathode side. Purging the anode side of the cell is carried out by a solenoid valve. Figure 3.8 Basic schematic of hydrogen supply system. dP MFC H 2 Cylin d e r Ma n u al Valv e Fil ter Pres s u re G au ge Ma s s Flow Con tro lle r Dif f e re n tia l Pres s u re Se n s o r H u m id ity an d T e m p e ra tu re Se n s o r Sole n o id Valv e Chec k Valv e T P T h e rm o cou p le Pres s u re Se n s o r Fu e l Ce ll U n it S R H , T 72 The experimental uncertainty analysis was performed for all measured data used for the validation of the transient model based on the specifications of the sensors in the fuel cell testing system. The errors of the sensors and data for the validation are shown in Table 3.2. Table 3.2 Errors of the sensors and data for validation. Name Model Error Local current sensor (A) LEM HX03-P ?0.03 Cell current sensor (A) LEM LAH 25-NP ?0.04 Cell voltage sensor (mV) AD 5B31-01 ?15 Cathode differential pressure sensor (Pa) Dwyer 616-0 ?1.6 Anode differential pressure sensor (Pa) Dwyer 616-00 ?0.8 Air humidity sensor (%) E+E Elektronik EE08-PFT-2-V11-E-602-T22 ?2 Air flow rate (sccm) AALBORG GFC17 Air 1000 ?15 H2 flow rate (sccm) AALBORG GFC17 H2 500 ?7.5 3.3 Instrumentation A schematic diagram of the segmented cell instrumentation is shown in Figure 3.9. Currents of ten anode segments by the aluminum screw terminals are separately measured using ten current sensors. Outputs of the sensors are voltages that are connected via DAQ board. 73 Figure 3.9 Diagram of cell instrumentation. Due to the low voltage output of the fuel cell that is less than 1V, the DC electronic load with an additional DC power was offset to enable the electronic load to draw high currents in the low voltage range. In our experiments, the output voltage of the DC power supply is fixed at 4V and the current of DC electronic load is controlled by the analog output of the DAQ board in the scale of 0-10V, which is corresponding to the 0-10A of the controlled fuel cell current. A digital camera is used to periodically capture images of the liquid water distribution along the GFCs, which is triggered by the LabVIEW program. All the sensors and transmitters are connected to the DAQ boards based on National Instruments (NI) hardware as shown in Figure 3.10. + - + - DC b o o st p o w e r su p p ly +- DC e lect ron ic loa d 1 0 Curr e n t se n so rs to DAQ Fue l ce ll cu rr e n t se n so r Fue l Cell Unit P ote nti al m e a s u r e m e n t 74 Figure 3.10 Schematic diagram of the DAQ system. S e g m e n t Cu r r e n t sen sor # 1 S e g m e n t Cu r r e n t sen sor # 2 S e g m e n t Cu r r e n t sen sor # 10 Fu e l c e l l c u r r e n t sen sor Fu e l c e l l p o te n tial D r y A i r fl o w Wet Ai r fl o w H 2 fl o w D r y A i r fl o w set Wet Ai r fl o w set H y d r o g e n fl o w set E - Load o u tp u t c u r r e n t set A i r i n l e t te m p e r atu r e A i r i n l e t r e l ativ e h u m i d i ty A i r i n l e t p r e ssur e A i r o u tlet te m p e r atu r e A i r o u tlet r e l ativ e h u m i d i ty H 2 i n l e t te m p e r atu r e H 2 i n l e t p r e ssur e H 2 o u tlet te m p e r atu r e H 2 o u tlet r e l ativ e h u m i d i ty Cat h o d e p r e ssur e d r o p A n o d e p r e ssur e d r o p A I (N I PCI - 6229 ) A O (N I PCI - 6229 ) A I (N I US B - 6210) 75 In addition the mass flow controllers and the DC electronic load are controlled by the analog output signals of the DAQ board. All the AI and AO connecting the NI board are isolated by the corresponding input and output blocks to protect the AI and AO modules. The fuel cell and testing station including the fuel cell, air supply system, hydrogen supply system, control, measurement and data acquisition system are shown in Figure 3.11. (a) Anode part of the fuel cell. (b) Cathode part of the fuel cell. 76 (c) The GDL and membrane. (d) Assembled fuel cell. 77 (e) The test station for the experiments. Figure 3.11 The fuel cell and experiment setup pictures. In this chapter, the design of a single fuel cell and its test station is introduced. A segmented anode gas channel plate and a transparent cathode gas channel plate are applied to measure the current distribution along gas flow channels and observe the liquid water distribution in the channels. The experiments were performed and the experimental data were compared with the simulation results for model validation in the next chapter. 78 Chapter 4 Model Integration, Validation and Comparison Simulations of the model explained in the Chapter 2 are reported in this chapter. Experimental data I-V curves for a cell were collected using the test station and compared with the results of the simulations. 4.1 Development of the Transient Model The model is developed using embedded M-code functions in MATLAB/Simulink as shown in Figure 4.1. The input variables are average current density (I_den_avg), cathode operating conditions that include flow rate (sccm_in_c_GFC), temperature (T_in_c_GFC) and RH (RH_in_c_GFC) at the inlet of cathode GFCs, and back pressure (p_back_c_GFC) at the outlet of the cathode GFCs. The anode operating conditions are the same input variables as those on the cathode side, and the value of the state variables (State_Variables) that are updated with a time interval by solving the ordinary differential equations (ODE). For the segmented model of the fuel cell, each state variable is a n element vector, where the n is the segmentation number along channel direction. Thus, the state variables can be considered as a n-by-56 state matrix. 79 The initial conditions, upper and lower limits of all the state variables are constants defined in a M-file and set in the integrator block as shown in Figure 4.1. When the model starts to run, all constants for this model are initialized before the model runs. Figure 4.1 Subsystem of the two-phase transient model. Other output variables such as cell voltage (V_cell), activation overpotential at cathode side (V_op_ca), liquid water saturation in GFCs (s_c_GFC and s_aGFC,) pressures and pressure drops along GFCs (p_c_GFC, p_a_GFC, Dp_c_GFC and Dp_a_GFC) and temperature in GFCs 80 (T_c_GFC and T_a_GFC) are set in the model for the purpose of the validation as shown in Figure 4.1. The FC_Update_Model embedded some functions as shown in Figure 4.2, where the V_cell_fun is a function of cell voltage introduced in Chapter 2, the Half_Cell_fun is a function to update the state variables in half a fuel cell (CP, GFC, GDL and CL), and the MEM_fun is a function to update the state variables in membrane. The Half_Cell_fun are used twice separately to model the transient performance at cathode and anode sides, respectively. Figure 4.2 Function structures of FC_Update_Model. V _Cel l _f un Cath od e Cata l y s t La y er S tat es (C O2 ,, s , T , ? ) A v erage Cur r en t Den s i ty ( I de n, av g ) A no de Cata l y s t La y er S tat es (C H 2 ,, s , T , ? ) Cu r r e n t De n s i ty ( I de n, v ec ) Cel l V ol tag e( V cel l ) Hal f _Cel l _f un Cath od e O pe r ati ng Con di ti on s Cur r en t Dens i ty ( I de n, v ec ) Cath od e S i de S tat es ( G F P , G F C, G DL, CL) ME M_ f un Hal f _Cel l _f un P E M S tat es F l ux f r om CL to ME M Der i v ati v e of S tat es F l ux f r om CL to ME M at c ath od e s i de Der i v ati v e of S tat es A no de O pe r ati ng Con di ti on s Cur r en t Dens i ty ( I de n, v ec ) Cath od e S i de S tat es ( G F P , G F C, G DL, CL) P E M S tat es F l ux f r om CL to ME M at an od e s i de P E M S tat es Der i v ati v e of S tat es F l ux f r om CL to ME M 81 The operating conditions for the simulation are set as constant over time. The input and output variables are written into data files with a specific sampling time. The final transient model of a fuel cell is shown in Figure 4.3. Figure 4.3 Final transient model of the fuel cell in MATLAB/Simulink. 4.2 Boundary Conditions of the Transient Model For the one dimensional model in chapter 2, the boundary conditions at z = 0 and z = length of unit are introduced for different layers (CP, GFC, GDL, CL and membrane) in a single cell. At the inlet of the GFC at cathode and anode sides, the operating conditions (flow rate, temperature and RH of inlet gas flow) are converted into the inlet boundary conditions of the gas flow in the channels, which include the gas species flux, liquid flux and heat flux calculated by: 82 G F CinG F CiniG F Cini UCj ,,,,, ? (4-1) 0,, ?GFCinlj (4-2) ?? G F CinG F CiniipG F Cin TjCq ,,,,,? (4-3) The species concentrations at the inlet of GFC can be obtained by gas state equation, and the inlet gas velocity is assumed to be proportional with the pressure difference between GFC inlet and the first segment as: G F Cin iG F CinG F Cini RT xpC , ,,, ? (4-4) ? ?)1(,,, G F CG F CinG F CinG F Cin ppkU ?? (4-5) Thus, if the pressures at the inlet of GFC are known, the species flux and heat flux at the inlet of GFC can be obtained by above equations explicitly. In our experiments, since the flow rates at the inlet are known and fixed, the inlet pressure can be calculated by substituting Eq. (4-5) into following equation: G F CinG F Cin G F CinG F CG F Cin URT pnN ,,,, ? (4-6) The inlet pressure is the positive root by solving a quadratic equation. At the outlet of the GFC, the gas species flux, liquid flux and heat flux are calculated by the similar equations as: 83 G F Co u tgG F CiG F Co u ti Ue n dCj ,,,,, )(? (4-7) G F Co u tllG F Co u tl Uj ,,,, ?? (4-8) ?? )(,,,, e n dTjCq G F CG F Co utiipG F Co ut? (4-9) The superficial gas velocity and liquid velocity is calculated by following equations by considering the liquid water effect as: ? ?b ac kG F CG F CG F Co utG F Co utg pendpendskU ?? )()( 5.1,,, (4-10) G F Co utgG F Cl G F CgG F CG F CG F CtpG F Co utl Uend end ends endsKU ,,,, 2 ,,, )( )( )(1 )( ??? ??? ? ??? ? ??? ? ? ?? ?? (4-11) The boundary conditions other layers (CP, GDL, CL and membrane) are assumed that the species flux and heat flux that are equal to 0 at inlet and outlet. 4.3 Comparison between Simulation and Experimental Results Experimental data are collected for the designed cell using the test station described above. I- V curve at the terminal is measured for two cases. The first one is obtained by sweeping a current from 0 to a maximum value of the cell for different air humidity. For the second case, the I-V curve was obtained for different gas flow rates. The data collected are compared with those simulated. 84 4.3.1 Measurement of I-V curves by current sweeping I-V characteristics of a cell are affected by amount of humidity on the air and flow rates of air and hydrogen. Thus, the RH of air and flow rates of reactant gases are changed to obtain I-V curves, where the cell current with a slope of 0.05A min-1 is applied on the cell by a electric load. The current density distribution, voltage and liquid water images in GFCs are recorded at every second. The RH of the air is set by controlling the ratio of flow rates between dry and wet air. The operating conditions for study of the effects of RH are summarized in Table 4.1. Table 4.1 Operating conditions for RH experiments. Operating conditions Cathode Anode Inlet flow rate (sccm) 133 56 Inlet temperature (?C) 23 23 RH (%) 20, 40, 60, 80 and 100 0 Back pressure (atm) 1 1 The I-V curves of these five cases for different RH of air are shown in Figure 4.4 to evaluate the air humidity effect on the performance of the fuel cell. The results show the RH significantly affects performance of the cell. The high RH decreases proton conductivity and increases the output voltage. However, improvements of the cell performance were not proportional to the amount of RH. Even though RH was increased from 80 % to 100%, the increase of the cell voltage is the smallest. 85 Figure 4.4 I-V curves of the fuel cell for different RH of air at cathode GFC inlet. Two types of graphs are shown in Figure 4.5. The first type is I-V curves of 10 segments of the cell along the GFCs at different air relative humilities and the second one is the histograms of limiting current of all segments. It was observed in I-V curves of Figure 4.5(a) that the currents of the segments near the outlet of GFCs are higher than that near the inlet because of high water vapor concentration at the outlet, which was reconfirmed by the histogram of limiting currents of the segments as shown Figure 4.5(a). When the reactants near the outlet of the GFCs are fully humidified, removal of 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 1 . 6 1 . 8 2 2 . 2 0 100 200 300 400 500 600 700 800 900 1000 C e l l C u r e n t ( A ) C e l l V o l t a g e ( m V ) R H = 2 0 % R H = 4 0 % R H = 6 0 % R H = 8 0 % R H = 1 0 0 % Cell Current (A) 86 water from inside of the cell gets difficult and more water is transported to membrane. As a result, membranes are more hydrated and the proton conductivity is improved. When the RH increases, the current difference among segment becomes less because of the higher inlet vapor concentration, which help hydrate the membrane near the inlet of GFCs in increase the cell voltage as shown in Figure 4.5(b), (c), (d) and (e). However, the limiting currents of segment #9 and #10 near the outlet do not increase, and even decrease slightly for the cases of RH = 60%, 80% and 100%, which means the variation of inlet RH does not affect the segments near the outlet greatly. In the I-V curves of segments for high RHs as shown in Figure 4.5(b), (c), (d) and (e), the I-V curves near the outlet of GFCs can cross the I-V curves near the inlet of GFCs. This cross-over effect is caused by that much more liquid water is accumulated near the outlet of GFCs and increases the mass transport voltage loss. The decrease of the limiting current is also caused by condensations of liquid water inside of the cell components as shown in histograms of Figure 4.5(c), (d) and (e), which is explained by that the liquid water in CLs covers the reaction site on catalysts and increases the resistance of reactant gas transport. (a) Current distribution along GFC for RH = 20% 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S e g m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) S e g # 1 S e g # 2 S e g # 3 S e g # 4 S e g # 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 87 (b) Current distribution along GFC for RH = 40% (c) Current distribution along GFC for RH = 60% (d) Current distribution along GFC for RH = 80% 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S e g m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) S e g # 1 S e g # 2 S e g # 3 S e g # 4 S e g # 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S e g m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) S e g # 1 S e g # 2 S e g # 3 S e g # 4 S e g # 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S e g m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) S e g # 1 S e g # 2 S e g # 3 S e g # 4 S e g # 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 88 (e) Current distribution along GFC for RH = 100% Figure 4.5 Current distribution along GFCs for different air RH. Liquid water distributions in the GFC as a function of different currents and RH = 20%, 60% and 100% are captured at several time instants and shown in Figure 4.6. Encircled red line shows the areas where droplets of liquid water are formed and attached in GFC. For all cases of RH experiment, the liquid water is firstly observed in the outlet area of the GFCs. The case of RH = 100% is earliest one that the liquid water droplets are observed in GFCs, while the case of RH = 20% is the latest one. The wet area of the case of RH = 100% is the largest one which almost cover 50% area in GFCs. As shown in the I-V curves and histograms of Figure 4.5, the voltages corresponding to the wet segments drop rapidly with the appearance and increase of the wet area. The histograms of limiting currents also shows that the case of RH=20% has the largest limiting current values at the segments near the GFC outlet for smallest wet area as shown in Figure 4.6. 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S e g m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) S e g # 1 S e g # 2 S e g # 3 S e g # 4 S e g # 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 89 Figure 4.6 Liquid water distribution along GFCs for different inlet RH of air. I = 0 .8 A I = 1 .2 A I = I l i m i t L i q u i d w a te r d i s tri b u ti o n i n G FC fo r RH = 2 0 % I = 0 .8 A I = 1 .2 A I = 0 .8 A I = 1 .2 A I = 1 .6 A I = 1 .6 A I = 2 .0 A I = I l i m i t I = I l i m i t L i q u i d w a te r d i s tri b u ti o n i n G FC fo r RH = 6 0 % L i q u i d w a te r d i s tri b u ti o n i n G FC fo r RH = 1 0 0 % I = 0 .8 A I = 1 .2 A I = I l i m i t L i q u i d w a te r d i s tri b u ti o n i n G FC fo r RH = 2 0 % I = 0 .8 A I = 1 .2 A I = 0 .8 A I = 1 .2 A I = 1 .6 A I = 1 .6 A I = 2 .0 A I = I l i m i t I = I l i m i t L i q u i d w a te r d i s tri b u ti o n i n G FC fo r RH = 6 0 % L i q u i d w a te r d i s tri b u ti o n i n G FC fo r RH = 1 0 0 % 90 Effects of flow rates of air and hydrogen on performance of a cell are experimentally studied. For this study, the RH for air is set to be 100% and others are according to Table 4.1. The I-V curves for different cases of flow rates are shown in Figure 4.7. The results shown that change of the flow rates does not significantly affect voltages and limiting currents, even thought the flow rate is at the highest (Air/H2 = 200/84 sccm). Figure 4.7 I-V curves for the cases of different flow rates of air and hydrogen. Effect of the flow rates on currents in the 10 segments are shown in Figure 4.8. For low flow rates of gases as shown in Figure 4.8 (a), the voltages of segments from #3-10 drop rapidly. Conversely, the voltage drops take place in the segments of #7-10 when flow rates are high, as 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 1 . 6 1 . 8 2 2 . 2 0 100 200 300 400 500 600 700 800 900 1000 C e l l C u r r e n t ( A ) C e l l V o l t a g e ( m V ) A i r / H 2 = 6 6 / 2 8 s c c m A i r / H 2 = 1 0 0 / 4 2 s c c m A i r / H 2 = 1 3 3 / 5 6 s c c m A i r / H 2 = 1 6 6 / 7 0 s c c m A i r / H 2 = 2 0 0 / 8 4 s c c m 91 shown in Figure 4.8 (e). The liquid water effect (voltage rapid drop) on the performance of fuel cell is more apparently observed in the cases of low flow rate of reactant gases. When the flow rate increases, liquid water residing in the channels is easily removed and as a result resistance against reactant flow is reduced. The limit of current of each segment is also shown in the histograms of Figure 4.8, where the limit current is almost evenly distributed, especially for the segment #5-10. Thus, the flow rates of reactant gases have little effect on the current distribution along the channels. (a) Air/H2 = 66/28 sccm (b) Air/H2 = 100/42 sccm 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S g e m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) S e g # 1 S e g # 2 S e g # 3 S e g # 4 S e g # 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S g e m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) 1 S e g # 2 e g # 3 e g # 4 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 92 (c) Air/H2 = 133/56 sccm (d) Air/H2 = 166/70 sccm (e) Air/H2 = 200/84sccm Figure 4.8 I-V curves of 10 segments for the cases of different flow rates. 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S g e m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) S e g # 1 S e g # 2 S e g # 3 S e g # 4 S e g # 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S g e m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) 1 2 3 4 S e g # 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 100 200 300 400 500 600 700 800 900 1000 S g e m e n t C u r r e n t ( A ) C e ll V o lt a g e ( m V ) S e g # 1 S e g # 2 S e g # 3 S e g # 4 S e g # 5 S e g # 6 S e g # 7 S e g # 8 S e g # 9 S e g # 1 0 1 2 3 4 5 6 7 8 9 10 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 S e g m e n t N o . L im it C u rr e n t (A ) 93 These images of liquid water distribution in GFCs for three cases of different flow rates are also captured shown in Figure 4.9, where the area circled by the red line is the wet area by liquid water. For the case of low flow rate (Air/H2 = 66/28), the wet area increases rapidly, and the ratio of the area reaches about 70% of that at limit cell current. For other two cases shown in Figure 4.9, even though more liquid water was observed at lower flow rate (Air/H2 = 133/56 sccm), the difference of wet areas was smaller, and the wet areas are almost the same as that at the limiting currents. As a result, an increase of air flow rate decreases the wet area in GFCs and improves performance of the fuel cell. This relationship needs to be considered for optimization of flow rates. 94 Figure 4.9 Liquid water distribution in GFCs for different gas flow rates. I = 0 .8 A I = 1 .2 A I = I l i m i t L i q u i d w a te r d i s tri b u ti o n i n G FC fo r Ai r/H 2 =6 6 /2 8 sccm I = 0 .8 A I = 1 .2 A I = 1 .6 A I = I l i m i t L i q u i d w a te r d i s tri b u ti o n i n G FC fo r Ai r/H 2 =2 0 0 /8 4 sccm I = 1 .6 A I = 2 .0 A I = 2 .0 A L i q u i d w a te r d i s tri b u ti o n i n G FC fo r Ai r/H 2 = 133 / 56 sccm 95 On the other hand, liquid water accumulated in GFCs increases pressure drop along the channels. Effects of different RH and flow rates on pressure drops along GFCs are shown in Figure 4.10 and Figure 4.11, as the current increases. It was observed in Figure 4.10 that the pressure drop was not significantly affected by RH, but propositional when the RH is relatively high like 60%, 80% and 100%. On the other hand, the pressure drop tends to follow the current density. When the current density becomes high, the pressure drop increases more rapid than that at the low current density. Figure 4.10 Effects of RH on pressure drop along GFCs at cathode side. 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 1 . 6 1 . 8 2 2 . 2 2 . 4 2 . 6 50 5 2 . 5 55 5 7 . 5 60 6 2 . 5 65 6 7 . 5 70 7 2 . 5 75 7 7 . 5 80 C e l l C u r r e n t ( A ) P r e s s u r e D r o p a lo n g C a t h o d e G F C ( P a ) R H = 2 0 % R H = 4 0 % R H = 6 0 % R H = 8 0 % R H = 1 0 0 % 96 When flow rates were increased, the pressure drop tends to follow them linearly but with a high slope at high current density as shown in Figure 4.11. The gas flow rate has more influence on pressure drop along GFCs than the RH does. Figure 4.11 Effects of air flow rates on pressure drop along GFCs. From above analysis, the increase of pressure drop is mainly caused by the accumulation of liquid water in GFCs. For high cell current, more liquid water is observed in the channels, which reduces section area of GFCs and increases the transport resistance of the gas flow in the channels. 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 1 . 6 1 . 8 2 2 . 2 20 40 60 80 100 120 140 C e l l C u r r e n t ( A ) P r e s s u r e D r o p a l o n g C a t h o d e G F C ( P a ) A i r / H 2 = 6 6 / 2 8 s c c m A i r / H 2 = 1 0 0 / 4 2 s c c m A i r / H 2 = 1 3 3 / 5 6 s c c m A i r / H 2 = 1 6 6 / 7 0 s c c m A i r / H 2 = 2 0 0 / 8 4 s c c m 97 4.3.2 Validation of the model and analysis Comparisons between simulations and experiments are carried out with respect to static and dynamic behavior of a cell. The parameters used for the model is listed in Appendix A, which reflect parameters of the designed cell that includes geometric size of the cell, and material properties of CP, GDL and membrane. Operation conditions are the same as those used in the experiments list in Table 4.1. The model developed is non-isothermal, segmented, and transient. Because of limitation of computational time, the comparison is performed only for a simplified model that considers isothermal, non-segmented and transient case. In order to reduce the computational time, simulation was performed using a high performance computer cluster (HPCC) that is available in Auburn University. The HPCC allows for parallel computation, where the original codes are converted to S-functions. The model is used to calculate I-V curves at different RH values, which is compared with those obtained from experiments at RH = 20%, 60% and 100%, as shown in Figure 4.12. For the case of RH = 20%, the difference of I-V curves at low current density is higher than that at the high current density. For the cases of RH = 40% and 100%, the I-V curves show similar trends, but large errors existed in low and high current range. The voltage difference in the low current density is dominantly caused by the cathode over-potential while in the high current the liquid water accumulated in the cell mainly affects the voltage drop. 98 Figure 4.12 I-V curves of model predication and experiments for different air RH. The liquid water saturation in GFCs at cathode side at different RHs was calculated and plotted in Figure 4.13. It can be recognized that a high RH firstly shows the liquid water in the channels. The amount of water is the largest, which is comparable to that in Figure 4.6. At a RH = 20%, water is not likely condensed and transported to the GFCs. Since the liquid water amount in GFCs cannot be accurately measured by the current test station, qualitative comparison was impossible. The comparison of Figure 4.6 and Figure 4.13 shows that the liquid water amount in GFC estimated from the circled of images has a similar increase trend to that obtained from the simulation results. 0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 0 100 200 300 400 500 600 700 800 900 1000 C u r r e n t D e n s i t y ( A / m 2 ) C e l l V o l t a g e ( m V ) R H = 2 0 % , M o d e l R H = 2 0 % , E x p R H = 6 0 % , M o d e l R H = 6 0 % , E x p R H = 1 0 0 % , M o d e l R H = 1 0 0 % , E x p 99 Figure 4.13 Liquid water saturation in GFCs at cathode side of model predication. To validate transient behavior of the cell, a multi-step current is applied and the responses are plotted in Figure 4.14. All the operating conditions are listed in Table 4.2. Table 4.2 Operating conditions for dynamic experiments. Operating conditions Cathode Anode Inlet flow rate (sccm) 133 56 Inlet temperature (?C) 23 23 RH (%) 100 0 Back pressure (atm) 1 1 0 500 1000 1500 2000 0 0 . 0 5 0 . 1 0 . 1 5 0 . 2 0 . 2 5 0 . 3 0 . 3 5 0 . 4 0 . 4 5 0 . 5 C u r r e n t D e n s i t y ( A / c m 2 ) L iq u id W a t e r S a t u r a t io n in G F C a t C a t h o d e S id e ( s c a , G F C ) R H = 2 0 % R H = 6 0 % R H = 8 0 % 100 The current applied has two steps. The first step is from 0A to 1.5A at 13sec and 1A at 26sec, and the second step is from 1A to 2A at 37sec and 2A to 0A at 52sec. Figure 4.14 The cell current of the fuel cell used in the experiments and simulation. The simulation and experiment response of the terminal voltage were plotted in Figure 4.15. Firstly, the overshoots for both cases show the similar dynamics. The response time for the step currents takes about 0.6sec that is relatively short. This overshoot can be explained by that the concentration of reactant gases response at reaction sites for the sudden change of current. When the cell current was suddenly increased, the consumption rates of reactant gases were also suddenly increased that caused the low concentration of reactant gases at the reaction sites. The low concentrations of reactant gases can decrease the OCV and increase the activation overpotentials. So the cell voltage is also suddenly reduced. When the reactant gases in CLs were 0 5 10 15 20 25 30 35 40 45 50 55 60 65 -1 - 0 . 5 0 0 . 5 1 1 . 5 2 2 . 5 3 T i m e ( s e c ) C e ll C u r r e n t ( A ) 101 replenished by the supplied gas in GFCs in a short time, the concentration of reactant gases at reaction site was increased and then reached steady state. The cell voltage also was increased and then keeps stable in a short time. For the sudden decrease of cell current, the inverse process occurred and caused the overshoot of decrease of the cell voltage. Figure 4.15 The cell voltage responses of experiment and simulation result. After the short overshoot, the terminal voltage tends to increase relatively slowly in experiment. The increase of the voltage is thought to be caused by decrease of the electrical or ionic resistance. The increase of liquid water amount as byproduct in cells with time can decrease the contact resistance between layers with relatively slow rate, and decrease the ohmic loss. Since the 10 20 30 40 50 60 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 T i m e ( s e c ) C e ll V o lt a g e ( m V ) E x p M o d e l 102 contact resistance is assumed to be constant in the current model of the single fuel cell, the ohmic loss due to electrical resistance does not change when the cell current is in steady state. The increase of water amount also can increase the conductivity of membrane in cells. When phase transition rate between dissolved water in the ionomers and water vapor in void space in CL is relatively small, then the amount of water in the membrane increases with time, which can increases the ionic conductivity in the membrane. As a result, the voltage loss due to ionic resistance in membrane decreases. Because of the phase transition parameters used in current model, phase changes in CL take place faster than that in the experiments and as a result the water as byproduct transports from the ionomers to void space in CL rapidly and cause the water content in the membrane and the ohmic loss of cell voltage stable in as short time. From above analysis in this chapter, the simulation result of the model shows that the I-V curves of model match the trend of the results based on the RH experimental data. The dynamic response of model simulation can predict the short time overshoot in the experimental data. 103 Chapter 5 Modeling of Fuel Delivery System In this chapter, the FDS including the two supply lines and two recirculation lines is modeled. All components in the FDS were modeled, and then integrated with the anode model of the fuel cell stack. For the purpose of the modeling, the FDS shown in Figure 1.7 is converted into the figure as shown in the Figure 5.1, where three manifolds are added before the ejector (ejector manifold), the inlet of fuel cell stack (supply manifold) and the outlet of the fuel cell stack (outlet manifold). These manifolds are modeled as control volumes to express the dynamic characteristics in these volumes at corresponding position. As shown in Figure 5.1, the components of the FDS include three manifolds (ejector, supply, and return manifolds), an ejector, a blower, a pressure regulator, a flow control valve, and a purge valve. The pressure and water dynamics were considered in the modeling of supply and return manifolds. The rotator inertia of the blower motor was also considered in the dynamic model of the blower. A static ejector model is presented to determine the flow rates at its two inlets and one outlet. Other components such as the flow control valve, low-pressure regulator, and ejector were built as static models. Based on the relationship of the flow rates of each component, the FDS model can be modeled by connecting all components. 104 Figure 5.1 Schematic diagram of FDS. To simply the modeling process for the purpose of control, several assumptions were made for the purpose of modeling as: 1) the outlet pressure of high-pressure regulator is stable; 2) there is no pressure drop along the pipe connections; 3) spatial variations are neglected in manifolds; 4) no contaminant gases are in the hydrogen supplied by the tank; 5) no gas crosses the membrane and no gas leaks in the fuel cell; 6) the ideal gas law applies to all control volumes; 7) all manifolds work in isothermal conditions; 8) liquid water is not re-circulated by the blower or ejector; and, 9) gaseous and liquid water are in equilibrium in all control volumes. The assumption 1 is based on that the change of pressure at the outlet of high pressure regulator is very small. The assumption 2, 3 and 6 are standard assumptions for lumped parameter models. The assumption 4 is made because of the small amount of contaminant gases S Hy dr o g en T a nk Hi g h P r es s ur e Regul a to r Ejecto r M a nifo l d F l o w C o ntr o l V a l v e Lo w P r es s ur e R eg ul a to r Ejecto r Bl o w er Supply M a nifo l d R eturn M a nifo l d P urg e V a l v e F uel C el l Sta ck 105 in supplied hydrogen. The gas crossover and leakage effects are neglected since they do not significantly affect the fuel cell performance for commercial membrane and gasketings. The assumption 7 is from the fact that manifolds in the FDS do not change rapidly in a large range when cooling system is applied on fuel cell stacks. The assumption 8 is based on that the liquid water can only be purged out by high speed gas flow when the purge valve is open. The assumption 9 is used to neglect the transient process of phase change between the water vapor and liquid water. Because of the back diffusion effect across the membrane, water can be transported into the anode of the fuel cell and condense in the CL, GDL, and GFCs. In addition, the liquid water can also flow into the return manifold by the entrainment effect of the gas flow, be discharged by purging operation, and condensed by the saturated water vapor in the return and supply manifolds. If large current is drawn out by the fuel cell stack, large amount of water is generated and exists in the fuel cell stack and manifolds, which may block the gas transport path in the stack and FDS. Thus, the gas-liquid two-phase phenomenon should be considered in the models of the supply and return manifolds and the purge valve. 5.1 Manifolds In this section, the three manifolds (ejector manifold, supply manifold and return manifold) are modeled separately. The two-phase (gas and liquid phases) condition is considered in the supply and return manifolds, while single-phase (gas phase) condition is considered in the ejector manifold. 106 As shown in Figure 5.1, the ejector manifold represents the pipe volume which connects the flow control valve and primary inlet of the ejector. Because pure hydrogen is supplied from the flow control valve, only the hydrogen pressure dynamic is considered in the model, ? ? pejf c vem emHem WWV TRdtdp ,2 ?? (5-1) where pem is the hydrogen pressure in ejector manifold, Wfcv is the mass flow rate of flow control valve and Wej,p is the mass flow rate of primary inlet of ejector. Since the hydrogen exhaust from a fuel cell stack can contain water vapor and can be re- circulated into the supply manifold by the recirculation pumps, water vapor is condensed if its partial pressure is higher than the saturation pressure in the manifold. Thus, the hydrogen pressure and water activity dynamic equations can be used to describe the dynamic change of the amount of hydrogen and water in the manifold as [39]: ? ? o u tsmHblHo u tejHl p rHsm smHsmH WWWWV TRdtdp ,,,,,,, 222222 ???? (5-2) ? ? o u tsmwblvo u tejvsmsms a t smOHsmw WWWVTp TRdtda ,,,,,, )( 2 ??? (5-3) where pH2,sm is the hydrogen partial pressure in supply manifold; WH2,lpr, WH2,lpr, and WH2,lpr are the mass flow rates of the low-pressure regulator, ejector outlet and blower. The water activity aw describes the total amount of water in gas and liquid phases in a controlled volume [39], and psat is the saturation pressure, which can be written as a function of temperature [39]. Wv,ej,out and Wv,bl t the vapor mass flow rate of the ejector outlet and blower as the total flow rate entering the 107 manifold, and Ww,sm,out is the water outflow rate of the supplied manifold involving liquid water and water vapor. Since the inflow of the return manifold is from the fuel cell stack, water vapor and liquid water can also exist in this manifold and the dynamic model is: ? ?p u r g eHblHsejHinrmH rm rmHrmH WWWWV TRdtdp ,,,,,,, 2222 22 ???? (5-4) ? ?p u r g ewblvsejvinrmw rmrms a t rmOHrmw WWWWVTp TRdtda ,,,,,,, )( 2 ???? (5-5) In the above equations, the water inflow rate of the return manifold Ww,rm,in and the water flow rate of purge valve Ww,purge are the sum of water flow rates in gas and liquid phases. The gas constant and specific heat capacity of the gas in the supply and return manifold can be obtained by averaging the values of species based on the mass fraction as: ?? ?? i iigi iig yccyRR , (5-6) where Rg represents the gas constant of mixed gas, cg refers the specific heat capacity (constant volume or constant pressure) of the mixed gas in manifolds, and yi is the mass fraction of species i obtained from the partial pressures. Then, the density and specific heat ratio of the mixed gas can be obtained by: gv gpg gg c ca n dTR p , ,, ?? ?? (5-7) 108 where p is the total pressure of mixed gas, and the ?g is the specific heat ratio of the mixed gas in manifolds. 5.2 Blower In Figure 5.1, a blower is used as another recirculation pump to actively change the mass flow rate of recirculation. Comparing the ejector, the blower speed is controlled by a voltage signal and consumes power from the system. The model of blower consists of two parts, a static blower model and a dynamic electric motor model. For the static blower model, the blower flow can be expressed as the function of pressures at inlet and outlet of the blower, inlet gas temperature and speed of the blower motor. Different curve fitting methods can be applied on the flow rate calculation based on different forms of the curve fitting function [70]. Jensen & Kristensen [70] presented a model by using the dimensionless head parameter to reduce to number of variables for the curve fitting, which is also used in the compressor modeling for the air supply at cathode side by Pukrushpan [36, 71]. For the blower, or compressor model, the dimensionless head parameter is defined as [70], 2 1 , 2 1 1, , bl rm sm strmp bl U p pTc rmg rmg ? ? ? ? ? ? ? ? ? ? ??? ? ? ??? ? ? ? ? ? ? (5-8) where Vbl is the tip velocity of the blower blade, which is given by: 2 bcbl bl DU ?? (5-9) 109 where Dbl is the diameter of the blower blade, and ?bc is the corrected angle velocity (rad s-1) of the blower rotator defined as: refrm blbc TT /?? ? (5-10) where Tref is the reference temperature at 288 K. Another important dimensionless parameter defined in the Jensen & Kristensen method [70] is the scaled blower flow rate ?bl given by: blblrm bcbl UD W 24??? ? (5-11) where Wbc is the corrected mass flow rates defined as: refrm refrmblbc pp TTWW //? (5-12) where Wbl is the mass flow rate of the blower, and pref is the reference pressure (101325 Pa). The Jensen & Kristensen method [70] used the following forms of function to express the relationship between the scaled blower flow rate, blower efficiency, and the dimensionless head parameter [70]: 3,2,1,, 213 21 ?????? iMaaaw h e r ea aa bliiibl blbl ??? (5-13) 3,2,1,, 3 213221 ??????? iMb Mbbbw h e r ebbb bli bliiiblblbl ??? (5-14) 110 where a and b are function parameters, and Mbl is the inlet Mach number defined by: rmrmgrmg blbl TRUM ,,? ? (5-15) When the blower dimensionless head parameter ?bl is calculated by Eq. (5-13) from the gas states at the inlet and outlet of the blower, the output corrected blower flow is given by: blblrmblbc UDW 24???? (5-16) The parameters in Eqs (5-13) and (5-14) can be obtained by the curve fitting method from the experimental data. For the current blower model, these parameters are given in Table 5.1. Table 5.1 Blower map function parameters. A Value b Value a11 -1.598?10-3 b11 -7923.8 a12 2.663?10-2 b12 1.502?104 a21 -3.062?10-2 b13 0.2144 a22 -0.1740 b21 24.91 a31 14.55 b22 -821.5 a32 -15.73 b23 -4.093?10-2 b31 -4.929?10-2 b32 0.8529 b33 1.715?10-2 111 A blower map as shown in Figure 5.2 presents the relationship of the flow rate and pressure ratio. To plot the blower map of above model, the inlet pressure, temperature and RH are set to 1.4Bar, 70?C and 100%. For each blower speed range from 6000 RPM to 18000 RPM with 2000 RPM step, the plot of mass flow rate and pressure ratio curve can be drawn based on the calculation of above model as shown in Figure 5.2, where the pressure ratio is defined as that of outlet pressure to the inlet pressure. Figure 5.2 Blower map of the mass flow rate and pressure ratio. From the blower map, the mass flow rate of the blower can be obtained by the pressure ratio and the blower angular speed if the inlet gas state is known. The angular speed is updated by the 0 . 5 1 1 . 5 2 2 . 5 3 3 . 5 4 x 1 0 -3 1 1 . 0 2 1 . 0 4 1 . 0 6 1 . 0 8 1 . 1 1 . 1 2 1 . 1 4 1 . 1 6 1 . 1 8 M a s s f l o w r a t e ( K g s - 1 ) P r e s s u r e r a t i o 6000 8000 10000 12000 14000 16000 18000 112 dynamic part of the blower based on the motor model, which is the only state variable of the blower updated by the following dynamic equations, ? ? blbmblbl Jdtd ??? ?? 1 (5-17) bl rm sm blbl rmrmp bl Wp pTc rmgrmg ?? ? ? ? ?? ? ? ? ?????????? ? 1, , 1 , ? ? ??? (5-18) ? ? blvblbmtbmbm kuRk ??? ?? (5-19) where Jbl is the rotational inertia of the rotator, ?bm, kt, kv and Rbm are motor parameters, and ubl is the control voltage of the blower motor. The species (hydrogen and water vapor) flow rates of the blower can be calculated by the total and the species mass fraction in the return manifold. 5.3 Ejector Study of ejectors has been taken for many years, especial for the applications of chimerical industry, oil plant, airplane jet propulsion and refrigeration [72-75]. Recent research [38, 76-79] introduced the ejector as the anodic hydrogen recirculation pumps, which allows utilize the high pressure gas from the hydrogen tank to entrain the relatively low pressure hydrogen exhausted into the supply line as a recirculation pump. The ejector has the characteristics of simple structure, no movable parts, and no power consumption which increases its reliability and ease of maintenance on the application of FDS of fuel cell stack. 113 The basic structure of an ejector is shown in Figure 5.3. In the FDS, the high-pressure pure hydrogen of the primary flow is chocked at the throat of the primary inlet and forms supersonic flow in ejector. When the primary flow expanded passing through the throat, its velocity increased and a low pressure zone is formed. Then, the relatively low-pressure hydrogen at the secondary inlet is entrained by high-speed primary flow and mixed in the mixing chamber. The shock wave occurs in the ejector depending on the outlet pressure. Figure 5.3 Basic structure of an ejector. Reviews of the mathematical modeling on ejectors have been published by Bartosiewicz [80] et al. and He et al. [81]. The models in the publications can be divided into two types, thermodynamic model and CFD model. The thermodynamic model is based on the thermodynamic analysis and solved in one dimension along the flow direction. The constant- pressure mixing model or constant-area mixing model are applied in the thermodynamic model as the mixing assumption [77, 82, 83]. CFD method is applied on the ejector modeling by Prim ary Flow Secon dar y Flow Mixed Flow Suction Cham ber Mixin g Cham ber Diffu ser 1 2 114 utilizing the powerful computation of computers on the solution of two dimensional or three dimensional problem [84-86]. In following section, a mathematical model of ejector is presented based on the shock circle by Zhu et al.[77], and the CFD model is used to validate the shock circle model and its parameters. 5.3.1 Mathematical modeling The primary flow rate at the primary inlet throat can be calculated by the convergent nozzle equation based on chocked flow and un-chocked flow conditions, which depends on the ratio of inlet pressure to back pressure [77]: ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ? ? ? ?? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ??? ? ? ? ? ? ? ? ? ? ? ?? ? 1 12 1 11 1 1 12 1 2 2 2 2 2 2 2 2 2 2 2 2 22 2 1 2 1 1 2 1 2 1 2 ?? ? ? ?? ? ? ?em rm ? ? em rm ? ? em rm em? ?e j , p t , e jem ? ?em rm ? ? ?em? ?e j , p t , e jem e j , p ?p p if p p ?p p ?R ?? ?p ?p p if ??R ?? ?p W (5-20) where At,ej is the throat area, and ?ej,p is the isentropic coefficient of primary flow. When the pressure ratio is less than a critical value (critical pressure ratio), the primary flow is choked, and the Mach number at the nozzle throat (section 1 in Figure 5.3) is 1. When the pressure ratio is greater than the critical value, the throat Mach number can be calculated by: ?? ? ? ? ?? ? ? ? ??????????? ? 112 2 2 2 1 H H em rm H t p pM ? ? ? (5-21) 115 It is assumed that the primary flow and secondary flow mix at section 2 in the mixing chamber as shown in Figure 5.1, where the static pressure of primary flow is equal to that of the secondary inlet pressure prm. Applying the isentropic equations from section 1 to section 2 for the primary flow, we have: ?? ? ? ? ?? ? ? ? ??????????? ? 112 2 2 2 1 2, H H rm em H p p pM ? ? ? (5-22)) 2 2,2, 2 11 2 pHpem MTT ??? ? (5-23) 2,2,2, 22 pHHpp TRMU ?? (5-24) )1(4 1 2 2 2, 2,e x p ,2, 2 2 2 2 )1(2 )1(2 ? ? ??? ? ??? ? ?? ??? HH tH pH p tejtp MMMMDD ? ? ? ? ? (5-25) where Mp,2, Tp,2, Up,2 and Dp,2 are the Mach number, temperature, gas velocity, and diameter of primary flow at section 2 In Zhu?s model [77], the velocity distribution of secondary flow at section 2 is assumed to be an exponential function along the direction of the radius. The mass flow rate of secondary flow can be calculated by integrating the velocity distribution function and expressed as: ? ?? ?? ?? ? 212 2,2,22,22,, ?? ???? vv pvpppsejs nn RnRRRRUW ?? (5-26) 116 where R2 and Rp,2 are the radius of section 2 and primary flow at section 2, nv is the exponent of the velocity function, and s? is the average density of the secondary flow, which is equal to the density of mixed gas in the return manifold. Zhu [77] also found that the exponent of velocity function nv is a variable related to the inlet pressures and diameter ratio of throat to mixing chamber, which is given by: 1 6 6 8.04 5 6.0)05.0/e x p (103 9 4.1 4 ???? ? Dpvn ?? (5-27) where 1.1,8.0, / ejpejsp PP?? (the pressure is in bar), and ?D=Dm/Dt. Thus, the average velocity of secondary flow at the section 2 Us,2 can be calculated by: ))(2)(1( )(2 2,2 2,2,22,2, pvv pvpps RRnn RnRRUU ??? ??? (5-28) The above calculation is based on the assumption of critical working mode of ejector, which means that the entrain ration defined as the mass flow rate ratio of secondary flow to primary flow, are not affected by the backpressure as shown in Figure 5.4. When the backpressure is higher than the critical pressure pc,b, the ejector will work in a subcritical mode, in which the entrain ratio will drop rapidly. If the backpressure increases further to another higher critical value p0,b, the ejector will work in back flow mode, in which the entrain ratio becomes 0 and no secondary flow is entrained. The back flow mode of the ejector should be avoided in real applications. However, the ejector in the FDS can work in a subcritical mode, and the calculation of secondary flow in this mode should be considered in the ejector model. 117 Figure 5.4 Working mode of ejector. The critical backpressure pc,b can be obtained using the equations of mixing and expansion processes in mixing and diffuser chambers, respectively. The mass, energy and momentum conservation equations are used to solve the states of mixtures of gases in the mixing chamber: ejsejpmmg mmm WWTR AVp ,,, ?? (5-29) ? ? rmrmpejsemHpejpmmmpejsejp TcWTcWUTcWW ,,,,2,,, ,21 2 ???????? ?? (5-30) ? ? ? ? ? ? 2,2,2,,2,2,2,,,, sssejspppejpejsejpm i xmmm ApVWApVWWWVAp ?????? ? (5-31) p back En t rai n rat i o p c , b p 0, b C r itical m o d e Su b - cr itical m o d e B ac k f lo w m o d e 118 where pm, Um and Tm are the pressure, velocity, and temperature of the mixed flow. The gas constant and specific heat capacity of the mixed flow are calculated by the mass fraction of primary and secondary flows, which are given by: ejsejp ejsrmg ejsejp ejpHmg WW WRWW WRR ,, ,, ,, ,, 2 ???? (5-32) ejsejp ejsrmp ejsejp ejpHpmp WW WcWW Wcc ,, ,, ,, ,,, 2 ???? (5-33) In Eq. (5-30), the gas states in the ejector manifold and return manifold are considered the stagnation states of primary and secondary flow. In Eq. (5-31), the ?mix is a coefficient that accounts for the frictional loss in the mixing process. Ap,2 and As,2 are the section areas of primary and secondary flows at section 2. The sum of Ap,2 and As,2 is equal to the mixing chamber section area Am,ej. The pressure of primary flow at section 2 is equal to prm, and the average pressure of secondary flow at section 2 is determined by: 12 2,,2, , , 2 11 ???? ? ??? ? ? ??? rmg rmg srmgrms Mp p ? ?? (5-34) where the Mach number of secondary flow at section 2 is rmrmgrmgss TRUM ,,2,2, / ?? . By solving the conservation equations, we can obtain the states of mixed flow and use them for the isentropic process in the diffuser chamber and the stagnation pressure of the mixed flow (the critical backpressure) is: 119 12,, , , 2 11 ??? ? ??? ? ? ??? mg mg mmgmbc Mp p ? ?? (5-35) where ?g,m is the specific heat ratio of the mixed flow, which is calculated by Eq. (5-7) where the specific heat capacity of the mixed flow is from Eq. (5-33). If the backpressure is equal to the critical pressure of back flow p0,b, the mass flow rate of the secondary flow is 0. In this situation, the pressure of the secondary flow at section 2 is equal to its stagnation pressure prm. Substituting Ws,ej=0 and ps,2=prm into Eqs. (5-29), (5-30) and (5-31), the stagnation pressure of the mixed flow p0,b, can be calculated by the conservation equations. If the mass flow rate of the secondary flow in critical mode by Eq. (5-26) is expressed as Wc,s,ej, the mass flow rate of the secondary flow can be derived from the working mode and the assumption that the entrain ratio drops linearly with increasing back pressure in the subcritical mode as shown in Figure 5.4: ? ? ? ? ? ? ? ? ?? ? ? ? ? bcsmejsc bsmbc bcb smb ejsc bsm ejs ppW ppp pp ppW pp W ,,, ,0, ,,0 ,0 ,, ,0 , 0 (5-36) where psm is the supply manifold pressure that is the backpressure of ejector. In the above ejector model, the stagnation pressure of the primary flow should be the highest among the inlet pressures and backpressure or the ejector will work in back flow mode. The ejector is working in the critical or sub-critical mode for real applications, and the mass flow rates of primary and secondary flows are a function of gas states of the ejector manifold, return, and supply manifolds. The species flow rates at the outlet of the ejector are derived by: 120 rmHsejpejo u tejH yWWW ,,,,, 22 ?? (5-37) )1( ,,,, 2 rmHsejo u tejv yWW ?? (5-38) 5.3.2 CFD validation To validate and optimize some parameters in the ejector model shown in equations of section 5.3.1, the CFD simulation or experiment can be performed for different cases in which the gas states at first and secondary inlets change separately. In this section, a CFD model is built in the commercial software FLUENT 12.0.16 (ANSYS, Inc). FLUENT is a commercial CFD software package which utilizes the finite-volume method to convert the conservation equation (PDEs) to algebraic equations that can be solved numerically. The mesh of the ejector model was made of 77,500 quadrilateral cells. A pressure- based, absolute, steady and axi-symmetric solver is used to solve the non-linear equations. Since the supersonic flow was thought to be compressible and turbulent, the RNG k-? model was selected to govern the turbulent flow. The new wall treatment was set as the ?standard wall function?, which gave reasonably accurate results for the wall bounded with very high Reynolds number flow [86].The species model with the option of species transport was selected for the calculation of species (hydrogen and water vapor) transport in the ejector. A grid sensitivity analysis for the ejector mesh was performed based on the cases with grids of 41,500, 77,500 and 127,000 cells. The grid with 77,500 cells was found to provide relatively better grid independent results. 121 The mixture of hydrogen-vapor is used as the working fluid of the model treated as an ideal gas for the relatively low gas pressure at the inlets of the ejector. The mixing law and ideal-gas- mixing law are applied on the computation of the specific heat and thermal conductivity of the gas mixtures, respectively. The boundary conditions of primary and secondary inlets were set as pressure-inlet, while the outlet of ejector was set as pressure-outlet. The RH variation of secondary inlet is set by the mole fraction of water vapor. Different cases are simulated based on the variation of the primary inlet pressure, secondary inlet pressure and RH of secondary inlet, and the mass flow rates were compared with the results of the thermodynamic model discussed in previous section. The calculation were considered as be convergent when the all types of the calculation residual are reduced lower than the specified value (in our cases, less than 10-3). The ejector geometry was drawn and meshed as a 2D axi-symmetric model in Gambit software based on a real design of a small ejector applied on the FDS of the fuel cell stack. The problem setup of all cases was finished in the FLUENT software on the personal computer. The calculation of each case was submitted as a job to the SGI Altix supercomputer of Alabama Supercomputer Authority (ASA) based on parallel computation. Normally, the calculation of each case in our study was no more than 8000 iterations. Fourteen cases are submitted to the supercomputers by changing the primary inlet pressure, secondary inlet pressure and secondary inlet RH. As simulation results of one case, the static pressure and Mach number are shown in Figure 5.5 and Figure 5.6, in which the boundary conditions are set as, primary inlet pressure pp,in=3Bar, secondary inlet pressure ps,in=1.42Bar, and RH at the secondary inlet RHs,in=100%. The temperatures of primary and secondary inlets were set to be 293K and 345K, respectively. 122 Figure 5.5 Contours of static pressure in an ejector. Figure 5.6 Contours of Mach number in an ejector. 123 In Figure 5.5, the high pressure of primary flow was decreased rapidly after it passed the primary throat, and a low pressure core is formed at the position after the primary throat. This suction area of the low pressure core can entrain the secondary flow into ejector. After the mixture process, the tow inlet flows have the same pressure at the outlet of the ejector. A pressure plot at the center line is given as shown in Figure 5.7. The pressure was decreased at the throat outlet with a following shock wave which made the vibration curve of the pressure at the centerline. After the mixing process, the pressure of the gas is steady and was increased to the outlet pressure. Figure 5.7 Pressure variation at the centerline of the ejector. In the Figure 5.6, a supersonic flow is formed after the throat of the ejector, where the Mach number is larger than one. A plot of the Mach number at the centerline is shown in Figure 5.8, 124 where the Mach number is about 1 at the outlet of the throat, and increased to supersonic flow. The shockwave is observed in Figure 5.6, and expressed a wavy variation of the Mach number in Figure 5.8. Figure 5.8 Mach number variation along the centerline of the ejector. Since the calculation of mass flow rates of primary and secondary inlets and the outlet are the goal of the modeling work of FDS, these flow rates are calculated in FLUENT by the integral of the gas properties on corresponding surfaces. The detail results of the CFD prediction and the calculation results of the thermodynamic model in previous section for each case are shown in Table 5.2. The inlet temperatures of the primary flow and secondary flow were 293 K and 345 K, respectively, and the outlet pressure is set to be 1.5 Bar. 125 Table 5.2 Comparison of ejector MATLAB model and CFD model pp,in (Bar) ps,in (Bar) RHs,in (%) Wp (kg s-1) (CFD) Ws (kg s-1) (CFD) Wp (kg s-1) (Model) Ws (kg s-1) (Model) 2.5 1.42 1 0.0011982 0.0011951 0.0012228 0.0011068 2.5 1.45 1 0.0011982 0.0014595 0.0012198 0.0010032 2.5 1.48 1 0.0011977 0.0016877 0.001216 0.00090324 3 1.42 1 0.0014675 0.0016662 0.0014725 0.0021797 3 1.45 1 0.0014698 0.00189 0.0014725 0.0021274 3 1.48 1 0.0014672 0.0019737 0.0014725 0.0020696 3.5 1.42 1 0.0017217 0.0020134 0.0017179 0.0025588 3.5 1.45 1 0.0017217 0.0021367 0.0017179 0.0025536 3.5 1.48 1 0.0017169 0.0021318 0.0017179 0.002545 3.5 1.42 0 0.0017122 0.001094 0.0017179 0.00089513 3.5 1.42 0.2 0.0017146 0.0013385 0.0017179 0.0012279 3.5 1.42 0.4 0.0017173 0.0015502 0.0017179 0.0015606 3.5 1.42 0.6 0.001719 0.0017231 0.0017179 0.0018933 3.5 1.42 0.8 0.0017198 0.0018546 0.0017179 0.002226 Comparison between the thermodynamic model and the CFD prediction results based on Table 5.2 are shown in Figure 5.9. In Figure 5.10, the results of secondary flow rates for the thermodynamic model and CFD precondition are compared based on the data in Table 5.2. The simulation result of the secondary flow rates of thermodynamic model has a larger relative error with the CFD simulation results than that of the primary flow rates. 126 Figure 5.9 Primary flow rate comparison between MATLAB and CFD models. Figure 5.10 Secondary flow rate comparison between the MATLAB and CFD models. 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 1 . 6 1 . 8 2 x 1 0 -3 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 1 . 4 1 . 6 1 . 8 2 x 1 0 -3 W p , C F D W p , M o d e 3 0 0 . 5 1 1 . 5 2 2 . 5 3 x 1 0 -3 0 0 . 5 1 1 . 5 2 2 . 5 3 x 1 0 -3 W s , C F D ( k g s - 1 ) W s , M o d e l ( k g s - 1 ) 127 5.4 Pressure Regulator The fuel supply line of low pressure regulator completely supplies the fuel at a low load, and partially supplies the fuel at medium and high load. The low pressure regulator can automatically adjust its outlet pressure to a setting pressure by manipulating the mass flow rate through an embedded nozzle in the regulator as shown in the Figure 5.11. Figure 5.11 Pressure and flow rate variation curve. When the outlet pressure of the regulator, plpr,out is higher than the setting pressure of the regulator, plpr,set, the flow rate through the regulator, Wlpr is 0 or very small in the lockup regime. When the outlet pressure is equal to the setting pressure, the mass flow rate though the regulator reaches the minimum controllable flow Wlpr,min. The pressure of the regulator drops from the setting pressure level as the flow rate increases. The maximum flow rate, Wlpr,max is obtained when the outlet pressure is as low as to make the plug of regulator open completely, where the mass flow rate is not changed for the choking condition as outlet pressure decreases further. 128 In this study, a static low pressure regulator (LPR) was used due to short response time. The scaled mass flow rate ?lpr and the pressure drop ?lpr are defined as, max,lpr lprlpr WW?? and r ef o utlp rs etlp rlp r p pp ,, ??? (5-39) where Wlpr is the flow rate of the low pressure regulator, Wlpr,max is the maximum flow rate, plpr,set is the setting pressure of the regulator, plpr,out is the outlet pressure, and pref is the reference pressure. The scaled mass flow rate is fitted by a polynomial, ),1m i n ( 432231 ccccl pr ???? ???? (5-40) The constants above also can be obtained using a curve fitting method with the experimental data. The parameters for the model are given in Table 5.3. Table 5.3 Constants for the LPR model. Parameter c1 c2 c3 c4 value -116.1 29.77 3.30 0.077 Since the state of hydrogen at the outlet of high-pressure regulator is assumed to be stable, the mass flow rate of the regulator is a function of its outlet pressure, psm. 5.5 Flow Control Valve The duty control solenoid valve in Figure 1.7 is simplified as a flow control valve shown in Figure 5.1 for the high working frequency of the valve. It is assumed that the control valve can 129 response for the control signal in a very short time, and mass flow rate through it can be changed linearly by the control signal as ma x,fcvfcvfcv WuW ? (5-41) where ufrv is the control input signal of the valve that varies from 0 to 1, and Wfcv,max is the maximum mass flow rate. 5.6 Purge Valve The liquid water that accumulates in the return manifold is removed periodically by the purge valve. The flow rate through the nozzle of the purge valve is governed by the nozzle equation as Eq. (5-20), where the outlet pressure of the nozzle is atmosphere pressure. The species flow rates of hydrogen and vapor are, rmHp u rg ep u rg eH yWW ,, 22 ? (5-42) ? ?rmHp u r g ep u r g ev yWW ,, 21 ?? (5-43) When the water activity in the return manifold is greater than one, the liquid water in the return manifold will be purged out by the high speed gas flow through the purge valve. The liquid water and the total water mass flow rate are obtained by, p u r g evrmwp u r g el WaW ,,, )1( ?? (5-44) p u rg elp u rg evp u rg ew WWW ,,, ?? (5-45) 130 The integrated FDS is shown in the Figure 5.12, where the manifolds, blower and fuel cell stack are dynamic models (red blocks), and others are static ones (blue blocks). The state variables in the dynamic models are written after the block names. Figure 5.12 Diagram of integrated model of FDS. There is only one state variable, hydrogen pressure (pem) in the ejector manifold. Two state variables, hydrogen pressure and water activity (pH2 and aw) exist in the supply manifold and return manifold, respectively. The angular velocity (?bl) is the variable state of the blower for hydrogen recirculation. The fuel cell stack model will modeled in next chapter based on the 1D transient model of a single fuel cell. L o w P r e ssu r e R e g u l a t o r F l o w C o n t r o l V a l v e Ej e c t o r B l o w e r ( ? bl ) P u r g e V a l v e S u p p l y M a n i f o l d ( p H 2,sm , a w ,sm ) F u e l C e l l S t a c k ( p H 2,gf c , a w ,gf c , p H 2,gdl , a w ,gdl ) R e t u r n M a n i f o l d ( p H 2,rm , a w ,rm ) W f c v W l pr p em W H 2,e j ,ou t W v ,e j ,ou t W p ,e j , W H 2,s,e j , u f cv p s m y H 2,sm W H 2,bl W v ,bl p rm y H 2,rm W H 2,pu rge W w ,pu rge Ej e c t o r M a n i f o l d ( p em ) W v ,s,e j , W H 2,sm ,ou t W w ,sm ,ou t W H 2,rm ,i n W w ,rm ,i n SR H2 u bl Pu rge S t art I st 131 Chapter 6 Fuel Delivery System Control Now, with the FDS and fuel cell models ready, the state feed-back controller should be designed to meet the control requirement of anode side of the fuel cell. In this section, the FDS with the fuel cell model was linearized at a chosen operating condition. Based on this linear state-space model, a state feed-back controller was designed based on the optimal control theory, and compared with the PI control. Finally, an observer was designed to estimate the immeasurable states based on the real working situation of the fuel cell. 6.1 Anode Model for Controller Design For the control of FDS, the model of anode side of fuel cell should be integrated into the FDS. The modeled fuel cell system in Chapter 3 should be simplified to design the state feed- back controller. For the control of fuel delivery system, the temperature of anode side does not change apparently because of the thermal management system. The transient phase change process (condensation and evaporation) between liquid water and water vapor is neglected. Since the water byproduct is generated at cathode catalyst layer, much water will exist in this layer to hydrate the membrane. Thus, following simplifications are made for the anode model for the controller design as: 132 1. The temperature of the anode side is assumed to be equal to the stack temperature. 2. The liquid water and water vapor are in equilibrium state. 3. The catalyst near the membrane at cathode side is assumed to be fully humidified. 4. Dynamics of the CL at anode side is neglected for its small thickness. 5. Single-layer lumped GDL model is considered. The role of the FDS is mainly to supply the sufficient hydrogen to stack, whereby the inlet fuel pressure is controlled at a given reference and liquid water present in channels is regularly removed. As the FDS mainly interacts with the anode side of the fuel cell, only behaviors of half a cell are considered in this study. A schematic diagram for mass transport of hydrogen and water in the half a cell is depicted in Figure 6.1. The hydrogen supplied to the GFCs diffuses through the GDL (GDL) and reaches the CL (CL). Water is transported across the membrane between the anode and cathode, the amount of which is balanced by electro-osmotic force and back diffusion. The water vapor transported may be condensed in the channels or porous GDL and then becomes a two-phase (water vapor and liquid water) flow of water that affects transport of reactant gas. Figure 6.1 Schematic diagram for mass transport in half a cell of anodic side. G a s f low c ha nne l Ga s dif fusion lay e r Anode c a tal y st lay e r Membr a ne Humidi fie d hy dr o ge n Hy dr og e n a nd wa ter Hy dr og e n tra nsport V a por tr a nsport L iqui d wa ter tra nsport Dissol ve d wa ter tra nsport in m e mbra ne C a thode c a tal y st lay e r 133 6.1.1 Gas flow channels The GFC in the fuel cell stack are regarded as one control volume without consideration of the spatial differences. The dynamics of hydrogen partial pressure and water activity in the GFC of single cell are described by: ? ? outgf cHac tgdlgf cHingf cHgf c stHgf cH WAjWV TRdtdp ,,,,,,, 22222 ??? (6-1) ? ? ou tgf cwac tgd lgf cwingf cwgf csts at stOHgf cw WAjWVTp TRdtda ,,,,,,, )( 2 ??? (6-2) where pH2,gfc and aw,gfc are the hydrogen partial pressure and water activity in the GFC, Tst is the stack temperature, and Vgfc is the total volume of the anode GFC in the single cell. The inlet mass flow rates of hydrogen and water vapor are calculated by the pressure difference between the supply manifold and GFCs as: ? ?g f csmsmHing f cing f cing f cH ppAkW ?? ,,,,, 22 ? (6-3) ? ?g f csmsmving f cing f cing f cv ppAkW ?? ,,,,, ? (6-4) where kgfc,in is the flow coefficient at the inlet of the GFC, Agfc,in is the inlet area of the GFC of a single cell, and the total pressure in the GFC, pgfc is calculated by the hydrogen partial pressure and water activity. In a working fuel cell stack, the pressure drop along the GFCs depends on geometrical dimensions and operating conditions that include the cross-section shape of the channels, fluidic conditions of the channel wall that determines contact angles and surface roughness, the channel 134 patterns such as serpentine, inter-digitized or straight, gas properties and amount of liquid water and flow pattern. For design of a controller, this complex model is simplified by considering, the psm-prm as the total pressure drop along the GFCs. The pressure drops are divided into two parts, inlet pressure drop psm-pgfc and outlet pressure drop pgfc-prm. The relationship between pressure drops and flow rate is further simplified using a constant inlet flow coefficient kgfc,in in Eq. (6-3). This flow coefficient is determined by the channel conditions such as cross section, channel wall condition and patterns. Effects of liquid water effect are neglected in the inlet flow coefficient because little water is expected to be presented in the inlet part of the GFCs. At the outlet, a similar outlet flow coefficient is applied, but effects of liquid water are considered using the relationship given between flow rate and pressure drop described in the following sections. If no purging takes place, no liquid water flows to the stack from the supply manifold. Otherwise, the liquid water entering the GFCs of a single cell is: ing f cvsmwing f cl WaW ,,,,, )1( ?? (6-5) If the water activity in the supply manifold is less than 1, no liquid water exists in the manifold and the Wl,gfc,in is 0. Then, the total water mass flow rate entering a single cell Ww,gfc,in in Eq. (6-2) is equal to the sum of Wv,gfc,in and Wl,gfc,in. The total inlet flow rates of multiple cells of a stack are equal to the total mass flow rate leaving the supply manifold: ing f cHcello u tsmH WNW ,,,, 22 ?? , ing fcwcello u tsmw WNW ,,,, ?? (6-6) Similarly, the mass flow rates of gas species at the outlet of GFCs are: 135 ? ?rmg f cg f cHoutg f coutg f coutg f cH ppAkW ?? ,,,,, 22 ? (6-7) ? ?rmg f cg f cvo u tg f co u tg f co u tg f cv ppAkW ?? ,,,,, ? (6-8) where kgfc,out is the outlet flow coefficient that is not a constant. In the model, the amount of liquid water is considered as the dominant variable that affects the relationship of outlet pressure drop and flow rate. Increase of liquid water in the channels leads to an increase of the flow resistance and as a result the flow coefficient becomes smaller. The flow coefficient at the outlet of the GFC can be described as: 5.1,, )1( gf cingf cou tgf c skk ?? (6-9) where sgfc is the liquid water saturation, which is defined as the liquid volume fraction in the control volume, and calculated by: satl wsat as ??? ? ?? )1( (6-10) where ?sat is the density of saturated vapor calculated by the saturation pressure and temperature in the control volume. The liquid water in the outlet of the fuel cell stack is entrained by the viscous force of the gas flow in the channels. The mass flow rate, Wl,gfc,out, is expressed as: outg f cglg f cgg f cg f coutg f clg f ctpoutg f cl Us sACW ,,, 2 ,,,, 1 ? ?? ?? ? ??? ? ? ?? (6-11) 136 where Ctp,gfc is the correction factor, ?g,gfc and ?l are the viscosity of gas in GFC and liquid water, and Ug,gfc,out is the superficial velocity of gas at the outlet of the channel. The viscosity of the mixed gas of hydrogen and vapor is obtained by the average value given by the semi-empirical formula proposed by Wilke [46], where the superficial gas velocity is given by: )( ,,, rmg f co u tg f co u tg f cg ppkU ?? (6-12) Then, the water mass flow rate of a single cell at its outlet is the sum of the flow rates of water vapor and liquid water. Then, the total mass flow rates exhausted from all cells of the stack are: o u tg f cHc e l linrmH WNW ,,,, 22 ?? , o u tg fcwcellinrmw WNW ,,,, ?? (6-13) Now, unknown variables in the dynamic equations, Eqs (6-1) and (6-2), are the species fluxes between the GFC and GDL that are derived in the following section. 6.1.2 Gas diffusion layers For the model of GDL, the volume of pores in the GDL is considered as an isothermal control volume. Then, the dynamics for the GDL are described as follows: ? ? clgd lHgd lgf cHgd lgd lgd l stHgd lH jjs TRdtdp ,,,,, 2222 )1( ??? ?? (6-14) ? ? clgd llclgd lvgd lgf clgd lgf cvgd lgd lsts at stOHgd lw jjjjTp TRdtda ,,,,,,,,, )( 2 ???? ?? (6-15) 137 where ?gdl is the porosity of the GDL, ?gdl is the thickness of the layer and j is the mass flux through the layer. The hydrogen and water vapor fluxes from the GFC to the GDL are calculated using the Fick?s law as: )( ],[],[,,,],[ g d lg f cg d lg f cmg d lg f c CChj ??? (6-16) where C[ ],gfc and C[ ],gdl refer to the concentrations of gas species in the GFC and GDL and are obtained by the hydrogen partial pressure and water activity. The mass transfer coefficient between GFC and GDL is given by: gdlmg fcm gdlmg fcmgdlg fcm hh hhh ,, ,,,, ??? (6-17) where the mass transfer coefficients in GFC and GDL are: gf ch gf cOHHgf c ac t c ov e rAgdlgf cgf cm dDShA rAh , ,,,, 22)1( ??? , gd l gd lOHHgd lm Dh ? ,, 222 ?? (6-18) where the Agfc,gdl is the area of gas contact interface between GFC and GDL for single cell, Aact is the active area of the fuel cell, Shgfc is the Sherwood number of the channels, dh,gfc is the channel hydraulic diameter, and rA,cover is the liquid water cover ration on the gas contact interface between GFC and GDL, that is defined and calculated as: )]1(,1m i n [ c ov, ,,, c ov g f cerg dlg f c g dlg f cl i q ui derA sKAAr ??? (6-19) 138 where Aliquid,gfc,gdl is the liquid water cover area on the interface between the GFC and GDL. The binary diffusion coefficient DH2-H2O in Eq. (6-18) is [7]: ? ? ? ??? ? ??? ? ??? ?? ??? ? ??? ? ? ? ? ? CLandGDLf o rpTsK G F Cf o rpTK D OHH OHH OHH 5.1 5.25.1 5.1 )1( 22 22 22 ? (6-20) where KH2-H2O is a constant, ? is the porosity of layers, s is the liquid saturation, T is temperature (K) and p is the pressure (Pa). Since the CL is much thinner than that of the GDL, the characteristic diffusion time in the CL given by ?2/D is also much shorter than that in the CLs. Thus, the CL is considered as the boundary for both the GDL and membrane. According to the mass conservation principle, the hydrogen flux from the GDL to the CL is equal to that consumed in the CL as: d enOHclg d lH IF Mj ?? 2 2 2 ,, (6-21) where Iden =Ist/(Ncell?Aact) is the area current density of the fuel cell stack. The vapor transfer from the GDL to the CL is given by: )( ,,,,,, clvg d lvclg d lmclg d lv CChj ??? (6-22) where Cv is the vapor concentration obtained from the water activity. Liquid water transport at the interface between the GDL and GFC is driven by the capillary pressure in pores. The liquid water flux is determined by the difference between the critical pressure in the GFC and the capillary pressure in GDL: 139 ? ? g d lcc r i tcg d lc a plg d lg f cl pphj ,,,,,. ??? (6-23) where hl,cap,gdl is the liquid water transfer coefficient in GDL, and pc,gdl is the capillary pressure in the GDL written as: lg dl g dlrg dllg dlc a pl KKh ??? ,,, 2?? (6-24) ? ?32 263.112.2417.1c o s SSS Kp cwc ??? ??? (6-25) where K is the permeability, Kr is the relative permeability, ?w is water surface tension, ?c is contact angle, S is the reduce liquid saturation. The Kr and S can be written as: 3)( imr ssK ?? , im imsssS ???1 (6-26) where sim is the immobile liquid saturation. Likewise, the mass flux of liquid from the GDL to the CL is determined by the capillary pressure difference between these two layers, given by: )( ,,,,,, clcg d lcclc a plclg d ll pphj ?? (6-27) where pc,cl is the capillary pressure in the CL, and hl,cap,cl is the liquid water transfer coefficient at the boundary of GDL near the CL, where it is assumed that the capillary pressure is continuous as: 140 clcclgdlc pp ,,, ? (6-28) where pc,gdl,cl is the capillary pressure of GDL at the boundary near the CL. Now, if the capillary pressure pc,cl is known, the reduced liquid saturation of GDL near the CL, Sgdl,cl is calculated from pc,gdl,cl by using the inverse function of Eq. (6-25), and as a result the liquid water transfer coefficient at the boundary of GDL near CL, hl,cap,cl, is obtained using Eqs. (6-24) and (6-26). Now, the mass flow rates of gas species and liquid water at the both boundaries of GDL in Eqs. (6-14) and (6-15) are calculated by using Eq. (6-16), (6-21), (6-22), (6-23) and (6-27). The unknown variables in these equations are the water vapor concentration, Cv,cl, and capillary pressure, pc,cl, in the CL that are a function of the water activity in the CL. The water activity in the CL can be calculated by the water balance in the GDL, CL and membrane described in following section. 6.1.3 Membrane Because of the thin thickness of the CL, the boundary between the GDL and the membrane, the water flux across the CL is assumed to be continuous and described as follows: ? ? caanpe mpe m pe mOHde nOHdr agclgd llclgd lv EW DMIFMnjj ??? ? ? ???? 22 ,,,, (6-29) wherr ?, ndrag and D? are calculated by equations of membrane model depicted in section 2.7 of Chapter 2. Hence, the terms on the left are the total water mass flux from the GDL to the CL. The terms on the right are the total water flux from the CL to the membrane, where the first one represents the water flux driven by the electro-osmotic drag force, and the second one is that driven by the 141 gradient of the water content in membrane. In fact, water is mostly generated in the cathode CL. Thus, the water content of the membrane near the cathode CL ?ca in Eq. (6-29) is assumed to be the liquid maximum value ?lmax with aw=3 as shown in Figure 2.11. Hence, the right hand side of Eq. (6-29) is a function of water activity in the anode CL when the current density Iden and stack temperature Tst are known. Then, the water flux across the membrane and the water activity in the anode CL can be solved by substituting Eqs. (6-22) and (6-27) into left hand side of the Eq. (6-29). All components along with the integrated fuel cell model are shown in Figure 6.2. The whole fuel cell system includes four state variables, and the input variables are the gas states in the supply manifold, return manifold and stack current. Figure 6.2 The diagram of integrated model of fuel cell stack. G a s f l o w c h a n n e l s ( p H 2,gf c , a w ,gf c ) p s m y H 2,sm p rm j H 2,gf c ,gdl j w ,gf c ,gdl W H 2,sm ,ou t W w ,sm ,ou t W H 2,rm ,i n W w ,rm ,i n SR H2 I st SR H2 C H 2,gf c a w ,gf c G as d i f f u si o n l a y e r ( p H 2,gdl , a w ,gdl ) a w ,gdl M e mb r a n e ? ca a w ,c l 142 6.2 Analysis of Integrated System of FDS and Stack For the purpose of controller design, it is assumed that FDS works in three modes: low, medium and high currents. I is defined that the low, medium, and high current working conditions refer to the ranges of 0-6000A m-2, 6000-8000A m-2, and above 8000A m-2, respectively. In the low mode, the flow control valve and ejector are shut down, and the low- pressure regulator supplies the hydrogen. The blower is the pump to re-circulate the exhausted hydrogen. For medium and high current of the stack, the flow control valve and ejector start work to supply and re-circulate more hydrogen. In low mode, the anode gas pressure is passively stabilized near its setting value by the low- pressure regulator, and the blower is used to adjust the recirculation mass flow rate to control the total flow rate of hydrogen fed to the stack. In medium and high modes, the anode gas pressure and hydrogen flow rates supplied to stack are expected to be controlled by the flow control valve and blower simultaneously. In the FDS model, there are 6 state variables including ejector manifold pressure, hydrogen pressure and water activity in supply manifold, hydrogen pressure and water activity in return manifold, and angle velocity of blower rotator. Four additional state variables exist in the anode model of the fuel cell stack, which are the hydrogen pressures and water activity in the GFCs and GDL, respectively. The control inputs of this integrated system involve the control signal of the flow control valve and blower control voltage. The stack current and purging operation are considered to be disturbance inputs. The FDS with the fuel cell model is simulated using MATLAB with a simple PI controller to find its steady state values with the parameters shown in Table 6.1. It is found that the water activities in supply and return manifolds cannot be stable because of the water condensation in 143 the manifolds. Thus, it is assumed that the water activities are assumed to be constant 1 for steady state analysis. Table 6.1 Model parameters. Symbol Value Symbol Value Wfcv,max 2.4?10-3 kg s-1 Tst 353 K Vem 2.5?10-3 m3 ?gdl 250 ?m Tem 293 K ?gdl 0.6 At,ej 8.04?10-6 m2 ?c,gdl 110? Am,ej 4.07?10-5 m2 Kgdl 1.76?10-11 m2 ?p 0.64 ?pem 2?103 kg m-3 ?s 0.9 ?pem 60 ?m ?exp 0.60 EWpem 1.1 kg mole-1 ?mix 0.90 pc,crit 800 Pa pset,lpr 1.5?105 Pa Jbl 2.6?10-3 kg m2 Wlpr,max 1.75?10-3 kg s-1 kt 0.15 N m A-1 Vsm 4?10-3 m3 kv 0.15 V rad-1 s Tsm 318 K Rbm 0.82 ohm Vrm 4?10-3 m3 ?bm 0.9 Trm 338 K Dbl 0.15 m Ncell 381 ?l 986 kg m-3 Aact 576 cm2 ?w 6.25?10-2 N m-1 Agfc,in 8 mm2 KD0,H2-H2O 1.7232?10-3 Agfc,gdl 288 cm2 Ctp,gfc 0.5 sim 0.01 Kcover 1.5 Shgfc 2.78 At,purge 5?10-6 m2 kgfc,in 0.002 m s-1 Pa-1 ?purge 0.81 Vgfc 2.8?10-5 m3 144 As one of the control objectives, the supply manifold pressure psm (Bar) for medium and high current mode is a function of current density as: )6 0 0 0(10249.1 6 ???? ? d e nsm Ip (6-30) The reference SRH2, instead of the flow rate, is set to constant 1.5 for all three modes. Ten operating points including the steady state variables and control signals are obtained for different current densities (from 1000Am-2 to. 10000A m-2 by step of 1000A m-2) as shown in Figure 6.3. The control signal the of blower ubl (0 to 350V) is normalized to ubl,n (0 to 1). 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 - 0 . 2 0 0 . 2 0 . 4 0 . 6 0 . 8 I d e n ( A m - 2 ) u u f c v u b l , n 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 1 2 3 4 5 I d e n ( A m - 2 ) p e m ( B a r) 145 Figure 6.3 Control and state variables at steady state for different current densities. As shown in Figure 6.3, the control voltage of blower motor and the angle velocity of the blower increase with current density in low mode, while they decrease in the medium mode and then increase in high mode. Since the ejector starts the recirculation in medium and high mode, the power that the blower consumed on the recirculation is reduced. The direct calculation for 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 1 . 1 1 . 1 5 1 . 2 1 . 2 5 1 . 3 1 . 3 5 1 . 4 1 . 4 3 I d e n ( A m - 2 ) p H 2 ( B a r) S u p p l y m a n i f o l d R e t u r n m a n i f o l d G a s f l o w c h a n n e l G a s d i f f u s i o n l a y e r 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 . 9 9 7 0 . 9 9 8 0 . 9 9 9 1 1 . 0 0 1 I d e n ( A m - 2 ) a w G a s f l o w c h a n n e l G a s d i f f u s i o n l a y e r 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 2 . 5 5 7 . 5 10 1 2 . 5 15 I d e n ( A m - 2 ) ? ( k R P M ) 146 control signals from the current density can form a static feed-forward (SFF) control based on steady state analysis. 6.3 Design of State Feed-back Control with an Observer Because of the shutdown of the supply line of the flow control valve and ejector recirculation at the low current mode, the FDS is operated with one supply line and one recirculation loop, which is considered to be a SISO system with the disturbance signals. A PI controller (KP,SR = 3, KI,SR=6) is applied on the system to control the SRH2 at the low mode of the FDS. At the medium and high modes, the integrated system of the FDS with the fuel cell stack can be written in the form of state equations as: ?? ??? ? ? ? ),,( ),,( ),,( wuxhz wuxgy wuxfx? (6-31) where x is the vector of state variables, u is the vector of control inputs, w is the vector of disturbance inputs, y is the vector of measurable outputs and z is the vector of control objectives. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? st blf c v Hsm T g dlwg dlHg f cwg f cHblrmwrmHsmwsmHem Iw uuu SRpz apapapappx , , ,,,,,,,,, 2 2222 ,,,,,,,, ? (6-32) The hydrogen pressure in supply manifold is considered as another control objective that should follow the pressure at cathode side. From the simulation results of an output feed-back controller (PI controller), it is found that water activities in supply and return manifolds are not 147 convergent and increase slowly with time. The water activities in the GFCs and the GDL are almost 1 and not measurable and observable. Thus, the dynamics of water activities is neglected in the controller design, and the values of these 4 water activities are assumed to be constant 1. In addition, the hydrogen pressure in the GDL is also not observable, and is removed from the state vector. Now, the current state vector becomes: ? ?TblrmHg f cHsmHem ppppx ,,,,, ,,, 222 ?? (6-33) where the pressures are in the units of Bars, and the angle velocity of blower ?bl is in units of kRPM. The measurable output vector y in Eq. (6-31) for the controller design in medium and high modes is chosen as: To utsmblrmsmem Qpppy ),,,,( ,?? (6-34) where the pressures pem, psm and prm are in the unit of Bar, and Qsm,out is the flow rate at the outlet of the supply manifold in the units of SLPM, which is also the flow rate entering the stack Qst,in. The blower angle velocity at steady state is shown in Figure 6.4, and other values of y at steady state for different stack current densities are shown in the Figure 6.3. Hence, the steady state values of the vector y can be obtained by the interpolation of these curves, and are considered to be an additional output of the static feed-forward block. In the vector z, the SRH2 at the inlet of the fuel cell is not measurable directly, which should be estimated by the Qsm,out as: 148 sm sms a tsm r e a c t e dH o u tsmH p TppQ QSR )(101 5 , , 2 2 ???? (6-35) where psm (Bar) is the total pressure in supply manifold, QH2,reacted (SLPM) is the rate of hydrogen consumption, which can be obtained by the equation of WH2,reacted=Ncell?MH2?Ist/(2F). Figure 6.4 Measurable output at steady state for different current density. A block diagram for a state feed-back controller designed is depicted in Figure 6.5, where the cathode pressure and SRH2 are the state variables to track. u* and y*is the steady state value that is calculated by the feed forward block, and z* is the reference values of the psm and SRH2. The SR estimator block outputs the psm that is measured by sensor directly, and SRH2, which is 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 1 . 4 1 . 4 5 1 . 5 1 . 5 5 1 . 6 I d e n ( A m - 2 ) p s m a n d p rm ( B a r) S u p p l y m a n i f o l d R e t u r n m a n i f o l d 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 500 1000 1500 2000 2500 3000 I d e n ( A m - 2 ) Q s t, in ( S L P M ) 149 obtained by Eq. (6-35). A state observer is used to estimate the perturbation of the state variables of the system. The K and L are the controller gain and observer gain obtained by the Linear- Quadratic-Gaussian (LQG) design. Figure 6.5 State feedback control with integral and observer. The state equation above for the FDS and stack is nonlinear and needs to be linearized for design of controllers and observers. Operating points for the systems are chosen at two currents, Iden=7000 A m-2 and Iden = 9000 A m-2, that represents the medium and high current mode. Thus, two linear equations are obtained for each of the modes. At these operating points, the linearized system of FDS is written in the form of state equations as: ?? ??? ??? ??? ??? wHuHxGz wDuDxCy wBuBxAx wu wu wu ???? ???? ???? ? (6-36) where ?(?)=(?)op?(?) is a perturbation describing the difference between the state variables and their the steady state values at the operating points, and the A, Bu, Bw, C, Du, Dw, G, Hu and Hw 150 are the system matrices. For the medium and high modes, the two linear dynamic models are obtained to calculate the controller and observer gains. The Linear-Quadratic-Integral control is applied in the state feedback control as shown in Figure 6.5, where the controller form is as: ];[ ixxKu ?? (6-37) where the xi is the integrator output. The controller gain K to minimize the cost function is optimized using the method of Linear Quadratic Regulator (LQR). The cost function is: ? ?? ? ??? 0 dtuRuxQxzQzJ TiiTizT ???? (6-38) where Qz, QI and R are the weighting matrices for the control objectives, error integral and control inputs. The same weighing matrices of the cost function are used for the controller design in the systems of medium and high current modes, which are given as: ?? ??? ??? ??? ??? ])101,101([ ])101,101([ ])101,101([ 42 810 68 d i a gQ d i a gQ d i a gQ R i z (6-39) where the matrices are all in diagonal form. The control gains K to minimize the cost function are obtained based on the linearized systems by the MATLAB command as K = lqi(ss(A,Bu,G,Hu),Q, R), where Q =blkdiag( CTQzC, Qi) is a 5-by-5 block diagonal matrix. The calculation results for the medium and high current modes are: 151 ?????? ??? ????? 77831169480 9 3 7 1 8096513315476 2 9 554 1 6 9 40 1694883 1 7 79 0 4 614732871 5 8 77871000727 . . . . . . . . . . . . . .K m i d ?????? ??? ???? 16964627442930357168478695715300720 46274961637371840961893222903097 . . . . . . . . . . . . .K h i g h (6-40) A linear observer is used to estimate the state variable by the measurable outputs as shown in Figure 6.5. The inputs of the observer are perturbations ?u and ?y=y-y*, and the output of the observer is estimated states perturbation ? . The structure of the observer is: )?(?? uDxCyLBuxA dt xd uu ????? ???? (6-41) where A, Bu, C, and Du are the linear system matrices in Eq. (6-36) for medium and high currents. To minimum the noise effect on the measurement signals and state estimation, the continuous Kalman filter design method is used to obtain the observer gain matrix L. The process and measurement noises covariance for different current are given as Qn = 10 and Rn = diag([1?10-4, 1?10-4, 1?10-4, 1?10-4, 1]). Then, the observer gains L (5-by-5 matrix) for medium and high modes are calculated by the MATLAB command [kest,L,P] = kalman(ss(A,[Bu,Bw],C, [Du,Dw]),Qn,Rn). Another disturbance signal affecting the responses of the controllers other than stack current is the purging operation that is not considered in the discussion of controller design. To prevent accumulation of impurities and contaminant species which lead to output voltage loss, the hydrogen compartment purging is necessary for the long-term operation. The purge time and purge stop time are the effective parameters in purge process controlling [87], which is a dynamic process determined by the concentrations of impurities and liquid water 152 volume in the system. In the current system with stable SRH2, the concentration of the impurities and the liquid water volume are assumed to be proportional to the amount of hydrogen consumed, which is calculated by the integration of the stack current. Thus, the control strategy of the purging process of the FDS is dependent on the integration of the stack current as: ?? ??? ?? o pe n c r i t t t st t EdtI s h u t d o w n ofl e n g t ht i m ef o ro p e n sita f t e rv a l v ep u r g et h eS h u t d o w n if v a l v ep u r g e O p e n t h e (6-42) where Ecrit is a constant that should be determined experimentally. The integration of the stack current is reset by the last shutdown of the purge valve as shown in Figure 6.6. In the current research, the integration of current density instead of stack current is used to control the purge valve. The critical integration of the current density Ecrit=5000A?sec m-2 and the opening time topen= 1sec are used in the model simulation. Figure 6.6 The state of purge valve. t o p en t c l o se t sh u t d o w n t st a rtt st a rt t sh u t d o w n O 1 c r i t t t st EdtI s h u td o w n ?? Pur ge Stat e 0?? t t st s h u td o w n dtIr e s e t t 153 For different current modes, the controller structure and gains will be switched based on the current request. The track performance of the pressure of supplied fuel and the stabilization of SR of hydrogen of the FDS will be investigated with variations of disturbance signals including stack current and purging operation in the next section. 154 Chapter 7 Analysis of Simulation Results of FDS with Feed-back Controllers The FDS with designed controllers (SFF, PI and SFB) were run in the MATLAB/Simulink environment under the step change of the stack current density input for different currents to study the pressure and SR responses. The integrated system of SFF control and FDS model is shown in Figure 7.1. The SFF block is implemented based on the steady state analysis of control signals of FDS, which is directly calculated from the stack current by interpolating the control signal curves shown in Figure 6.3. Figure 7.1 SFF control of FDS. FDS w ith An o d e Mod el (10 state vari ab le s) z = [ p sm , S R H2 ] T P u rg e Con trol O pe n or c l os e I st u * SFF 155 A decentralized PI controller (ufcv-psm loop: KP,psm=40; KI,psm=80; and ubl,n-SRH2 loop: KP,SR = 3, KI,SR=6) is designed to compare with the SFB control in medium and high current load as shown in Figure 6.2. In low current mode, the FDS is controlled by a SISO PI control with gains (Kp,SR=3 and KI,SR = 6 for ubl,n-SRH2 loop). The SFF block and SR estimator block are used to obtain the reference control signals and estimated control objectives. Figure 7.2 PI control of FDS. The SFB control system of the FDS is shown in Figure 7.3, where the SFB control is implemented based on the block diagram shown in Figure 6.5. The SFF block is also used in the SFB control to obtain a reference control signals for SFB control. The control gain K and observer gain L for medium and high loads in the SFB block are tuned by the stack current value. When the stack current is in the range of low mode, the PI controller is switch on to connect the FDS automatically. Thus, the FDS is controlled by three sets of control parameters with two different control structures for three ranges of current loads of fuel cells tack. FDS w ith An o d e Mod el (10 state vari ab le s) z = [ p sm , S R H2 ] T y = [ p em , p sm , p rm , ? bl , Q sm , ou t ] T SFF P u rg e Con trol O p e n o r c l o s e u = [ u f cv , u bl ] T I st z * = [ p ca , S R H2 * ] P I Cont r ol l er u * S R E s ti m ato r + _ K p K I ? + + + + e ? u z est 156 Figure 7.3 SFB control of FDS. In following sections, the system responses of the single step and multi-step changes of the stack currents were simulated and investigated for the comparison of these three different controllers under purging operations. 7.1 Response of Step Change Current In the case of step change responses, the decentralized PI controller (ufcv-psm loop: KP,psm=40; KI,psm=80; and ubl,n-SRH2 loop: KP,SR = 3, KI,SR=6) and static feed forward controller (SFF) are compared with the SFB controller in the medium and high currents. For the medium current mode, a step change of the current density from 7000A m-2 to 7100A m-2 occurs at 10th second. The responses of psm and SRH2 for this step change are shown in Figure 7.4, where the reference value of psm is calculated using Eq. (6-30). The psm response of SFB controller reaches the steady value in about 1 second with small overshoots, while the PI and SFF controllers require about 2.5 seconds to stabilize the psm. The SRH2 response in Figure 7.4 shows that the SFB and PI FDS w ith An o d e Mod el (10 state vari ab le s) z = [ p sm , S R H2 ] T y = [ p em , p sm , p rm , ? bl , Q sm , ou t ] T SFF SFB (M i d d l e a n d h i g h m o d e ) PI (l o w m o d e ) P u rg e Con trol O p e n o r c l o s e u = [ u f cv , u bl ] T I st z * = [ p ca , S R H2 * ] Chan ge K an d L y* u * 157 controllers both enable the SRH2 to return to 1.5 in about 1.2 second, while the SFF controller can stabilize the SRH2 with a steady state error. Figure 7.4 Step change response of medium mode controller. Figure 7.5 shows the responses of psm and SRH2 in the high current mode for the step change of current from 9000A m-2 to 9100A m-2 at 10th second. The settle time of the SFB is only about 2.5 seconds, while that of PI controls is more than 3 seconds for the response of psm. In fact, the SFF could not reach the desired steady state without the feedback signal. In the response figure of SRH2, the response curve of SFF is still not convergent to the reference value, while the settle time of SFB control (about 0.5 second) is smaller than that of PI control. 5 7 . 5 10 1 2 . 5 15 1 7 . 5 20 1 . 4 9 1 8 1 . 4 9 1 9 1 . 4 9 2 1 . 4 9 2 1 1 . 4 9 2 2 1 . 4 9 2 3 T i m e ( s e c ) p s m ( B a r) S F B PI S F F R e f 5 7 . 5 10 1 2 . 5 15 1 7 . 5 20 1 . 4 8 1 . 4 9 1 . 5 1 . 5 1 T i m e ( s e c ) SR H 2 S F B PI S F F 158 Figure 7.5 Step change response of high mode controller. 7.2 Response of Multistep Change Current In the working fuel cell stack applied on the vehicles, the stack current and cathode pressure vary continuously in the three modes along with periodic purging operations. To mimic real operations, a multi-step stack current density is applied to the FDS with the fuel cell stack and at the same time a multi-step noise is added to the cathode pressure (pca) that calculated by Eq. (6-30) as the dotted line in Figure 7.6. The solid line shown in second figure of Figure 7.6 represents the real reference pressure of the anode supply manifold. In addition, the purging valve is controlled by on-line calculations as shown in Figure 7.6 depending on the variations of 5 7 . 5 10 1 2 . 5 15 1 7 . 5 20 1 . 4 9 5 8 1 . 4 9 5 9 1 . 4 9 6 1 . 4 9 6 1 1 . 4 9 6 2 1 . 4 9 6 3 T i m e ( s e c ) p s m ( B a r) S F B PI S F F R e f 5 7 . 5 10 1 2 . 5 15 1 7 . 5 20 1 . 4 8 5 1 . 4 9 1 . 4 9 5 1 . 5 T i m e ( s e c ) SR H 2 S F B PI S F F 159 the stack current. The results show that the time interval between sequent purging operations gets shorter in the high current range, but larger in the low current range. Figure 7.6 Multi-step changes of stack current, cathode pressure and purging state. 0 10 20 30 40 50 60 2000 3000 4000 5000 6000 7000 8000 9000 10000 T i m e ( s e c ) I d e n ( A m -2 ) 0 10 20 30 40 50 60 1 . 4 8 9 1 . 4 9 1 1 . 4 9 3 1 . 4 9 5 1 . 4 9 7 1 . 4 9 9 1 . 5 0 1 1 . 5 0 3 1 . 5 0 5 p c a ( B a r) E q . ( 5 5 ) R e f 0 10 20 30 40 50 60 0 1 T i m e ( s e c ) S ta te o f P u rg e V a lv e Eq.(6-30) 160 The system initial values of the hydrogen partial pressures and angle velocity are given by the interpolation of curves of the steady state from the current density, while the initial value of pem for low mode for the given multi-step change of current density is 1.4 bar. All initial values of water activities are set to be 10 to observe the purging effect on the liquid water. The responses of psm under the change of stack current and purging operating for the three controllers are shown in Figure 7.7, where the PI gains are the same as those in the previous case. When the load current is low, the PI gains for SRH2 control are the same as that used in the SFB control for the SISO system. Thus, the psm at the beginning and end of the stack current curve is stabilized by the low-pressure regulator as shown in Figure 7.7. When the load current changes between the medium and high modes, the two controllers and observers are automatically switched. The actual value of the pressure, psm, tracks the reference curve by the controllers. The settle time of the SFB control is the smallest, while the overshoot of the PI control is the largest among the three controllers. The step response of the SFF shows that the control cannot follow the reference value of the supply manifold pressure at steady state. At purging, it takes about 1 second for the SFF control to reach the steady state, but with a constant error, while SFB and PI controls are able to reach the reference value as shown in Figure 7.7, while purging is operating at the 17th second. The settle time of the SFB control is about 0.7 second after the start of the purge valve, and 0.3 second after the shutdown of the purge valve, while a longer settle time of the PI control is observed, as shown in Figure 7.7. 161 Figure 7.7 The supply manifold pressure response for the multi-step current change. With respect to response of SRH2 shown in Figure 7.8, where performances of three controllers are different. The SFF control has the largest steady state error and overshoot among others, while the SFB and PI controls have almost the same settle time. The settle time of the SFB and PI controls are comparable, but the SFB control has a slightly higher overshoot than the PI control, while the SFB control converges more rapidly to the reference, especially when purging is in operation as shown in Figure 7.8. The SFB control allows for stabilizing the SRH2 to a reference within about 0.3 second after purging, while it takes about 1 second for the PI control to reach the reference. 0 10 20 30 40 50 60 1 . 4 8 5 1 . 4 9 1 . 4 9 5 1 . 5 1 . 5 0 5 1 . 5 1 T i m e ( s e c ) p s m ( B a r) S F B PI S F F R e f 17 1 7 . 5 18 1 8 . 5 19 1 9 . 5 20 1 . 4 9 4 1 . 4 9 5 1 . 4 9 6 1 . 4 9 7 1 . 4 9 8 1 . 4 9 9 1 . 5 1 . 5 0 1 T i m e ( s e c ) p s m ( B a r) S F B PI S F F R e f 162 Figure 7.8 The SRH2 response of multi-step current change. As shown in Figure 7.9, water activities in the supply manifold rapidly drop to about 1 at the first purging for three controllers. Since the water activity is 1 when the water vapor is saturated and there is no liquid water, the Figure 7.9 indicates that there is little liquid water in the manifold. Then, the water activities change in the range of 1 to 2 periodically with the purging process when the SFB and PI control are applied. In contrast, the SFF control shows a highest peak value of the water activity for the SFF control after the first purging. As a result, the amount 0 10 20 30 40 50 60 0 . 8 1 1 . 2 1 . 4 1 . 6 1 . 8 2 2 . 2 2 . 4 T i m e ( s e c ) SR H 2 S F B PI S F F 17 1 7 . 5 18 1 8 . 5 19 1 9 . 5 20 1 . 4 6 1 . 4 8 1 . 5 1 . 5 2 1 . 5 4 1 . 5 6 1 . 5 8 T i m e ( s e c ) SR H 2 S F B PI S F F 163 of liquid water in the supply manifold is very small under the current purging strategy at the normal working conditions of a fuel cell stack. It should be noted that the water activity in the supply manifold frequently becomes less than 1 during purging, which affects the estimation of the SRH2 as shown in Eq. (6-35), where the supply manifold water activity is assumed to be greater than 1. The error of the estimation is shown in Figure 7.9, where the estimated SRH2 gets lower than the real value when purging. 0 10 20 30 40 50 60 0 1 2 3 4 5 6 7 8 9 10 11 12 T i m e ( s e c ) a w ,s m S F B PI S F F 11 12 13 14 15 16 17 18 19 1 . 4 2 1 . 4 4 1 . 4 6 1 . 4 8 1 . 5 1 . 5 2 1 . 5 4 T i m e ( s e c ) SR H 2 S F B S F B e s t i m a t o r 164 Figure 7.9 The water activity responses in manifolds of multi-step change of stack current. The water activity in the return manifold varies like a triangle wave as plotted in Figure 7.9, where the range of water activity varies from 9 to 15. The peak values of the water activities in the return manifold decrease with time. Thus, the liquid water amount in the return manifold is larger than that in supply manifold, and purging can restrict and reduce the amount of liquid water under the current purging control strategy. As shown in Figure 7.10, the water activities in the GFC decreases to 1 in about 40 seconds because of the purging that decreases the water activity in the GFC in the first 40 seconds. However, the purging does not directly affect the water activity in the GDL. After the water activity of the GFC has reached to 1, the water activity in the GDL drops from 10 to 1 in about 12 seconds. In the last seconds, the water activities in the GFCs and GDL becomes near the constant 1, which indicates that a small amount of liquid water may appear in the anode of the fuel cell under the normal working conditions of current variations and purging. 0 10 20 30 40 50 60 8 10 12 14 16 18 T i m e ( s e c ) a w ,r m S F B PI S F F 165 Figure 7.10 The water activity responses in fuel cells of multi-step change of stack current. From the above analysis, the SFF control always leads to steady state error and the largest amount of liquid water accumulated in the supply and return manifolds. The SFB control can make the anode supply manifold pressure track the cathode pressure with the smallest response time and overshoot among the three controllers, and stabilize the SRH2 at the reference value with shortest settle time and an acceptable overshoot under the multi-step change of stack current and periodic purging. Since the track performance of cathode pressure is more important than the stabilization of SRH2, the SFB control outperforms the other controls on the rejection of the disturbance signals of stack current and purging. The SFB control discussed above can be implemented in a microcontroller as shown in Figure 7.11. The stack current, cathode pressure and the measurement output vector y of the plant are inputted to a micro-controller by corresponding sensors and data acquisition system from the fuel cell stack and FDS respectively. The blocks of SFF, SFB, PI and purge control can be implemented by programming in the controller to determine the values of control signals and purge valve state that are input signals of FDS. 0 10 20 30 40 50 60 0 1 2 4 6 8 10 12 T i m e ( s e c ) a w ,g fc a n d a w ,g d l a w , g f c , S FB a w , g f c , P I a w , g f c , S FF a w , g d l , S FB a w , g d l , P I a w , g d l , S FF a w , g d l a w , g f c 166 Figure 7.11 Diagram of the control system for a real FDS and fuel cell stack. z = [ p sm , S R H2 ] T y = [ p em , p sm , p rm , ? bl , Q sm , ou t ] T SFF SFB (M i d d l e a n d h i g h m o d e ) PI (l o w m o d e ) P u rg e Con trol O pe n or c l os e u = [ u f cv , u bl ] T I st z* Chan ge K an d L Fue l Ce ll Stack SR H2 * = 1.5 p ca FDS ( Pla nt) S e n so rs F DS c on tr ol s i n m i c r o - c on tr ol l er 167 Chapter 8 Conclusion In my research work, control strategies for a new hybridized FDS were proposed. The strategies were designed based on a one dimensional segmented and layered transient model of a single cell that considers two-phase and non-isothermal effects. The model was control-oriented and capable of presenting static and dynamic behaviors of physical and electrochemical properties that includes concentrations of species, phase changes, pressure drops, temperature distributions, as well as current distributions along channels. The occupation of liquid droplet in gas flow channels was originally considered in the mass and heat transfer for the transient model based on wall condensation and evaporation. Performance of the model for a cell was compared with a single cell that was designed with respect to I-V characteristics and responses of step currents. The design of combination of segmentation and water visualization was firstly applied in the single fuel cell experiments to investigate the current and liquid water distributions along GFCs. The model has shown a fairly good match with the experiments, which can be used for analysis of the FDS system. A dynamic stack model at the anode side was derived from the model above, and models for components of the FDS was also developed that includes two supply lines and two recirculation loops, blowers, ejectors and other valves. The FDS model was then connected to the dynamic stack model that considered the liquid water effect in manifolds, GFC, GDL and CL. The two- 168 phase effect was firstly involved in the FDS modeling and the simulation. Based on the integrated model, a control strategy was proposed, which consists of a PI controller for low current mode and a new state feedback controller with integral and an observer for medium and high current mode in addition to purging controls. The SFB was compared with two classic controllers and its performance was analyzed. The major findings of the research work are summarized as follows: 1. A control-oriented one dimensional two-phase transient model of a single fuel cell including 15 layers across membrane and 53 state variables was developed that can use the total cell current as the system input, and be utilized in fuel cell stack modeling and integrated with controllers of BOP. 2. Experiment results shown that that the change of air humidity in gas flow channels (maximum voltage change between the cases of RH = 20% and RH = 100% was about 0.3V by comparing the I-V curves) had larger influence on the cell performance than the change of air flow rates did (maximum voltage change between the cases of Air/H2 = 66/28 sccm and Air/H2 = 200/84 sccm was about 0.1V by comparing the I-V curves). 3. The static and transient performance including I-V curves (voltage errors between the experiment and simulation were less than 0.1V), liquid water amount in GFCs and cell voltage response (voltage errors between the experiment and simulation were less than 0.18V) were predicted by the simulation results of the fuel cell model. 4. The FDS with the controller was likely to be unstable without purging process since the liquid water caused by condensation and transport was found in the anode side of fuel cell stack and accumulated in gas flow channels and manifolds. 169 5. Evaluations of three control strategies (SFF, PI and SFB) indicated that the SFB control with integral and observer shown the best performance with respect to dynamically stabilize the stoichiometric ratio of hydrogen (SFB settle time of SRH2 is the shortest ones that are 1.2 seconds and 0.5 second for step change of stack current in medium and high modes, respectively) and follow the cathode pressure (SFB used 2.5 second to stabilize the fuel pressure for step change of stack current in medium and high loads) under the disturbance of purging operating and step change of stack currents (SFB settle time is less than 0.3 second to stabilize the pressure and SRH2 for the typical purging operating). 6. Analysis shows that the water activities in GFC and GDL are regulated at about constant one, which means the liquid water amount was limited in a small range, for the dynamic purging strategy dependent on the time integral of stack current load change, which effectively prevented the water flooding at anode side and in the supply and return manifolds. The presented model of fuel cell systems and model-based control of the FDS can accelerate the design process and flexibly tune the system and control parameters. The SFB control with integral and observer used fewer realizable sensors than the full state feed-back controls and can be optimized based on the importance of the control variables. The simple mathematical structure of the SFB control can be easily implemented in the microcontroller of the BOP control systems to fulfill the control of hydrogen pressure and stoichiomitric ratio simultaneously. In the future, the one dimensional two-phase transient model and the anode stack model will be optimized. The FDS system modeling and control will include effects of temperature changes in manifolds and fuel cells, and automatic tuning of control parameters dependent upon the load current using advanced control algorithms. 170 References [1] DOE, " Hydrogen and fuel cells program plan (2010 Draft)," http://www1.eere.energy.gov/hydrogenandfuelcells/pdfs/program_plan2010.pdf, (2011). [2] Wang, Y., Chen, K. S., Mishler, J., Cho, S. C. and Cordobes, X., "A review of polymer electrolyte membrane fuel cells: Technology, applications, and needs on fundamental research", J Power Sources, 88, pp. 981-1007, (2011). [3] James, B. D., Kalinoski, J. A. and Baum, K. N., "Mass production cost estimation for direct H2 PEM fuel cell systems for automotive applications: 2010 update," http://www1.eere.energy.gov/hydrogenandfuelcells/pdfs/mass_production_cost_estimatio n_report.pdf, (2010). [4] Rodatz, P., Tsukada, A., Mladek, M. and Guzzella, L., "Efficiency improvements by pulsed hydrogen supply in PEM fuel cell systems", Proc. 15th triennial world congress, Barcelona, Spain, (2002). [5] Li, X. and Sabir, I., "Review of bipolar plates in PEM fuel cells: flow-field designs", J Power Sources, 30, pp. 359-371, (2005). [6] Cindrella, L., Kannan, A. M., Lin, J. F., Saminathan, K., Ho, Y., Lin, C. W. and Wertz, J., "Gas diffusion layer for proton exchange membrane fuel cells-A review", J Power Sources, 194, pp. 146-160, (2009). [7] Weber, A. Z. and Newman, J., "Modeling transport in polymer-electrolyte fuel cells", Chem Rev, 104, pp. 4679-4726, (2004). [8] Yao, K. Z., Karan, K., McAuley, K. B., Oosthuizen, P., Peppley, B. and Xie, T., "A review of mathematical models for hydrogen and direct methanol polymer electrolyte membrane fuel cells", Fuel Cells, 4(1-2), pp. 3-29, (2004). [9] Faghri, A. and Guo, Z., "Challenges and opportunities of thermal management issues related to fuel cell technology and modeling", Int J Heat Mass Transfer, 48, pp. 3891- 3920, (2005). [10] Djilali, N., "Computational modeling of polymer electrolyte membrane (PEM) fuel cells: Challenges and opportunities", Energy, 32, pp. 269-280, (2007). 171 [11] Haraldsson, K., "Evaluating PEM fuel cell system models", J Power Sources, 126, pp. 88-97, (2004). [12] Wang, C. Y., "Fundamental models for fuel cell engineering", Chem Rev, 104, pp. 4727- 4766, (2004). [13] Cheddie, D. and Munroe, N., "Review and comparison of approaches to proton exchange membrane fuel cell modeling", J Power Sources, 147, pp. 72-84, (2005). [14] B?y?ko?lu, A., "Reviewof proton exchange membrane fuel cell models", Int J Hydrogen Energ, 30, pp. 1181-1212, (2005). [15] Ma, L., Ingham, D. B., Pourkashanian, M. and Carcadea, E., "Review of the computational fluid dynamics modeling of fuel cells", J Fuel Cell Sci Tech, 2, pp. 246- 257, (2005). [16] Siegel, C., "Review of computational heat and mass transfer modeling in polymer- electrolyte-membrane (PEM) fuel cells", Energy, 33, pp. 1331-1352, (2008). [17] Berning, T. and Djilali, N., "Three-diemensional computational analysis of transport phenomena in a PEM fuel cell-a parametric study", J Power Sources, 124, pp. 440-452, (2003). [18] Chen, C. H., Jung, S. P. and Yen, S. C., "Flow distribution in the manifold of PEM fuel cell stack", J Power Sources, 173, pp. 249-263, (2007). [19] Dutta, S., Shimpalee, S. and Van Zee, J. W., "Numerical prediction of mass-exchange between cathode and anode channels in a PEM fuel cell", Int J Heat Mass Transfer, 44, pp. 2029-2042, (2001). [20] Gerteisen, D., Heilmann, T. and Ziegler, C., "Modeling the phenomena of dehydration and flooding of a polymer electrolyte membrane fuel cell", J Power Sources, 187, pp. 165-181, (2009). [21] He, W., Yi, J. S. and Van Nguyen, T., "Two-phase flow model of the cathode of PEM Fuel Cells using interdigitated flow fields", AIChE J, 46(10), pp. 2053-2064, (2000). [22] Le, A. D. and Zhou, B., "A numerical investigation on multi-phase transport phenomena in a proton exchange membrane fuel cell stack", J Power Sources, 195(16), pp. 5278- 5291, (2010). [23] Maharudrayya, S., Jayanti, S. and Deshpande, A. P., "Pressure losses in laminar flow through serpentine channels in fuel cell stacks", J Power Sources, 138, pp. 1-13, (2004). [24] Maher, A. R., Sadiq, A. B., Haroun, A. K. and Shahad, A. J., "Modeling optimizes PEM fuel cell performance using three-dimensional multi-phase computational fluid dynamics model", Energ Convers Manage, 48, pp. 3102-3119, (2007). 172 [25] Promislow, K., Chang, P., Haas, H. and Wetton, B., "Two-phase unit cell model for slow transients in polymer electrolyte membrane fuel cells", J Electrochem Soc, 177(7), pp. A494-504, (2008). [26] Sivertsen, B. R. and Djilali, N., "CFD-based modelling of proton exchange membrane fuel cells", J Power Sources, 141, pp. 65-78, (2005). [27] Wang, Y., Basu, S. and Wang, C. Y., "Modeling two-phase flow in PEM fuel cell channels", J Power Sources, 179, pp. 603-617, (2008). [28] Wu, H., Berg, P. and Li, X., "Non-isothermal transient modeling of water transport in PEM fuel cells", J Power Sources, 165, pp. 232-243, (2007). [29] Wu, H., Li, X. and Berg, P., "On the modeling of water transport in polymer electrolyte membrane fuel cells", Electrochim Acta, 54, pp. 6913-6927, (2009). [30] Pei, P., Ouyang, M., Feng, W., Lu, L., Huang, H. and Zhang, J., "Hydrogen pressure drop characteristics in a fuel cell stack", Int J Hydrogen Energ, 31, pp. 371-377, (2006). [31] Chang, P. A. C., St-Pierre, J., Stumper, J. and Wetton, B., "Flow distribution in proton exchange membrane fuel cell stacks", J Power Sources, 162, pp. 340-355, (2006). [32] Baschuk, J. J. and Li, X., "Modelling of polymer electrolyte membrane fuel cell stacks based on a hydraulic network approach", Int J Energ Res, 28, pp. 697-724, (2004). [33] Yi, J. S., Yang, D. and King, C., "Water management along the flow channels of PEM fuel cells", AIChE J, 50(10), pp. 2594-2603, (2004). [34] Rodatz, P., Buchi, F., Onder, C. and Guzzella, L., "Operational aspects of a large PEFC stack under practical conditions", J Power Sources, 128, pp. 208-217, (2004). [35] Argyropoulos, P., Scott, K. and Taama, W. M., "Pressure drop modeling for liquid feed direct methonal fuel cells Part 1. Model development", Chem Eng J, 73, pp. 217-227, (1999). [36] Pukrushpan, J. T., Peng, H. and Stefanopoulou, A. G., "Control-oriented modeling and analysis for automotive fuel cell systems", J Dyn Syst-T ASME, 126, pp. 14-25, (2004). [37] Bao, C., Ouyang, M. and Yi, B., "Modeling and control of air stream and hydrogen flow with recirculation in a PEM fuel cell system II: Linear and adaptive nonlinear control", Int J Hydrogen Energ, 31, pp. 1897-1913, (2006). [38] Karnik, A. Y., Sun, J. and Buckland, J. H., "Control analysis of an ejector based fuel cell anode recirculation system", Proc. 2006 American Control Conference, Minneapolis, Minnesota, USA, (2006). 173 [39] He, J., Choe, S.-Y. and Hong, C.-O., "Analysis and control of a hybrid fuel delivery system for a polymer electrolyte membrane fuel cell", J Power Sources, 185, pp. 973- 984, (2008). [40] Baschuk, J. J. and Li, X., "Modeling of polymer electrolyte membrane fuel cells with variable degrees of water flooding", J. Power Sources, 86, pp. 181-196, (2000). [41] Weber, A. Z., "Gas-crossover and membrane-pinhole effects in polymer-electrolyte fuel cells", J Electrochem Soc, 155(6), pp. B521-B531, (2008). [42] Springer, T. E., Zawodzinski, T. A. and Gottesfeld, S., "Polymer electrolyte fuel cell model", J Electrochem Soc, 138, pp. 2334-2342, (1991). [43] Carey, V. P., Liquid-Vapor Phase-Change Phenomena, 2nd ed: Taylor & Francis, (2008). [44] Kays, W., Crawford, M. and Weigand, B., Convective Heat and Mass Transfer, 4th ed, (2005). [45] Zong, Y., Zhou, B. and Sobiesiak, A., "Water and thermal management in a single PEM fuel cell with non-uniform stack temperature", J Power Sources, 161, pp. 143-159, (2006). [46] Wilke, C. R., "A viscosity equation for gas mixtures", J Chem Phys, 18(4), pp. 517-519, (1950). [47] White, F. M., Viscous Fluid Flow, 3rd ed: Mcgraw-Hill, (2006). [48] Shah, A. A., Ralph, T. R. and Walsh, F. C., "Modeling and simulation of the degradation of perfluorinated ion-exchange membrane in PEM fuel cells", J Electrochem Soc, 156(4), pp. B465-B484, (2009). [49] Hinatsu, J. T., Mizuhata, M. and Takenaka, H., "Water uptake of perfluorosulfonic acid membranes from liquid water and water vapor", J Electrochem Soc, 141(6), pp. 1493- 1498, (1994). [50] Ge, S., Li, X., Yi, B. and Hsing, I.-M., "Absorption, desorption, and transport of water in polymer electrolyte membranes for fuel cells", J Electrochem Soc, 156(2), pp. A1149- A1157, (2005). [51] Weber, A. Z. and Newman, J., "Transport in polymer-electrolyte membranes II. mathematical model", J Electrochem Soc, 151(2), pp. A331-A325, (2004). [52] Motupally, S., Becker, A. J. and Weidner, J. W., "Diffusion of water in Nafion 115 membranes", J Electrochem Soc, 147(9), pp. 3171-3177, (2000). [53] Mench, M. M., Wang, C. Y. and Ishikawa, M., "In situ current distribution measurements in polymer electrolyte fuel cells", J Electrochem Soc, 150(8), pp. A1052-A1059, (2003). 174 [54] Yoon, Y.-G., Lee, W.-Y., Yang, T.-H., Park, G.-G. and Kim, C.-S., "Current distribution in a single cell of PEMFC", J Power Sources, 118, pp. 193-199, (2003). [55] Strickland, D. G., Litster, S. and Santiago, J. G., "Current distribution in polymer electrolyte membrane fuel cell with active water management", J Power Sources, 174, pp. 272-281, (2007). [56] P?rez, L. C., Brand?o, L., Sousa, J. M. and Mendes, A., "Segmented polymer electrolyte membrane fuel cells-A review", Renew Sust Energ Rev, 15, pp. 169-185, (2011). [57] T?ber, K., P?cza, D. and Hebling, C., "Visualization of water buildup in the cathode of a transparent PEM fuel cell", J Power Sources, 124, pp. 403-414, (2003). [58] Spernjak, D., Prasad, A. K. and Advani, S. G., "Experimental investigation of liquid water formation and transport in a transparent single-serpentine PEM fuel cell", J Power Sources, 170, pp. 334-344, (2007). [59] Ous, T. and Arcoumanis, C., "Visualisation of water droplets during the operation of PEM fuel cells", J Power Sources, 173, pp. 137-148, (2007). [60] Hussaini, I. S. and Wang, C.-Y., "Visualization and quantification of cathode channel flooding in PEM fuel cells", J Power Sources, 187, pp. 444-451, (2009). [61] Satijaa, R., Jacobson, D. L., Arifb, M. and Werne, S. A., "In situ neutron imaging technique for evaluation of water management systems in operating PEM fuel cells", J Power Sources, 129, pp. 238-245, (2004). [62] Trabolda, T. A., Owejana, J. P., Jacobson, D. L., Arifb, M. and Huffmanc, P. R., "In situ investigation of water transport in an operating PEM fuel cell using neutron radiography: Part 1-Experimental method and serpentine flow field results", Int J Heat Mass Transfer, 49, pp. 4712-4720, (2006). [63] Hartniga, C., Mankeb, I., Kardjilovb, N., Hilgerb, A., Gr?nerbela, M., Kaczerowskia, J., Banhartb, J. and Lehnerta, W., "Combined neutron radiography and locally resolved current density measurements of operating PEM fuel cells", J Power Sources, 176, pp. 452-459, (2008). [64] Hickner, M. A., Siegel, N. P., Chen, K. S., Hussey, D. S., Jacobson, D. L. and Arif, M., "In situ high-resolution neutron radiography of cross-sectional liquid water profiles in proton exchange membrane fuel cells", J Electrochem Soc, 155(4), pp. B427-B434, (2008). [65] T?tzke, C., Manke, I., Hilger, A., Choinka, G., Kardjilov, N., Arlt, T., Mark?tter, H., Schr?der, A., Wippermann, K., Stolten, D., Hartnig, C., Kr?ger, P., Kuhn, R. and Banhart, J., "Large area high resolution neutron imaging detector for fuel cell research", J Power Sources, doi:10.1016/j.jpowsour.2011.01.049, (2011). 175 [66] Sinha, P. K., Halleck, P. and Wang1, C.-Y., "Quanti?cation of liquid water saturation in a PEM fuel cell diffusion medium using X-ray microtomography", Electrochem Solid St, 9(7), pp. A344-A348, (2006). [67] Manke, I., Hartnig, C., Gr?nerbel, M., Lehnert, W., Kardjilov, N., Haibel, A., Hilger, A., Banhart, J. and Riesemeier, H., "Investigation of water evolution and transport in fuel cells with high resolution synchrotron x-ray radiography", Appl Phys Lett, 90(17), pp. 174501-174505, (2007). [68] Manke, I., Hartnig, C., Kardjilov, N., Riesemeier, H., Goebbels, J., Kuhn, R., Kr?ger, P. and Banhart, J., "In situ synchrotron X-ray radiography investigations of water transport in PEM fuel cells", Fuel Cells, 10(1), pp. 26-34, (2009). [69] Kr?ger, P., Mark?tter, H., Hau?mann, J., Klages, M., Arlt, T., Banhart, J., Hartnig, C., Manke, I. and Scholta, J., "Synchrotron X-ray tomography for investigations of water distribution in polymer electrolyte membrane fuel cells", J Power Sources, doi:10.1016/j.jpowsour.2010.09.042, (2010). [70] Moraal, P. and Kolmanovsky, I., "Turbocharger modeling for automotive control applications", International Congress and Exposition, Detroit, Michigan, USA, (1999). [71] Pukrushpan, J. T., Modeling and control of fuel cell systems and fuel processors, PhD Dissertation, Mechanical Engineering, The University of Michigan, (2003). [72] Huang, B. J., Jiang, C. B. and Hu, F. L., "Ejector performance characteristics and design analysis of jet refrigeration system", J Eng Gas Turb Power, 107(3), pp. 792-802, (1984). [73] Munday, J. T. and Bagster, D. F., "A new ejector theory applied to steam jet refrigeration", Ind Eng Chem Proc DD, 16(4), pp. 442-449, (1977). [74] Aphornratana, S. and Eames, I. W., "A small capacity steam-ejector refrigerator: experimental investigation of a system using ejector with movable primary nozzle", Int J Refrig, 20, pp. 352-358, (1997). [75] Rogdakis, E. D. and Alexis, G. K., "Investigation of ejector design at optimum operating condition", Energ Convers Manage, 41, pp. 1841-1849, (2000). [76] Marsano, F., Magistri, L. and Massardo, A. F., "Ejector performance influence on a solid oxide fuel cell anodic recirculation system", J Power Sources, 129, pp. 216-228, (2004). [77] Zhu, Y. and Li, Y., "New theoretical model for convergent nozzle ejector in the proton exchange membrane fuel cell system", J Power Sources, 191, pp. 510-519, (2009). [78] Ferrari, M. L., Traverso, A., Magistri, L. and Massardo, A. F., "Influence of the anodic recirculation transient behaviour on the SOFC hybrid system performance", J Power Sources, 149, pp. 22-32, (2005). 176 [79] Kim, Y.-B., "Improving dynamic performance of proton-exchange membrane fuel cell system using time delay control", J Power Sources, 195, pp. 6329-6341, (2010). [80] Bartosiewicz, Y., Aidoun, Z., Desevaux, P. and Mercadier, Y., "Numerical and experimental investigations on supersonic ejectors", Int J Heat Fluid Flow, 26, pp. 56-70, (2005). [81] He, S., Li, Y. and Wang, R. Z., "Progress of mathematical modeling on ejectors", Renew Sust Energ Rev, 13, pp. 1760-1780, (2009). [82] Huang, B. J., Chang, J. M., Wang, C. P. and Petrenko, V. A., "A 1-D analysis of ejector performance", Int J Refrig, 22, pp. 354-364, (1999). [83] Zhu, Y., Cai, W., Wen, C. and Li, Y., "Shock circle model for ejector performance evaluation", Energ Convers Manage, 48, pp. 2533-2541, (2007). [84] Pianthong, K., Seehanam, W., Behnia, M., Sriveerakul, T. and Aphornratanac, S., "Investigation and improvement of ejector refrigeration system using computational fluid dynamics technique", Energ Convers Manage, 48, pp. 2556-2564, (2007). [85] Rusly, E., Aye, L., Charters, W. W. S. and Ooib, A., "CFD analysis of ejector in a combined ejector cooling system", Int J Refrig, 28, pp. 1092-1101, (2005). [86] Sriveerakul, T., Aphornratana, S. and Chunnanond, K., "Performance prediction of steam ejector using computational fluid dynamics: Part 1. Validation of the CFD results", Int J Therm Sci, 46, pp. 812-822, (2007). [87] Mokmeli, A. and Asghari, S., "An investigation into the effect of anode purging on the fuel cell performance", Int J Hydrogen Energ, 35, pp. 9276-9282, (2010). 177 Appendix A Parameters of the transient model of the fuel cell Symbol Name Value cw,l Specific heat capacity of liquid water 4.1813?103 J kg-1 K-1 cCP Specific heat capacity of CP 200 J kg-1 K-1 cGDL Specific heat capacity of GDL 150 J kg-1 K-1 cCL Specific heat capacity of CL 160 J kg-1 K-1 cPEM Specific heat capacity of membrane 160 J kg-1 K-1 dpore,GDL Average pore diameter in GDL 1.2?10-5 m dpore,CL Average pore diameter in CL 5?10-7 m kads Adsorption rate 80 s-1 kdes Desorption rate 50 s-1 kcon Condensation rate 100 s-1 kevp Evaporation rate 9.8717?10-4 s-1 Pa-1 kH2 Thermal conductivity of hydrogen 0.2040 W m-1 K-1 kO2 Thermal conductivity of oxygen 0.0296 W m-1 K-1 kv Thermal conductivity of vapor 0.0237 W m-1 K-1 kN2 Thermal conductivity of nitrogen 0.0293 W m-1 K-1 kin,ca,GFC Inlet flow coefficient of cathode GFC 2.84?10-2 m s-1 Pa-1 178 kin,an,GFC Inlet flow coefficient of anode GFC 3.82?10-2 m s-1 Pa-1 kl Thermal conductivity of liquid water 0.66 W m-1 K-1 kr Water release rate 4?106 mol m3 s-1 Pa-1 ku Water uptake rate 4?106 mol m3 s-1 Pa-1 kw,con Wall condensation rate 1?10-3 s-1 kw,evp Wall evaporation rate 1?10-8 s-1 Pa-1 kCP Thermal conductivity of CP 25 W m-1 K-1 kGDL Thermal conductivity of GDL 1.67 W m-1 K-1 kCL Thermal conductivity of CL. 0.67 W m-1 K-1 kPEM Thermal conductivity of membrane 0.67 W m-1 K-1 ksp,ca,GFC Flow coefficient of gas phase in GFC at cathode side 0.0656 ksp,an,GFC Flow coefficient of gas phase in GFC at anode side 0.0309 ns1 Power parameter in Eq. (2-24) 1.5 (assumed) ns2 Power parameter in Eq. (2-33) 5.5 (assumed) ns3 Power parameter in Eq. (2-43) 2/3 (assumed) pc,crit Critical capillary pressure in GFC 0 pa sim,GDL Immobilized liquid saturation in GDL 0.2 sim,CL Immobilized liquid saturation in CL 0.01 Aact Active area of fuel cell 140 c m2 Cp,H2 H2 specific heat capacity at constant pressure 28.82 J mol-1 K-1 Cp,O2 O2 specific heat capacity at constant pressure 29.38 J mol-1K-1 Cp,v Vapor specific heat capacity at constant pressure 36.93 J mol-1K-1 179 Cp,N2 N2 specific heat capacity at constant pressure 29.14 J mol-1K-1 Ctp,GFC Two phase correction factor for liquid flow rate 0.2 DH2-H2O Binary diffusion constant of H2-H2O (307.1K, 1atm) 9.15?10-5 m2 s-1 DO2-H2O Binary diffusion constant of O2-H2O (293.15K, 1atm) 2.40?10-5 m2 s-1 DO2-N2 Binary diffusion constant of O2-N2 (293.15K, 1atm) 2.19?10-5 m2 s-1 DN2-H2O Binary diffusion constant of N2-H2O (307.5K, 1atm) 2.56?10-5 m2 s-1 EWPEM Equivalent weight of membrane 1.1 kg HGFC Height of GFC 0.6 mm KGDL Permeability of GDL 8.76?10-12 m2 KCL Permeability of CL 4.15?10-13 m2 Kl,cov,GFC Liquid water cover coefficient 2.08 MH2 H2 molecular weight 2.016?10-3 kg mol-1 MO2 O2 molecular weight 32?10-3 kg mol-1 MH2O H2O molecular weight 18.015?10-3 kg mol-1 MN2 N2 molecular weight 28.013?10-3 kg mol-1 NGFC Number of GFCs 5 NuGFC Nusselt number of GFC 3.21 ShGFC Sherwood number of GFC 1.98 WGFC Width of GFC 1 mm Wrib,GFC Rib width between GFCs 1 mm ?l Liquid water density 0.986?103 kg m-3 ?CP density of CP 3000 kg m-3 180 ?GDL GDL density 1500 kg m-3 ?CL CL density 1800 kg m-3 ?PEM membrane density 2000 kg m-3 ?w Water surface tension 0.0625 N m-1 ?CP Thickness of CP 3 mm ?GDL Thickness of GDL 240 ?m ?CL Thickness of CL 10 ?m ?PEM Thickness of membrane 50 ?m ?GDL Porosity of GDL 0.8 ?CL Porosity of CL 0.3 ?i,CL Volume fraction of the ionomers in CL 0.3 ?c,GDL Contact angle in GDL 110? ?c,CL Contact angle in CL 95? ?c,i Contact angle in the ionomer channels 90.02? ?San Enthalpy change at anode reaction -0.104 J mol-1 K-1 ?Sca Enthalpy change at cathode reaction 326.36 J mol-1 K-1