Quasi-three dimensional dynamic model of a proton exchange membrane fuel cell for system and controls development

A ﬁrst principles dynamic model of the physical, chemical, and electrochemical processes at work in a proton exchange membrane fuel cell has been developed. The model solves the dynamic equations that govern the physics, chemistry and electrochemistry for time scales greater than about 10ms. The dynamic equations are solved for a typical but simpliﬁed quasi-three dimensional geometric representation of a single cell repeat unit of a fuel cell stack. The current approach captures spatial and temporal variations in the important physics of heat transfer and water transport in a manner that is simple enough to make the model amenable to PEMFC system simulations and controls development. Comparisons of model results to experimental data indicate that the model can well predict steady state voltage–current relationships as well as the oxygen, water, and nitrogen spatial distribution within the fuel cell. In addition, the model gives dynamic insight into the distribution of current, water ﬂux, species mole fractions, and temperatures within the fuel cell. Finally, a control system test is demonstrated using the simpliﬁed dynamic model.


Background and introduction
Due to quiet operation, high efficiency and low pollutant emissions performance characteristics, proton exchange membrane fuel cells (PEMFC) are receiving increased attention for applications in the automobile, portable power and backup power industries. PEMFC utilize air and hydrogen to generate electrical power that can be used to meet a wide variety of electrical load demands. Use of PEMFC in automotive applications especially, and in many other applications, will require a significant capability to follow dynamic load changes. In addition, to ensure safe and efficient operation and to develop control strategies for PEMFC systems it is important to understand both the steady state and dynamic performance characteristics of the PEMFC.
Modeling has been used extensively to garner insight into PEMFC performance. Most of the PEMFC modeling efforts in the literature address steady state performance through multidimensional simulation. More recently, reduced order (bulk * Corresponding author. Tel.: +1 949 338 3953; fax: +1 949 824 7423. E-mail address: jb@nfcrc.uci.edu (J. Brouwer). or one-dimensional) transient fuel cell modeling and multidimensional fuel cell component transient modeling efforts have been presented in the literature. The current dynamic model is novel in that it attempts to resolve some of the geometric features of a PEMFC in a simplified manner sufficient to make it amenable to use in PEMFC system simulation and controls development. The geometric simplification employed in the current work is a 2 + 1D reduction (quasi-3D) with one dimension (1D) resolved through the membrane electrode assembly and two dimensions (2D) resolved in the fuel cell repeat unit crosssection accounting for gas and thermal transport in lumped serpentine channels.
An excellent review of steady state and multi-dimensional PEMFC modeling has been developed presented by Weber and Newman [1]. Weber and Newman characterize these models by affiliation as being derived from the significant early contributions of Springer et al. [2], Bernardi and Verbrugge [3,4], and those using computational fluid dynamics or other approaches. Since publication of the Weber and Newman paper, steady state PEMFC modeling has continued to be an active area of research as characterized by recent contributions from many groups including Lum and McGuirk [5], and Shimpalee et al. 0378 [6], Pasaogullari and Wang [7], Guvelioglu and Stenger [8], Gorgun et al. [9], Nguyen et al. [10], and Freunberger et al. [11], among others. In addition, Burt et al. [12] and Campanari and Iora [13] have developed similar multi-dimensional models for solid oxide fuel cells. While many additional groups are making significant contributions to steady state fuel cell models, discussion of these contributions is limited herein since the current paper addresses transient PEMFC modeling. Amphlett et al. [14,15], Wohr et al. [16], and van Bussel et al. [17] were some of the first researchers to develop and use dynamic PEMFC models in the late 1990s. These dynamic models were generally bulk models that accounted for the dynamic physical, chemical and electrochemical processes in each single cell component (e.g., gas diffusion layer, membrane). These transient models were able to predict fuel cell voltage and heat losses as a function of time due to various changes imposed on the PEMFC. Subsequent reduced order (bulk or one-dimensional) dynamic PEMFC models have been developed by many groups. The bulk dynamic model of Yuyao and Choe [18] elucidated the mechanisms of PEMFC dehydration. Yerramalla et al. [19] developed one of the first bulk dynamic PEMFC models for control system development in Simulink ® . Ceraolo et al. [20] and Pathapati et al. [21] each developed dynamic bulk PEMFC models accounting for the primary processes that govern membrane electrode assembly behavior in Simulink ® . Pukrushpan et al. [22] were the first to demonstrate control system development based on the dynamic understanding of water transport in a bulk PEMFC model. Xue et al. [23] developed a dynamic PEMFC model that could predict the effects of temperature, gas flow, and capacitance on system transient behavior. Lemes et al. [24] used a bulk dynamic PEMFC model in a hardware-in-the-loop testing effort. In addi-tion, several groups have developed bulk dynamic models for the related direct methanol fuel cell, including Sundmacher et al. [25], Kulikovsky [26] and Xu et al. [27]. While each of these models makes a significant contribution, none resolve the planar distribution of species, temperature, membrane water content, etc. within the PEMFC and some use simplified approaches for modeling the dynamics of water transport.
More recently, multi-dimensional dynamic modeling of PEM fuel cells has been explored. The focus of these multidimensional dynamic models has been on resolving internal cell transient behavior. For example, Um and Wang [28] have investigated the interaction between mass transport and electrochemical kinetics at the cell level. The more recent contribution of Um et al. [29] presents a transient, multi-dimensional model that accounts simultaneously for electrochemical kinetics, current distribution, hydrodynamics, and multi-component transport. Chen et al. [30] investigated the effect of water transport dynamics on cell transient capability, and Yan et al. [31] have investigated the effects of flow distribution and gas diffusion layer (GDL) morphology on the transient performance during startup operation. Berning and Djilali [32] presented a computational fluid dynamics multiphase model of a PEM fuel cell that accounted for three-dimensional (3D) transport processes with phase change and heat transfer in the gas diffusion layers, the gas flow channels, and the cooling channels. Yu et al. [33] developed a dynamic water and thermal management model for a Ballard PEMFC stack. Shimpalee et al. [34,35] present a threedimensional numerical simulation of the transient response of a PEM fuel cell (PEMFC) subjected to a variable load. Transient responses of the cell are predicted based on local distributions of the current density and gas mole fractions. The predictions show transients in the current density that overshoot the final state when the cell voltage is abruptly changed from 0.7 to 0.5 V for fixed initial flow rates in excess of stoichiometric to normal [34]. Similar transient predictions are made for normal to minimal flow conditions in [35]. Wang and Wang [36] also developed a three-dimensional transient PEMFC model, which when applied to a single channel PEMFC showed that the time for fuel cells to reach steady state was on the order of 10 s due to the effect of water accumulation in the membrane. Wang and Wang observed overshoot and undershoot of the current densities during step changes of some operating conditions. These multi-dimensional transient investigations have typically utilized finite difference methods with detailed resolution of geometrical features in multiple cells and are quite computational intensive. In addition, many are developed using commercially available computational fluid dynamics packages that make them not amenable or suitable for system level transient study and controls development.
The objective of the current work is to develop a model that can capture both the dynamic response characteristics and some effects of cell geometry in a PEMFC model that is simplified enough to be useful in system level simulation and control system development. The present model is a balance between the high geometric resolution of prior steady state and dynamic cell and stack models and the low to zero geometric resolution of previous reduced order dynamic models. The approach simplifies the geometry by only considering a quasi-three dimensional representation of a single fuel cell repeat unit of the fuel cell stack. In addition, only the physics, chemistry and electrochemistry extant in the fuel cell that contribute to dynamic performance of a PEM fuel cell at time scales greater than about 10 ms are resolved. All other processes are assumed to be in equilibrium or quasi-steady state. Finally, the model is developed in the Simulink ® framework to allow development and testing of control strategies and techniques.
At the timescales of interest to the current effort (i.e., greater than 10 ms) the dynamic voltage response of a PEMFC depends primarily upon: (1) the amount of current drawn from the fuel cell, (2) localized species mole fractions, (3) local temperatures, (4) pressures at local triple phase boundaries, and (5) water content of the membrane. Therefore, in order to garner detailed insights into dynamic fuel cell performance characteristics, it is necessary to understand and resolve the local behavior of species mole fractions, temperatures, pressures, and water transport within the fuel cell during transient operation. In addition, it is desirable to develop such a detailed dynamic model in a framework that allows the development and testing of control strategies for safety and enhanced performance.
By discretizing the PEMFC in two dimensions and resolving heat and mass transfer and electrochemical processes in the third dimension (quasi-three dimensional) it is possible to predict the dynamic distribution of current, temperature, species mole fractions, and membrane water content in the fuel cell dynamically. The local resolution provides insight into the local performance of the fuel cell, local thermal stresses, local hydrogen and oxygen starvation, as well as potential flooding conditions and locations during transient operation. In addition, the local resolution enables more accurate prediction of overall fuel cell dynamic performance compared to bulk dynamic models, which is beneficial for control and system design. To evaluate the model and verify model accuracy, experimental fuel cell data were acquired and compared to results from the current model. Once verified, the modeling technique as described herein can be applied to other PEMFC geometries and used to evaluate system configuration and control strategies. To demonstrate the utility of the model for control design, a robust power and current-based fuel control is developed and demonstrated using the model.

Model development
The current model was developed to dynamically resolve local states and operating features of the fuel cell, but yet be simple enough for simulation on a personal computer (∼2 GHZ) and incorporation into systems level models for control system development and testing. Since the model is to be utilized for controls development it was developed in Matlab-Simulink ® . Simulink ® is a preferred framework for controls development, but imposes a major limitation to the resolution that can be captured in the model. As a result, neither all geometric features of the particular PEM fuel cell stack nor all of the physical and chemical processes occurring in the fuel cell stack can be captured without the model becoming too computationally intensive. Thus, many simplifying assumptions and geometric simplifications have been made to enable solution of the dynamic equations that govern the PEMFC.
Only processes that affect the larger timescales (>10 ms) associated with PEMFC performance are resolved dynamically. An order of magnitude analysis was accomplished in the current effort to determine whether processes must be dynamically resolved or can be assumed to be in equilibrium for the timescales of interest (>10 ms). The results of this order of magnitude analysis are presented in Table 1. Note that the complete kinetic analysis of electrochemical reactions would require the determination of the rate constants for each elementary reaction step. Hence a governing equation for electrochemical reaction rate is not presented in Table 1. However, the electrochemical time scale is on the order of 10 −3 s as reported by [36]. From the order of magnitude analysis, it has been determined that dynamics of electrochemical reactions, and charging or discharging of the charge double layer can be ignored since only time scales greater than 10 ms are of interest. On the other hand, species transport and diffusion, water accumulation in the membrane and heat transfer effects should be captured. Heat effects should be captured when investigating transients on the order of tens of seconds since heat generation in the membrane can have time scales in this range.
The current model was developed in the Simulink ® framework using the numerical method of lines as the basic solution technique. Each of the components that comprise the repeat unit of a PEMFC is discretized in space using control volumes. In this fashion the mechanisms of system energy and species conservation can be described in terms of ordinary differential equations that can be solved in Simulink to determine local temperatures and species mole fractions in each control volume. From local temperatures and species mole fractions of each control volume, the reaction rate, molar capacity, heat transfer, and thermodynamic properties can be evaluated throughout the quasi-three dimensional cell. In addition, by assuming quasi-steady electrochemistry (the kinetics of which only affects dynamic performance at timescales less than about 10 ms), the fuel cell voltage and current generation can be determined. Of course the losses associated with the electrochemical kinetics are accounted for through a bulk activation polarization term. In this fashion the model predicts local temperatures, species mole fractions, and flows as well as the extent of heat transfer and local electrical work generated throughout the cell.

Assumptions
The following is a list of the major assumptions made: 1. Control volumes are characterized by a single lumped temperature, pressure, and set of species mole fractions condition. 2. All gases are ideal gases [11]. 3. A uniform gas pressure is assumed. The pressure drop along the gas flow channels and GDL are neglected [11]. 4. One-dimensional fully developed laminar flow occurs along the stream-wise direction [11]. 5. Parallel diffusive fluxes in the gas diffusion layer and membrane are ignored. Convective transport inside the flow channel by far dominates parallel diffusive fluxes (Pe = 10 10 ) [11]. 6. The solid gas diffusion layer (GDL) and membrane electrode assembly (MEA) have a lumped temperature. (The respective Biot number was found to be much less than 0.1.) 7. Each cell in the stack is assumed to operate identically, so that a single cell simulation is taken as representative and used to calculate full stack performance [11,12]. 8. All electrodes are good conductors for which an equipotential electrode surface is assumed [13]. 9. Quasi-steady electrochemistry is assumed, since the electrochemistry is rapid (occurring at time scales on the order of 10 −3 s) [37]. 10. A single activation polarization equation is used to capture the effects of all physical and chemical processes that polarize the charge transfer process. 11. All reactants generate their ideal number of electrons, and no fuel or oxidant crosses the electrolyte [11].

Discretization
To solve for local states of the fuel cell, the fuel cell is discretized into control volumes quasi-three dimensionally. The fuel cell is discretized perpendicular to the flow through one repeat unit of the fuel cell stack as well as two dimensionally in the flow plane to capture the superficial area of the cell and the serpentine flow pattern as well. Flow perpendicular discretization is accomplished in a fashion similar to Yuyao and Choe [18] and Freunberger et al. [11]. The primary components of the fuel cell (bipolar plate, anode flow, anode GDL, MEA, cathode GDL, cathode flow, cathode plate and coolant flow) are discretized into perpendicular control volumes ( Fig. 1) for each of the nodes used to discretize the fuel cell in the flow plane. In this fashion the fuel cell is discretized in the vertical direction into eight control volumes, using five types of control volumes: (1) solid plate, (2) bulk gas, (3) GDL, (4) MEA, and (5) coolant.
Each of the eight control volumes in the flow perpendicular direction, were chosen because the operating voltage and local current production of a fuel cell is governed by the localized species mole fractions, temperature, and pressure at the triple phase boundary and the extent of membrane hydration, which can only be resolved with understanding of local heat transfer, electrochemistry, water transport, etc. in this perpendicular direction. In order to resolve the local species mole fractions at the triple phase boundary, the anode and cathode gas and GDL control volumes were desired. The gas control volume represents the bulk flow of the gases. The gas from the bulk flow then diffuses through the GDL to the triple phase boundary, between the GDL and the MEA. Diffusion fluxes and osmotic water transport to the GDL control volume is determined for the boundary between the electrolyte and GDL control volume. In this fashion the species calculated to diffuse through the GDL control volume represent the species mole fractions at the triple phase boundary of the fuel cell.
It is possible to model a PEMFC without resolving the GDL or the triple phase boundary by assuming rapid diffusion between the triple phase boundary and the bulk flow. In this case, the concentration gradient between the bulk flow and the triple phase  (5) MEA Water content Eqs (11) and (12)  Coolant Temperature Eq. (2) boundary as it impacts cell voltage can be somewhat accounted for using a bulk concentration polarization term [38]. Use of this assumption would make it possible to simplify the model by removing the GDL control volume. However, removing the GDL control volume makes it quite difficult to accurately model water transport and membrane hydration, which significantly affects conductivity and overall cell performance. Research has indicated that there can be significant water gradients between the bulk gas phase, triple phase boundaries and the electrolyte due to water formation at the cathode triple phase boundary, and osmotic drag [7]. As a result, to more accurately model the extent of membrane hydration it is beneficial to discretize the fuel cell utilizing both GDL and MEA control volumes. This makes it possible to predict water movement in the membrane, accounting for water formation, osmotic drag, diffusion, and water storage within the fuel cell. The solid plate, bulk gases, coolant and a lumped GDL and MEA control volume are needed to resolve the perpendicular temperature profile of the fuel cell (see Fig. 3). While a separate GDL and MEA control volume is essential to resolve the water transport and concentration gradients within the fuel cell, the GDL and MEA can be lumped together into one perpendicular control volume to resolve the perpendicular temperature profile. This is justified because the applicable Biot number of the GDL and MEA is much less than 0.1. The perpendicular discretization is summarized in Table 2.
One set of eight control volumes in the perpendicular direction represents a single node of the cross-sectional area of the fuel cell. The fuel cell cross-sectional area can then be resolved by solution of local performance characteristics in each node as it represents a discrete part of the fuel cell serpentine flow path. Individual nodes are therefore connected to one another in the primary flow direction in a two dimensional extension of the approach introduced by Yi and Nguyen [39] and used also by Freunberger et al. [11]. Appropriate accounting for heat and mass transport between and amongst adjacent nodes in the flow direction and across flow channels is made. The flow parameters (flow rate, species mole fractions, and temperature) in the bulk gas and coolant control volumes are then passed from node to node in the stream-wise direction. Since the flow paths through the typical PEMFC are not straight (e.g., serpentine), the current discretization results in two-dimensional resolution of the flow field and fuel cell performance, and three-dimensional resolution of the temperature distribution within the fuel cell repeat unit. Fig. 1 demonstrates the PEMFC discretization in the two directions for a two by two node representation of a fuel cell that contains a total of 32 control volumes.
The particular experimental fuel cell that we desire to simulate contains seven small channels for oxidant and fuel transport, each of which deliver gases and remove products in a seven pass serpentine fashion as shown in Fig. 4. To simulate this particular geometric configuration of the fuel cell, the seven individual small channels are lumped together, into single flow channels in each of the nodes used in the model. In addition, the serpentine path is captured by the use of five discrete nodes along the width and seven discrete nodes along the length of fuel cell cross-section (Fig. 4).
Thus the current cell model contains a five by seven nodal representation of the flow field and fuel cell performance characteristics for a total of 35 nodes, and 280 control volumes. The discretized model geometry is presented diagrammatically in Fig. 4. Cross-channel flow effects can be important in PEMFC operation. In the current model, however, no cross-channel flow is allowed between the serpentine paths. Note that the current model flow paths already lump several actual flow channels of the experiment into one larger channel. In effect, this accounts for some of the cross-channel flow effects amongst the channels that would otherwise have comprised the actual flow paths of the experiment, while resolving the primary overall direction of flow.
The model was discretized in a manner similar to Freunberger et al. [11] with eight control volumes in the flow perpendicular direction and many nodes along the flow direction. The current model could benefit from a more detailed discretization in the flow direction, except for the goal of assuring amenability to system simulation and controls development. Detailed CFD models for transient or steady state analysis typically contain on the order of 10,000 states, while Freunberger et al. discretized their model into approximately 2000 states. Using 2000 states in the current model for system control design is unrealistic as this will increase the computational burden on Simulink by a factor of 10 compared to the present model, resulting in convergence difficulty, and long simulation time. Hence a simplified discretization is used to retain practicality as a control and system design tool.

Conservation equations
Each control volume temperature, species mole fractions, molar flow rate, and water content is solved from conservation first principles. Energy conservation and heat transfer equations are used to solve for temperatures, species conservation and mass transport equations are used to solve for species mole fractions, molar flow rates and water content. The following describes the set of dynamic conservation equations that is solved in each control volume. Subsequent sections describe the flux of species, reactions, and heat transfer between and amongst control volumes as needed to solve each of the conservation equations.

Energy conservation
Temperatures throughout the model are determined from solution of the dynamic energy conservation equation. The temperature of each solid plate control volume is found from solving the following ordinary differential equation (ODE): (1) whereQ in represents heat transfer to the control volume. In a similar fashion each of the bulk gas and coolant control volume temperatures is determined by solution of the dynamic energy conservation equation: whereṄ in h in represent the enthalpy flux into the control volume, N out h out the enthalpy flux out of the control volume,Q in the heat transfer to the control volume, andQ latent represent the latent heat of liquid water evaporating into the bulk gas or vice versa. The temperature of the control volume is taken to be the exit temperature of the control volume, and the inlet temperature of the control volume is that of the upstream (from the "flow" perspective) node exit temperature. Note that the "upstream" node can be different for each of the solutions for anode gas, cathode gas, and coolant streams since flows may traverse through the cell in different manners. This feature must be carefully accounted for in the model. The molar flow rate of each species is determined from species conservation. It is important to mention that the enthalpy flux in the gas control volumes is that of the bulk gas flow, as well as the gas that diffuses to-from the GDL control volume. The enthalpy of the fluid or gas of Eq. (2) is determined as: where ៝ X are the component mass fractions and C P are the temperature dependent specific heats at constant pressure for each species. The number of moles in the control volume of Eq. (2) can be determined from the ideal gas law: The temperature of each lumped GDL and MEA control volume is found by combining the gas and solid control volume conservation equations. In addition, irreversibilities associated with the electrochemical reactions are modeled as heat generated in this control volume as follows: (5) where H is the enthalpy of formation of water, I the fuel cell current, and V is the fuel cell voltage. From each conservation of energy ODE it is possible to directly determine the temperature of each control volume, which is assumed to be the same as the control volume exit temperature.

Species conservation
As energy conservation is used to determine the temperature of each control volume, species mole numbers at each gas control volume exit is determined from species mass conservation. Specifically, the exit mole number of each bulk gas control volume is found from the following species conservation equation: The exit molar flow rate is determined from the total species conservation equation: In Eqs. (6) and (7), ៝ Φ is the species diffusion flux from adjacent GDL control volume. In order to solve for the exit mole fraction, and molar flow rate, the control volume mole fractions are assumed to be that of the control volume exit condition. We apply this often dubbed "perfectly stirred" assumption to all gas and liquid control volumes in the model. Local electrochemical reactions rates for each species (៝ R) are modeled in the GDL control volume. The species mole number of each GDL control volume is determined as follows: where Ψ H 2 O is the water diffusion flux from the adjacent MEA control volume, and Θ H 2 O is the amount of water osmotic flux through the MEA (between the two GDL control volumes) within each node. From species mole numbers, the species mole fraction and concentrations within GDL control volumes can readily be determined: Species conservation is further used to calculate the amount of water in the MEA from which membrane hydration (λ) can be determined by: where V m is the membrane dry volume, EW the membrane dry equivalent weight, and ρ m is the membrane dry density. This local membrane water content, which significantly impacts membrane ionic conductivity, is critically important to resolve to accurately determine fuel cell voltage model. From species mass conservation ordinary differential equations it is possible to determine the species mole number in each of the bulk gas streams, each GDL, and in the MEA control volume of a single node. In addition, with the perfectly stirred assumption these ODEs define exit molar flow rates in bulk gas control volumes. The species mole number can then be used to determine mole fractions, concentrations, and membrane water content as appropriate. In general, the current approach and set of assumptions allows solution of dynamic equations that can account for all species in each major subcomponent of the fuel cell repeat unit. In the present model four species are considered, namely: hydrogen, oxygen, nitrogen, and water since the modeled fuel cell operates on pure hydrogen and air.

Heat transfer
As described, energy conservation equations are used to determine temperatures throughout the fuel cell (Table 2). Specifically Eq. (1) is used to determine solid plate temperatures, Eq. (2) is used to determine bulk gas and coolant temperatures, and Eq. (5) is used to determine temperatures of MEA and GDL control volumes. In each of the three equations, the extent of heat transfer between adjacent control volumes needs to be defined.
Heat is generated in the fuel cell due to irreversibilities in the electrochemical reactions. Recall that the impact of these irreversibilities on energy conservation is accounted for in the solution of the bulk GDL and MEA thermal control volume as the difference between the total energy available from the global electrochemical reaction, and the amount of energy exiting the fuel cell as electricity ( H(1/nF ) − (VI/1000)). This concept can be visualized in Fig. 2. In the perpendicular direction (see Fig. 3) the heat generated in the GDL and MEA layer, is transferred in part to the bulk gas through convection heat transfer and by diffusion of species from the GDL to the bulk gas. In addition, some heat from the electrolyte and GDL is transferred to the solid plate through conduction. The heat in the solid plate is then transferred to the Fig. 2. V-I curve illustrating the amount of heat generated in the fuel cell. coolant channel, through convection. The electrode (solid plate) is in thermal contact with the GDL, coolant channel, as well as the bulk gas. As a result heat can transfer between the bulk gas and the solid plate electrode. The heat transfer modeled in the perpendicular direction is illustrated in Fig. 3. Note in addition, that heat is transferred between and amongst adjacent nodes as indicated in the two nodes of Fig. 3. The amount of heat generated at each location in the fuel cell can vary depending on the amount of current generated locally in the fuel cell. This variation, combined with the effects of flow oriented convective heat transfer and variations in gas and coolant flow rates and properties can lead to significant temperature gradients in the MEA and solid plates of the fuel cell. These gradients are somewhat "smoothed" by in-plane conduction heat transfer. Thus, in addition to capturing heat transfer in the perpendicular direction as shown in Fig. 3, conduction heat transfer between adjacent solid plates and natural convection at the edge of each plate is captured in the flow plane as shown in Fig. 4. Conduction heat transfer is not modeled in the MEA, because this layer is very thin.
Throughout the model convection heat transfer between solid and gas nodes are determined from Newton' law of cooling as follows: The convection coefficient (h) is determined from the Nusselt number (Nu D ) provided by Incropera and Dewitt [40] as listed in Table 3: where k f is the fluid conduction heat transfer coefficient, and D H is the hydraulic diameter. Conduction heat transfer throughout the model is determined from Fourier's law: The thermodynamic quantities used in the model are presented in Table 3.   17) and the cathode half reaction is:

Species diffusion osmotic drag and reactions
with a global reaction as shown in Eq. (13). From Faraday's law the reaction rate of both half reactions is directly proportional to the current as follows: Species diffusion is captured in the perpendicular direction between gas, GDL, and MEA control volumes. Species transport from the gas channel to the GDL accounts for the convection driven by a concentration gradient and diffusion in the GDL. The mass transport coefficient (៝ g m ) at the gas channel and GDL interface is obtained based on the Reynolds analogy between heat and mass transfer: where Sh is the Sherwood number, ៝ D m the diffusion coefficient, and D H is the hydraulic diameter of the gas flow channel. The diffusion coefficients for species are functions of temperature and pressure and are modified via the Bruggeman correlation to account for the effects of porosity and tortuosity in the GDL as follows: where ៝ D o is the species diffusion coefficient at standard pressure and temperature, ៝ D eff m the effective species diffusion coefficient, and ε is the GDL porosity. The species diffusion flux between the GDL and bulk gasses is then as: where ៝ R dif is the total diffusion resistance of each species and t gdl is the thickness of the GDL.

Water transport
Since water content in the membrane strongly affects ionic conductivity, the current dynamic model captures the details of water behavior in the MEA. Two types of water molecule transport from anode or cathode GDL to the MEA are considered: (1) the electro-osmotic drag, and (2) diffusion due to a concentration gradient between control volumes.The rate of water molecule transport via osmotic drag from anode-electrolyte interface to cathode-electrolyte interface is proportional to current density and the electro-osmotic drag coefficient as follows: The osmotic drag coefficient (n d ), is calculated from the membrane water content (λ), which depends on the water activity (a), as follows [41]: The water diffusion due to the concentration gradient between the two GDL and MEA control volumes is calculated by: (27) where D w is the diffusion coefficient of water in the electrolyte and t mea is the thickness of the MEA. The water concentration in the membrane is calculated from the membrane water content: The diffusion coefficient of water in the electrolyte (D w ) is calculated from the empirical equation [41]: where

Electrochemical model
By resolving species mole fractions, molar flow rates, and temperatures dynamically throughout the fuel cell and by assuming the fuel cell voltage is in quasi-equilibrium with the dynamic state of the fuel cell (Assumption 9), it is possible to determine the fuel cell operating voltage, and current generation distribution throughout the fuel cell. It is important to note that because the fuel cell electrodes are good conductors the voltage difference across any one cell is assumed to be constant for all parts of the cell (equipotential Assumption 8). As a result each nodal current must be determined such that all node voltages are equivalent. Before explaining the solution procedure, the voltage-current relationships are explained.
Each nodal voltage is determined by subtracting locally calculated activation, and ohmic polarization from the Nernst voltage: The Nernst voltage at each node is determined based on the MEA temperature, and the local species mole fractions in the GDL control volumes, which are representative of concentrations at the triple phase boundary: The activation polarization is modeled from the Tafel equation based on the local GDL control volume states: where a c is the activation polarization coefficient, and i o is the current exchange density. The ohmic polarization is modeled as determined in [2] for Nafion 117 based on the electrolyte control volume temperature and hydration as follows: The voltage is determined at each node from Eqs. (24)- (27), but to satisfy the equipotential assumption, each nodal voltage must be equivalent. In addition, the fuel cell must satisfy Ohm's law for the external circuit: That is, the sum of all the nodal currents multiplied by the external resistance must be equal to each nodal voltage, or equivalently the cell voltage. As a result, each nodal current is iterated until all the nodal voltages are equivalent and ohms law is satisfied. It is important to note that the amount of current at each node effects species mole fractions, and temperatures throughout the fuel cell, which in turn affect the voltage. Fig. 6 illustrates the algebraic loop that must be solved to resolve the fuel cell voltage and current distribution. While this solution strategy is being executed for each time-step, the external resistance can be manipulated, using a feedback loop, to control the fuel cell current, power, or voltage.

Experimental setup
A unit cell of a PEM fuel cell with an active area of 25 cm 2 was used to validate the simulation results. The unit cell contained 7-channel serpentine flow channels for gas delivery to both the cathode and anode of the cell with a counterflow direction. This basic configuration is presented schematically in Fig. 4. The channel depth is 1 mm and width is 0.7 mm, and the channel turns are 7 mm wide near the cell edge. The MEA used in this study is based on Nafion 112. The detailed specifications of the experimental cell configuration are presented in Table 3. Fig. 7 shows a schematic diagram of experimental setup, which consists of a gas supply unit, a unit fuel cell, an electronic load device, various sensors, a personal-computer (PC) based data acquisition system, gas sampling equipment, and a gas chromatograph (HP 5890 Series II).
High-purity hydrogen at the anode and dry air at the cathode are used as reactant gases and are humidified by passing each gas stream through a bubble type external humidifier. The flow rates of hydrogen and air were controlled and measured by two mass flow controllers. The humidification temperatures of the reactant gases were controlled by adjusting the humidifier temperature.
In order to prevent water from condensing in the sample line for species measurement, a line heater with controller was used to maintain the temperature of the gas sample above 100 • C. Species sampling points on the cathode side of the fuel cell are shown in Fig. 8. The size of each sample hole is 0.7 mm, which is the same as the width of the flow passage. The sampling points are evenly located along the gas channel spaced at 1/8th of the total channel length. To minimize the effects of gas sampling on the cell performance, the gas sampling rate was less than 5% of the total flow rate in the cathode channel. A cartridge heater was used to maintain the cell temperature as a constant during the experiment, which was controlled by a PID temperature controller.

Steady state validation and analysis
Comparisons between the simulation and the experiment were made in terms of the cell polarization curve as well as oxygen, nitrogen, and water vapor molar concentration along the length of the cathode gas channel. The polarization curve comparison result is shown in Fig. 9. The model well simulates the overall cell polarization. Note that this good comparison was achieved by selecting a single point (labeled "reference point") for establishing the constants in the cell polarization Eqs. (34) and (35) and holding them constant for all other conditions simulated. Relative humidity and temperature of both anode and cathode inlets are held at 96% and 70 • C, respectively, and the unit cell is also maintained at a constant operating temperature of  Comparison of oxygen, nitrogen, and water vapor mole fraction between experiment and simulation at different sampling locations is shown in Fig. 10. The relative humidity and temperature of both anode and cathode is 50% and 70 • C, respectively, in this case. The operating current is 15 A. Oxygen mole fraction at the cathode inlet is 0.167 and monotonically decreases along the cathode flow channel due to electrochemical reaction in the cathode compartment. Along the flow channel some oxygen is consumed at the MEA and two H 2 O molecules are produced for oxygen molecule consumed. Thus, the total molar flow rate increases and since water vapor is not condensed in the channel due to low relative humidity the oxygen and nitrogen mole fraction along the cathode channel is decreases considerably. If there were no H 2 O transport between the GDL layers by electro-osmotic drag or by back diffusion, the theoretical value of oxygen mole fraction at the cathode exit would be 0.113. Note that this value is close to the simulation result at the cathode outlet. This indicates that net H 2 O transferred from the anode GDL to cathode GDL is very small compared to the amount of H 2 O generated in the cathode due to electrochemical reaction. As a result, the magnitude of H 2 O transport by osmotic drag is almost equivalent to that of back diffusion from cathode to anode in the current experiment. The simulation result of oxygen, nitrogen, and water vapor molar concentration along the cathode channel well predicts the experimental data.
One of the main features of the current model is the ability to predict the local characteristics of current, temperature, and species concentrations throughout a PEM fuel cell without a complicated computational fluid dynamics (CFD) simulation. Because of the simplified resolution, the model can be integrated into system level models. To demonstrate the model capability to resolve current density distribution, membrane water content, and water flux providing insight into potential operating issues, a high current density operation (1 A cm −2 at 0.4 V) was simulated. The fuel cell operated at 70 • C, with the inlet cathode and anode gas humidified to 60% at 70 • C. The current Fig. 11. Distribution of (a) current density and (b) membrane water content at pressure of 1.1 bar, cell temperature of 70 • C, and relative humidity of 60%. density and water content distribution are presented in Fig. 11. The local current density is fairly uniform but increases slightly along the anode gas stream channel. As Fig. 11b indicates, the membrane water content is also quite uniform, so that the MEA is adequately hydrated over the entire active area even though the relative humidity of both the anode and cathode inlet flows is low. Thus, the local ohmic loss in the membrane, which is strongly dependent upon the water content, is almost constant and the high current density observed near the anode outlet is due to high oxygen concentration (affecting the Nernst potential) near the cathode inlet.
Local distribution of hydrogen is shown in Fig. 12. Hydrogen mole fraction along the length of the anode channel is almost  constant, even though the overall hydrogen utilization is 0.52. This is due an almost constant net water flux between the anode to the cathode. About 28% of the H 2 O that entered via the anode gas stream was transferred to the cathode side due to high electro-osmotic drag, even though back diffusion transported some water back to the anode GDL. As a result, even though more than half of the hydrogen is consumed, H 2 O mole fraction in the anode gas stream increases only slightly from 0.167 to 0.25. Fig. 13 shows the water flux magnitude of the electro-osmotic drag process, the diffusion process as driven by the concentration gradient, and the net flux from the anode to the MEA. As the water flux due to the electro-osmotic drag is proportional to the protonic flux, the water flux due to electroosmotic drag gradually increases along the anode channel length. But the water flux due to back diffusion peaks at the anode inlet and decreases along the anode channel length. As a result, the net flux between anode GDL and cathode GDL changes sign at a distance of about one-fourth of the channel length, which shows the internal circulation of water that well hydrates the membrane.
Water molar concentrations at various locations in the cell components (perpendicular direction) are shown in Fig. 14 for several specific nodes (numbers 1, 18, and 35). Note that the water concentration in the membrane is the fictitious vapor concentration that is in thermodynamic equilibrium with membrane water content. One-to-one correspondence between the fictitious vapor concentration and water concentration based on water content can be obtained from Eqs. (27) and (29). The water concentration gradient is significant across the membrane because the mass diffusivity of water in the membrane is two orders of magnitude smaller than that in the other regions. Therefore, the majority of the water that is produced at the cathode membrane interface is transported through the cathode GDL into the cathode gas channel. At the current temperature of operation, the water vapor is already saturated in the cathode GDL by the middle of the channel length. Although water content in both the anode and cathode compartments increases, the water molar concentration along the cathode channel length (from node 35 to node 1) increases more significantly than it does along the anode channel length (from node 1 to node 35). On the anode side the water vapor is transferred from the GDL to the channel near the anode inlet and then the direction of water transfer is reversed due to the substantial hydrogen reduction along the length of the anode channel and back diffusion. But on the cathode side water vapor is always transferred from the GDL to the gas channel.
Oxygen mole fraction is substantially decreased along the air flow direction from the cathode inlet due to cathodic electrochemical reactions that consume oxygen and produce H 2 O. Thus, in this high current density case with good membrane hydration, oxygen mole fraction ends up being the major factor that determines the local current density. Relative humidity in the cathode gas stream is above 100% near the cathode outlet, which indicates that water condensation would likely occur under these operating conditions causing flooding in the cathode channel. The current model indicates that the majority of the water produced on the cathode is transported through the cathode GDL into the cathode gas channel. Cathode flooding by liquid water would significantly hinder the access of oxygen to the cathode catalyst layer and reduces the number of electrochemically active sites, resulting in a substantial decrease in local current. The current model does not contain a detailed understanding of water condensation physics or droplet formation, but, it can indicate the likely location of channel flooding as near the point of water saturation. No mechanism for affecting the cell performance due to flooding was included in the current model. As a result, the predicted local current density near the cathode outlet (beyond the point where water saturation occurred) may be higher than observed current densities.

Transient comparison and controls
While the model can be used for steady state analysis, the main feature of the model is transient capabilities resolving a simplified geometry in such a way it can be integrated into system level models. Results showing system level modeling the quasi-three dimensional PEM model are not presented herein because the focus of the paper is to present the quasi-three dimensional model and validate it. The models voltage response to an instantaneous 10-15 A increase was compared to the single cell experiment (Fig. 15). The fuel cell operated at 70 • C with gas at a 50% relative humidity. To avoid hydrogen and oxygen starvation the utilization of each was 0.4748 and 0.2891 respectively. Following the current increase the voltage dropped from 0.69 to 0.58 V instantaneously. The instantaneous decrease in voltage is a result of increased ohmic and activation polarization. Following the instantaneous decrease the voltage slowly increased to 0.62 V in about 75 s. The increase in voltage is a result of increased membrane hydration which results in a decrease in the cells internal resistance. As the current increased, the rate of water formation increases at the cathode which results in an increase in membrane hydration. The simulated transient was on the same order of magnitude as expected from the simple transient analysis presented in Table 1. Fig. 16 shows the membrane hydration during the voltage undershoot and at steady state following the current increase, illustrating the increase in membrane hydration that causes the slow voltage transient. Overall the simulated response represents very well the experimental response. There is a very slight offset in the voltage steady state response and a slightly larger undershoot in the simulated transient, but overall the simulated responses well match the experiment.
During load changes, transients in membrane hydration are impossible to control exactly. Therefore, robust control strategies will have to be implemented in fuel cell systems to reject disturbances such as small local variations in membrane hydration. The model developed herein is a powerful tool that can be used for such controls development. To demonstrate this feature of the current model a power controller supplemented by current-based fuel control is implemented using the cell model in Simulink ® . The power controller is shown in Fig. 17. Current is manipulated based on power using a proportional and integral feedback together with feed forward on the power demand. Integral feedback makes it possible to track power with zero steady state error.
Controlling the fuel flow rate in proportion to current makes it possible to operate the fuel cell at a desired high utilization, increasing the efficiency of the fuel cell. Utilization is the ratio of hydrogen consumed in the fuel cell to the amount of hydrogen From Eq. (19), it is known that the amount of hydrogen consumed is proportional to the current. Hence utilization can be represented as:  Therefore, to operate the fuel cell at constant utilization, the flow rate can be controlled by: Eq. (37) is a simple rearrangement of Eq. (36) for a fixed utilization. Current-based fuel control, as shown in Eq. (37), together with the robust power controller shown in Fig. 17 were implemented to control a power demand increase using the quasi-three dimensional PEM model. A controlled power increase from 7 to 10 W is demonstrated as shown in Fig. 18. The set-point power demand was tracked almost exactly by the model. During the transient the utilization remained nearly constant, with only a slight variation in the fuel cell.

Conclusions
A quasi-three dimensional dynamic model of a PEM fuel cell has been developed in Matlab-Simulink ® . The model is based upon the dynamic equations that govern the first principles of mass and energy conservation, mass and heat transfer, and electrochemical reaction with special attention to the details of water transport within the cell. In addition, the fuel cell model captures some of the geometric features of the fuel cell as these dynamic equations are solved in a simple characteristic geometric configuration. The fuel cell is discretized into (1) 5-7 control volumes through the membrane electrode assembly, and (2) 35 nodes in the stream-wise direction. Five perpendicular control volumes are used to solve the dynamic species and mass conservation equations and six perpendicular control volumes are used to solve the dynamic energy balance and to capture the details of MEA/GDL behavior, such as water transport, which is critical to accurately determine polarization losses and performance. The dynamic conservation equations, primary heat transfer equations and equations of state are solved in each control volume of the node, and each of these nodes is interconnected with adjacent nodes (in the stream-wise direction) to simulate local cell performance characteristics throughout the cell cross-section.
Experiments were conducted to measure the performance of a single cell PEMFC under laboratory operating conditions for comparison to simulated results. A comparison of simulated and observed polarization curves shows that the dynamic model well predicts overall cell performance. In addition, since the model is discretized it can predict distributions of local current density, membrane water concentration, and species concentrations at both the anode and the cathode. Comparison of the simulated and observed oxygen water and nitrogen concentrations in the cathode compartment shows that the model well predicts species concentrations in the PEM fuel cell. The ability of the model to predict dynamic variations of performance characteristics in time and space based on first principles can be useful in the design of a fuel cell stack and for optimizing and controlling the dynamic cell performance for a variety of operating conditions. Also this model is shown to be a useful tool for investigating the effects of inlet conditions, operating parameters, and overall load and temperature conditions on both overall and local performance characteristics. Finally, such a dynamic model may be useful for the development of control strategies for enhancing overall PEM system performance.
The current model can capture the local current density, membrane water content, and species distributions as they depend upon operating conditions, which can be utilized to improve cell design and optimize performance for various operating conditions. This model can be used to analyze the characteristics of the individual fuel cell transients and provide the basis for designing control strategies for integrated fuel cell systems. The first principles approach provides fundamental insights into local performance characteristics and allows capturing of the important physics, chemistry, and electrochemistry that governs the dynamic performance of a PEM fuel cell. The simplified geometry approach provides insight into spatial variations and challenges while at the same time being simple enough to include in simulations for the development of control systems.