PEMFC system simulation in MATLAB-Simulink® environment

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 8 0 4 5 e8 0 5 2 Available at www.sciencedirect.com journal homepage: www.elsevier.com/locate/he PEMFC system simulation in MATLAB-Simulink
environment Fabio Musio a,*,1, Fausto Tacchi b, Luca Omati b, Paola Gallo Stampino b, Giovanni Dotelli b, Stefano Limonta a, Davide Brivio c, Paolo Grassini c a Genport srl, Via Garcia Lorca 29, 23871 Lomagna (LC), Italy Politecnico di Milano, Dipartimento di Chimica, Materiali e Ingegneria Chimica “G. Natta”, Piazza Leonardo da Vinci 32, 20133 Milano, Italy c SAATI S.p.A., Via Quasimodo, 33 20025 Legnano (MI), Italy b article info abstract Article history: Modelling and simulation of PEM fuel cell stack operation is developed in Simulink
Received 21 May 2010 environment and validated through experimental data. The present work is the starting Received in revised form point for the development of a user friendly and versatile tool aimed at controlling and 6 December 2010 optimizing the operation of a PEMFC stack; in addition, it could be of help in stack and BOP Accepted 18 January 2011 components design, for instance feeding and humidification systems, cooling circuit, Available online 24 March 2011 temperature control logic and electrical interface. The constitutive equations used to model the FC stack operation are the fundamental equations of electrochemistry. First, the Keywords: model is used to describe the behaviour of a single cell under steady-state conditions upon PEMFC varying variables such as temperature, pressure and relative humidity of reactants; then, it Stack modelling is applied to simulate the operation of a stack configuration, including also fluid-dynamics Matlab Simulink aspects, thermal and kinetic behaviour of feed systems. In particular, thermal control Thermal modelling modelling is based on a simplified approach where different heat removal mechanisms are Steady-state conditions accounted for in a separate way. In its present state, the simulation tool so developed allows a feasible investigation of some process variables influence on the FC stack performances. The stack modelling is tested against benchmark results obtained from a 300W 20-cell air-cooled stack under variable operative conditions. MEAs based on Nafion 112 and Carbon cloth GDLs developed ad hoc are assembled into each cell of the stack. Although the model is quite simple, these preliminary results point out that it may be an adequate tool to set design targets and support further steps of optimization. Copyright ª 2010, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. 1. Introduction Several PEM Fuel Cell models are available in the recent literature: a large number of them investigate the fundamental physical interactions outlining the cell behaviour [1e4]; a lot of works propose the dynamic modelling without a validation by experimental data. Only few authors [2,3,5] reported such validations for a more appropriate physical approach, instead of an empirical model. In the present work it was developed a user friendly and versatile simulation tool able to take into account electrochemical and thermodynamic aspects necessary for describing * Corresponding author. Ing. Fabio Musio, Via Garcia Lorca 29, 23871 Lomagna (LC), Italy. E-mail addresses: [email protected], [email protected] (F. Musio). 1 URL http://www.genport.it. 0360-3199/$ e see front matter Copyright ª 2010, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2011.01.093 8046 i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 8 0 4 5 e8 0 5 2 the single cell behaviour as well as system process variables linked to the cooling and gas feeding circuits. The ultimate goal is to have a simulation tool which may be of help in setting design targets and consequently supporting further steps of optimization in the realization of a portable FC stack device. Modelling has been first applied to describe the operation of a single cell under steady-state conditions and then to simulate the operation a 20-cell stack. The outputs were tested against benchmark results obtained from a 300 W 20cell air-cooled stack under variable operative conditions. This work mainly aims at investigating the control and operation optimization process of medium size FC stacks and has the ambition of developing in the long period a useful tool for the design of FC systems, including Balance of Plant (BOP) components, such as feed and humidification systems, cooling apparatus and temperature controls. A generalized steady state fuel cell model is created and it is tested and compared with experimental data. The electrochemical approach estimates the cell voltage for particular operating conditions, setting parameters such as temperature, pressure and relative humidity. The model of the stack system is based upon a set of fundamental equations accounting for fluid-dynamics, thermal dynamics and kinetic behaviour of feed systems, although here are still in a simplified way, in order to investigate the influence of process variables for design optimization of fuel cell systems. Thermal control modelling is also developed for an air cooling system, trying to separate the various heat removal phenomena. The proposed modelling approach is implemented in Matlab-SIMULINK environment; this approach can be found in only few works [6,7]; this environment is highly flexible and can be easily adjusted for different stack systems and operating conditions. 2. Model development 2.1. General The theoretical cell potential at standard conditions of pressure and temperature is given by Nernst equation when the fuel cell is not connected to the electrical circuit; when the fuel cell has the circuit closed, there are different kinds of voltage losses caused by several phenomena. In a fuel cell there are different kinds of voltage losses, namely activation, concentration and ohmic losses. Thus, the voltage of a fuel cell can be expressed as: Vcell ¼ ENernst hact hohmic hconc (1) where ENernst is the Nernst voltage, hact is the activation overvoltage, hohmic represents the ohmic losses and hconc the concentration losses. 2.2. Thermodynamic potential E¼ DG nF (2) The Nernst equation gives the open circuit cell potential; at standard condition of 25  C and pressure of 1 atm, the theoretical potential can be calculated as: Er Trif ; Prif
¼ (3) The reversible Nernst potential is a function of pressure and temperature [8]: Er T; Prif
2.3. Overvoltages ¼ Er Trif ; Prif
DS T nF Trif
pH2 p0:5 RT O2 þ ln pH2 O nF ! (4) As expressed in (1), there are three kinds of losses (or overvoltages) which make the voltage decrease. 2.3.1. Activation overvoltages In order for a reaction to activate, it is necessary to outstrip the limit of activation energy. The activation losses represent the difference from the equilibrium that is needed to get the reaction started, and this is mainly due to sluggish electrode kinetics; these losses happen at both anode and cathode side. The activation losses can be expressed as: DVact ¼   RT i ln anF i0 (5) where a is the transfer coefficient of the cathode side or anode side, i0 is the effective exchange current density at anode and cathode that represents the rate constant of the chemical reaction. It happens that some small amounts of hydrogen diffuse from anode to cathode whereas some electrons may find a “shortcut” through the membrane. These losses appear insignificant in fuel cell operation, but when the fuel cell is in the open circuit mode, or when there is a very low current density, these losses may have a dramatic effect on cell potential. To simulate this phenomenon, the term iloss is added to the current i in (6) with a fixed value of 3 mA cm 2 [8].So in a fuel cell the activation losses are: DVact ¼ DVact;c þ DVact;a ¼ 2.3.2.     RT i RT i þ ln ln ac F i0;c aa F i0;a (6) Ohmic overvoltages These losses exist because of the resistance in the electrolyte and to the flow of electrons through the electrically conductive fuel cell components, such as bipolar plates and gas diffusion layers (contact resistance). The losses can be expressed with Ohm’s law: DVohm ¼ irohmic The maximum amount of electrical energy generated in a fuel cell corresponds to Gibbs free energy, then the theoretical potential of the fuel cell is: DG 237340 Jmol 1 ¼ 1:229Volts ¼ nF 2  96485 Asmol 1 (7) The following expression for Nafion membrane resistivity is used in the model [1]: i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 8 0 4 5 e8 0 5 2 rohmic h i
T 2 2:5 181:6 1 þ 0:03i þ 0:062 303 i
  ¼ ½l 0:634 3iŠexp 4:18 T T303 (8) where l is the water content of the membrane. A fixed value of resistivity of 5 mU cm was defined to represent contact resistance. 2.3.3. Concentration polarisation occurs when a reactant is rapidly consumed at the electrode by the electrochemical reaction so that concentration gradients are established. The relationship that describes this phenomenon is given by: DVconc ¼   RT iL ,ln iL i nF (9) iL is the limiting current density, that represents the maximum current delivered by the fuel cell: it occurs when the current density becomes so large that the reactant concentration on the surface of the cell falls to zero. The current density averaged over the electrode surface can be expressed as [9]: iL ¼ nFhm " Cin C  out  in ln CCout # approximates the behaviour of the high current field of the polarisation curve. The gain factor K used is 15 and this makes the curve decrease earlier than how it would decrease just through the theoretical formulation. Thus the final formulation for the concentration losses used in the simulation is: hconc ¼ Khconc Concentration overvoltages (10) where n is the number of electrons transferred per mol of reactant consumed, F is Faraday’s constant, hm is the mass transfer coefficient, Cin is the concentration of the reactant at the channel inlet and Cout at the channel outlet [10]. The model was built based upon Eq. (9) multiplied by an empirical factor (called gain factor from now on) that 8047 (11) The value K ¼ 15 was chosen once for all. 2.4. Fuel cell stacks Connecting many fuel cells in series (stack), it is possible to increase the total voltage of the system and, consequently, the power delivered. The stack potential can be found as the sum of individual cell voltages or as a product of average individual cell potential and the number of cells in the stack: Vstack ¼ Ncell X Vi ¼ Vcell ,Ncell (12) i¼1 A pressure drop is determined through each bipolar plate; this pressure drop can be approximated by the equation for incompressible flow in pipes and conduits with sufficient accuracy as long as the pressure drop is less than 30% of the inlet pressure: DP ¼ f L v2 X v2 KL r r þ 2 DH 2 (13) where f is the friction factor, L and DH are the channel length and the hydraulic diameter respectively, both expressed in m, r is the fluid density (kgm 3), v is the average velocity of Fig. 1 e Scheme of the model developed in Simulink environment in order to simulate single cell behaviour under steadystate conditions. 8048 i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 8 0 4 5 e8 0 5 2 the gas in m s turns [10]. 3. 1 and KL is the local resistance in the sharp Single cell modeling The model built in Simulink is obtained from Eq. (1) as it is shown in Fig. 1. For each type of overvoltage there is a subsystem in which the input parameters are inserted and the relative values of voltage are obtained as output according to Eqs. (4), (6), (7) and (11). There is another subsystem, which describes the water uptake of the cell depending on the relative humidity of the inlet gases. The water content of the membrane can be expressed as [11]: l ¼ 0:043 þ 17:18a 39:85a2 þ 36a3 (14) where a is the water activity defined as: a¼ p psat In the ohmic losses region, i.e. between about 0.2 and 0.7 A/cm2, the agreement between experimental and theoretical data is very good (Fig. 2); it should be remarked that this current density range is representative of common steadystate operating conditions. At high current density (over 0.8 A/cm2) the comparison between model and experimental data is less satisfactory because of the influence on losses of porous media microstructure, of cell water uptake and of lack of precise reactant gases distribution at the cathode side. These microstructure features are not properly accounted for in the present modeling. 4. Stack modelling In order to be able to simulate the behaviour of the operating system, the stack model has as a core the single cell model validated before. The model developed to simulate the stack comprises a series of single cells connected as they are in the (15) where p is the water partial pressure and psat is the saturation pressure at a given temperature. The saturation pressure can be expressed, through an empirical formulation, as [11]:
log10 psat ¼ 2:1794 þ 0:02953  T þ 9:1837  10 þ 1:4454  10 7  T3 5  T2 (16)  where T is the temperature expressed in C. The molar fraction of water vapour in the incoming gas stream is simply the ratio of the saturation pressure and to total pressure: xH2 O;in ¼ psat pin (17) The ratio of nitrogen and oxygen in dry air is known to be 79/21, so the inlet oxygen fraction can be found with the following expression: xO2 ;in ¼ 1 xH2 O;in 1 þ ð79=21Þ (18) The model gives as output a steady-state polarisation curve, depending on temperature, inlet flow rate and relative humidity. 3.1. Validation of the single cell model The validation process started from experimental data carried out “ad hoc” in laboratory. The cell for the validation has an active area of 22 cm2, with the following configuration: - Catalyst coated membrane (CCM) Nafion
112 membrane with both electrode Pt/C 0.5/0.5 mg/cm2 - Carbon cloth gas diffusion layer (GDL) treated with PTFE 10% wt and coated with a micro-porous layer (MPL) The bipolar plates have a single channel serpentine to feed the membrane with the gases; both air and hydrogen are humidified by bubbling the gases through water. Fig. 2 e Experimental and simulated polarization curves of the single cell at different cell temperatures: a) 60  C, b) 50  C and c) 40  C. i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 8 0 4 5 e8 0 5 2 stack, a pressure drop block (for the losses into the bipolar plates), thermal management and cooling system blocks, and a feed system block. The feeding system for the reactant gases has a “Z-shape” configuration [8] that permits a more uniform flow distribution; moreover, exhaust oxygen leaves on the stack side where hydrogen is fed and vice versa for the exhaust hydrogen (Fig. 3). So, the inlets and the outlets are on opposite sides of the stack and the flows in inlet and outlet manifolds are parallel to each other. This means also that as for gas feedings the cells are connected in parallel. However, the pressure drop across each cell is not exactly the same depending on the flow distribution which in turn involves a flow network problem consisting of N-1 loops where N is the number of cells in the stack. In any case, the Z-shape configuration of the stack flows enables a much more uniform flow distribution. In Fig. 4, an example of 6-cell stack model is represented: each cell is treated independently and the pressure drop subsystem is inserted at the beginning of the series of cell in order to simulate the pressure drop across the stack for both anode and cathode losses; in fact, a specific pressure drop block is inserted between each couple of cells to properly simulate the pressure drop across them. On the left, there is the cathode feed subsystem of the stack and on the right the anode feed subsystem, while on the bottom there is the thermal management subsystem. 4.1. Thermal management subsystem The performance of a fuel cell is greatly dependant on stack temperature. High temperature generally shifts the polarisation curve upwards, but it may cause membrane dehydration and stack degradation. Generally, in a PEM fuel cell system, it is attempted to keep the stack temperature below 80  C. The system temperature can be related to the heat absorbed by the stack with a first order differential equation: Ct dT ¼ Q_ stack dt (19) where Ct is the thermal capacitance (J  C 1), T is the stack temperature expressed in  C and Q_ stack is the rate of heat absorption by the stack in J s 1. Considering the different methods of dissipation, Eq. (19) can be further modified: Ct dT ¼ Ptot dt Pelec Q_ cool Q_ loss (20) where Ptot is the total power delivered into the stack, Pelec is the power consumption by electrical load, Q_ cool is the heat flow 8049 rate of the cooling system and Q_ loss is the heat flow rate through the stack surface expressed in W. The heat of exhaust gases is taken into account indirectly by considering as a reference for the overall chemical reaction ambient temperature, i.e. 25  C. The total power input to the stack is directly related to the amount of hydrogen consumed; since hydrogen consumption is dependant on the stack current and the number of cells, the following expression can be used in determining the total power input: _ H2 ;used DH ¼ Ptot ¼ m I  ncell DH 2F (21) where DH is the combustion enthalpy of hydrogen _ H2 ;used is the rate of hydrogen cons(285.5 kJ mol 1), and m umption expressed in mol s 1. The electrical power output is given by: Pelec ¼ Vstack  I (22) The rate of removal of heat by forced air could be expressed as a function of two different components: external cooling on the surface of the stack and the heat removal by convective flux into the cooling channels. In fact, in the stack used for the validation, the bipolar plates have a gas side, with the channels used to feed the cell, and a cooling side, with parallel channels for cooling purposes. When two cooling sides are attached, cooling channels inlets appear on the lateral surface of the stack. The external cooling consists of a flux of air generated by a fan that blows directly on the lateral surface of the stack (surface 1, Fig. 5) parallel to the inlets air cooling channels and flows parallel to the top and bottom stack surfaces (surfaces 2 and 3, Fig. 5). End plate surfaces (5 and 6, Fig. 5) are thermally separated from graphite plates by a polyethylene layer, which is considered in the present model as a thermal insulator; so, heat transfer through these surfaces is neglected.The rate of heat removed from a body blown by a transversal flux can be calculated by the following expression: Q_ ¼ hAðT Tamb Þ (23) where h is the coefficient of convection heat transfer  expressed in W m-2 C 1, A is the surface area in m2, T and Tamb are the temperature of the body considered and the ambient temperature respectively, both in  C. For the forced heat transfer into the channels the temperature at the output of the channel can be considered lower than at the inlet, but this value is not easily available; for this reason the following equation was used: Fig. 3 e Z-shape flow configuration for hydrogen or air feeding adopted in the simulation of the 20-cell stack and reproducing the experimental feeding systems. 8050 i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 8 0 4 5 e8 0 5 2 Fig. 4 e Schematic view of the modelling of a 6-cell stack as developed in Simulink environment, highlighting the cathodic (left side) and anodic (right side) feed sub-systems as well as the thermal management (block located below the six cells); the pressure drop block has been depicted only once on the left upper side for brevity; in fact, sub-systems blocks are inserted between each couple of cells. T Q_ ¼ hA  LMTD ¼ hA
Tair;in h T ð ln T ð T Tair;in Þ Tair;out i Tair;out Þ
(24) Tair,out is the temperature of air at the exit of the cooling channels, dependent on the internal temperature of the stack as: Tair;out ¼ KT (25) where K is an estimation of the quality of the cooling system: the higher is the value of K, the better will be the efficiency of the heat transferred. Natural heat transfer is calculated as the sum of the natural convection and the radiation to the surrounding: ðT Q_ ¼ Tamb Þ Rth (26) Ct in Eq. (20) is the product of the mass of the stack and the specific heat of the materials which compose the stack. The mass of the stack was approximated as the weight of the bipolar plates in polymer graphite composite. For this, the thermal capacitance has a value of Ct ¼ 978 J  C 1. For the forced heat transfer a hypothetic fan was used to generate the flux of air; the fan is modelled as a variable flow of air with a velocity range of 0e8 m s 1. 4.2. Fig. 5 e Schematic view of the fan cooling air flows with respect to stack surfaces. Air cooling channels (inlets on surface 1) are parallel to the air flow. Stack model validation The stack used for the validation is composed by a series of 20 fuel cells, built and tested at the “Zentrum fur Brennstoffezellen Technik (ZBT)” of Duisburg (Germany) for Saati S.p.A.. The membranes have an active area of 50 cm2, with the following configuration: 8051 i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 8 0 4 5 e8 0 5 2 80 20 18 70 14 50 12 40 10 8 30 Voltage [V] Stack current [A] 16 60 6 20 4 10 2 52 3 49 4 46 5 43 6 40 7 37 8 34 9 32 0 29 1 26 2 23 3 20 4 17 5 14 6 11 7 88 59 0 30 1 0 Time [s] I [A] Air [l/min] H2 [l/min] Temperature U [V] Fig. 6 e 20-cell stack experimental input (current and cathode flow) and outputs (voltage and anode flow) data used to validate the model. Flow data as well as temperature are in arbitrary units. - Catalyst coated membrane (CCM) Nafion 112 membrane with electrode Pt/C 0.3/0.6 mg/cm2 - Carbon clothes gas diffusion layer (GDL) treated with PTFE 10% wt and with micro-porous layer (MPL) Overvoltage losses were calculated taking into proper account the Pt/C loads, which are different from the one used in single cell validations [10]. Experimental tests on the 20-cell stack were used to validate the system model so far developed. The starting set point of the stack is 12V and 70  C; then, by varying the inlet flow of the cathode, voltage change is monitored; finally it is carried out a polarisation curve at constant temperature. The input data are test current and the inlet flow at the cathode; the inlet flow at the anode depends on the current demand with a fixed stoichiometry set at 2.3. In Fig. 6 relevant input (current) and output (voltage) experimental data are shown; cathode and anode flow as well as temperature data are reported in the same figure in arbitrary units in order to highlight how cathode flow changes affects stack performances. As shown in Fig. 7a, the simulation constantly underestimates the stack voltage by about 1V; for a 20-cell stack this means a difference of 0.05V per cell. In Fig. 7b, it is possible to see some fluctuations of the temperature in experimental data, after 400 s, due to the starting of the polarisation curve test that makes the temperature move away from its equilibrium value. As there was a lack of information about how the cooling system of ZBT works, it was decided to operate the model in such a way that the temperature stays closer as much as possible to the set point, i.e. 70  C. However, the simulated and experimental cooling systems should have been reasonably comparable because the gap between experimental and simulated temperature is always below 4  C. Studying the importance of the various aspects of the cooling system [10], it has been noticed that natural convection and thermal radiation are almost constant being the temperature held constant; they are negligible if compared to the heat generated during the test: this makes clear that a cooling system is absolutely necessary to remove the excess heat generated by the electrochemical reaction. According to the model, the air flowing into the channels, and so cooling the stack in its internal parts, seems to Fig. 7 e Experimental and simulation of the steady-state behaviour of a 20-cell stack operation; a) voltage profile (upper) and b) stack temperature (lower). 8052 i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 6 ( 2 0 1 1 ) 8 0 4 5 e8 0 5 2 contribute to the heat removal by an order of magnitude higher than the air surrounding the stack externally. 5. Conclusions The main goal of the work was to construct a valid and modular model of a PEM fuel cell stack that could be useful for various applications, such as designing and developing single system components in order to improve the configuration for a specific use. At first, a generalized steady-state fuel cell model was generated, and it was run and compared with experimental data. The results of the validation looked satisfactory; the polarisation curves approximated finely the trends of the experimental data, certainly at the most typical operating conditions. In the end, the fuel cell steady state model was used as a basis to develop an analogous steady-state stack model where feeding cooling systems were incorporated. In particular, attention was mainly concentrated on the pressure drop of reactant gases in the flow channels and on the behaviour of the air cooling of the item on test. A stack operative situation was simulated with the same inputs adopted in the experimental tests, in terms of inlet flows, inlet pressures and current demand. No variable parameters are used in the simulation of the stack in order to fit the experimental data. The only parameter introduced to fit simulation with experiments was the gain factor K used in the modelling of the concentration losses of the single cell. Once chosen the best value at the level of single cell model validation, the parameter was no longer varied and any subsequent simulation on the stack was performed without any variation of parameters. Comparing the simulated and the real voltage, it results that the model slightly underestimates the experimental data, probably due to some simplifying assumptions. The Simulink approach here adopted resulted valid to describe a large range of stack operations under steady-state conditions; the modelling was quite straightforward and not so much demanding in terms of computing time, so it may be a useful tool in the design of BOP systems, in particular for what concerns the cooling system. references [1] Mann RF, Amphlett JC, Hooper MAI, Jensen HM, Peppley BA, Roberge PR. Development and application of a generalized electrochemical model for a PEM fuel cell. J Power Sources 2000;86:173e80. [2] Ceraolo M, Miulli C, Ponzio A. Modeling static and dynamic behaviour of PEMFC on the basis of electro chemical description. J Power Sources 2003;113:131e44. [3] Wagner N. Characterization of membrane electrode assemblies in polymer electrolyte fuel cells using AC impedance spectroscopy. J Appl Electrochem 2002;32:859e63. [4] Zawodzinski Jr TA, Lopez C, Jestel R, Valerio J, Gottesfeld S. A comparative study of water uptake by andtransport through ionomeric fuel cell membranes. J Electrochem Soc 1993;140: 1981e5. [5] Amphlett JC, Mann RF, Peppley BA. A model predicting transient responses of proton exchange membrane fuel cells. J Power Sources 1996;61:183e8. [6] Khan MJ, Iqbal MT. Modeling and analysis of electrochemical, thermal and reactant flow dynamics for a PEM fuel cell system. Fuel Cell 2005;4:463e75. [7] Khan MJ, Iqbal MT. Dynamic modelling and simulation of a fuel cell generator. Fuel Cell 2005;1:97e104. [8] Barbir F. PEM fuel cells theory and practice. Oxford UK: Elsevier Academic Press; 2005. [9] Spiegel C. PEM fuel cell modelling and simulation using MATLAB. Oxford UK: Elsevier Academic Press; 2008. [10] Musio F. Tacchi F, Model and simulation of PEM fuel cells systems, Master degree thesis, Politecnico di Milano; 2008. [11] Berning T, Djilali N. Three-dimensional computational analysis of transport phenomena in a PEM fuel cella parametric study. J Power Sources 2003;124:440e52.