Daming Zhou

M Babak Nahid-Mobarakeh

M Hubert Razik

M Mahesh Krishnamurthy

M Demba Diallo

M Abdellatif Miraoui

M Damien Paire

M Fei Gao

M Alexandre Ravey

Modeling and Multi-Dimensional Analysis of A Proton Exchange Membrane Fuel Cell By

Why we need a PEMFC model

I would like to give the most sincere thanks to Prof. Dr. Fei Gao, who are respectively my supervisor in FEMTO-ST and FCLAB. I am grateful for his guidance and support in my research work and the preparation of this dissertation. I would also like to express my deepest appreciation to Prof. Dr. Alexandre Ravey and Prof. Dr. Abdellatif Miraoui, who are also my thesis supervisors, for their insightful discussion about this work and constructive comments on the dissertation. I would also thank the colleagues in FEMTO-ST and FCLAB, for their working on the experiments corresponding to this work. Many thanks for the friends in Belfort.

I would like to give my sincere thanks to Miss MengLi Yin for her meticulous care and companionship in Strasbourg and Belfort.

Finally, I would also like to express my deepest gratitude to Miss JiaQi Shang, she gives me her companionship and a very precious journey in Meerane, Germany. This is a really unforgettable memory during my life. Thank you again for accompanying me, my forever dear Meerane, forever dear JiaQi, forever dear Jiji Zha.

INTRODUCTION

During the last few decades, electric vehicles (EVs), and most recently hybrid electric vehicles powered by fuel cells (FCHEVs) have a fast growing interest due to environment pollution and energy crisis. Compared to conventional thermal machine, fuel cell power generation system shows some significant advantages, such as high conversion efficiency, reduced greenhouse gas emissions, and fast fueling time. With these advantages, fuel cells have been widely considered as a more suitable energy device for long-range EVs. Among different fuel cell types, proton exchange membrane fuel cell (PEMFC) has the aforementioned advantages. In addition, it can operate at lower temperature and lower pressure with higher power density compared to other types of fuel cells.

Nevertheless, before mass commercialization of PEMFC, there are still many research works need to be done: 1) Design of appropriate control strategies and auxiliaries to achieve PEMFC optimal working modes (air compressor, humidifier, cooling circuit, power converters, etc.);

2) Fuel cell lifespan should be increased in order to meet the requirements of different applications, 5000 hours are required for transportation applications, and 100,000 hours are required for stationary;

3) Fuel cell vehicles are currently too expensive to compete with hybrids and conventional gasoline and diesel vehicles. The production costs of the PEMFC stack and hydrogen storage needs to be reduced;

4) The solutions of producing, transporting, and dispensing hydrogen need to widespread implementation.

To achieve these goals, it is important to develop advance real-time control and monitoring methods to optimize the fuel cell operation. However, it is very difficult to observe the internal variables and state of fuel cell during its operation. Since PEMFC incorporates different control variables in different physical domains with different time constant ranges, such as voltage transient due to double layer capacitance, gas pressure variation due to the volume of gas manifold, water content variation due to the water absorption in membrane and temperature variation due to the cell thermal capacity, all these dynamic phenomena in different physical domains are indeed intercoupled between each other and the variation of one would influence another.

A good understand of how these parameters impact the fuel cell performance would be very useful for fuel cell system design and control development. Therefore, an accurate multi-physical PEMFC model could greatly help the system control strategy development and the parameters sensitivity analyses. The main objective of this thesis is to present a dynamic multi-dimensional multi-physical PEMFC (electrochemical, fluidic, and thermal) model, which can be used for control coupling analysis, spatial distribution of physical quantities prediction, real-time control implementation and prognostic. The main contributions of this thesis can be summarized as follows:

1) A multi-dimensional dynamic modeling approach for a PEMFC is developed. The proposed PEMFC model covers multi-physical domains for electrochemical, fluidic, and thermal features;

2) A variable coupling analysis of fuel cell dynamic behaviors is presented and discussed based on the developed dynamic PEMFC model. This coupling analysis can help engineers to design and optimize the fuel cell control strategies, especially for the water and thermal management in fuel cell systems;

3) A 2-D modeling approach for PEMFC for real-time control implementation is developed. The practical feasibilities of the modeling approach for advanced real-time control of PEMFC systems have been experimentally demonstrated; 4) Two novel approaches for PEMFC performance degradation prediction are proposed. These prediction methods have been experimentally validated and their strong capabilities on forecasting the future trend of PEMFC degradation voltage under different fuel cell operation conditions have been demonstrated. This thesis is organized as follows: Preface presents a brief introduction of proton exchange membrane fuel cell, and gives an overview on PEMFC multi-dimensional modeling approaches.

Chapter I proposes a dynamic multi-physical model of a proton exchange membrane fuel cell which considers electrical, fluidic and thermal domains. In addition, an innovative 2-dimensional modeling approach who considers in particular the fuel cell flow field geometric form is presented, in order to fully consider the characteristics of reactant gas convection in the serpentine gas pipeline and diffusion phenomena through the gas diffusion layer (GDL).

Based on the PEMFC dynamic multi-physical model developed in the previous chapter, a variable coupling analysis of fuel cell dynamic behaviors is presented and discussed in the first part of Chapter II. The analyses of dynamic phenomena step responses are conducted using the relative gain array (RGA) for various control input variables.

Chapter III proposes a novel real-time modeling approach based on the 2-D PEMFC multi-physical model developed in the Chapter I. In this approach, differential equations for reactant gas convection and diffusion phenomena in serpentine channels are transformed into tridiagonal systems of equations, in order to use an efficient numerical solver tridiagonal matrix algorithm (TDMA). In addition, a three levels bisection algorithm has been developed to solve spatial physical quantity distributions of electrochemical domain.

In the first part of Chapter IV, A multi-physical aging model has been proposed in order to predict the output voltage degradation of PEMFC. In the proposed aging model, three most important aging phenomena during PEMFC operation are considered. In addition, particle filter and extrapolation approach are used to estimate the aging parameters. In the second part of Chapter IV, an innovative approach for PEMFC aging prediction based on a combination of model-based and data-driven methods is presented. The proposed hybrid prognostic method is able to capture both the fade trend and non-linear features observed in the fuel cell voltage degradation data.

A conclusion is given at last.

As another energy power source widely used in the FCHEV powertrain, the lithium-ion batteries have advantages of high energy density and long cycle life. In order to accurately estimate state of charge (SOC) of battery during the FCHEVs operation, a novel approach for battery SOC estimation based on multi-models data fusion technique is presented Appendix A.

The last focus of this thesis is energy management strategy for FCHEVs, since it directly affects the efficiency and performance of energy storages in FCHEVs. For example, by using a suitable energy distribution controller, the fuel cell system can be maintained in a high efficiency region and thus saving hydrogen consumption. In Appendix B, an on-line adaptive energy management control is proposed based on extremum seeking method and fractional-order calculus, in order to improve both the performance and durability of PEMFC used in the FCHEVs. Nowadays, the research on fuel cells is a very active field, since fuel cells have been considered as one of the most attractive green energy generation device [1]. Different from lithium-ion battery [2] [3], fuel cells are energy conversion devices which require a continuous source of fuel and oxidant (usually air or pure oxygen) to sustain the chemical reaction, in order to converts chemical energy into electricity, whereas the lithium-ion battery is energy storage device.

There are many different fuel cell types. For example, a single proton exchange membrane fuel cell generally consists of two separate electrodes cathode and anode, as well as an electrolyte. The gaseous fuels and oxygen respectively transport in the gas pipeline of cathode and anode, and diffuse from the gas diffusion layer to the catalyst layer. In the anode catalyst layer, the fuel is oxidized to generate electrons and protons (positively charged hydrogen ions), these protons flow from the anode side to the cathode side through the electrolyte (membrane), and are further reduced by absorbing electrons and protons to produce water. At the same time, electrons are drawn from the anode electrode to the cathode electrode through an external electrical circuit, producing directly the electricity. The complete electrochemical reaction in the fuel cell system can be written as the following [4]:

Based on their functionalities, a single proton exchange membrane fuel cell in a fuel cell stack can be separated into eight functional layers: one cooling channel, one membrane, two catalyst layers, two gas diffusion layers and two gas supply channels [5].

It should be noted that, electrical potential of an individual cell is relatively low, this voltage is generally around 0.7V. In order to provide sufficient electrical power to meet the requirement of applications, a fuel cell stack is commonly composed (placed in series) of a dozen or even a hundred individual cells [6].

Different types of Fuel Cells

Fuel Cells can be mainly categorized by electrolyte types. The electrolyte types determine the catalysts types and chemical reactions inside the fuel cell, as well as the operation temperature. The commonly used fuel cell types are listed in the following subsections.

Proton Exchange Membrane Fuel Cell (PEMFC)

Proton exchange membrane fuel cell is considered as a more suitable energy device for mobile applications, such as hybrid electric vehicles or portable power supply, since it can operate at lower temperature and lower pressure with higher power density compared to other types of fuel cells [7]. In the proton exchange membrane fuel cell Nafion (copolymer fluoropolymer) are usually used materials for exchange membrane, platinum is widely used as the catalyst for PEMFC. The operation temperature of proton exchange membrane fuel cell system is around 80 , its chemical reactions can be expressed as follows

The electrochemical half-reaction occurs at the anode electrode:

The electrochemical half-reaction occurs at the cathode electrode:

Alkaline Fuel Cell (AFC)

The alkaline fuel cell (or hydrogen-oxygen fuel cell) is one of the most developed fuel cell, which is used in Apollo space program to provide the source of electrical energy and drinking water [8]. The design of alkaline fuel cell is similar to the PEMFC, but the electrolyte used in the alkaline fuel cell is generally a porous matrix saturated with an aqueous alkaline solution, for example the potassium hydroxide (KOH) or Sodium hydroxide (NaOH). Its operation temperature is similar to PEMFC, the materials of catalyst required for the electrodes can be selected from a number of relatively inexpensive chemicals, but with a lower current density. The chemical reactions for alkaline fuel cell system can be expressed as follows:

The electrochemical half-reaction occurs at the anode electrode:

The electrochemical half-reaction occurs at the cathode electrode: 0.2.3 Phosphoric Acid Fuel Cell (PAFC)

In the phosphoric acid fuel cell, the phosphoric acid is used as a non-conductive liquid acid (electrolyte) to transfer positive hydrogen ions from anode to cathode through an external electrical circuit [9]. Since the phosphoric acid fuel cell has a simple and stable structure, it is generally used in the stationary applications. The phosphoric acid fuel cell system can operate efficiently in the temperature range from 150 to 200 , its chemical reactions can be expressed as follows:

The electrochemical half-reaction occurs at the anode electrode:

The electrochemical half-reaction occurs at the cathode electrode: 0.2.4 Molten Carbonate Fuel Cell (MCFC)

The Molten carbonate fuel cell is quite different from the previous fuel cell types: it uses an electrolyte composed of molten carbonate salts suspended in a porous ceramic matrix and chemically inert solid electrolyte of alumina beta [10]. The molten carbonate salts in this type produce the migrate ion from the cathode to the anode, then the hydrogen with carbonate ions from the electrolyte to produce water, carbon dioxide, and electrons. The operating temperature of molten carbonate fuel cells is above 650 , its chemical reactions can be expressed as follows:

The electrochemical half-reaction occurs at the anode electrode:

The electrochemical half-reaction occurs at the cathode electrode: 0.2.5 Solid Oxide Fuel Cell (SOFC)

The solid oxide fuel cell uses a solid material as the electrolyte, such as dense ceramic materials yttrium-stabilized zirconia (YSZ), which separates gases from the anode and the cathode, blocks electrons and conducts oxygen ions from the cathode to the anode [11]. Since the solid oxide fuel cell requires high operating temperatures (from 800 to 1000 ), it is generally used for stationary applications. The chemical reactions for solid oxide fuel cell system can be expressed as follows:

The electrochemical half-reaction occurs at the anode electrode:

The electrochemical half-reaction occurs at the cathode electrode: 0.2.6 Direct Methanol Fuel Cell (DMFC)

The methanol fueled fuel cell is derived based on the proton exchange membrane fuel cell. It directly uses methanol ( , in a liquid form) as the fuel. The main advantage of methanol fueled fuel cell is that the methanol is a relatively stable liquid, which is easy and low cost for transportation [12]. The efficiency of methanol fueled fuel cell is lower than other fuel cell types, it is generally used for portable applications, where the energy density is more important than efficiency. During methanol fueled fuel cell operation, the methanol is firstly converted to carbon dioxide and hydrogen at the anode, and the remaining steps of the reaction are similar to the PEMFC, its chemical reactions can be expressed as follows:

The electrochemical half-reaction occurs at the anode electrode:

The electrochemical half-reaction occurs at the cathode electrode:

Why We Need A Proton Exchange Membrane Fuel Cell Model

Nowadays, research on proton exchange membrane fuel cells (PEMFC) has made major advances in sustainability, cost and compactness [13] [14], compared to other types of fuel cell, the PEMFC can provide higher power density for transport and portable applications with relatively short start-up time and lower operation temperature and pressure [15].

Nevertheless, before mass commercialization of PEMFC, one of the major challenges in PEMFC research is the development of appropriate control strategy for PEMFC stack and system auxiliaries (i.e., air compressor, cooling circuit, power converter), in order to maintain the optimal operation conditions of fuel cell system [16]. In addition, the PEMFC stack is a very compact device since it incorporates different phenomena in different physical domains. During fuel cell operation, these dynamic phenomena are indeed inter-coupled between each other and the variation of one would influence another. In practice, it is very difficult to observe the internal variables and state of fuel cell during its operation.

In order to get a good understand of how these parameters quantitatively impact the fuel cell performance, and further help engineers to design and optimize the fuel cell control strategies, one possible solution for this problem is using a model-based control method. This brings up the need of an accurate and precise PEMFC dynamic model, which at least considers the following issues:

1) Dynamic behavior should be considered for transient state control;

2) The developed model should have a capability to provide multi-dimensional behaviors, which is very useful to give insights into the interaction effects of parameters on the fuel cell spatial performance;

3) The practical feasibilities of advanced real-time control of PEMFC systems should be considered, in order to effectively perform quantitative analysis of fuel cell performance and make fast control decisions.

4) Consideration of major internal physical aging phenomena of fuel cell for degradation prediction, for example the fuel cell ohmic losses, reaction activity losses, and reactants mass transfer losses.

The above issues are detailed presented in the following subsections.

Dynamic Behaviors

The control of the fuel cell system (air compressor, cooling circuit, power converter, etc.) is a very complicate work because it incorporates different control variables in different physical domains [17]. During fuel cell operation, different dynamic phenomena within different time constant ranges, such as voltage transient due to double layer capacitance, gas pressure variation due to the volume of gas manifold, water content variation due to the water absorption in membrane and temperature variation due to the cell thermal capacity, can be clearly observed during load transient [18]. These dynamic phenomena in different physical domains are indeed inter-coupled between each other and the variation of one would influence another. This inter-coupling effect is especially important between the dynamic phenomena which have similar transient time constants. Thus, all these dynamic phenomena should be considered in the developed PEMFC model.

Spatial Distribution of Physical Quantities

Compared with one-dimensional models, a multi-dimensional PEMFC model has a capability to provide spatial distribution of physical quantities, which is very useful for spatial non-uniformity and control coupling analysis [19].

For example, the one-dimensional modeling of fuel cell bipolar plate flow field are too simplified and do not represent accurately the pressure distribution characteristics, since the fuel cell gas supply pipeline is assumed to be single and straight. In reality, cathode and anode gas supply channels may be of different patterns like single serpentine, parallel serpentine or inter-digital channels. Therefore, a comprehensive representation of non-homogeneous gas phenomenon by fully taking the geometric form of the fuel cell pipeline into consideration is particularly useful to achieve highly accurate spatial distribution information for two-dimensional model of PEMFC.

Real-Time Applications

Different from the common modeling approach, a real-time oriented fuel cell model has more restrictions: the accuracy and computational efficiency of a real-time fuel cell model are both crucial for model based control process [20]. A sophisticated fuel cell model can provide comprehensive physical quantities for model-based control design and optimization. While the high performance computation of a fuel cell model ensures the model-based controller can be efficiently implemented in real-time applications with a low cost of computations.

Degradation Prediction

It is meaningful to develop a multi-physical aging model for degradation prediction of fuel cell performance [21]. This multi-physical aging model considers the real physical aging phenomena during the PEMFC degradation process. Although the model-based methods need large computations and complex physical model, it can predict not only the system degradation trend (fuel cell output voltage decay over time), but also the information about the internal physical parameters during the degradation process. A common drawback of these models is that the different fuel cell dynamic phenomena, especially the ones with similar time constant, are not considered simultaneously or over-simplified, thus make them unsuitable for dynamic variable coupling analysis.

Compare to one-dimensional models, multi-dimensional models have a capability to provide local phenomena and spatial distribution physical variables, which is very useful for spatial non-uniformity and control coupling analysis. However, calculating such complex physical quantities leads to higher computational requirements. Many PEMFC multi-dimensional models have been previously proposed in the literature [5]- [8].

Y. Shan et al. [5] propose a dynamic two-dimensional PEMFC model, which considers the fluidic and thermal behaviors. In order to obtain the temperature dynamic distribution along the gas channel direction and through-plane direction, a numerical solver is used based on SIMPLE algorithm. In addition to the thermal dynamic behaviors, the proposed 2-D model can also predict the current density and oxygen concentration dynamic distribution. B. Sivertsen et al. [8] introduce a comprehensive non-isothermal three-dimensional PEMFC model. In the presented model, the fluid transport inside the fuel cell gas channels and gas diffusion layer, as well as the thermal behaviors are developed and solved based on the framework of a CFD code. This CFD computational model can accurately predict the cathode over-potential distribution.

A common drawback of these works is that the presented fuel cell bipolar plate flow field (gas channels) models are over-simplified (or not even considered). Thus they cannot describe accurately the non-linear and non-uniform pressure distribution characteristics. On the other hand, as a commonly used modeling technique for multidimensional model, the CFD models [7]- [8] are not suitable for real-time model-based controller since the computational burdens are too heavy.

Multi-Physical PEMFC Model

In this section, the presented PEMFC dynamic model is based on a developed multiphysical PEMFC model in the previous works [9]- [10]. The PEMFC stack level and the single cell level are shown in the Figure 1.1. As shown in Figure 1.1, one fuel cell stack level can be separated into eight cell layers, which consist of:

1) Cathode cooling channel layer;

2) Cathode gas supply channel layer;

3) Cathode gas diffusion layer; The advantage of this cell layer structure is that each modeling layer can be described separately by its own physical equations and the boundary conditions. Each cell layer is considered as a control volume in the presented model. For each cell layer, the modeling of different physical domain is presented in the following subsections.

Electrical Domain Modeling

As an electricity-converting device, the PEMFC converts fuel energy into electricity through electrochemical reactions. Therefore, the electric domain is included in the proposed PEMFC model.

The total output voltage of a single-cell can be calculated by the following equation:

where is the single fuel cell thermodynamic voltage (V), is the ohmic voltage drop (V), is the cell activation voltage drop (V).

The thermodynamic voltage is calculated from the following Nernst equation [11]:

where is the catalyst layer temperature (K), is the Faraday constant (C/mol), is the ideal gas constant, is the oxygen pressure (atm) at the interface of cathode catalyst layer, is the hydrogen pressure (atm) at the interface of anode catalyst layer (please refer to the last paragraph of this section).

The membrane resistance (Ω) is calculated by [12]:

where is the section surface of membrane (m 2 ), is the membrane thickness (m). is the resistivity of membrane (Ω m) which can be calculated by the following equation [12]:

Thus, the cell ohmic voltage drop can be calculated from Ohm's law equation [13]:

The electrochemical activation voltage drop of single cell can be calculated by Butler-Volmer equation [14]:

where is the stack current (A), is the charge transfer coefficient, is the electrons number. The exchange current density (A/m 2 ) can be calculated by an empirical equation [15]:

where and are empirical parameters, is the oxygen activation energy on the electrode catalyst interface. It should be noted that, the in anode side for fuel cells of PEMFC type can be neglected due to the easy electrochemical process.

It can be seen from the equation 1.6 that, the term activation voltage drop is in an implicit form. In order to explicitly calculate this non-linear implicit Butler-Volmer equation to obtain the value of , an iterative solver is used. For a high stack current, the Butler-Volmer equation can be written as the well-known Tafel equation [16]:

The dynamic behavior of activation losses voltage in the electrical domain due to the "double layer effect" can be expressed by:

Fluidic Domain Modeling

It should also be noted that, since the and used in equation 1.2 are reactant gas pressure at the catalyst layer interface instead of the gas supply channels,

another fuel cell over-potential term due to pressures drop through the GDL, well known as "concentration losses", has been already implicitly considered in the proposed fluidic model. Therefore, the fluidic behaviors inside fuel cell, such as reactants convection in the channels and diffusion through gas diffusion layer, have a great impact on fuel cell performance. In this section, the fluidic domain modeling is presented.

Cooling Channels

In the gas channel, the Reynolds number can be calculated by the following equation [17]:

where is the channel hydraulic diameter (m), is the mean fluid velocity in the channel (Pa s). The mean velocity (m/s) of gas can be calculated based on the following equation:

where is the fluid mass flow (kg/s), A is the total section of cooling channels (m 2 ). The fluid density (kg/m 3 ) can be calculated by ideal gas equation of state:

The gas pressure drop of serpentine channel depends on the surface friction losses of straight pipeline, which can be modeled by the Darcy-Weisbach equation [18]: with where , and are respectively the pressure at inlet, outlet, and center of cooling channel, is the total length of straight channel (m), is the Darcy friction factor which can be obtained from the empirical equation [19]:

The gas pressure dynamic response in the fuel cell is generally due to the channels volume. Thus, the dynamic behaviors of fluid in the cooling channel can be given by the mass balance equation:

where is the gas molar mass (kg/mol), is the volume of the cool channels (m 3 ), is the cooling channel temperature (K), is the gas pressure in the cool channels (Pa) and is the fluid mass flow rate (kg/s) entering or leaving the channels.

Gas Supply Channels

The total pressure of the center of the gas supply channels can be calculated by: where , , are respectively the oxygen, nitrogen, and vapor pressure in the center of cathode supply channel. , are respectively the hydrogen and vapor pressure in the center of anode supply channel. The gas pressure in the center of channels is defined as follow:

where is the gas pressure at the channel inlet, and is the gas pressure at the channel outlet. Thus, the dynamic behaviors of fluid in the center of gas supply channels can be also written based on the mass balance equation:

where is the gas mass flow at the inlet of channel, is the gas mass flow at the outlet of channel, is the volume of the gas supply channels (m 3 ) and is the temperature of gas supply channel (K).

Gas Diffusion Layer (GDL)

To obtain the reactant gas pressure at interface of the catalyst layer, the diffusion phenomenon in the gas diffusion layer can be described by modified Fick's law [20] [21]:

where is pressure of specie x in gas diffusion layer (Pa), is thickness of gas diffusion layer (m), is the gas diffusion layer temperature (K), is the gas molar flow rate of specie x (mol/s), is the gas diffusion layer area (m 2 ), is the gas diffusion coefficient (m 2 /s) between the species x and y can be calculated from [22]:

where is the total pressure of species (atm), is the critical temperature of species (K), is the critical pressure of species (atm), and is the molar mass of species (kg/mol), is the porosity of the GDL and is the GDL tortuosity. The coefficients and depend on whether one of the species is a polar gas or not and are determined accordingly, which are given as follows [22]:

For pair of gases contains no polar gas:

For pair of gases contains polar gas:

Catalyst Layers

As mentioned before, the reactant gas mass flow rate through the GDL to the catalyst layer is directly proportional to the fuel cell stack current [23] [24]. Thus, the oxygen mass flow (kg/s)at the cathode side can be expressed by: the hydrogen mass flow (kg/s)at the anode side can be expressed by: and the mass flow of produced water (kg/s) at the cathode side due to the electrochemical reaction can be calculated by:

Dynamic Membrane Water Content

Because the membrane ionic conductivity is highly dependent on water content in polymer membrane [25], a more detailed knowledge of transient behavior of the would give a more accurate value of Ohmic losses , as shown in the equation 1.3. Moreover, the dynamic phenomena of plays an important role on the dynamic performance of PEM fuel cell due to its relatively long transient time (up to some minutes) [26].

The dynamics of the membrane water content is generally influenced by two water flow effects in the membrane: the electro-osmotic drag flow due to proton conduction from the anode to the cathode; the water back diffusion flow caused by the concentration gradient between anode and cathode side. The membrane water content is defined as the relationship of the number of water molecules per charged site (sulphonate site) [27]:

where is the water activity factor, which can be obtained based on the water local vapor partial pressure (pa), and the local vapor saturation pressure (pa):

where the local vapor saturation pressure is calculated by: where is the vapor temperature (K). Thus, the dynamics of the water content in the membrane can be obtained by considering the water molar flows balance at two sides (i.e. anode and cathode) of membrane, and the mass conservation of water, as shown in Figure 1.2. Thus, the dynamics of the water content can be described by: where is the membrane dry density (kg/m 3 ), is the molecular mass of membrane (kg/mol), and represents different water molar flow (mol/s) entering or leaving the membrane due to electro-osmotic drag and water back diffusion flow.

Gas supply channel

Gas

The water molar entering or leaving the membrane due to electro-osmotic drag can be described by: where the is the coefficient of electro-osmotic drag for maximum hydration conditions. The water molar entering (from cathode to membrane) or leaving (from membrane to anode) the membrane due to back diffusion flow can be described by [28]:

where is the dry density of the membrane (kg/m 3 ), is the equivalent mass of the membrane (kg/mol). The boundary water content at anode and cathode side can be expressed as a function of water activity which can be deduced from the water vapor partial pressure equation 1.28. The membrane water diffusion coefficient (m 2 /s) can be calculated from the empirical equations [28]:

As shown in Figure 1.2, in the case of non-humidified hydrogen supply and anode deadend mode operation, the anode side water accumulation is only caused by the water diffusion from the membrane to the anode , and the cathode side water accumulation depends on three factors: comes from humidified air supply at cathode , the produced water at cathode side during electrochemical reaction and the electro-osmotic drag flow from the anode to the cathode .

Under the dead-end mode operation (no water accumulation at anode side), the water molar flow entering into the membrane from anode due to electro-osmotic drag is equal to the water molar flow leaving out the membrane to anode due to back-diffusion (dashed portion as shown in figure 1.2). Thus, the dynamics of the membrane water content can thus be simplified by: with anode side water content :

Thermal Domain Modeling

In addition to the fluidic phenomena, the effect of temperature on the fuel cell performance should also be considered in the fuel cell modeling. For example, the heat transfer changes the gas convection and diffusion behaviors, and further influences the electrochemical quantities inside the fuel cell.

The fuel cell temperature transient behavior in the thermal domain is due to the heat generation and thermal conduction and convection phenomena. Like dynamic behavior of membrane water content, dynamics of fuel cell temperature is an important phenomenon due to its relatively long transient time. This dynamic behavior can be generally described as follows:

where is the mean layer volume density of stack (kg/m 3 ), is the layer volume of stack (m 3 ), is the layer thermal capacity (J/kg K), is the temperature (K) of each control volume and stands for the different types of heat (J) flows entering or leaving the layer respectively: conduction, convective flow, forced convection and internal heat sources.

The heat flows due to conduction can be expressed according to Fourier's law [29]:

where is the control volume thermal conductivity (W/m K), is the section of the control volume in heat transfer direction (m 2 ), and is the control volume thickness (m).

The convective heat flow due to the mass transfer entering or leaving the control volume can be calculated by where is the mass flow rate (kg/s).

The heat transfers by forced convection can be written according to Newton's cooling law:

where is the coolant temperature and is the contact area (m 2 ), is the forced convection heat transfer coefficient (W/m 2 K), which can be calculated by [29]:

where is the fluid thermal conductivity, is the Nusselt number of the fluid, which can be calculated by the empirical equations [30]:

Where is the coolant thermal capacity (J/kg K).

At last, the linear expression of the heat sources as a function of the cell temperature can be obtained:

1

.3 Multi-Dimensional Modeling Considerations

Compared with 1-D models [31]- [39], a 2-D PEMFC model has a capability to provide two-dimensional behaviors, which is very useful for spatial non-uniformity and control coupling analysis. This analysis can give detailed and valuable spatial physical quantity information under different fuel cell operation conditions by taking multiple spatial dimensions into consideration. For example to prevent local "hotspot" on electrode due to non-homogeneous distribution of reactants, and can be further employed in a model-based real-time controller.

Many PEMFC 2-D models have been previously proposed in the literature [40]- [45].

However, a common drawback of these works is that the presented fuel cell bipolar plate flow field (gas channels) models are over-simplified (or not even considered). Thus they cannot accurately describe the non-linear and non-uniform pressure distribution characteristics.

For example, as shown in the upper part of figure 1.3, the gas pressure prediction results of a model without the consideration of channel geometric form, could lead to an inaccurate gas diffusion phenomenon in the serpentine pipeline, which would further impact the accuracy of electrode current density analysis. In these models, the gas supply channel is assumed to be straight and single. In fact, the gas supply pipeline in In order to take the geometric form of the fuel cell pipeline into consideration, a 2-D modeling of fluidic domain is developed firstly, followed by a 2-D electric modeling.

Two-Dimensional Approach in Fluidic Model

To accurately model the reactant gas pressure distribution on the electrodes surface, a comprehensive modeling of gas convection-diffusion phenomena in the gas supply channel and GDL is presented hereafter by precisely considering the fuel cell gas channel geometric form in this section.

Gas Supply Channels

Different from single and straight channel assumption in the previously developed PEMFC model, the geometric patterns of gas supply channels (both anode and cathode sides) are now considered in the improved model, such as single serpentine, parallel serpentine channels, with the consideration of sharp and curved U-bends (channel angles). Under the same inlet air supply conditions, the gas pressure distribution on the surface of GDL depends highly on the flow field form. Thus, a detailed representation of non-homogeneous gas pressure distribution by considering the flow field form can be very useful to achieve accurate modeling results.

In this thesis, the geometric form of gas supply channels is taken from Ballard NEXA fuel cell as shown in figure 1 where is the Darcy friction factor, is the fluid density in the channel (kg/m 3 ), is the total length of straight channel (m), is the mean fluid velocity in the channel (m/s), is the hydraulic diameter of the channel (m), is the total length of elongated section of the bends (m), the Kays friction factor can be given by [START_REF] Maharudrayya | Pressure losses in laminar flow through serpentine channels in fuel cell stacks[END_REF]:

where is the pipeline width (m), is the pipeline thickness (m). is the Reynolds number of the fluid in the channels, channel aspect ratio is defined by:

The excess bend loss coefficient of i th bend is given by [START_REF] Maharudrayya | Pressure losses in laminar flow through serpentine channels in fuel cell stacks[END_REF]:

where is the spacer length between two neighboring duct, and is the curvature ratio of the bend, is given by [START_REF] Maharudrayya | Pressure losses in laminar flow through serpentine channels in fuel cell stacks[END_REF]: where is the mean radius of bends, is the duct hydraulic diameter which is calculated by: where the pipeline width , pipeline thickness , and spacer length between two neighboring duct can be seen clearly in the figure 1.7.

From figure 1.7, the rectangle pane is represented as the cross sectional area of pipeline.

The anode channel curved U-bends are the same as those in the cathode channel.

Gas Diffusion Layer

The modeling equations for reactant gas diffusion phenomena in the GDL are described in this section. Based on the gas supply channel geometry, the GDL, which is directly adjacent to the gas channels layer, can be divided into two sections, denoted as "fluid adjacent volume" and "solid adjacent volume", as illustrated in figure 1.8. From figure 1.8, the thin line GDL control volume is adjacent to the channel fluidic section, denoted as "fluid adjacent volumes" of GDL. The gas flows to a "fluid adjacent volumes" come from the gas channel pipeline and the adjacent volumes. In contrast, the bold line GDL control volumes in figure 1.8 are adjacent to the channel solid section, denoted as "solid adjacent volume" of GDL. The gas flows to a "solid adjacent volume" come only from the adjacent volumes in gas diffusion layer.

In the previous section, the gas convection phenomenon in the serpentine pipeline has been well developed. As shown in the figure 1.8, the convective gas flow direction in the gas pipeline is marked by arrow 1 in the gas pipeline A, and by arrow 3 in the gas pipeline B. Then, the gas diffusion phenomena through the GDL can be divided into six categories:

1) Diffusion from pipeline to "fluid adjacent volume" (z-axis, marked with arrows 2 and 4);

2) Diffusion between two adjacent "fluid adjacent volumes" (y-axis, marked with arrows 5 and 6);

3) Diffusion between adjacent "fluid adjacent volume" and "solid adjacent volume" (x-axis, marked with arrow 7). It should be noted that, when the stack current increases, the diffusion from gas pipeline to "fluid adjacent volume", and the diffusion from volumes to catalyst layer become more important due to the increase of mass flow from gas channels to the catalyst layer, which is proportional to the fuel cell current. All these diffusion phenomena in the gas diffusion layer, except the above-mentioned current-driven ones, can be modeled by the modified Fick's diffusion equation 1.22.

In fact, the benefit of parallel pipeline is its lower pressure drop due to shorter single channel length. However its main drawback is water droplet accumulation during fuel cell operation. The single serpentine can improve water removal while introduce a larger pressure drop. This parallel serpentine flow field combines thus the advantages of both patterns, in order to achieve better performance of fuel cell system.

Non-Uniform Control Volume Consideration

In order to fully describe the gas flow in the serpentine pipeline gas channels, and further accurately obtain the two-dimensional physical quantity distribution both in fluidic As shown in the figure 1.9, the two surface dimensions are considered for the proposed 2-D model, the non-uniform control volume distribution of each side is based on the geometry form of channels (i.e. three-parallel serpentine pattern at cathode side and single-parallel serpentine pattern at anode side, denoted by black mesh in the figure 1.9.

Two-Dimensional Approach in Electric Model

The characteristics of reactant gas convection and diffusion in the pipeline and gas diffusion layer are fully described in the previous section. However, the spatial physical quantity distribution on the surface of electrode (e.g. the current density distribution) cannot be directly obtained using the non-uniform mesh grid. In order to unify the mesh segments distribution in homogenous material such as electrode and electrolyte, the non-uniform mesh grid of gas channels layer of each side are then linearly converted to uniform mesh grid for GDL and membrane layers denoted by red mesh in the figure 1.10. [10]).

Based on the developed multi-physical model, a 2-D modeling approach for a protonexchange-membrane fuel cell (PEMFC) is then proposed. The proposed model covers multi-physical domains for both fluidic and electric features. In order to accurately describe the characteristics of reactant gas convection and diffusion in the pipeline and gas diffusion layer, the gas pressure drop in the serpentine pipelines is comprehensively analyzed, especially for the reactant gas pressure drop due to the pipeline U-bends, followed by a comprehensive description of gas diffusion layer 2-D modeling by fully considering the geometric form of flow field [START_REF] Zhou | Development of a multi-physical multidimensional modeling of proton exchange membrane fuel cell[END_REF]. The experimental validation and coupling effects analysis are presented in the following chapter.

Chapter II: Experimental Validation and Coupling Analysis

This chapter presents the experimental test and simulation in order to validate the

Literature Review

Although the coupling of dynamic phenomena in a fuel cell has an important influence on the control design of fuel cell system, very little information has yet been published in the literature on the analysis of interaction between fuel cell dynamic variables in different physical domains.

Zhao et al. [1] present a decoupling control strategy for the strong coupling between mass flow and pressure in centrifugal compressor system, which is a key component for supplying compressed air to the fuel cell cathode channel, but the dynamic phenomena coupling in fuel cell system have not been considered. Cheah et al. [2] give a detailed analysis of coupling effects of electro-osmotic drag, water diffusion and interfacial water transport. Carnes et al. [3] present a analysis of coupling effect of the water transport within the proton exchange membrane and the partially saturated gas diffusion electrodes. Cao et al. [4] propose a single neuron adaptive proportionalintegral derivative (PID) feedback controller for a solid oxide fuel cell, which combines the advantages of robust control and PID control, in order to automatically adjust control parameters when system encounters uncertainties and disturbances. However, they only considered the coupling effect of dynamic water transport in fluidic domain without taking the thermal domain dynamics into account.

Model Experimental Validation and Discussion

Dynamic Membrane Water Content Results and Discussions

The same current profiles (figure 2

Effect of the Gas Supply Serpentine Channels and Discussions

The major improvement of gas supply channels modeling by taking the channel geometric form into consideration can give a more accurate pressure distribution prediction in the fuel cell, which can lead to a more precise prediction of fuel cell From figure 2.8 (a), it can be seen clearly from the figure that, the predicted pressure drop is much higher with the assumption of straight gas channel form. By taking the channel geometric form into consideration, the obtained pressure drop in the channel can be differed as much as 2.94 % of total channel pressure compared to the "straight channel" assumption. It has to be noted that, this pressure prediction error would not lead to a significant error in the fuel cell voltage value. However, it could lead to a wrong gas pressure distribution pattern on the surface of electrode for the electrode current density analysis for example. Thus, in order to get a more accurate gas channel pressure results, the channel geometric form has to be taken into account in the fuel cell model. The comparison of simulation results of pressure drop of two modeling approaches with short step current is also presented in figure 2.9.

2-D Model Simulation Results and Discussions

Model Grid Independence Analysis

The model grid independence analysis determines the minimum mesh grid number needed for a multi-dimensional model. When increasing the mesh number, if the changes of model outputs are less than a pre-defined acceptable error, this mesh number can be considered to meet the independence criteria. is because the oxygen supply excess rate is higher, and further lead to a more uniform oxygen pressure distribution on the surface of electrode.

Under the same current, the current density distributions on the surface of electrode are also shown in the figure 2.12.

It can be seen from the figure 2.12 that, the current density distribution is similar to that of oxygen pressure. The oxygen pressure at catalyst layer is higher at channel inlet than outlet, which leads to a higher current density at air inlet. As shown in the figure 2.12 (a), under the condition τ=1.68, the maximum difference of current density is about 45% between fuel cell inlet and outlet. It can be also observed from the figure 2.12 (b) that, under a higher τ condition τ=4.94, the current density distribution gradient is less significant compared to figure 2.12 (a). In this case, the maximum difference of current density is only about 15% between fuel cell inlet and outlet.

Dynamic Phenomena Coupling Analysis

The dynamic fuel cell model has been presented for the purpose of fuel cell control performance optimization. The fuel cell dynamic behaviors in different physical domains are indeed inter-coupled between each other and the variation of one would influence another. The dynamic variable coupling analysis should be performed for fuel cell system in order to develop an optimized control algorithm.

The water and thermal management is very important for the fuel cell performance and efficiency. For the thermal domain, the fuel cell temperature dynamic behavior is the most significant dynamic in the fuel cell stack. The time constants of fuel cell system in thermal domain can be relatively long, due to large thermal capacities and volumes of cell components (bipolar plates, membrane, etc.) [5]. It should also be noted that, as shown and discussed in the previous section, a nonuniform distribution of current density in the fuel cell electrode can be observed due to the reactant gas pressure variation in the gas supply channels. A more homogenous distribution of current density can maintain the stability and improve the long-term performance of the fuel cell system. The current density distribution is simultaneously affected by the dynamic phenomena in both thermal and fluid domains. Therefore, it should also be carefully considered for coupling analysis in order to provide insights into the variable interaction among three different physical domains.

In this section, detailed expressions of time constant for temperature and membrane water content are given at first, followed by analyses of step responses for various fuel cell system input variables. The corresponding dynamic variable coupling analysis is introduced and discussed at last.

Expressions of Time Constant for Temperature and Membrane Water Contents

The expression of time constant of temperature in thermal domain can be deduced from the general first order dynamic form [5]:

Thus, in the thermal domain, the fuel cell temperature transient response time can be estimated to be 497s in the case of studied 1.2 kW Ballard NEXA fuel cell stack.

For the time constant of membrane water content, the left hand side of equation In order to understand the effects of fuel cell operational parameter coupling from a control point of view among three different physical domains, the fuel cell system can be considered as a multi input and multi output (MIMO) system, where four possible control input variables are listed as follows:

1) The coolant inlet temperature (controlled by heat exchanger and bypass circuit);

2) The gas channel inlet temperature (controlled by inlet/outlet gas heat exchanger);

3) The gas supply channel inlet water vapor pressure (controlled by gas humidifier);

4) The inlet air molar flow rate (controlled by air compressor).

And four controlled output variables in the proposed MIMO system are listed as follows:

1) The membrane electrode assembly (MEA) temperature ;

2) The bipolar plate temperature ;

3) The membrane water content ;

4) The uniformity coefficient of current density distribution on the electrode .

where the uniformity coefficient is the current density ratio between the highest and the lowest value on the same electrode. This coefficient is proposed to describe the degree of uniformity of current density distribution.

Thus, the non-linear state space equations of this MIMO system can then be expressed as follows:

As the above fuel cell MIMO system state space shown, the vector form representations of the manipulated input variables is , and the controlled output variables is . , , and can be derived using the formula of the physical modeling equations presented in chapter II. Since there is a complicated non-linear mathematical relationship between input and output variables, examining the dynamic responses of controllable outputs after step changes of inputs are particularly useful for having insight on the possible variable coupling. The first operating point is set to between 1/3 and 1/2 rated power point, which corresponds to the high efficiency operating region of fuel cell system. Figure 2.13 shows the dynamic response of the MIMO system outputs after step changes of the coolant inlet temperature. An increase of the coolant inlet temperature is set at 300 s, which lead to an increase of temperature both in the bipolar plate and MEA. It further results in an increase of the saturation vapor pressure. Therefore, the water activity decreases, as well as the membrane water content at cathode side at 300 s, as the green line shows at the bottom of the

Rated power Operating Region

For a well-designed fuel cell hybrid powertrain with a proper energy management strategy, the operating range of the fuel cell system would be either in its "max efficiency" zone, or in its "rated power" zone. Thus, a coupling analysis is also performed for the operating point in "rated power" region (i.e. for the studied NEXA stack, it corresponds to a fuel cell current of 42 A), where the fuel cell system power is close to its rated value. The simulated dynamic responses of the MIMO system outputs after step changes are shown in the figures 2.17 -2.20.

It can be seen from figures 2.17 -2.20 that, for the operating point in the "rated power" region, the four variation ranges of membrane water content are reduced compare to that of figures 2.13 -2.16 as shown in the previous section. That is because when the fuel cell system operates at "max power" region, the dynamic behavior of membrane water is changed and lead to an insensitive . Based on the presented analyses of step responses, detailed analyses of coupling effects between inputs and outputs are further given in the following section. where represents the control inputs , and respectively. represents the system outputs , and , respectively.

In order to give a clear implications of interactions among control loops, the RGA analysis rules are expressed based on the values of :

1) , a change of produce an opposite direction effect on when other control loops are closed, the sign of open-loop gain between and is changed. In this case, the selection of control pairing -should be avoided;

2) , a change of can not affect and therefore should not be used to control . In this case, there is no interaction between control pairing -and other control loops (i.e., coupling does not exist); since the open-loop gain between and is not affected by the interaction from the other control loops. In this case, there is no coupling between variable pairing -and other control loops;

5)

, the effect of other control loops enhance the open-loop gain between and , the larger is above the unity, the greater will be this effect. In this case, variable pairing -is recommended when is not very large, a multivariable control design [9] could be used to achieve optimized control performance; 6) , a value of greater than 3 indicates that the system is difficult to control due to strong coupling effects from other control loops, and is also sensitive to input uncertainty (e.g., cause by neglected actuator dynamics) [10].

In this case, a robust decoupling control strategy [4] could be used to achieve optimized control performance. , the formula of this four inputs four outputs system can be described as follows:

It should be noted that, the transfer function is highly non-linear. In order to linearize the non-linear system equation 2.8, two common operating points of fuel cell system in typical fuel cell hybrid powertrain are selected for the analysis. The first operating point is set to between 1/3 and 1/2 rated power point, which corresponds to the high efficiency operating region of fuel cell system.

On the other hand, the magnitude range difference of each physical parameter is very large. For example, the variation range of control input is from 297.15 K to 303.15 K, the variation range of control input is from 2814 Pa to 4700 Pa. In order to analyze the coupling effect of different physical parameters in the same RGA, their numerical variation ranges are normalized firstly prior to the RGA analysis. In addition to the high efficiency operating region, the second operating point is set to around its rated power point, which corresponds to the high power operating region of fuel cell system. In this case, table 2.2 presents the calculated corresponding RGA values between different input/output variables. From table 2.2, the similar coupling effects can be observed, and the same results of coupling analysis can be obtained. Thus it can be concluded that, there exist the similar coupling effects between thermal and fluidic domains in the proposed MIMO system for both fuel cell system typical operation points.

Coupling Analysis of Sub-System

From the analyses in the previous section, it could be possible to separate the proposed MIMO system into two control sub-systems by minimizing control coupling effects 61 between each other. The linearized formula of first possible sub-system can be described as follows: existence of coupling effect. In this case, a robust decoupling control strategy [4] is recommended to achieve the optimal control objectives for output variables. For example, in [4], a single neuron adaptive PID feedback controller is proposed, which can eliminate the interference derived from the coupling effect.

On the other hand, as the previously analyzed, the variable could almost be independently controlled by input variable , which can be considered as the second one input one output control sub-system: Furthermore, by conducting the same analysis for the rated power operating point, similar sub-systems separation can also be obtained for the proposed fuel cell MIMO system.

Based on the above analysis, it can be concluded that, different fuel cell operational parameters coupling can be observed among different physical domains during fuel cell operation. When a coupling effect exists, special attention should be paid for control system design, in order to achieve an optimized control performance for fuel cell systems.

Conclusion

Although the coupling of dynamic phenomena in a fuel cell has an important influence on the control design of fuel cell system, very little information has yet been published in the literature on the analysis of interaction between fuel cell dynamic variables in different physical domains. The first part of this chapter investigates in particular the coupling effect between the dynamic behaviors during fuel cell transient operation, based on the proposed improved dynamic multi-physical proton exchange membrane fuel cell model, which can be found in the previous works [11] [12].

Chapter III: PEMFC 2-D model numeric solver development for Real-Time Control Implementation

As mentioned before, the quantitative analyses in Chapter II can provide us useful information for modeling assumptions (simplification), which can be used to simplify the 2-D modeling. For example, since the tortuosity ι and porosity ε can be considered as insensitive parameters for a lower fuel cell stack current, the diffusion mass transport can be reasonably neglected. In this case, there is no total pressure gradient in the gas diffusion layer. During the development of fuel cell model with these simplifications, the computational complexity can be effectively reduced while maintaining a high accuracy.

The main objective of this paper is to present a novel 2-D PEMFC modeling approach based on a numerical solver tridiagonal matrix algorithm (TDMA) for real-time control implementation. The proposed PEMFC model covers multi-physical domains in both fluidic and electrochemical. The major contributions of this paper can be summarized as follows:

1) A novel non-uniform control volumes mesh grid is defined in fluidic domain modeling based on channel geometric patterns, in order to thoroughly describe the gas flow characteristics by taking the fuel cell flow field geometric form into consideration. In addition, the differential equations of reactant gas convection and diffusion phenomena are transformed into tridiagonal systems of equations, which can be efficiently solved by tridiagonal matrix algorithm.

2) An implicit iterative solver has been developed to solve spatial physical quantity distributions for electrochemical domain. This original iterative solver algorithm is composed by three interactive computational loops and uses a robust convergence method for real-time computation.

3 drawback of these models is the computational complexity of mathematical operations.

For example, as mentioned before, the commonly used complex computational fluid dynamic (CFD) models are not suitable for real-time simulation, since the computational burdens are too heavy.

Tridiagonal Matrix Algorithm for Real-Time Simulation

The tridiagonal matrix algorithm is an efficient numerical solver, which can be applied iteratively for solving multi-dimensional problems [7]- [9]. The tridiagonal matrix algorithm uses a special form of Gaussian elimination, in order to solve a set of equations for tridiagonal system in a backward substitution. Therefore, it can reduce the computational time and memory usage considerably [9].

Tridiagonal Matrix Algorithm

A The TDMA is an efficient numerical solution for solving tridiagonal matrices. The TDMA consists of two steps: a forward elimination procedure and a backward substitution procedure. A tridiagonal system can be written as the following equations in the tridiagonal matrix form:

where denoted the non-zero inputs of tridiagonal system.

For the first line:

Divide both sides of the equation 3.

Modeling Hypotheses

In order to be able to simulate the model in real-time while keeping the accurate spatial non-homogeneous effect prediction and model accuracy, some assumptions are used when modeling the PEMFC stack.

1) The two-phase flow of water is ignored, but the liquid water saturation, vapor transportation and pressure gradient is considered in the proposed model;

2) The gas flow in the channel and diffusion through the GDL is considered in steady state, since the transient time constant of fluid is relatively short (microsecond or millisecond);

3) The diffusion mass transport due to total pressure gradient is neglected;

4) The activation losses of the anode side are neglected due to the fast electrode kinetic of hydrogen gas in PEMFC;

5) The Ohm losses are only determined by resistance of Nafion membrane, the layer contact resistance, electrode resistance are negligible;

6) The reactants are considered as ideal gases;

7) The geometric characteristics of each layer remain unchanged.

Solve Reactant Gas Convection 2-D Model Using Tridiagonal Matrix Algorithm

The control volume partitions of cathode parallel serpentine channel are shown in the figure 3.2 [10]. Based on the gas supply channel geometry, the gas supply channel can be divided into two sections, denoted as "straight volume" and "bend volume", as illustrated in figure 3.2. From figure 3.2, the thin dotted line channel control volume is in the straight section, denoted as "straight volumes" of gas channel. In contrast, the bold dotted line channel control volume is in the curved section, denoted as "bend volumes" of gas channel. The direction of gas convection flow is marked with arrows. where is the gas mass flow at gas supply channel inlet Based on the assumption 2), the fluid behaviors considered in strady state. Thus, the discretized form of equation 1.13 can be expressed as the following form:

. .

Channel inlet

1 q 2 q 0 q 1 i q  i q 1 2 i i+1 i+2 i+3 i+4 2i+2 2i+3 1 i q  2 i q  3 i q  4 i q  2 2 i q  2 1 i q  2 3 i q  Channel outlet
Thus, the fluid mass flow can be written as:

The equation 3.11 can be written as:

where is the pressure at gas supply channel inlet. Thus, for equations of all the control volume:

the following equations can be obtained:

where is the distance between two adjacent control volumes and . Equation 3.15 can then be written as the tridiagonal matrix:

where and is considered as boundary conditions, since they are known.

Equation 3.16 describes the reactants flow behaviors in a tridiagonal matrix form, which can be directly employed the TDMA. After solving equation 3.16 using TDMA, the pressure distribution due to the mechanical losses suffered in straight channels ( , ) can be obtained. Then the pressure distribution in the serpentine channels can be further obtained based on the equation 1.43.

Solve Reactant Gas Diffusion 2-D Model Using Tridiagonal Matrix Algorithm

In order to clearly show the calculation procedures of TDMA, the "fluid adjacent volumes" are denoted as "F1-F6", the "solid adjacent volumes" are denoted as "S1-S3", as illustrated in figure 3.4 [10]. Take the "solid adjacent volumes S2" as an example. Due to the diffusion phenomena, the reactant gas molar flow (mol/s) entering the "S2" including three parts: diffusion from "F3" to "S2" (x-axis, marked with arrow 1); diffusion from "F4" to "S2" (x-axis, marked with arrow 2); diffusion from "S1" to "S2" (y-axis, marked with arrow 3). The reactant gas molar flow (mol/s) leaving the "S2" including two parts:

diffusion from "S2" to "S3" (y-axis, marked with arrow 4); diffusion from "S2" to catalyst layer (z-axis, marked with arrow 5).

In order to compute all the pressure of "solid adjacent volumes S2, S3, S3...", the above equation can be extended to a set of equations:

Thus, equation 3.22 describes the reactants diffusion behaviors through GDL in a tridiagonal matrix form, which can be directly employed the TDMA.

Implicit Iterative Solver

In the proposed 2-D PEMFC model, the calculations of fluidic gas channel model is based on the non-uniform control volume (black mesh in the figure 1.10), while the calculations of electrode/electrolyte related physical quantities is based on the unified control volume (red mesh in the figure 1.10).

It should be noted that, since the activation loss appears in an implicit form in the Butler-Volmer equation 2.6, an iterative solving method should be developed to calculate [11]. In addition, this iterative algorithm should also solve the current of each segment (current density distribution) and cell potential (fuel cell output voltage), which cannot be calculated explicitly a priori in a 2-D modeling approach. A detailed schematic diagram of the proposed iterative algorithm is presented in the following figures 3.5. It should be noted that, the iterative algorithm used in this thesis is completely independent of commercial platform, and can be easily implemented to any embedded controller of PEMFC systems.

Real-Time Performance Comparison

Comparison with Newton's Method

In order to show the advantages of the proposed 2-D model in terms of computation time, its performance is compared with a recently published research [12] in Table 3.2. The pseudo-2D fuel cell model in [12] uses Newton's method to solve the ordinary differential equations, which describe reactant transport only along the straight channel direction. It can be seen from Table 3.2 that, the proposed full-2D real-time modeling approach can achieve a faster computation speed with a more comprehensive twodimensional consideration and a slightly larger mesh number. It is also important to mention that, the model in [12] uses a commercial software solver (Matlab/Simulink), while the developed iterative solver is completely independent of commercial platform, and can be easily implemented to any embedded controller of PEMFC systems.

Comparison with Gaussian Elimination Method

As mentioned before, as the proposed approach uses a special form of elimination and solves the tridiagonal system equations in backward substitution, the arithmetic complexity of tridiagonal equations system can thus be exponentially reduced compared to the classical Gaussian elimination method. This advantage is particularly evident for a larger mesh number. To show further this advantage, another Gaussian elimination based 2-D real-time modeling approach [13] has been used for performance comparison.

It can be seen from the model computation time of different mesh numbers in Table 3.3 that, the model CPU execution time in [13] is faster than that in this paper for a low mesh number (N<64). However, the CPU execution time in [13] significantly increases with increasing of mesh number. It can be observed that, when the mesh number N≥328, the CPU execution time could exceed 100 ms. And as mentioned in [13], for a large mesh number, the CPU execution time could exceed 1 second, thus making this model not suitable for real-time model-based controller applications. Although the model CPU execution time in this paper is higher than that in [13] for a low mesh number, the execution time quasi-linearly increases with mesh number when the mesh number N is larger than 128. And it can be observed that from Table 3.3, when the mesh number N=600, the execution time is 50.429 ms. Such short CPU execution time further demonstrates the effectiveness of the proposed modeling approach, especially for large mesh numbers. However, the first physical-based method cannot be used in the on-line prognostic applications, since this single-step prognostic method makes predictions for a fixed horizon only. In order to improve the proposed method, iteration-based prediction methods are particularly useful to achieve a better performance in the on-line prognostic applications.

Conclusion

For this purpose, an iteration-based prediction algorithm for performance degradation of PEMFC is proposed in the second part of this chapter. A novel approach using the moving window method is applied, in order to dynamically retrain the models during the forecasting process with new data inputs for iterative data training at next prediction step. In addition, since the proposed hybrid prognostic method is based on a combination of model-based and data-driven prognostic methods, it is able to capture both the fade trend and non-linear features observed in the fuel cell voltage degradation data.

Literature Review

As an advanced concept, the Prognostic and Health Management (PHM) of PEMFC is designed to minimize maintenance costs while increasing operational availability and utilization of PEMFC [1]. As a key process of PHM, the prognostics use a set of monitoring data from actual life cycles of PEMFC system, in order to indicate the future degradation trends as well as its current state of health [2]. Based on the degradation trends, the PEMFC system impending faults and remaining useful life (RUL) can be predicted, the advent of failure can be further forecasted [3].

Generally speaking, there are two main approaches for degradation prediction: modelbased methods [4]- [6] and data-driven methods [7]- [12]. Model-based methods consider the real physical aging phenomena during the PEMFC degradation process.

Although model-based methods need large computations and complex physical model, it can predict not only the system degradation trend (fuel cell output voltage decay over time), but also the information about the internal physical parameters during the degradation process. However, as mentioned before, the complexity of physical degradation phenomena makes it very difficult to build a reliable physical model of degradation. Several PEMFC model-based aging prediction methods can be found in the literature.

Jouin et al. [4] present a prognostics framework to provide RUL predictions based on three voltage degradation empirical models: a linear, an exponential and a log-linear model. However, the empirical voltage degradation models are too simplified, the fuel cell operating conditions, such as operating current and temperature are also not taken into account.

Bressel et al. [5] have proposed an empirical model of degradation. Based on Extended Kalman Filter (EKF) method, the proposed prognostic algorithm is able to estimate the PEMFC state of health and to predict its RUL under a variable load profile. However, the electrochemical kinetics in the proposed aging model is described by Tafel equation.

Under a variable load profile, this aging model can lead to a large error of electrochemical activation loss, especially for small current values.

Chen et al. [6] have developed a PEMFC lifetime quick evaluation method by taking the various changes of the operating condition into consideration. This method can achieve the RUL prediction in real-time applications. However, the fuel cell performance degradation is described by a linear aging model, which cannot truly reflect the nonlinear aging trend. Indeed, the fuel cell is a multi-physical system (electrochemical, fluidic, and thermal). Therefore, it is necessary to take into account aging process in different physical domains.

In the proposed model-based method, a multi-physical aging model is developed to consider major internal physical aging phenomena of fuel cells, including Ohmic losses, reaction activity losses, and reactants mass transfer losses. In addition, the Butler-Volmer equation is used in the proposed aging model to accurately calculate the electrochemical activation loss for all possible values of PEMFC stack current.

On the other hand, the data-driven methods use pattern recognition, statistical or machine learning techniques, such as Artificial Neural Network (ANN) [7], Support Vector Machine (SVM) [8], and adaptive neuro-fuzzy inference system (ANFIS) [9], to track and predict the system non-linear characteristics. Compared with model-based approach, data-driven method can be easily implemented without specific knowledge of physical degradation.

Marra et al. [10] proposed a neural network estimator of solid oxide fuel cells for its diagnostic. The proposed degradation estimator is trained based on a set of experimental data, in which includes stack current, temperature and reactant gas mass flow rate.

Mao et al. [11] investigated sensor selection algorithms for prognostic of PEMFC. With the identified optimal sensor, an adaptive neuro-fuzzy inference system (ANFIS) is used to predict the performance of fuel cell system.

Ibrahim et al. [12] proposed a data-driven approach for PEMFC prognostic based on Wavelet Transform technique. The prediction process consists of decomposition and reconstruction. This approach is able to predict the future power and estimate RUL under static and dynamic operating conditions of PEMFC system.

However, a common drawback of these data-driven methods is that the degradation trend and non-linear behaviors cannot be simultaneously captured, thus make them inaccurate for prognostic. The proposed hybrid prognostic method is able to capture both the fade trend and non-linear features observed in the fuel cell voltage degradation data. In addition, a novel approach using the moving window technique is applied in order to iteratively update the prediction process when the newly measured data become available. This iteration-based prediction method is particularly useful to achieve a more accurate prediction when the initial training samples are limited.

Degradation Prediction Based on Multiphysical Aging Model with Particle Filter Approach

In this section, a multi-physical aging model of a PEMFC is developed at first. The presented aging model considers different physical aging phenomena including fuel cell ohmic losses, reaction activity losses and reactants mass transfer losses. The proposed aging model is then initialized by fitting the PEMFC polarization curve at the beginning of lifetime.

During the prediction process, the aging dataset is then divided into two parts: learning and prediction phases. The particle filter framework is used to study the degradation characteristics and update the aging parameters during the learning phase. The suitable fitting curve functions are then selected to satisfy the degradation trends of trained aging parameters, and used to further extrapolate the future values of aging parameters in the prediction phase. By using these extrapolated aging parameters, the prediction results are thus obtained from the proposed aging model. Three experimental validations with different aging testing profiles have been performed. The results demonstrate the robustness and advantages of the proposed prediction method.

Description of Aging phenomena

In this section, the presented aging model is developed based on a previously developed multi-physical PEMFC model [13] [14]. In this model, each single cell is also divided into seven different element layers, and each layer modeling covers three physical domains: electrical, fluidic and thermal domain. It should be noted that, the time constant of aging process is much higher than the other physical dynamics in PEMFC system (thermal dynamic, membrane water dynamic, etc.). Thus, all the physical dynamics are removed in the proposed aging model compared to the original multi- e  e  Cations poisoning with the presence of Fe 3+ , Cu 2+ etc. Mechanical degradation resulted from the excessive stresses due to humidity/temperature cycling operating profile or even platinum particle precipitation in certain region of the membrane (Xo plane).

Input: Hydrogen

Input: Oxygen H  H  H  H  2 H 2 H O
Carbon support in catalyst layer:

Membrane:

Bipolar plates:

Corrosion of the carbon support in the catalyst layer or the gas diffusion layer due to gas starvation especially in the start-stops.

Lumped parameter of ohmic losses:

, ( )

ohmic aging R t

Aging Parameter Membrane Resistance

The first aging phenomenon considered in the proposed aging model is described by ohmic losses of fuel cell, including layer contact resistance, electrode resistance and membrane resistance [16]. More explanations and causes of fuel cell ohmic losses are estimation accuracy when non-linear battery models are used, since it focus on approximating the solution for non-linear and non-Gaussian system [19]- [22].

In this section, the state space representation of the proposed model is presented at first.

Then, the Bayesian estimation method is introduced, followed by a detailed description of particle filter algorithm.

State Space Model for Aging

In order to give a clear structure of the proposed aging model, the proposed non-linear model can be written in the form of state-space representation as follows:

where represents the state variable of system at time k, the observed variables of system, the input variables of system, the non-linear system observation equation can be derived using the formula of the electrical and fluidic domain modeling equations presented in chapter I. and the process and observation noises respectively. The system state includes three aging parameters:

The schematic diagram of non-linear function is shown in the figure

It can be seen from the figure 4.2 that, three aging coefficients are considered as the system state variables. The initial values of the three aging parameters , , and are obtained by fitting the polarization curve of PEMFC at the beginning of lifetime (described hereafter in section 4.2.3.1). Other modeling parameters (geometrical values or physical properties of the Ballard Nexa 1.2kW fuel cell stack) can be obtained in the previous work [13] [14]. Thus, the relationship between and can then be expressed by:

( ) k  ( ) k  ohmic V 2 ,
It should be noted that the system state variables follow a first order Markov process, and each observation are statistically independent.

Particle Filter Framework

From the previous analysis of Bayesian estimation based on the sequential importance sampling, a detailed particle filter algorithm can be described as follows:

1) Initialization: the set of particle samples is obtained from the initial distribution at time t=0:

, where is the number of particles.

2) Importance sampling: the proposal distribution is generated by:

3) Weights calculation: the state estimation is optimal with the assumption of . It means that probability density function only depends on and , then each sample particle can be obtained. Each sample particles weights can be calculated from equation 4.12 and 4.13:

The importance weight at time t=0 is set to .

4) Normalize the weights: the set of new particle sample is obtained by normalized weights:

5) Re-sampling: the effective sample size is defined by: , the resampling step is performed when is smaller than the given threshold .

6) State prediction: calculate the state by the equation:

Prediction Methodology Implement

In this section, the proposed model-based PEMFC performance degradation prediction method is given in detail. The initialization of three aging parameters is performed based on the initial polarization curve of PEMFC. Then, the proposed prediction method is divided into two stages: learning and prediction phases. During the learning phase, the particle filter estimator studies the system non-linear behavior, and at the same time updates the corresponding aging parameters. Based on these updated aging parameters, suitable fitting curve functions are then selected to extrapolate the values of aging parameters in the prediction phase. At last, by using these extrapolated aging parameters, the proposed aging model predicts the PEMFC degradation voltage in prediction phase.

Initialization of Aging Parameters

The performance of the proposed prediction method depends strongly on the initial estimation of aging parameters. In order to accurately estimate the initial values of aging parameters, these values are determined based on the reliable physical equations and experimental calibrations.

Initialization of Aging Parameter Ohmic Resistance

At the beginning of PEMFC lifetime, the ohmic losses are considered only in the membrane, the electrode resistance and layer contact resistance are negligible. Thus, the initial value of ohmic resistance can be calculated by the equation:

where is the thickness of the membrane (m), is the membrane section surface (m 2 ), is the membrane resistivity and is highly dependent on water content in polymer membrane, where . More detailed information about ohmic losses and water content equations can be found in [13]. With the membrane properties data given in [13], the initial value of first lumped aging parameter is around 2.1e-03 (Ω).

Initialization of Aging Parameter Exchange Current Density

The initial value of exchange current density is calculated by the following equation:

where is the oxygen pressure (atm) at the interface of cathode catalyst layer, and are two empirical parameters need to be identified through fuel cell experimental tests, is the oxygen activation energy at electrode platinum interface (J/mol). With , , the initial value of second lumped aging parameter is around 0.1154 (A/m 2 ).

Initialization of Aging Parameter Diffusion Coefficient

3) Initialization of : The initial value of diffusion coefficient between species and is calculated by the equation [13]:

where is the total pressure of species (atm), is the critical temperature of species (K), is the critical pressure of species (atm), and is the molar mass of species (kg/mol). The coefficients and depend on whether one of the species is a polar gas or not and are determined accordingly, is the porosity of the GDL and is the GDL tortuosity. With the characteristic data of gas diffusion layer given in [13], the initial value of third lumped aging parameter is around 6.7e-06 (m 2 /s).

Genetic Algorithm

In order to demonstrate the modeling accuracy and identify the initial values of aging parameters, the proposed multi-physical aging model is used to fit the experimental measurement of polarization curve at the beginning of PEMFC lifetime. As a commonly used strategy, the Genetic Algorithm (GA) is particularly suitable for such multiparametric and non-linear system. The purpose of GA is to find optimal solution for objective function, defined as follow:

where are the estimated parameter values at time in the generation , is the measured output value at time , and is the predicted output value.

The flow diagram of aging parameters initialization using GA algorithm is shown in the figure 4.3. The main idea of GA is to generate a population of solutions and then to improve it using techniques of natural evolution, such as inheritance, mutation, selection, and crossover method. This generational process is repeated until an appropriate solution can satisfy (optimize) the objective function. More detailed content about GA can be found in [23] [24].

As shown in the figure 4.3, the inputs of GA include fuel cell stack geometry parameters and physical properties, as well as operating condition of auxiliary components. In addition, the numeric variation range of above each aging parameter has been set to ±5% of their calculated initial values for the GA tuning. The particle filter in the previous learning phase estimates and updates the state variables at every sampling step, thus the non-linear behaviors of voltage degradation before the learning time h can be fully captured. In order to well represent the captured aging behavior, the multiply exponential function is commonly used in the literature [25]- [28] to describe the degradation trend. Furthermore, linear, exponential, power, and Fourier series functions [29] can also well express the PEMFC aging process under constant current and temperature condition. Since the degradation trends of different aging parameters represent different physical aging progresses and degradation mechanisms, the above-mentioned fitting curve functions should be properly selected for each aging parameter extrapolation, in order to further achieve accurate prediction voltages in the prediction phase.

Experimental Results and Discussions

In this section, a 400 hours experimental degradation voltage dataset is fully investigated to perform the fitting curve function selections at first. Then, based on the selected fitting functions, the proposed prediction method is then applied to 3 other 400 hours PEMFC experimental degradation tests under different operation conditions. The prediction performance evaluation of the proposed approach with shorter duration of learning phase is also shown and discussed at last. Taking the first aging parameter as an example. In the learning phase, the

exchange current density. (c). diffusion coefficient. (operating current 44A at 40°C)

The prediction results show that, for the third experimental aging test, the Ohmic resistance increases by more than 76%, the exchange current density and diffusion coefficient decrease by 15% and 40%, respectively. It worth also to mentioned that, when the initial training samples are limited, iterationbased prediction strategies [30] are particularly useful to achieve a better performance of the proposed prediction method. The purpose of using an iteration-based prediction strategy is to provide more accurate short-term prediction results and to dynamically update the prediction result with the newest dataset (newly measured data) for iterative data training at next prediction step.

Conclusion

In this section, a novel approach for PEMFC output voltage prediction based on a multiphysical aging model with particle filter and data extrapolation approach was proposed (please refer to previous work [13] [14]). The proposed multi-physical aging model fully considers the three most important aging phenomena during PEMFC operation: ohmic losses, reaction activity losses, and reactants mass transfer losses.

The proposed prediction method is divided into two stages: learning phase and prediction phase. During the learning phase, the particle filter is applied to study the non-linear aging behavior and update the proposed aging parameters. Then, different fitting curve functions are used to represent the trend of aging parameters in the learning phase, and further extrapolate the future values of aging parameters in the prediction phase. At last, by using the extrapolated aging parameters, the proposed aging model predicts fuel cell voltages in the prediction phase.

The experimental validations show that, in order to fully and accurately represent the non-linear trends of aging parameters and further achieve a better performance of proposed prediction method, a relatively large amount of data should be learned by particle filter in the learning phase, and suitable fitting curve functions should be used to extrapolate the values of aging parameters in the prediction phase. It is also important to note that, even when the acquired training data is extremely limited, by knowing a prior the parameter-function pairs, the proposed method is still able to provide acceptable results, especially at immediate short prediction time range.

Initialization of Model-Based Approach

The detailed steps of initialization for model-based diagnostic approach are described as follows:

1 2) Then, the training time point is set to and thus the prediction phase is long. It means that, 0 h to h, the state variable is trained and updated using particle filter, the remaining degradation data from to 400 hours is use to evaluate the prediction results.

3) The initialization of state variables is performed by fitting the experimental degradation data into the initial model. From the known proposal distribution, the particle filter uses a set of random sample particles based on the Monte Carlo sequence, in order to estimate the posterior probability density of non-linear system. More detailed content about particle filter can be found in the previous section. where the sample step .

Data-Driven Prognostic Method

The proposed model-based prognostic method is essentially based on an exponential regression model, it cannot accurately catch the non-linear and uncertain behaviors during the aging process. Compared with model-based prognostic method, the datadriven approach can well describe the local non-linear characteristics of degradation voltage. Especially for the short-range and medium-range prediction time, the datadriven prognostic method can give a better representation of uncertainties in the degradation process. In this thesis, a data-driven prognostic method has been chosen as the second prediction method.

As a data-driven approach, the artificial neural networks can be effectively used to implicitly indicate the complex non-linear characterization between system inputs and outputs. The aging process of PEMFC can be considered as a non-linear autoregressive time series model. The actual fuel cell voltage (model output) can be considered as a variable determined by an unknown non-linear process from the previous voltage values.

In this case, the non-linear autoregressive neural network (NARNN) model can then be used to track the non-linear characteristics of PEMFC degradation.

In this section, a NARNN model for fuel cell degradation is presented at first. An optimal training strategy is also proposed in order to achieve good prediction performance of the developed NARNN model.

Non-linear Autoregressive Neural Network Model

A NARNN model is suitable to describe the non-linear dynamic in a wide variety of system and have been extensively implemented in various applications [31]- [33]. The general formulation of NARNN model can be expressed as follow:

where

Hybrid Prognostic Approaches

As an integration process of multiple prognostic results from different methods, a suitable fusion approach is critical to demonstrate their advantages [36]- [42]. It increases the process reliability and robustness by combining the complementary information from different prognostic methods in intelligent ways. Therefore, a hybrid prognostic approach can be expected to provide more accurate and robust prediction results compared with single model based method.

In order to combine the advantages of different prediction methods in the proposed hybrid prognostic approach, a good understanding about the specific characteristics of each prediction method is important. The model-based method can effectively forecast the aging trend for a long-range prediction time, however it cannot accurately describe the local non-linear characteristics of aging. In contrast, the data-driven method can provide a good non-linear representation of uncertainties in the degradation process for the short-range and medium-range prediction time. But its long-term prediction result is not accurate due to predicted data fluctuation. Therefore, in order to effectively taking the advantages of each method into account, there are two critical issues that need to be considered when using a hybrid prognostic approach: 1) the range selection of training and prediction time, and 2) weight factors determination of each corresponding method.

As an efficient computational strategy, the moving window method is considered to be a good solution for the above-mentioned issues. In this section, a moving window method is presented firstly in order to take different prediction time range into account.

A weight factors adjustment method is then proposed using the results from the moving window method. The proposed method can dynamically adjust the weights vector at each step of the moving window.

Moving Window Method

The primary purpose of using a moving window method is to update and add the newest dataset for iterative data training, and provide dynamic weight factors to further improve the prediction accuracy. In order to perform iterative training during the forecasting process, the model input dataset are updated continuously by moving window approach, as shown in the figure 4 By applying only the data-driven prognostic method with the same moving window coefficient, the prediction results from data-driven method are illustrated in figure 4 where the weight coefficients is a array, the bias coefficients is a vector. The second hidden layer has four neurons, thus the output of the second hidden layer is a vector, which can be expressed as follow:

where the weight coefficients is a array, the bias coefficients is a vector. The output of the NARNN is:

where the weight coefficients is a array, the bias coefficients is a vector.

Similar to the previous model-based prediction method, for each prediction step, the It should be noted that, the time interval between 2 data points is around 1.5 hours in this test case. Thus, a moving window coefficient 100 to 700 corresponds to a forecast time of 150 to 1050 hours respectively. It can be concluded from the figure 4.29 that, by using the proposed hybrid prognostic method, the accuracy improvement decreases when the forecast time length increases. This is mainly because the errors of data-driven method increase significantly with the forecast time. In contrast, it can be seen that the errors of model-based method are not visibly affected by the forecast time.

Thus, in order to maintain simultaneously high precision and relatively long forecast time for on-line prognostic applications, a compromise between prediction accuracy and forecasting time length should be appropriately defined based on the specific application requirement. From a decision making point of view, the prognostic of fuel cell system must be performed with a long enough forecasting time to allow reaction.

For a stationary fuel cell application [12], around 160 hours long prediction (one week) is considered appropriate. It should be noted that, the proposed hybrid prognostic method has a rapid execution time for all the training, evaluation and prediction parts.

With different moving window coefficient from 100 to 700, the forecasting process is capable for on-line prognostic with a time step level in the order of several seconds.

Conclusion

In this section, an innovative approach based on a combination of different prognostic methods has been proposed for the application of PEMFC performance degradation prediction (please refer to previous work [43]). The proposed hybrid prognostic method considers two basic prediction approaches: model-based approach and data-driven approach. Model-based approach can be efficiently used to forecast the fuel cell aging trend in a long-range forecast time. Data-driven approach can accurately describe the local non-linear characteristics of degraded voltage for the short-range and medium-range prediction time.

By combining the advantages of the above two common prediction methods (i.e.

model-based approach and data-driven approach), the proposed method in this thesis can simultaneously and accurately capture the long-term fuel cell voltage degradation trend, as well as the non-linear voltage variation characteristics. In addition, a novel approach using the moving window technique is applied in order to iteratively update the parameters during the prediction process. Furthermore, the prediction performance evaluation of the proposed hybrid prognostic approach with different moving window length is further shown and discussed.

Three experimental validations with three different PEMFC stacks and different aging test

PREFACE 0. 1

 1 Introduction to Fuel Cell: What is Fuel Cell

Figure 0. 1 Structure

 1 Figure 0.1 Structure

X.

  Wang et al. [6] present a three-dimensional non-isothermal PEMFC model. This model uses two different water transport equation to describe the water two-phase transportation during the fuel cell operation. Based on the developed model, a parameter sensitivity analysis is performed to show effects of different parameters on the fuel cell polarization curve. S. Um et al. [7] develop a multi-dimensional transient model for PMFC. The proposed model simultaneously considers the electrochemical kinetics and hydrodynamics. In order to predict not only the experimental polarization curves, but also detailed distribution of electrochemical and fluidic features, the conservation equations are numerically solved using finite volume based computational fluid dynamics (CFD) approach.

Figure 1 . 1

 11 Figure 1.1 PEMFC stack level and the single cell level

Figure 1 . 2

 12 Figure 1.2 Dynamic water flow behaviors in fuel cell membrane.

the anode and cathode sides have different geometric patterns, as shown in the figure 1 Figure 1 . 3

 113 figure 1.4.

Figure 1 . 4

 14 Figure 1.4 Actual geometry form of gas channel of NEXA PEMFC: the left one is cathode air channel, the right one is anode hydrogen gas channel

Figure 1 . 5

 15 Figure 1.5 Structure of a single cell of fuel cell stack.

. 6 .Figure 1 . 6

 616 Figure 1.6 Three-parallel serpentine channel for cathode air supply; and single serpentine channel for anode hydrogen supply.

Figure 1 . 7

 17 Figure 1.7 Schematic diagram of pipeline width , pipeline thickness and spacer length between two neighboring duct in the channel curved U-bends.

4 )Figure 1 . 8

 418 Figure 1.8 Gas diffusion phenomena in the GDL.

  domain and electrochemical domain, the control volume mesh grid definition of gas channels in the proposed 2-D PEMFC model is considered in a non-uniform manner. It means that, the geometry form of each control volume follows the channel geometric patterns. The 2-D channel model can be then implemented by defined control volumes with the physical equations presented in the previous section. The control volume 2-D mesh grid at both cathode and anode sides are depicted in the following figure 1.9.

Figure 1 . 9

 19 Figure 1.9 Control volumes 2-D mesh grid at cathode/anode sides.

Figure 1 .

 1 Figure 1.10 Control volumes 2-D mesh grid at cathode/anode sides, and uniform segments for electrochemical calculation

proposed 2 -

 2 D model with a commercial Ballard NEXA 1.2 kW PEMFC stack. With the developed 2-D model, the spatial physical quantity information can be accurately observed and analyzed by taking the multiple spatial dimensions into consideration. These spatial results are very useful to help to quantitatively analyze the coupling effects in different physical domains, and study the influences of model parameters on the fuel cell spatial performance.

AFigure 2 . 1

 21 Figure 2.1 Experiment platform: 1.2 kW Ballard NEXA 47 cells PEMFC stack, Ballard control system board, and measurements

Figure 2 . 2

 22 Figure 2.2 Experimental validation of polarization curves.

Figure 2 . 3

 23 Figure 2.3 Experimental validation under different stack current profile: (a) Long current step. (b) Short current step.

Figure 2 . 4

 24 Figure 2.4 Experimental validation of stack voltage with different step current profile: (a) Long current step. (b) Short current step.

Figure 2 . 5

 25 Figure 2.5 Experimental validation of stack temperature with different step current profile: (a) Long current step. (b) Short current step.

Figure 2 . 6

 26 Figure 2.6 Dynamic responses after a current step: (a) Water content. (b) Ohmic resistance.

Figure 2 . 7

 27 Figure 2.7 Dynamic responses after current step variation: (a) Water content. (b) Ohmic resistance.Figure 2.7 (a) shows that the dynamic responses to quick current step variation. From

Figure 2 .

 2 figure 2.7 (a), the anode water content decreases sharply at a high current step, due to the high water electro-osmotic drag flow. It can be concluded from the figure that, the dynamic membrane water content is directly related to the stack current variations. The corresponding dynamic behaviors of Ohmic resistance are shown in figure 2.7 (b).

  voltage. A comparison of gas pressure drop in the channel between the developed flow field model and the model using straight channels assumption (as in the most of literature) is given hereafter, in order to highlight the importance of channel geometric form on the pressure modeling accuracy.Comparison of simulation results of pressure drop in the channels of two modeling approaches with different step current profile are shown in figure2.8 and figure 2.9.

Figure 2 . 8

 28 Figure 2.8 Comparison of pressure drop in the channels of two modeling approaches with long step current: (a) Cathode channels. (b) Anode channels.

Figure 2 . 9

 29 Figure 2.9 Comparison of pressure drop in the channels of two modeling approaches with short step current: (a) Cathode channel. (b) Anode channels.

2. 3 . 2

 32 Results and Discussions The proposed 2-D model thoroughly considers the geometric form of parallel serpentine flow field in the fuel cell. The oxygen pressure distribution in the cathode three-parallel serpentine channel, and hydrogen pressure distribution in the anode single-parallel serpentine channel are shown in figure 2.10 (non-uniform mesh grid distribution of each side).

Figure 2 .

 2 Figure 2.10 Gas pressure distribution in the parallel serpentine channels: (a) cathode side. (b) anode side.

Figure 2 .

 2 Figure 2.11 Oxygen pressure distribution on the surface of GDL under different oxygen stoichiometry: (a) τ=1.68. (b) τ=4.94

Figure 2 .

 2 Figure 2.12 Current density distribution on the electrodes under different oxygen stoichiometry: (a) τ=1.68. (b) τ=4.94

  For the fluid domain, The dynamic of the membrane water content is another significant dynamic phenomenon in the fuel cell stack, which has similar transient time constant compared to temperature dynamic in the thermal domain. It is thus necessary to analyze the variable coupling between these two fuel cell operational parameters.

Figure 2 .

 2 13 -2.16 show the simulated dynamic responses of controllable outputs after step changes of different inputs.

Figure 2 .

 2 Figure 2.13 The MIMO system outputs response after step change of .

Figure 2 .

 2 Figure 2.14 The MIMO system outputs response after step change of .

Figure 2 . 13 .

 213 The uniformity coefficient, presented by the purple line on the same figure, slightly increases at 300 s. It means that the distribution of current density on the electrode is less homogeneous. Similarly, when the coolant inlet temperature decreases at 1650 s, it results in an opposite effect on the four fuel cell variables.

Figure 2 .

 2 Figure 2.14 shows the dynamic response of the MIMO system outputs after step changes of the channel inlet air temperature. An increase in the bipolar plate temperature and MEA temperature can be observed when the channel inlet air temperature is increased at 300 s. With the MEA temperature increase, the membrane water content decrease, as the green line shows. At the same time the purple line indicates a slightly less homogenous distribution of current density with the increase of gas channel inlet temperature.

Figure 2 .

 2 Figure2.15 shows the dynamic response of the fuel cell MIMO system outputs after step changes of the gas supply channel inlet water vapor pressure. It can be seen from the figure 2.15 that membrane water content increases due to an increase of inlet water vapor pressure at 300 s. For the thermal domain, a step change of the gas supply channel inlet water vapor pressure has no significant effect on the temperature. For the electrical domain, an increase of channel inlet water vapor pressure makes the current density distribution slightly more dispersed.

Figure 2 .

 2 Figure2.16 shows the dynamic response of the MIMO system outputs after step changes of the air inlet mass flow rate. An increase of air inlet mass flow rate has weak effects on the thermal domain. Due to the increase of air flow rate, the water removal rate increases, and further results in a decrease of the membrane water content. For the electrical domain, an increase of air inlet mass flow rate makes the current density distribution more homogenous. That is because the oxygen is supplied under a higher flow ratio. For the same oxygen consumption rate (same current), the oxygen pressure through the cathode gas channel is thus more uniformly distributed, and further leads to more homogenous distribution of current density on the fuel cell electrode.

Figure 2 .

 2 Figure 2.15 The MIMO system outputs response after step change of .

Figure 2 .

 2 Figure 2.16 The MIMO system outputs response after step change of .

Figure 2 .

 2 Figure 2.17 The MIMO system outputs response after step change of .

Figure 2 .

 2 Figure 2.18 The MIMO system outputs response after step change of .

Figure 2 .

 2 Figure 2.19 The MIMO system outputs response after step change of .

3 )

 3 , there exist two effects: the effect of open-loop gain between and , and the effect from other control loops. The effect of open-loop gain from control pairing -increases with the value of . When, the two effects are equivalent, there is a strongest coupling effect between the two control variable pairs for a two inputs two outputs system. When , the openloop gain of control pairing -becomes the dominant effect. In this case, the most suitable control variable pairing should be -, and optimized control results could be achieved by using a decoupling control method[1] [8]; 4) , the relative gain indicates the most suitable control variable pairing -,

) 3 . 1

 31 The practical feasibilities of the proposed 2-D model in advanced real-time control of PEMFC systems have been experimentally demonstrated in a RT-LAB real-time simulator. The computing technologies presented in this paper are original for real-time PEMFC model and completely independent of commercial platform. It can be easily implemented to any embedded controller of PEMFC systems. Literature Review Different from the common modeling approach, a real-time oriented fuel cell model has more restrictions: the accuracy and computational efficiency of a real-time fuel cell model are both crucial for model based control process [1]. A sophisticated fuel cell model can provide comprehensive physical quantities for model-based control design and optimization. While the high performance computation of a fuel cell model ensures the model-based controller can be efficiently implemented in real-time applications with a low cost of computations. Several real-time control-oriented PEMFC models have been previously presented in the literature [2]-[4]. Jung et al. [2] present a PEMFC real-time model, which considers both the electrical and thermal dynamics. In order to reduce the computational burden, three optimization strategies are used: minimizing algebraic calculation, model separation and reducing the layer structure. Gao et al. [3] develop a cell-level dynamic PEMFC model, which covers electrochemical, fluidic and thermal domains. A top-down design approach is used to provide an efficient PEMFC model structure. By using VHDL-AMS language, the developed model can be used in the hardware-in-the-loop application. Colclasure et al. [4] describes a physical-based transient solid oxide fuel cell (SOFC) model, which considers the coupled interactions of multiple physics. In order to facilitate the real-time control applications, linear model reduction method is used. However, their models remain in 1-D. During the model-based control process, the spatial physical quantity distribution is neglected, such as gas pressure gradient in the channel, or current density distribution on the surface of electrode. Many PEMFC 2-D models have been previously proposed in the literature [5] [6]. However, a common

  2 by : with In order to forward eliminate , the equation 3.3 multiplied by and minus the second row of equation 3.1: Divide both sides of the equation 3.4 by : with Similarly, in order to forward eliminate , the equation 3.5 multiplied by and minus the second row of equation 3.1: Similar forward elimination procedure is repeated until the row, the tridiagonal matrix in equation 3.1 can be transformed into a upper triangular matrix: with The inputs of tridiagonal system can then be solved by backward substitution: Thus, the general solution for tridiagonal system equation 3.1 can be written as the following equations: with A detailed schematic diagram of the TDMA is presented in the following figure 3.1. It is worth specific mention that, by using this special form of Gaussian elimination, and solving the tridiagonal matrices in such backward substitution, compared with straightforward Gaussian elimination , the arithmetic complexity of TDMA exponential decays to operations . It means that, if total number of control volumes is 32 (the elements number of inputs in tridiagonal system), the model computation speed can be 1024 times faster compared with Gaussian elimination. Such fast solving speed allows to significantly reduce the computational time and memory usage.

ForFigure 3 . 1

 31 Figure 3.1 Schematic diagram of the TDMA.

Figure 3 . 2

 32 Figure 3.2 Control volume partitions of cathode gas supply channel based on the geometric form.

Figure 3 . 3 2 .

 332 Figure 3.3 Control volume mesh grid definition of the "gas pipeline A" marked in the figure 3.2.

A 2 -

 2 D multi-physical real-time model of proton exchange membrane fuel cell has been presented in this chapter. An innovative non-uniform 2-D mesh grid method is proposed for real-time simulation performance consideration. In order to efficiently calculate the physical quantities distribution in 2-D plane, an iterative solver is also developed in the model. The proposed original iterative solver algorithm is composed by three interactive computational loops and uses a robust convergence method for real-time computation. The proposed 2-D model has been tested in a RT-LAB real-time simulator and has been experimentally verified using a Ballard NEXA 1.2kW PEMFC stack. The experimental results demonstrate the practical feasibilities of the proposed 2-D model for advanced realtime control of PEMFC systems with a control loop time level on the order of milliseconds [10] [14]. Such short execution time of the proposed 2-D model makes control decisions and actions based on the predicted local phenomena and spatial distribution physical variables inside the fuel cells.Chapter IV: Degradation Prediction of Proton Exchange Membrane Fuel Cell StackTwo different degradation prediction methods for proton exchange membrane fuel cell (PEMFC) performance are proposed in this chapter. In the first part, a novel degradation prediction approach for proton exchange membrane fuel cell (PEMFC) performance is proposed based on a multi-physical aging model with particle filter and extrapolation approach. The proposed multi-physical aging model considers major internal physical aging phenomena of fuel cells, including fuel cell ohmic losses, reaction activity losses, and reactants mass transfer losses. By knowing a prior the parameter-function pairs, even when the acquired training data is extremely limited, this physical-based method is still able to provide acceptable results at immediate short prediction time range.

  physical model. In order to model the aging process over PEMFC operating lifetime, the proposed model uses time-variant modeling coefficients to describe three most important aging phenomena in different physical domains [15] [16], including fuel cell ohmic losses, reaction activity losses in electrical domain, and reactants mass transfer losses in fluidic domain, as shown in the figure 4.1.

2 H

 2 Loss of hydrophobicity of polytetrafluoroethylene (PTFE) treatment coating on the surface of GDL during long term operation impedes the water removal and reactant transfer in the pore structure. Physical deformation due to carbon losses or unbalanced force distribution from either the membrane or bipolar plates.Ostwald ripening or the dissolution of the nano platinum particles on surface of the catalyst layer due to potential cycling Irreversible adsorption on surface of catalyst layer with the presence of the pollutants from the air, reactants or other degraded components Loss of carbon support in the catalyst layer during gas starvation.Formation of resistive layer on surface of bipolar plates due to long term corrosionChemical degradation due to radical attack from hydroxyl (OH•) or hydro peroxyl (•OOH).

Figure 4 . 1

 41 Figure 4.1 Degradation mechanisms considered in the individual lumped aging parameters.

Figure 4 . 2

 42 Figure 4.2 Schematic diagram of proposed aging model.

Figure 4 . 3

 43 Figure 4.3 Flow diagram of GA for fitting the initial polarization curve.

Figure 4 .

 4 Figure 4.4 shows the polarization curve measured at PEMFC beginning of lifetime, and the GA fitting result of the proposed model. It can be seen from the figure 4.4 that, the identified aging model output fits very well to the experimental measurement over the entire current range by using the GA approach (the fitness value is 97.27%). The identified initial values of three aging parameters are, respectively, =2.081e-03 (Ω), =0.1149 (A/m 2 ), and =6.839e-06 (m 2 /s). These three initial values are further used for particle filter to study the degradation characteristics in the learning phase of the proposed aging prediction method.

Figure 4 . 4 Figure 4 . 5 Figure 4 . 6

 444546 Figure 4.4 Initial measurements of polarization curve and GA fitting results.

  particle filter studies the voltages degradation behaviors and updates the parameter values, as shown in the figure 4.7 (blue points in the learning phase). The upward tendency of in the learning phase can be well described by abovementioned linear, exponential, power, and Fourier series functions. These four functions also generate different extrapolated values of parameter in the prediction phase, as shown in the figure 4.7 (different colors points in the prediction phase).

Figure 4 . 7

 47 Figure 4.7 Fitting curves of ohmic resistance parameter using different functions.

Figure 4 . 8

 48 Figure 4.8 Prediction voltage results using extrapolated parameter of .

Figure 4 . 9

 49 Figure 4.9 Extrapolated result of aging parameters: (a). ohmic resistance. (b). exchange current density. (c). diffusion coefficient. (operating current 12A at 30°C).

Figure 4 .

 4 Figure 4.10 Prediction result (current 12A at 30°C, learning time is set to 250 h).

Figure 4 .

 4 12 clearly points out again that, 4.2.4.2.3 Third Experiment: PEMFC Operation under 44A at 40°C In the third experimental aging test, the PEMFC stack runs for 400 hours under a current of 44A at 40°C. By applying the proposed prediction method, the extrapolation results of the aging parameter and voltage degradation prediction result are shown in figure 4.13 and figure 4.14, respectively.

Figure 4 .

 4 Figure 4.13 Extrapolated result of aging parameters: (a). ohmic resistance. (b).

Figure 4 .

 4 14 clearly shows again that, the prediction curve can describe the trend of the non-linear experimental data set. The RMSE and MAPE of prediction result are 0.3168 and 0.0083, respectively.

)

  The degradation data and corresponding experimental conditions are presented at first. A PEMFC stack (commercial Ballard NEXA 1.2 kW PEM fuel cell stack) is used in a 400 hours experimental degradation test. It can be seen from the figure 4.19 that, during the 400 hours degradation testing (operating condition: current 35A at temperature 40°C), the fuel cell degradation voltage data is measured and recorded every 10 minutes. It should be noted that, before applying the proposed algorithm to the PEMFC stacks aging dataset, a data post-processing is performed to eliminate the unreliable measurement data points. Thus, the available degradation data consists of 91 data points (by downsampling the original measured data), which can be seen in the figure 4.19.

Figure 4 .

 4 Figure 4.19 Fuel cell stack experimental aging voltage under current 36A at 35 °C.

4. 3 . 1 . 3

 313 Prediction Phase of Model-Based Approach In the previous training phase, the particle filter estimates and updates the state variables at every sampling step. Thus, the voltage degradation behaviors of training phase have been fully captured. As mentioned before, in order to well represent the captured aging behavior in the prediction phase, all of the trained state variables during the training phase should be taken into account in the empirical model equation 4.25. In this case, the most efficient way of integrating all the state variables is to use their average value during the training phase (from 0 to hour). The average value of state prediction results is expressed as: where is the number of degradation data (blue points), which are used for particle filter learning during the training phase, as shown in the figure 4.19. Thus, the prediction voltage results in the prediction phase are calculated by the output equation of state space equation 4.25:

1  2  3 Figure 4 . 20 Figure 4 .

 1234204 Figure 4.20 Structure of the NARNN model for PEMFC aging prediction.

Figure 4 .

 4 Figure 4.24 Prediction results and identified parameters of model-based prediction method.

Figure 4 .

 4 Figure 4.24 (b)-(e) shows the identified model parameters at each prediction step, it can be seen that during forecasting process, the degradation model is dynamically retrained and the particle filter iteratively update the model parameters. As mentioned in section 4.3.1.3, the average value of each identified parameter at each step is further used to predict the model output in the prediction part, as shown in equation 4.28. These average values of model parameters are given in the table 4.3.

. 25 .

 25 figure 4.26.

Figure 4 .

 4 Figure 4.25 Prediction results of data-driven prediction method.

Figure 4 .

 4 Figure 4.26 Identified NARNN model for data-driven method.

  NARNN model is dynamically retrained and all the weight coefficients array and bias coefficients vector are iteratively updated. By applying now the proposed hybrid prognostic method, the comparison of prediction results is shown in figure4.27.

Figure 4 .

 4 Figure 4.27 Comparisons of three prediction method results.

  measured and post-processed under operation condition with a current density of and a temperature of . This experiment aims to analyse the influence of different moving window coefficient (length) on the performance of proposed algorithm for the long-term degradation prediction. By applying the proposed hybrid prognostic method, the coefficient of moving windows is respectively set to 100, 200, 300, 500 and 700. The obtained prediction results are shown in the following figure 4. 28.

Figure 4 .

 4 Figure 4.28 Prediction results of proposed hybrid prognostic method using different moving window coefficient (length).

Figure 4 .

 4 Figure 4.29 The accuracy improvement with different moving window coefficient .

  sampling the original measured data, the 1200 hours degradation voltage and power of PEMFC are shown in the figure 4.30.

Figure 4 .

 4 Figure 4.30 The third case 1200 hours aging test under a static current density of : fuel cell degradation voltage and power.

Figure 4 .

 4 Figure 4.31 Comparisons of three prediction method (voltage and absolute error).

Figure 4 .

 4 Figure 4.32 Comparisons of three prediction method results (power and absolute error).

  profiles have been performed to verify the accuracy and effectiveness of the proposed fuel cell performance degradation prediction method. The presented results can help engineers to appropriately choose the moving window length, in order to achieve simultaneously high prediction precision and relatively long forecast time for on-line prognostic, for example, the fuel cell hybrid electric vehicle (FCHEV) [] [].

  .3) are applied to the developed model. The model simulation results are shown in figure 2.6.Figure 2.6 (a) shows that the dynamic responses of water content at both cathode, anode sides and in the membrane after the long

	current step. From figure 2.6 (a), the cathode water content	increases at 100 s,
	because more water is produced at cathode side at high current. In contrast, the
	anode water content	decreases at 100 s, due to a higher water electro-osmotic
	drag flow (higher current) from anode to cathode. The green line represents thus
	membrane average water content , which depend on the boundary conditions at
	both cathode and anode sides. The dynamic of membrane Ohmic resistance, which is

highly depend on membrane water content based on equation 1.3, is also shown in figure 2.6 (b). It can be seen from the figure that, the transient time of membrane water content, thus the transient of Ohmic resistance, can last about 500 s.

Table 1 .1 Mesh Grid Independence Check for 2D Model Output

 1 The proposed 2-D model output voltage differences for different control volume numbers of uniform segments are shown in the table 1.1. It can be concluded from the table 1.1 that, when the mesh number is larger than 256, the model outputs difference is less than 0.01%, thus it can be consider that the model outputs are no longer affected by the change of mesh grid size.

	Mesh number	Model outputs difference (%) (to previous mesh number results)
	32	-
	48	0.314
	64	0.132
	128	0.015
	256	0.007

  It can also be seen that, temperature in thermal domain and membrane water content in fluidic domain have very similar dynamic time constant. It is thus necessary to analyze the variable coupling between these two fuel cell operational parameters.

	It can be seen clearly from equation 2.4 that	is determined by stack current . Thus,
	the system time constant value of water content	can be estimated to 108.75 s
	when stack current is 27 A, means that the response time of membrane water content
	variation can last about 435 s. This value is in agreement with the results shown in figure
	2.6.
	2.4.2 Analyses of Step Responses
	2.4.2.1 High Efficiency Operating Region
			1.34
	becomes zero in steady-state, and the steady-state membrane water content value
	(s) can be obtained:
	Then the water content dynamic in transient state can be obtained in a general first
	order dynamic form:
	And the expression of time constant of water content in membrane	(s) can then
	be obtained by the following equation:

  2.4.3.2 Coupling Analysis of the Proposed MIMO system (High Efficiency Operating Region)

	When the MEA temperature	, the bipolar plate temperature	, the
	membrane water content	, and the uniformity coefficient of current density
	distribution	on	the	electrode	are	viewed	as	outputs

Table 2 .1 Relative Gain Array of System among Thermal, Fluidic and Electrical Domain in Fuel Cell High Efficiency Operation Range

 2

	-2.6948	3.6892	0.0056	0.0000
	3.5376	-2.5309	-0.0083	0.0016
	0.1859	-0.1899	1.0369	-0.0329
	-0.0287	0.0316	-0.0343	1.0314
	Then by applying RGA method to the normalized parameters of MIMO system
	equation 2.8 in fuel cell high efficiency operation range, table 2.1 presents the
				of	,	,

calculated corresponding steady-state RGA values between different input/output variables. It can be seen from table 2.1 that, the RGA elements absolute value

Table 2 .2 Relative Gain Array in Fuel Cell High (Rated) Power Operation Range

 2

	-3.0861	4.0767	0.0094	0.0000
	4.0384	-3.0351	0.0135	-0.0167
	0.1227	-0.0966	0.9817	-0.0078
	-0.0750	0.0551	-0.0046	1.0245

Table 2 .3 Relative Gain Array of First Sub-System

 2

	-2.7098	3.7040	0.0057
	3.5451	-2.5362	-0.0088
	0.1647	-0.1678	1.0031

Table 2 .

 2

3 presents the corresponding RGA values for new three inputs three outputs subsystem in fuel cell high efficiency operation range. As analyzed in the previous section, large RGA elements indicate that the control design is critical and challenging due to

Table 3 .1 2-D Real-Time Model Benchmark Results (2.5GHz CPU)

 3

	Mesh number	Model CPU execution time	Model time step used	CPU occupation
	32	7.429 ms	10 ms	73.3%
	48	12.308 ms	15 ms	82.1%
	64	17.731 ms	20 ms	88.7%
	128	43.061 ms	50 ms	86.12%
	256	101.218 ms	120 ms	83.35%
	600	257.193 ms	300 ms	85.73%

Table 3 .2 2-D Real-time performance comparison with Newton's method [12].

 3

	Simulation	Simulation duration Simulation time cost	Ratio
	The proposed model with	800 s	297 s	37.12%
	N=32 (2.5 GHZ CPU)	1200 s	473s	39.41%
	Model in [12] with	800 s	431 s	53.87%
	N=25 (2.4 GHZ CPU)	1200 s	611 s	50.91%

Table 3 .3 2-D Real-time performance comparison with Gaussian method [13].

 3

	Mesh number	104	200	328
	Execution time	6.3 ms	39.5 ms	162 ms

Table 4 .1 Operation Conditions of NEXA Fuel Cell Stacks

 4 A at 40 °C (the dataset is shown in the figure4.5). The learning phase time is set to 250 h. This experiment aims to find appropriate fitting function for each aging parameter extrapolation. By applying the proposed prediction method, the particle filter studies the non-linear aging behaviors and updates the aging parameters in the learning phase. Then, an important step of the proposed method is to select appropriate fitting functions, which are used to represent the trends of aging parameters in the learning phase, and further extrapolate the future values of aging parameters in the

	Test in section	1 st test in	2 nd test in	3 rd test in
		4.1.4.1	section 4.1.4.2.1	section 4.1.4.2.2	section 4.1.4.2.3
	Stack type	Ballard NEXA 1.2 kW commercial PEM fuel cell stack
	Operation mode		Dead-end mode
	Air supply		Air blower + filter
	Cooling		Air fan cooled
	Active area			150
	Fuel supply		99.99% dry H2 @1.2 bar
	Operating hours			400 hours
	Air stoichiometry	2.0	4.2	2.2	2.0
	Stack temperature	40	30	35	40
	Current density	0.24	0.08	0.20	0.30
	4.2.4.1 Fitting Function Selection and Extrapolation Method
	In the first experimental aging test, the Ballard NEXA 1.2 kW PEMFC stack runs for 400
	hours under 36

  4.3.1.2 Training Phase of Model-Based ApproachIn the training phase, the purpose of particle filter is to estimate recursively the state variable based on the initialization model and experimental aging data in the training phase. It should be noted that, during the particle filter based estimating process,

	the variance selection of process noise	need to be chosen carefully. A large variance
	of	ensures a more extensive sampling distribution, since it provides sample particles
	with a great diversity. On the other hand, a smaller variance of	allows a sufficiently
	fast convergence. There already exists the observation noise in the experimental aging
	data, since the fuel cell degradation voltage data are measured directly by sensors.

  .22.Third case study: 600 W stationary PEM fuel cell stack 1200 hours aging test under static current load.4.3.4.1 First Case Study: Comparison of the Proposed Hybrid Model and Single Model MethodsIn the first case, a Ballard NEXA 1.2 kW commercial PEM fuel cell stack runs for 400 hours aging test under the working conditions of current 35A at temperature 40°C (the aging dataset is presented in figure4.19). This experiment aims to clearly show the accuracy and robustness of the proposed hybrid method. The detailed operation conditions of NEXA fuel cell stack are listed in table 4.2.

Table 4 .2 Operation conditions of the first case study: Ballard NEXA fuel cell stack

 4 By applying only the model-based method, with a moving window coefficient of 15, the prediction results and identified parameters are illustrated in figure4.24.

	First case study

Table 4 .3 Average values of identified parameters used in the prediction part at each step (model-based)

 4

	1-th	29.7299	-6.2637e-04	0.7831	3.9282e-03
	2-th	29.7294	-4.4274e-04	0.7827	3.9289e-03
	3-th	29.7297	-4.9277e-04	0.7835	3.7804e-03
	4-th	29.7305	-4.8636e-04	0.7833	4.1389e-03

Table 4 .4 RMSE of Prediction Results

 4 Second Case Study: Performance Evaluation with Different Moving Window Length In this second case, an 8 kW PM 200 PEMFC stack (96 cells) is operated for 10000 hours aging test. The operation conditions of PM 200 fuel cell stack are listed in table 4.5.

		Model-based method Data-driven method Hybrid prognostic method
	RMSE	0.3362	0.3665	0.2352
	MAE (V)	0.7949	0.8915	0.6302
	4.3.4.2

Table 4 .5 Operation conditions of the second case study: PM 200 fuel cell stack

 4 In this study, the total available aging data contains around 6500 data points, which are

	Second case study

  4.3.4.3 Third Case Study: Comparison of the Proposed Hybrid Method with Other MethodsIn the third case, a 600 W PEM fuel cell stack (5 cells) is operated for 1200 hours aging test under stationary condition with a static current density of . This third case experiment aims to show the prediction performance comparison between the proposed hybrid method and other methods. The operation conditions of this case experiment are listed in table 4.6.

Table 4 .6 Operation conditions of the third case study: 600W fuel cell stack

 4

	Third case study

Table 4 .8 Comparison of RMSEs between the proposed methods and ARIMA [12]

 4

		Model-based	Data-driven	Hybrid	ARIMA
		method	method	method	algorithm [12]
	Week 3	0.6184	0.5832	0.4778	0.4
	Week 4	0.6458	0.6377	0.5145	0.6
	Week 5	0.8026	0.8571	0.6633	0.7
	Week 6	0.8027	0.7523	0.6265	1.1
	Total four weeks	0.7198	0.7075	0.5705	0.7

ACKNOWLEDGMENTS PUBLICATION LIST

ACKNOWLEDGMENTS

This work was done while the author was with the FEMTO-ST (UMR CNRS 6174), Energy Department, Univ. Bourgogne Franche-Comte, UTBM and with the laboratory FCLAB in Belfort.

Relative Gain Array

In order to quantitatively analyze the interactions among multiple control loops, the relative gain array (RGA) method in control theory are carried out in this study. RGA method is a theory developed by Bristol in 1966 [6] [7] to quantify the degree of interaction of input and output variables. For a inputs outputs system, the matrix RGA is given by:

The open-loop gains of the -th controlled outputs are determined based on the response to a change of the -th manipulated inputs, and all other manipulated inputs remain constant: In the steady state, the gas molar flow entering the "S2" is equal to that of leaving:

where is the gas molar flow. Since the gas molar flow is determined by the pressure difference between two control volumes, thus the equation 3.17 can be expressed as:

All these diffusion phenomena in the gas diffusion layer, except the above-mentioned current-driven one , can be described by the modified Fick's diffusion law equation:

where is the gas molar flow rate of specie x, is the gas diffusion area, is the gas diffusion coefficient, is the ideal gas constant, is the temperature of GDL, is the center point distance between two control volumes. Thus, equation 3. 18

2-D modeling of electrochemical domain based on three levels bisection algorithm

Calculate thermodynamic

Voltage based on Eq. (1.1)

Figure 3.5 Schematic diagram of the proposed iterative solver (First level).

As shown in the figure 3.5, the proposed solver consists of three levels bisection algorithm.

The first level solver is the top level algorithm, which is used to compute cell potential The RT-LAB real-time simulator used in this thesis has a real-time processor operating at 2.5 GHz. For a real-time simulation, the model time step should be settled based on model complexity and computation performance.

Experimental Results

The real-time simulation bench mark results of the proposed model under different 2-D mesh number are shown in the table 3.1.

It can be seen from the table 3.1 that, the model execution time increases quasi-linearly with mesh number, since there are more iterations should be solved for a higher number of control volumes. It can be concluded from this benchmark results that, the proposed 2-D multi-physic PEMFC model in this thesis can be effectively used in real-time control implementations with a time step level on the order of milliseconds.

shown in the figure 4.1 (first group). In order to describe the fuel cell ohmic losses during the degradation process, the first lumped aging parameter ohmic resistance (Ω) has been proposed as following:

where the aging coefficient indicates the aging degree of ohmic resistance.

Aging Parameter Exchange Current Density

In order to describe the reaction activity losses during fuel cell operation, the second lumped aging parameter exchange current density (A/m 2 ) has been proposed

with the following form:

where the aging coefficient reflects the aging degree of exchange current density over time. As a parameter who represents electrode kinetics, the exchange current density determines how easily the reaction can occur on the electrodes. The reaction activity losses may result from degradation of active electrode surface area, Nafion in contact with active area, or loss of catalyst material [16]. More explanations and causes of reaction activity losses are shown in the figure 4.1 (second group).

Aging Parameter Gas Diffusion Coefficient

In order to describe the losses of reactants mass transfer during fuel cell aging, the third lumped aging parameter gas diffusion coefficient (m 2 /s) between the species x and y has been proposed under the following form:

where the aging coefficient reflects the aging degree of fuel cell gas diffusion layer material that influences the reactants mass transfer. The losses of reactants mass transfer may result from the corrosion of electrode supporting material, increase of tortuosity due to catalyst particle ripening, or more difficult water removal due to degradation of polymer material [16]. More explanations and causes of reactants mass transfer losses are shown in the figure 4.1 (third group).

Estimation Method

Compared with the Extended Kalman filter (EKF) which only focuses on linear systems and Gaussian noise [17], [18], the particle filter (PF), which use the Monte Carlo sequence for solving integration problem in Bayesian estimation, has significantly better

Bayesian Estimation

In the framework of Bayesian estimation, the state variable in the state space of aging model (equation 4.1-4.3) is estimated based on the observation sequence . Then, the optimal estimation of can be written in the form of the conditional expectation:

The main idea of Bayesian importance sampling is that, a set of random sample particles, which are sampled from the known proposal distribution , are used to approximate the posterior probability density function . Thus, can be rewritten as:

Then the Summation form of equation 4.7 can be expressed by: where are the particles sampling from , are the normalized form of importance weights , which can be expressed by: Bayesian importance sampling is an effective method using Monte Carlo sequence.

However, each step of importance weights calculation depends on all the previous observations, thus its computation is increasing with time. In order to avoid this deficiency, the sequential importance sampling is proposed. In this case, the importance weights are calculated recursively, and the proposal distribution can be rewritten as:

From equation 4.10, the importance weights at time -1 can be described as follows:

From equation 4.10 and 4.11, the importance weights at time can be described as follows:

the proposed prediction method shows a precise conformity with the validation data.

The RMSE and MAPE of prediction result are 0.1865 and 0.0052, respectively. Thus, it can be concluded from the above experimental validation that, although the acquired aging data is extremely restricted, based on the prior known parameterfunction pairs, the proposed method has the potential to provide fairly acceptable prediction results, especially at immediate short prediction time range.

Degradation Prediction Using a Moving Window Based Hybrid Prognostic Approach

The above physical-based prediction method cannot be used in the on-line prognostic applications, since this single-step prognostic method makes prediction for a fixed horizon. In order to improve the proposed method, an innovative robust prediction algorithm for PEMFC performance degradation is proposed based on a combination of model-based and data-driven prognostic method. In the proposed hybrid method, a novel approach using the moving window method is applied, in order to 1) train the developed models; 2) update the weight factors of each method and 3) further fuse the predicted results iteratively. In the proposed approach, both model-based and datadriven methods are simultaneously used to achieve a better accuracy.

In order to verify the proposed method, three experimental validations with different aging testing profiles have been performed. The results demonstrate that the proposed hybrid prognostic approach can achieve a higher accuracy than conventional prediction methods. In addition, in order to find the satisfactory trade-off between the prediction accuracy and forecast time for optimizing on-line prognostic (for example the dynamic operating conditions in fuel cell hybrid electric vehicles), the performance variation of proposed approach with different moving window length is further shown and discussed.

In this section, a model-based fuel cell voltage degradation prediction model using empirical equation and particle filter approach, and a data-driven prediction model using NARNN are presented, respectively. Then, a hybrid prognostic approach is further given based on combining the fusion approach and moving window techniques.

Model-Based Prediction Method

A model-based approach has been chosen as the first prediction method in this thesis.

In this section, an empirical fuel cell voltage degradation model is presented at first. In order to estimate the state variables in the model, a particle filter based identification algorithm is applied. The corresponding prognostic method is presented and discussed at last.

Based on regression analysis of experimental degradation data, it has been found that an empirical equation form can well describe the fuel cell voltage degradation process [26]- [28]:

where is the fuel cell output voltage at time , the coefficients and are related to the fuel cell internal impedance, the coefficients and are related to the fuel cell aging rate. The discrete-time state space of this voltage degradation model is depicted as follows:

where is the sampling step, the state variables need to be estimated can be described by , the observed variable represents the fuel cell voltage, and is the system stochastic normal distributed noises.

The equation 4.26 is clearly in a non-linear form. In order to accurately estimate the state variable , the particle filter approach is used, since it can effectively solve the Bayesian estimation problem of non-linear system based on the Monte Carlo sequence.

The purpose of using particle filter for non-linear model parameter identification is to capture the trend of fuel cell degradation voltage during the training phase, and further correctly represent the captured aging trend in the future (prediction phase). Therefore, this moving window method is an efficient strategy to dynamically retrain the models during the forecasting process with new data inputs. In addition, by iteratively evaluating the fitting ability of each method at each prediction step, the corresponding weight factor can also be adjusted dynamically. Moreover, by using the proposed moving window method, the forecast time can be easily changed by using different values of .

Weight Factors Calculation

As mentioned before, the data of gray shadow region are used to evaluate the model By assuming that both prediction processes (model-based and data-driven) are conditionally independent at each prediction step, the overall prediction results can be described using weighted average as:

where representing the total number of methods, is the predicting data, is the normalized weight factor, which can be calculated by: The proposed overall prognostic algorithm can then be illustrated in figure 4.23. In figure 4.23, is the moving window coefficient, is the number of predicted data, k is the actual prediction step, the total prediction step K (the threshold step in a specified length of voltage data) determines the retraining cycle time.

Experimental Results and Discussion

In order to experimentally validate the accuracy and robustness of the proposed hybrid prognostic method, experimental degradation datasets from three different types of fuel cell stacks are presented respectively in this section:

First case study: Ballard NEXA 1.2 kW commercial PEM fuel cell stack 400 hours aging test under the working conditions of current 35A at temperature 40°C;

Second case study: PM 200 8.0 kW fuel cell stack 10000 hours aging test working under stationary prime power application;

It can be seen from the figure 4.31 that, the proposed hybrid prognostic method can represent both the fade trend and non-linear features observed in the fuel cell voltage degradation data. The comparison of mean prediction error between the proposed methods and Adaptive Neuro-Fuzzy Inference System (ANFIS) algorithm [11] is shown in table 4.7. It can be seen that the proposed hybrid prognostic method can ensure a higher prediction accuracy than the both single prediction method, and could achieve nearly 21% improvement on the mean prediction error compared with the ANFIS algorithm. In addition, it is important to consider that the algorithm in [11] uses more degradation data (from 0 to 825 hours) for training the ANFIS, while the proposed moving window based prognostic approach requires fewer degradation data (only from 0 to 330 hours) for the first training, and can dynamically retrain the models during the forecasting process using the newly measured data. The above prediction results

and analysis demonstrate the robustness and effectiveness of the proposed hybrid prognostic method.

Table 4.7 Comparison of mean prediction error between the proposed methods and ANFIS algorithm [11]

Model-based method Data-driven method Hybrid prognostic method ANFIS algorithm [11] Mean error 0.0093 0.0089 0.0069 0.0087

In order to further show the prediction performance comparison of each moving window step between the proposed hybrid method and other method, the fuel cell power degradation dataset (shown as red curve in figure 4.30) is used to perform the proposed prediction methods. The forecast time is also set to one week (moving window horizon hours), the comparison of prediction power is shown in the figure 4. 32.

The comparison of Root Mean Square Errors (RMSEs) between the proposed methods and Auto-Regressive Integrated Moving Average (ARIMA) algorithm [12] is shown in table 4.8.

It can be seen from the table 4.8 that, the proposed hybrid method shows again a better prediction performance than single model method. Compared with the ARIMA algorithm in [12], the proposed hybrid method could achieve higher forecasting accuracy from the fourth week to the sixth week, the RMSE result of the proposed hybrid method is nearly 19% less than that of ARIMA algorithm for the total four weeks.

CONCLUSION

Before mass commercialization of proton exchange membrane fuel cell, the research on the design of appropriate control strategies and auxiliaries need to be done for achieving proton exchange membrane fuel cell (PEMFC) optimal working modes. An accurate mathematical PEMFC model can be used to observe the internal variables and state of fuel cell during its operation, and could further greatly help the system control strategy development.

A comprehensive multi-physical dynamic model for PEMFC is developed in chapter I. However, the proposed physical-based method is a single-step prognostic method, which cannot be directly used in the on-line prognostic applications. In order to achieve on-line prognostic, it is important to develop an iteration-based prediction strategy, which allows the prediction result to be dynamically updated with the newest dataset (newly measured data) for iterative data training at next prediction step.

For this purpose, in the second part of chapter IV, a moving window based prediction method is developed, in order to dynamically retrain the models during the forecasting process with new data inputs. In addition, this approach is based on a combination of model-based and data-driven prognostic methods. By combining the advantages of the above two common prediction methods, the proposed hybrid method can simultaneously and accurately capture the long-term fuel cell voltage degradation trend, as well as the non-linear voltage variation characteristics. Furthermore, the prediction performance evaluation of the proposed hybrid prognostic approach with different moving window length is further shown and discussed. Three experimental validations with three different PEMFC stacks and different aging test profiles have been performed to verify the accuracy and effectiveness of the proposed hybrid method. The presented results can help engineers to appropriately choose the moving window length, in order to achieve simultaneously high prediction precision and relatively long forecast time for on-line prognostic applications.