Component-based modelling of PEM fuel cells with bond graphs

A polymer electrolyte membrane (PEM) fuel cell is a power generation device that trans-forms chemical energy contained within hydrogen and oxygen gases into useful electricity. The performance of a PEMFC unit is governed by three interdependent physical phenomena: heat, mass, and charge transfer. When modelling such a multi-physical system it is ad-vantageous to use an approach capable of representing all the processes in a uniﬁed fashion. This paper presents a component-based model of PEMFCs developed using the bond graph (BG) technique in Modelica language. The basics of the BG method are outlined and a number of relevant publications are reviewed. Model assumptions and necessary equations for each fuel cell component are outlined. The overall model is constructed from a set of bond-graphic blocks within thermal, pneumatic and electrical domains. The model output was compared with the experimental data gathered from a two-cell stack and demonstrated a good accuracy in predicting system behaviour. In the future the designed model will be used for fuel cell reliability studies. the resistances each of the be in conﬁguration than


Introduction
Fuel cells (FCs) are a technology for electricity generation that promises to play an important role in the future energy landscape. These devices convert chemical energy stored within fuel directly into electricity through an electrochemical process. The principal fuel consumed in such conversion is pure hydrogen gas, although certain types of fuel cells can operate on hydrocarbons. Because the by-product of such energy generation is only water and heat, fuel cells are more environmentally friendly. The overall CO 2 emissions generated as a result of fuel cell operation can be minimised or even cut down to zero when the fuel is obtained with renewable electricity via water electrolysis. In addition to that, due to the absence of mechanical parts, fuel cells can be easily scaled to satisfy multiple applications including portable, automotive and stationary power generation [1].
Among the different types of fuel cells, polymer electrolyte membrane (PEM) FCs are the most suitable for automotive applications because the low operating temperatures of 60e80 C provide fast start-up times. However, despite the aforementioned advantages, widespread market adoption of PEMFC technology is hindered by high manufacturing costs and comparatively low durability and reliability [2].
Production and operation of more efficient, durable, and reliable PEM fuel cells remains an active area of research and development. Achieving such goals can be done via several routes. Chemists and materials engineers are working on improved catalysts with better electro-chemical properties, while new types of membranes are being developed and evaluated for enhanced water and ion transport capabilities. Additionally, efficient control and maintenance strategies need to be developed in order to maximise the fuel cell lifetime.
An appropriate control strategy designed to maintain desirable operating conditions ensures maximum efficiency with minimum degradation. Meanwhile, a good maintenance strategy maximises the system reliability and availability with minimum intervention into it's operation. To investigate the effectiveness of options for increased lifetime performance of PEMFC's, it is highly beneficial to create models of the process in question.
The physical interactions occurring within a fuel cell system are inherently interdependent and consist of phenomena in fluidic, thermal and electrical domains. Therefore the main objective of this paper is to present a model of the system that encompasses all of the dominant physical processes in a consistent fashion. In order to achieve this goal, a modelling technique called 'bond graphs' (BG) is employed in this paper to create a component-based model of the system. A brief outline of the core concepts of the BG method are provided in the following section. Next, a number of relevant publications are reviewed and the modelling approach used in this paper is established. The individual components of the fuel cell model are designed using a number of basic default as well as custom-created bond-graphic elements. The simulation results of the resulting model are compared to experimental measurements acquired from a 2-cell Pragma Pro-RD stack.

Bond graph method
The bond graph (BG) method [3,4] is designed to unify the representation of systems involving multiple physical domains. It uses equivalent analogies to link various phenomena. For example, electrical resistance in wires is analogous to pneumatic resistance in a pipe. At the core of the method is the principle of energy conservation, which states that in a closed system, total energy is never destroyed or lost, but instead it is converted from one form to another. Such energy conversion is performed by various components in engineering systems through storage, transfer, dissipation, and ultimately usage of it to perform useful work. The rate at which energy is transferred between components is power. In BGs

Nomenclature
Greek Water vapour 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 4 2 ( 2 0 1 7 ) 2 9 4 0 6 e2 9 4 2 1 this power transfer is denoted as a half-arrow (called power bond, or simply, bond) as shown in Fig. 1, where A and B are two different energy objects. The power flow is characterised by two power variables: effort (e) and flow (f). These variables are selected such that their product equates to the total power transferred: e Â f ¼ Power. The two variables have different meanings depending on the physical domain that the model is focused on. Some of the most commonly used analogies are provided in Table 1.
Bond-graphic elements located at the nodes of a bond graph represent different energy manipulation mechanisms. Sources of effort (Se) and flow (Sf) are active elements and provide input to the system. Such elements controlled by an external signal are called 'modulated' and denoted by a prefix 'm', e.g. mSe. Energy dissipation and storage phenomena are implemented via resistive (R), capacitive (C) or inductive (I) elements. Each of the BG elements carry unique constitutive equations with them as shown in Table 2. Detectors of effort (De) and flow (Df) are shown with a full arrow to emphasise that they do not participate in energy exchange, but rather simply act as sensors and measure corresponding power variables.
Multiple power bonds can meet at one of two junction types: 0-and 1-type, which enforce the laws of energy conservation within the system. In the electrical domain, 0-and 1-junctions correspond to Kirchhoff's current and voltage laws accordingly. Another junction structure called transformers (TF) act as energy transducers converting the transferred power from one physical domain to another. TF elements can only have two bonds connected. The corresponding equations are listed in Table 3.

Pseudo-bond graphs
The power variables listed in Table 1 are applicable to a wide array of physical problems. However, some engineering systems exhibit coupled energy phenomena, such as cooling systems where thermal energy transfer is associated with a moving volume of liquid or gas. Such thermo-fluidic systems may require a different set of power variables in order to represent them efficiently. Karnopp [5] proposed so called "pseudo-bond graphs" which have two separate, but coupled bonds for depiction of thermo-fluidic power transfer. One bond corresponds to pressure (P) and mass flow rate ( _ m), while the second bond transfers temperature (T) and enthalpy flow rate ( _ H) or heat flow rate ( _ Q). Fig. 2 illustrates the most commonly used depiction of pseudo bonds.
Such bonds are called "pseudo" because in contrast to the definition of power variables, the product of either of these pairs is not equal to power. Because of this, it is required to design some ad-hoc elements to link pseudo bond graphs to the rest of the graph that uses the true power variables outlined in Table 1. This may lead to an inconsistent representation of power within a model and care must be taken when employing both pseudo-and true bonds in one graph [6]. Nevertheless, either one of the modelling methods can provide results of equal accuracy. Pseudo bonds are widely used when modelling thermofluid and linear heat transfer problems because they offer a more intuitive way of representing such processes.

Literature review
Bond graph is a modelling technique for representation of multiple physical phenomena in a unified graphical notation.   SfH   Generally speaking bond graphs are a-causal and can be used to create hierarchical models of systems [3,4,7]. The method is moderately popular within the fuel cell modelling community due to its versatility with regards to the application area. Bond graphic models were demonstrated to be used in fuel cell system simulation [8e13], controller optimization [14e17], fault diagnostics [18] and hardware design [19]. Saisset et al. created a dynamic model of PEMFC with the aim of studying fuel cell interactions with the DC/AC converters [8]. The voltage dynamics at the anode and cathode electrodes was represented by two dissociated double-layer capacitors. The model also takes into account thermal dynamics that include heat conduction through the layers of the cell. Heat transfer to the environment through convection and radiation was also considered, although the coefficients were assumed constant. Additionally, the model does not include hydraulic losses in the gas channels and the hydration changes of the membrane. The analysis is also limited to nonhumidified fuel and oxidant. The model performance was only compared to a current-voltage characteristic, so it is not clear how well the model performs in thermal domain.
Mzoughi et al. adopted the model from Ref. [8] and extended it through the addition of diffusion dynamics at the cathode GDL and the catalyst layer [9]. Furthermore, the authors included a more detailed representation of the cooling process by setting the convection coefficients depend on the coolant flow rate. Similarly to [8], only the current-voltage performance of the model was evaluated. The difference with experimental results was explained by omission of changing membrane water content and pressure losses within channels.
In another paper by Mzoughi et al. the authors focused on the control-oriented aspects of fuel cells [14]. The authors employ similar approach as in Ref. [9] but simplify the thermal domain dynamics by neglecting heat transfer between individual layers of the cell. The model also assumes no pressure losses within channels and constant membrane hydration. The solenoid valve that controls hydrogen flow rate and the air compressor were modelled as transfer functions. The authors used iterative approach to tune the coefficients of the controllers and demonstrated the improvement of simulated voltage output when compared to open-loop configuration.
Peraza et al. designed a simple 0D, static, control oriented model of a single PEMFC [10]. The model uses block diagrams to compute the reversible potential of the fuel cell at constant temperature and pressure. The electrical part of the model is modelled using the bond graph method and considers activation, ohmic and concentration losses within the fuel cell. The coefficients for polarization curve simulation were obtained using curve fitting with an average relative error of 1.18%. The resulting model is simple, yet accurate and effective at demonstrating how structural information can be embedded within bond graph models.
Rabih et al. published a bond graph model of the hydraulic part for the gas supply into the PEMFC stack [11]. The model is composed of a sequence of resistive R-and C-elements. The model assumes laminar flow of pure H 2 and O 2 at constant temperature. The pressure losses are determined as a function of gas viscosity and channel geometry. This model was connected to electro-chemical and thermal model from Ref. [8] in order to study the effects of current, temperature and channel dimensions. The authors presented the results of simulations which show that pressure losses are significantly higher at cathode side of the fuel cell and they linearly increase with increasing current. Additionally, the pressure difference nonlinearly decreases with increasing the pipe diameter.
McCain and Stefanopoulou created a distributed model of a GDL consisting of 3 identical sections [12]. Liquid water condensation, evaporation and transport are implemented within the model and its effect on the lumped current density during flooding is considered. The model was experimentally validated using a 24-cell PEMFC stack. The authors emphasise that even though the model is not too complex it still contains 25 states, which is too high for control applications. This is why the authors employ a special model order reduction algorithm in order to eliminate model components that do not bear significant effect on model performance. As a result, the authors reduced the number of states to 19.
Bruun in her PhD thesis developed a 1D distributed model of a tubular SOFC for use in marine power plants [13]. The author first created a lumped parameter section of the cell, which was later used to create a model consisting of six identical sections. Each section incorporates thermal, fluid dynamics and electrochemical phenomena. The model also features auxiliary components such as ejectors, fuel reformers, valves, afterburners and heat exchangers. The author presented the results of simulating several model configurations to compare the performance of a lumped and distributed models, however, no experimental validation was presented.
Hung et al. developed a control oriented bond graph model of a single PEMFC [15]. The bond graph is separated into parts: thermofluid and electrochemical. The flows of binary gas mixtures are controlled by inlet and outlet valves, and chemical potentials of species are calculated within special chemical reactor blocks. The gas diffusion through the GDL and pressure losses within channels were neglected. Thermal dynamics were simplified to include only convective heat transfer between the gases and the body of the cell. Additionally, heat generation due to electrochemical conversion and voltage losses is also ignored. The model was linearised around a set operating point and transformed into state-space form. The authors did not present any experimental validation results.
Vijay et al. developed a control-oriented, 0D, bond graph model of a SOFC [16]. The model includes the representation of binary gas mixtures, thermal dynamics and gas entropy variation due to changes in molar compositions. The control strategy was designed to ensure 80 to 90% fuel utilization and air stoichiometry of more than 8. Secondary control objective was to maintain the cell temperature at acceptable levels was also included. The model was validated with data from literature and after performing simulations the authors concluded that the defined control objective was found to be conflicting and needs to be adjusted. The same team of researchers continued their research into optimal control strategy of SOFC and published another paper suggesting and algorithm to determine the optimal operating conditions for a given current load [17]. The authors also included a model of afterburner and heat exchangers. One drawback of the modelling approach employed in Refs. [16] and [17] is the fact that the model was developed such that chemical potentials of gas species were used not only for characterising electro-chemical reactions, but also for describing the flow rates within the fuel cell. Although this yields a unified power representation, the approach is counter-intuitive.
Yang et al. designed a bond graph model of a PEMFC suitable not only for simulation but also for fault diagnostics [18]. The diagnosis procedure relies on the structural and causal information from the bond graph to detect and identify a set of faults such as flooding or drying of the membrane, loss of catalytic activity and failures related to H 2 supply. The temperature of the FC was assumed constant, as thermal domain modelling was not included. Unfortunately, neither the model nor the diagnosis procedure were verified experimentally in this paper.
Robin et al. presented a 2D model of a PEMFC based on bond graph methodology to study the heterogeneous processes occurring on the surface of the cell layers [19]. The authors created a 2D grid of the MEA with each block representing a control. These control volumes were linked by various transfer mechanisms such as diffusion, conduction and convection. Simulation results were compared against a six-cell stack with an integrated printed circuit board containing and array of sensors for measuring the magnetic field. As a result, the authors were able to obtain the current density and temperature distribution maps of the MEA. The authors were interested in analysing the effects of operating conditions on the occurrence of heterogeneities. The resulting model can be used for fine tuning the channel geometries and material properties of the fuel cell for optimised performance.
From the conducted literature review it can be concluded that there is a great variety of published bond graph models of PEMFCs. However, no publications use bond graphs in an object-oriented modelling environment. Additionally, there seems to be a lack of experimental validation. The model proposed in this paper aims to extend the existing BG models by incorporating all the major physical phenomena occurring in the fuel cell system including the thermal and water management. The model takes into account the chemical composition of the inlet gases and computes the electrolyte ionic resistance based on the amount of water vapour present. The model is designed in an open-source modelling language called Modelica. This provides the opportunity to make full use of the object-oriented features of the language to create individual components of a fuel cell as separate objects. Such a modular architecture provides the possibility of creating 1D models of multi-cell stacks to study the performance of individual cells located in different parts of the stack.

Modelling approach
Bond graphs are well suited for hierarchical modelling, therefore it is possible to develop a set of independent physical components that can be arranged to create the overall system. Such a modular approach allows for models that are easy to read and understand, and can be modified very quickly.
The model hierarchy is as follows: basic bond graphic elements that describe energy storage and transfer mechanisms are at the very base. A set of BG elements describing pneumatic and heat transfer phenomena construct two bipolar plate (BPP) components for the anode and cathode sides of the fuel cell. A third component consisting of BG elements for electrochemical, transport and thermal phenomena is dedicated for representation of the membrane electrode assembly (MEA). The combination of these three components are encapsulated to form the fuel cell component. Additionally, cooling channels and the end plates are implemented as separate components.
The overall model architecture is shown in a word-bond graph in Fig. 3. The diagram shows blocks resembling physical components of the actual PEM fuel cell and bonds connecting them depict power flows between components. Anode/cathode inlet and outlet blocks correspond to mass flow controllers or valves and regulate the flow of matter in or out of the cell. This is shown by bonds labelled with P,T as efforts and _ m; _ H as flows. The source of electric current, denoted as mSf, represents the load demanded from the fuel cell. Electrochemical phenomena within the MEA component calculate the rates of reactants consumption _ m H2 and _ m O2 as well as the rate of product formation _ m H2O . Transport phenomena determine the diffusion flows through the MEA, while thermal effects evaluate heat flows _ Q between the bipolar plates and the membrane.

Model assumptions
The following assumptions are used for the purposes of this model: A.1 Gases obey the ideal gas law due to the fact that the temperatures and pressures are low compared to their critical values.
A.2 The thermal properties of gases and solids are considered constant due to the fact that heat capacities do not change significantly over normal operating temperatures of PEMFCs (60 + C to 80 + C).
A.3 Reactant gases are humidified. This means the anode feed consists of a two-component mixture of pure hydrogen and water vapour, while the cathode feed is a mixture of oxygen, nitrogen, and water vapour.
A.4 The gas mixtures are assumed to be at the same temperature and the energy required to create the homogeneous mixture (entropy of mixing) is neglected.
A.5 Air diffusion through the membrane is ignored. A.6 The water within the fuel cell is assumed to be in vapour phase at all times.

Model implementation
The Modelica modelling language [20] provides suitable tools for BG implementation. It is an open source object-oriented language for convenient representation and simulation of complex physical systems.
Modelica does not have built-in support for bond graphs, but it is possible to create libraries of components that have the necessary properties [21]. The BondLib package, implemented by Cellier [22] is an open-source bond graph library that contains all the essential building blocks for BG models, such as those mentioned in Tables 2 and 3. The library contains partial classes (within Interfaces package) upon which additional single-and multi-port bond-graphic components can be built as required.
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 4 2 ( 2 0 1 7 ) 2 9 4 0 6 e2 9 4 2 1 In order to simplify graphical representation of complex systems, an additional library called MultiBondLib is also available online [23]. The library extends the functionality of BondLib by introducing new types of bonds: multi-bonds (also called vector-bonds). Such multi-bonds can carry an arbitrary number of bonds of the same or related physical domains and represent multi-dimensional power transfer. This makes them ideal for representation of pseudo-bonds, where two coupled bonds are enclosed in a single two-dimensional multi-bond, as shown in Fig. 4. Multi-bonds from Multi-BondLib are well suited to represent sets of pseudo-bonds that carry all the power variables corresponding to a thermofluidic flow of a gas mixture. This means that each multibond carries 4N pseudo-power variables (P,T and _ m; _ H), where N is the number of components in the gas mixture.
The following sections describe the modelling procedure for each physical component of the PEM fuel cell. Firstly, a valve component implemented as a bond graphic element that computes the mass flow rates of gases into and out of the fuel cell, is described. Next, thermal and pneumatic phenomena are analysed and implemented within bipolar plates. A library of necessary components is created using the Modelica language.

Inlet/outlet valves
The total mass flow rate into and out of the volume depends on the pressure difference and is described by Equation (1) [4].
where u is the valve discharge coefficient, A is the orifice area, r mix is the density of gas mixture and P u and P d are upstream and downstream pressures respectively. The mass flows of individual components in the gas mixture are determined by mass fraction w i : where w i depends on the molar composition of the gas and molar masses of the individual chemical species M i : where N is the number of species in the mixture The corresponding enthalpy flow rate is calculated according to equation (4): where c p;mix ¼ P i x i c p;i is the specific heat capacity of the mixture. Equations (1)e(4) are governing equations of the directed bond graph element Rth depicted in Fig. 5. This thermo-hydraulic restrictor operates under assumption that the upstream pressure at the source is always higher than downstream pressure, thus determining the direction of gas flow (left to right). The block 'Gas in' calculates partial pressures of components in the mixture based on provided temperature, pressure and relative humidity. The calculated values are then used as inputs to a multi-port source of effort mSe, as shown in Fig. 5. The signal 'u' is the control signal that regulates the state of the valve.

Bipolar plates
Bipolar plates play multiple roles in fuel cells. They provide the structural backbone to the fuel cell stack and conduct the generated electrical current. The gas flow channels embedded in them are designed to distribute reactant gases equally to the electrodes, while aiding with removing excess liquid water. In addition to that, some bipolar plates can have cooling channels incorporated in them. In this model, each plate is approximated by two control volumes (CVs), each of which corresponds to the specific physical domain: CV1 Gas flow channels (pneumatic domain). CV2 Solid material of the plate (thermal domain).
Pneumatic domain. Optimal mass transport within the gas flow channels is important to provide the fuel cell with a steady supply of reactants at sufficient concentrations. According to the law of mass conservation, the masses of gases must be conserved within CV1. For the cathode this is described by Equations (5)e (7): where _ m H 2 O , _ m O 2 , and _ m N 2 are mass flow rates of water vapour, oxygen, and nitrogen respectively. Subscripts in and out are assigned to quantities entering and leaving the control volume. Subscript MEA correspond to the amounts transported between the channels and the MEA.
Similarly, the mass conservation equations for the anode are: Pressure losses within the channels arise due to friction with the channel walls. Because of this the mean flow velocity in all parallel channels can be approximated as follows: where A ch is the cross-sectional area, L ch is the length and N ch is the number of parallel channels. Hydraulic diameter D h for rectangular channels is: where a ch and b ch are channel depth and width. The product of frictional losses f with the Reynolds number Re for rectangular channels: fRe ¼ 24 À 1 À 1:3553a ch þ 1:9467a 2 ch À 1:7012a 3 ch þ 0:9564a 4 ch À 0:2537a 5 ch Á (12) where a ch ¼ b ch =a ch . The dynamic viscosity m mix is calculated using Sutherland's formula [1]. The total fluidic dynamics within the channels volume are governed by the ideal gas law: Thermal domain. Thermal dynamics play a crucial role in fuel cell performance. The temperature needs to be maintained in a specific range of approximately 60 Ce80 C in order to facilitate effective water management and optimal electrochemical activity. To achieve this, a temperature control system is required. Temperature dynamics involve heat radiation to the environment, conduction through the materials, and convection by gases and liquids within the cell.
Thermal energy carried by the gases participates in the overall energy balance of the volume as shown by Equation (14).
where, E is the total energy of the volume, _ H in and _ H out are total enthalpies of gases flowing in and out of the volume. Quantity _ Q chÀBPP describes the convective heat transfer between the solid part of the plate and the volume of gas flow channels. The thermal balance of the CV2 is governed by conductive heat transfer to the neighbouring layers ( _ Q left and _ Q right ), convective heat transfer to and from the cooling and gas flow channels ( _ Q chÀBPP , _ Q cool ) and energy exchange with the environment ( _ Q BPPÀamb ): The convective heat transfer between channels and the solid part of the plate is calculated according to Newton's cooling law: where A surf is the area of solid material in contact with the gas. T bulk and T surf are temperatures in the middle of the channel and at the surface of the material, the coefficient of forced convection h forced is calculated as: where k gas is thermal conductivity of the gas, Nusselt number Nu depends on the properties of the fluid characterised by the Prandtl number Pr ¼ c p m mix =k gas and channel geometry [24]: The heat flow between the surfaces of two solid materials can be expressed by Fouriers' law: where k is the material thermal conductivity, A surf is the contact area, DT is the temperature difference between the two layers, and Dx is the material thickness. The amount of heat lost due to thermal radiation and described by the following relation: where s is the Boltzmann constant, k air is thermal conductivity of air. The heat capacity of materials dictates how fast a component can heat up or cool down. The dynamics of which is described by the following equation: where m is the mass of the component, c p is specific heat capacity of the layer and _ Q is the heat flow through the layer. Equations 16 and 19e21 can be modelled with a combination of R-and C-elements in pseudo-bond graph context. Fig. 6 depicts a bond graph structure of the bipolar plate component. Multiports 1 and 2 is where the gas mixture enters and leaves the gas flow channels of the BPP. Three 0junctions 0 1 -0 3 enforce the mass conservation relations of Equations (5)e(7). Since the diffusion of nitrogen is not considered in this model, only water vapour and oxygen flows to and from the MEA component are modelled by multiports 3 and 4. The gas dynamics within the gas flow channels are enforced by Equation (13) inside thermo-hydraulic capacitive Cth-element [4]. The element has Nþ1 ports, where N is the number of chemical species in the gas mixture and another port for thermal energy exchange with its surroundings. The element receives the mass flow rates of matter and computes the resulting temperature and pressure change.
The bond graph structure around junctions 1 1 -1 4 is within thermal domain. Energy conservation is ensured by junctions 0 4 and 0 5 that represent Equations (14) and (15) respectively. Bond graphic elements R 1À2 compute _ Q chÀBPP the combined convective and conductive heat flow between the gas flow channels and the solid part of the plate. Elements R 3À5 calculate _ Q BPPÀamb e the amount of heat lost to the environment by conduction, convection and radiation. Ports 5 and 6 represent points of physical contact of the plate with neighbouring layers and transfer thermal energy via conduction ( _ Q left and _ Q right ). Port 7 models the convective heat exchange with the cooling channels _ Q cool . Tables 4 and 5 summarise the properties of the BPPs and dimensions of cooling channels.

Membrane electrode assembly
Combination of the gas diffusion layers (GDLs), catalyst layers (CLs) and polymer membrane (PEM) are modelled as a single block named MEA. According to Fig. 3, this block deals with electrochemical, transport and thermal phenomena occurring within the MEA.
Gas diffusion layers. The GDL is required in order to redistribute the reactants evenly to the catalyst sites. GDLs are usually made out of carbon fibre or paper.
As the reactants are being consumed in the electrochemical reaction at the CL, the concentration of gases across the GDL changes according to Fick's law: Where J i is the molar flux of species across the GDL, d GDL is the layer thickness, D ij is the binary diffusion coefficient of species in empty medium, c i (x) is the molar concentration of species i across the layer. Molar flux is defined as: Fig. 6 e Bond graph for a cathode bipolar plate component. The molar concentration of an ideal gas can be expressed using Equation (13) in the following form: Using Equations (23) and (24) Equation (22) can be rewritten in order to be consistent with the power variables used in the bond graph as: The diffusion coefficient depends on the temperature and pressure of the mixture and it can be computed using Formula 26.
where p ci;cj are the critical pressures and T ci;cj are the critical temperatures of the gases in the mixture. Constants a and b depend on the kind of gases involved [1]. The calculated value of D ij needs to be adjusted for the case when the process occurs within a medium. GDLs are often made out of carbon fibre paper, and, therefore consist of random fibrous media. Nam and Kaviany [25] suggest a percolation type correlation for the calculation of an effective diffusion coefficient that takes into account the amount of liquid water present inside the layer: 1 À 0:11 0:785 Â ð1 À sÞ 1:5 (27) where ε is the porosity of the layer and s is the liquid water saturation of the volume and calculated as s ¼ V l =V p , where V l is the volume occupied by liquid water and V p is the total volume of the pores. Equations (25)e (27) are constitutive equations for another bond-graphic restrictive elements (labelled Rd) that calculates the diffusion mass flow of gases through the GDL due to the pressure difference. In order to calculate the partial pressure of reactants at the catalyst surface, the volume of GDL pores is approximated as additional control volume, represented by previously mentioned multiport capacitive element Cth. The bond graph in Fig. 7 is a part of the MEA component and illustrates the model of the cathode GDL. The multi-ports 3 and 4 correspond to those in Fig. 6. Port 8 represents the mass flow of oxygen consumed in the reaction, while port 9 corresponds to the water vapour generated in the reaction. The amount of water transferred through the membrane is represented by port 10. The mass conservation is ensured by two 0-junctions.

Catalyst layer
Catalyst layers in PEM fuel cells consist of nano-particles of platinum alloy suspended on a thin layer of carbon support. The total surface area of catalyst particles is called the electroactive surface area (EASA) and needs to be as big as possible in order to provide enough reaction sites. At the anode electrode, the hydrogen oxidation reaction occurs, while the cathode electrode facilitates the oxygen reduction reaction. The amount of reactant gas consumed in each reaction is determined by the current load applied to the cell and is expressed as molar flow rate. It is, therefore, needed to transform power variables used in the pneumatic domain ( _ m and p) into chemical ones ( _ n and m). The molar flow rate is related to mass flow rate as where M i is molar weight of the gas. In order to calculate the chemical potential of species, changes on enthalpy Dh and entropy Ds of the reactions are needed. These thermodynamic variables are calculated with the following equations: where h 0 i and s 0 i are enthalpy and entropy values of individual species at reference state conditions . These values are summarised in Table 6 [26].
The chemical potential is then calculated as: Heat generated by each chemical reaction can be evaluated as: An additional element called the chemical transformer (labelled as cTF) is introduced to convert pneumatic variables  Fig. 7 e Bond graph representation of the cathode GDL.
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 4 2 ( 2 0 1 7 ) 2 9 4 0 6 e2 9 4 2 1 into chemical ones using Equations (28)e(31) as illustrated by Fig. 8. This block also acts as a source of heat generated in the reaction. It important to note that the energy is not conserved within cTF block which doesn't comply with definitions of bond graphs, however this does not impede the accuracy of the resulting model. The reaction rate depends on temperature, concentration of reactants and reaction rate coefficients [1,27]. However, at present these coefficients are not known, so in this model the reaction rate is simply governed by the electrical current drawn from the fuel cell according to Faraday's law: where z is the number of electrons participating in the reaction, F is the Faraday's constant.
Gibbs free energy DG is the amount of energy that can be extracted from the chemical species calculated as follows: The reversible potential of the cell is expressed as: The effect of partial pressure and temperature change on the reversible cell potential is found according to the Nernst equation:

Polymer electrolyte membrane
The primary function of the polymer electrolyte membrane is to conduct H þ -ions released during the H 2 oxidation reaction at the anode and transfer them over to the cathode. Conventionally, Nafion is used as an electrolyte material and it's ionic conductivity properties are strongly dependent on it's level of humidification. Since water is generated at the cathode electrode, the membrane is humidified non-uniformly along its thickness. The amount of water at each side of the membrane can be estimated using the empirical equation [1]: where k2fan; cag and a H 2 O;k is water activity equal to the ratio of water vapour pressure and saturation pressure P sat ðTÞ at the given temperature: The mean water content in the membrane is therefore: There is also a mass flow of water occurring through the membrane which is governed by two phenomena: electroosmotic drag and back diffusion. The former is caused by water molecules dragged by hydrogen ions as they travel through the membrane and the latter is caused by the difference of the partial pressure of water vapour at each side of the membrane.
The mass flow of water due to electro-osmotic drag is calculated according to: where i is the current density drawn from the fuel cell, A mem is the active surface area of the membrane and n d is the coefficient of electro-osmotic drag calculated using [28]: The back diffusion of water from the cathode to the anode depends on the difference of water content and physical properties of the membrane: where r mem is the density of the dry membrane, M mem is the equivalent weight of the polymer, d mem is membrane thickness and the effective diffusivity coefficient D w is calculated using: where diffusivity of water is calculated by Ref. [28]: 10 À10 ; l < 2 10 À10 ð1 þ 2ðl À 2ÞÞ; 2 l < 3 10 À10 ð3 À 1:67ðl À 3ÞÞ; 3 l < 4:5 1:25 Â 10 À10 ; 4:5 l Combining Equations (39) And (41), the overall water mass flow through the membrane is: These calculations are performed by a non-bondgraphic block called 'PEM' (see Fig. 10) that receives partial pressures of water from the anode and the cathode sides, temperature of the membrane, and electrical current density. It then computes the water content l and the total mass flow of water through the membrane.

Voltage output
The reversible potential of a fuel cell from Equation (35) is reduced by various voltage loss effects. In electro-chemistry such losses are often referred to as over-voltages There are three main types of voltage loss mechanisms in PEM fuel cells: activation (h act ), ohmic (h ohm ), and concentration (h con ) loss. The voltage produced by a cell is expressed with the following equation: Activation losses are predominant in lower currentdensities and are calculated using the well-known Tafel equation [1]: where a is the transfer coefficient and i 0 is the exchange current density. The charge transfer through the bipolar plates, GDL's and catalyst layers is hindered by their corresponding ohmic resistances, which remain relatively constant throughout the lifetime of a fuel cell. However, the membrane's ability to conduct ions largely depends on it's level of humidification. The ionic conductivity of Nafion membranes is estimated from the following empirical relation [1]: sðT; lÞ ¼ ð0:005193l À 0:00326Þexp 1268 where l is calculated using Equation (38). The ionic resistance is therefore: Ohmic losses occur due to the combination of ionic R ion and ohmic R el resistances of the electrolyte: The concentration over-potential h con occurs at higher current densities when the reactants are consumed faster than they are provided. 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 4 2 ( 2 0 1 7 ) 2 9 4 0 6 e2 9 4 2 1 Under such conditions, the concentration of fuel or oxidant at the catalyst layers approaches 0, what puts a limit on the fuel cell performance. This effect is described by the limiting current density i L : Each of the over-voltages is represented by a 2-port nonlinear thermo-electrical resistor denoted as RS. In contrast to ideal 1-port R resistive elements, where no thermal energy is dissipated, RS-elements calculate the amount of thermal energy generated according to: A generic RS-element is shown in Fig. 9 with electric port on the left and thermal port on the right. The generated heat is the flow variable of the pseudo-bond that leaves the electrical domain and enters the thermal domain.
Additionally, voltage dynamics are determined by cell double layer capacitive effect described by Equation (53): The bond graph of the MEA component (excluding the GDL part of Fig. 7) is depicted in Fig. 10. The modulated source of flow mSf provides the demanded load current signal, while detector of effort De outputs the generated voltage of the fuel cell. Three elements RS 1À3 calculate the voltage losses using Equations (46)e(50) as well as generated heat with Equation (52). The Nernst Equation (35) is implemented by four transformers TF 1À4 which apply the corresponding stoichiometry coefficients and three chemical transformers cTF 1e3 which compute the chemical potentials and molar flows of gases. The temperature of the MEA is determined by the element C 1 and two resistors R 1,2 that calculate the conductive heat transfer to the BPPs. The energy conservation is ensured by junctions 0 1 and 0 2 . Table 7 contains the parameters of the MEA component used in this model.

Complete model
The individual components discussed above, are put together as illustrated by Fig. 11. The model closely resembles the topology of a real fuel cell and corresponds to the diagram in Fig. 3. The ports shown in Fig. 11 are included to show how internal components are connected to external ones.
It is possible to encapsulate Fig. 11 and augment it with supporting systems such as cooling channels to create the overall fuel cell model of Fig. 12. The illustrated configuration consists of two cooling loops with liquid water at either side of the cell and two end plates supporting the cell. Such modular approach allows the creation of models with multiple fuel cells, each of which having different characteristics. In Fig. 12 inlet mass flows are regulated by Rth 1 and Rth 2 elements. Each cooling loop is comprised of a single Cth-element and elements Rc 1À4 calculate the heat transfer rates according to equation (16). End-plate components compute thermal losses to the environment, same as it was illustrated for bipolar plates in Fig. 6. The parameters for end-plates and cooling channels used in the model are provided in Tables 8 and 9.

Model validation and parameter identification
In order to check if the designed model is correct, simulation results are compared to experimental data collected from a test rig. The fuel cells used in this study were manufactured by Pragma Industries [29]. The graphite bipolar plates are 6 mmthick with integrated 7-fold parallel serpentine gas flow channels. The liquid coolant can flow through the channels  Fig. 11 e A single cell assembly.
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 4 2 ( 2 0 1 7 ) 2 9 4 0 6 e2 9 4 2 1 also embedded into the plates. The MEA has 100cm 2 active area and the polymer membrane is Nafion XL with 27.5mm thickness. The GDL is a Sigracet gas diffusion media 10 BCE has thickness of approximately 415mm. Mass flow controllers and humidifiers for both gas feeds maintained the appropriate levels of stoichiometry and relative humidity during the experiments.
First, a single cell was subject to evaluation by electrochemical impedance spectroscopy (EIS) in order to determine      the values of electrical resistance and double-layer capacitance. Frequency responses were obtained at 10, 25, and 50 A and the resulting Niquist plots are shown in Fig. 13. The value of the double layer capacitance can be estimated using the following formula [30]: where f co is the cut-off frequency of the response, R c is the activation resistance calculated from R c ¼ R f À R el . Table 10 summarises the findings from the EIS measurements. From these measurements, the estimated double layer capacitance C dl was found to be around 0.015F and the ohmic resistance of the membrane R el z 0.38Ucm 2 .
During the polarization curve test the cell operated at ambient pressure and at the temperature of 60 C and measurements were taken as shown in Fig. 14. Setting the corresponding model parameters to the obtained values and imposing the same current load, the simulated voltage output can be compared to the experimental one as shown in graph in Fig. 14. Fig. 14 shows that in the regions of activation losses (between 0 and 0.1A/cm 2 ) and ohmic losses (between 0.1 and 0.55 A/cm 2 ) the agreement between the experimental and modelled results is very good. The deviation between results increases in the region of concentration losses (from 0.55 A/ cm 2 and higher). This is due to the fact that the model doesn't fully account for the electrode porosity and the effects of liquid water formation within the fuel cell.
The simulations were performed using a OpenModelica's default Differential Algebraic System Solver (DASSL). This is a popular solver which is well suited for simulating engineering systems which are computationally stiff because they often contain phenomena occurring at different time scales [31]. Although, the total number states within the final 2-cell model is 33, on a computer with an Intel Core i5 2.2 GHz and 8 Gb RAM the model was compiled in 20s while simulating 140min of time as shown in Fig. 15 took 121s. If faster simulation performance is required, the order of the model can be reduced by lumping individual cells into a single block and simplifying the thermal domain dynamics by combining individual thermal masses of fuel cell layers into one.
Thermal characteristics of the fuel cell materials obtained from the literature and hardware specification are summarised in Tables 8 and 9 [26]. In order to validate them, another test was performed in which the stack was subject to the current load profile depicted at the top of Fig. 15. The voltage response shown in the middle of Fig. 15 also shows good correlation with the experimental measurements. The differences in the time-periods 229e244, 274e289 and 319e334 min are due to increased water formation at the high current loads. The measured and simulated temperatures of the coolant outlet are shown at the bottom of Fig. 15. The average temperature difference between the two results is 0.56 K signifying a good agreement. Due to the modular structure of the model, each cell within the stack can be analysed individually. For example, Fig. 16 shows the Fig. 15 e Temperature and voltage response of a 2-cell stack to a variable current load. 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 4 2 ( 2 0 1 7 ) 2 9 4 0 6 e2 9 4 2 1 calculated mean water content and area-specific resistances of each membrane of the stack. It can be seen that in current configuration cell 2 has lower ohmic resistance than cell 1 due to extra humidification at the cathode side of the stack.

Conclusions
The goal of this paper was to design an accurate and detailed bond-graphic model of a PEM fuel cell. A general description of the approach was provided and a number of published papers on the subject was reviewed. Although bond graph formalism may seem superfluous or even arcane to a beginner at first, mastering the approach provides a unique perspective on engineering systems. The unified representation of all components allows the modeller to gain insight not only into dynamical but also into structural and causal properties of the system. It is also clear that bond graphs are very well suited for modelling fuel cell devices as demonstrated by the reviewed publications and modelling work performed in this paper.
The designed model relies on a pseudo-bond graph representation of thermo-fluid phenomena and is implemented in Modelica modelling language. The unique feature of the model is it's hierarchical structure. Several blocks corresponding to their physical counterparts (mass flow controllers, bipolar plates, cooling channels, MEA) are created. Each block is an encapsulated model of the separate physical process. For example, bipolar plate sub-model exclusively computes the pressure within the gas flow channels and heat transfer inside the plate. This means that these physical processes can be modelled and analysed independently of other components in the system. As a result, the model topology makes the analysis of component interactions straightforward and clear.
Another innovative aspect of the proposed model is the adaptation of multi-bonds for the analysis of multicomponent mixtures. This enables a streamlined graphical portrayal of the process and is highly beneficial for bondgraphic implementation of such phenomena.
In order to perform model simulations, parameters were gathered from literature and hardware description. As it was demonstrated, experimental and simulation results show strong correlation between each other. Consequently, the designed model can be used as a tool for studying system behaviour in various operational scenarios, such as those demonstrated in Section Model Validation and Parameter Identification.
The model can be further improved by incorporating a more detailed description of diffusion mechanisms and implementation of the liquid water formation in the channels and GDLs. Thanks to the object oriented nature of the model, each individual block can be augmented with features enabling reliability assessment of the component under varying operational conditions. The addition of such features is the next step of the research.