Open-Source Dynamic Matlab / Simulink 1D Proton Exchange Membrane Fuel Cell Model

: This work presents an open-source, dynamic, 1D, proton exchange membrane fuel cell model suitable for real-time applications. It estimates the cell voltage based on activation, ohmic and concentration overpotentials and considers water transport through the membrane by means of osmosis, di ﬀ usion and hydraulic permeation. Simpliﬁed equations reduce the computational load to make it viable for real-time analysis, quick parameter studies and usage in complex systems like complete vehicle models. Two modes of operation for use with or without reference polarization curves allow for a ﬂexible application even without information about cell parameters. The program code is written in MATLAB and provided under the terms and conditions of the Creative Commons Attribution License (CC BY). It is designed to be used inside of a Simulink model, which allows this fuel cell model to be used in a wide variety of 1D simulation platforms by exporting the code as C / C ++ .


Introduction
Proton exchange membrane fuel cells (PEMFC) can contribute to achieving the goal of a sustainable energy supply and production.High efficiency and power density as well as zero emissions are beneficial for both stationary and mobile applications.However, designing the water management, cooling and media supply of a PEMFC-system is challenging.A model-based approach for the simulation of such a system can be a valuable tool in this matter.
Many PEMFC-models have been developed in the past, with varying objectives.Some models intend to deliver highly accurate results through means of computational fluid dynamics (CFD) [1][2][3][4] whilst others target a faster simulation by reducing the complexity [5,6].However, in terms of computational time, these models still operate in the range of min or hours.This makes them unsuitable for large system simulations like complete vehicles or even real-time evaluation.To account for this, simplified fuel cell models have been developed in recent history to reduce the required CPU-time at the expense of accuracy [7][8][9].Some of the aforementioned models have been made publicly and freely available, others remain closed-source.Additionally, depending on the chosen programming language as well as structure of model inputs and outputs, the compatibility with various simulation environments might be restricted.
After considering the previous research on fuel cell modeling, the motivation for creating the presented model was to combine the following three aspects.First, the compatibility with commonly used 1D simulations' environments.Since each software has its strengths, weaknesses and price, individual programs might not be available for everyone who would benefit from a fuel cell model.In particular, educational facilities often lack the funds to provide expensive licenses.Offering a model Energies 2019, 12, 3478 2 of 12 that can be used in multiple environments, and therefore can be utilized by a large audience, was a primary motivation.Second, there is a low computational demand.Since large systems (e.g., cars or trains completing a drive cycle) consist of many components, each individual model is required to be of reduced complexity in order to keep the overall simulation time low.However, the operating conditions inside of such a dynamic system are constantly shifting.This is why a compromise between speed and accuracy, but with a focus on speed, was another objective.Third, there is an open-source software license.Since open-source models can be further expanded upon and offer value for research, education and development, this represented the third requirement.
The model at hand was designed to be used either as a cell model or as part of a fuel cell stack.In its core, it represents a cell model with internal (optional) scaling for stack parameters like cell area and quantity.In order to ensure compatibility, the MATLAB-Simulink (version used: R2016b by The MathWorks, Inc.) environment was chosen.Furthermore, to allow for fast calculations and even real-time applications, simplified equations are used, even though the underlying Nernst-equation has recently been found to show considerable inaccuracies [10].Additionally, to further reduce the computational demand, no discrete model was used for the membrane electrode assembly (MEA).A compromise was made between simulation quality regarding water and gas transport and the speed of the calculations.
Possible use cases for the presented model are the creation of polarization curves or the cell performance estimation under varying conditions like in a moving vehicle. Figure 1 depicts the basic structure of the program code, which was designed to be run inside of a MATLAB-function-block as part of a Simulink model.Based on various inputs and physical parameters, the model estimates cell voltage as well as several other outputs in every time step.A complete list of the model inputs, outputs and parameters is provided as Supplementary Material.The calculations inside the code consist of simplified equations, presented in the following section.In particular, educational facilities often lack the funds to provide expensive licenses.Offering a model that can be used in multiple environments, and therefore can be utilized by a large audience, was a primary motivation.Second, there is a low computational demand.Since large systems (e.g., cars or trains completing a drive cycle) consist of many components, each individual model is required to be of reduced complexity in order to keep the overall simulation time low.However, the operating conditions inside of such a dynamic system are constantly shifting.This is why a compromise between speed and accuracy, but with a focus on speed, was another objective.Third, there is an open-source software license.Since open-source models can be further expanded upon and offer value for research, education and development, this represented the third requirement.The model at hand was designed to be used either as a cell model or as part of a fuel cell stack.In its core, it represents a cell model with internal (optional) scaling for stack parameters like cell area and quantity.In order to ensure compatibility, the MATLAB-Simulink (version used: R2016b by The MathWorks, Inc.) environment was chosen.Furthermore, to allow for fast calculations and even realtime applications, simplified equations are used, even though the underlying Nernst-equation has recently been found to show considerable inaccuracies [10].Additionally, to further reduce the computational demand, no discrete model was used for the membrane electrode assembly (MEA).A compromise was made between simulation quality regarding water and gas transport and the speed of the calculations.
Possible use cases for the presented model are the creation of polarization curves or the cell performance estimation under varying conditions like in a moving vehicle. Figure 1 depicts the basic structure of the program code, which was designed to be run inside of a MATLAB-function-block as part of a Simulink model.Based on various inputs and physical parameters, the model estimates cell voltage as well as several other outputs in every time step.A complete list of the model inputs, outputs and parameters is provided as supplementary material.The calculations inside the code consist of simplified equations, presented in the following section.

Mathematical Model
In this fuel cell model, the real cell voltage  is estimated by subtracting the voltage losses inside the cell from the ideal cell voltage  , .These are summarized as activation  , ohmic  and concentration  overpotentials [11]: Figure 1.Model structure and calculation order.Italic functions are implemented as a subfunction.

Mathematical Model
In this fuel cell model, the real cell voltage E cell is estimated by subtracting the voltage losses inside the cell from the ideal cell voltage E T,p .These are summarized as activation E act i , ohmic E ohm and concentration E con i overpotentials [11]: Energies 2019, 12, 3478 3 of 12

Ideal Voltage
The theoretical maximum voltage of a PEMFC under reference conditions E 25C,1atm can be calculated with the Gibbs free energy ∆G as well as the Faraday constant F and the number of electrons involved n [11]: Variations in temperature T and partial pressure of reactants p i can be accounted for by using the Nernst-Equation [11]: Neglecting influences of changing enthalpy ∆H and entropy ∆S, as well as assuming the product water to be in the liquid phase, the ideal cell voltage E T,P can be expressed as follows [11]:

Activation Overvoltage
Based on the Butler-Volmer equation, the activation overpotential E act,i can be estimated as a function of current density i, exchange current density i 0 as well as temperature T and charge transfer coefficient α i [12,13]: where R denotes the universal gas constant and F the Faraday constant.
The exchange current density i 0 of a platinum electrode as a function of partial pressure p i , temperature T, catalyst loading L i and specific area a i can be calculated on the basis of a reference value i 0,ref [11]: For simplicity, a constant value was used for the activation energy ∆G i , even though it can vary under real operating conditions [14].

Membrane Water Content
For estimating the water content λ inside of a Nafion-membrane, a function of water activity a is used [15]: As a simplification, ab-and desorption of water were neglected.Furthermore, the distribution of water along the membrane geometry was assumed to be uniform.
The water activity a can be expressed as [16] where RH denotes the relative humidity (for ideal gas properties) and s the liquid water volume fraction.For this model, it is assumed that liquid water is only present in the catalyst layer.For RH and s, a logarithmic average is used to account for nonhomogeneous distribution inside the channels (see Equation ( 17)).This can be disabled by suppling the model with equal values for input and output.
Energies 2019, 12, 3478 In case of a cold start at subzero temperatures, the behavior of frozen water inside of a Nafion-membrane is approximated using the following terms [16]: To estimate the water concentration C w , a simple proportional correlation between the membrane density ρ mem , equivalent weight and water content λ is used [11]:

Ohmic Overvoltage
Regarding ohmic resistances inside the cell, only the membrane resistance as the most influential factor is considered.Therefore, the overpotential can be calculated with the current density i, membrane conductivity σ mem and (wet) thickness δ mem .Since membrane thickness-at typical PEMFC operating conditions-is only marginally affected by swelling, the thickness is assumed to be constant [17,18]: To estimate the membrane conductivity σ mem , the following correlation is used [9,12,19]: It is mostly dependent on membrane water content λ and temperature T mem , whilst also being affected by the material properties of equivalent weight EW and membrane density ρ mem as well as the density of water ρ W (T). A simple arithmetic average of anode and cathode site values is used for further calculations.
The density of water ρ w (T) at 1[bar] as a function of temperature T can be approximated by [20]: T in • C.

Concentration Overvoltage
Using the Nernst-equation, concentration overvoltage E con i is described as a function of temperature T, current density i and limiting current density i L [11] where R is the universal gas constant, n the number of electrons involved and F the Faraday constant.Since this ideal equation often underestimates the real overvoltage, e.g., because of uneven gas concentration, a correction factor (see Supplementary Materials) is used to adjust the results [11,21].
Energies 2019, 12, 3478 5 of 12 For the limiting current density i L , a simplified expression including the Faraday F constant, number of electrons n, diffusion coefficient D i , gas concentration C i and diffusion distance (here: electrode thickness) δ e is used [11]: To calculate the diffusion coefficient D i , the model expects an external reference value D i,ref for the gas mixture.It therefore relies on the external calculation of gas concentration and diffusivity.The reference value will then be adjusted for electrode porosity and tortuosity τ as well as temperature T, pressure p and liquid water volume fraction s [9,22,23]: Porosity and tortuosity τ are used to approximately describe the geometry of the electrodes, while a value of 1 for the liquid water volume fraction s represents a fully flooded channel.As a simplification, it is assumed that the gas diffusivity in liquid water is 0. The approach of relying on D i,ref as a model input for the calculation of D i ensures compatibility with arbitrary gas mixtures and composition of external models for fluid mechanics-for instance, when used in a vehicle model and being connected to components for the media supply.
A simple logarithmic average is used to account for nonhomogeneous distribution of the gas concentration C lm .If desired, this can be disabled by supplying the same value for input and output [5]:

Cell and Stack Performance
What is labeled as the electrical efficiency η electric in this work relates to the lower heating value of hydrogen E LHV [11]: and, furthermore, is used to calculate the heat flow .
Q of the cell/stack based on the power delivered P stack .As a simplification, it is assumed that the product water fully evaporates before leaving the cell/stack [11]: .

Water Transport
In this model, the estimation of the flow of water j w from anode to cathode is divided in three categories: osmotic j osmo and diffusive j diff flow as well as hydraulic permeation j hyd .A positive value means an increase in water concentration: Three major simplifications are applied to reduce the complexity of water transport mechanisms.First, only the flow through the membrane is considered, whilst transport through the porous media of catalyst and gas diffusion layers is neglected.Second, it is assumed that liquid water is only present once the gas mixture is saturated.Finally, liquid and vapor phases are not directly considered.Instead, the chosen equations for diffusive and hydraulic flow are adjusted by several functions (Figure 2) created with the MATLAB curve fitting application and experimental data from Adachi et al. [24].As a result, three cases are described in the following sections: vapor-vapor (VVP), liquid-vapor (LVP) and liquid-liquid permeation (LLP).24)-(32), adjusted with data from Adachi et al. [24]); VVP: cathode 96%, anode 38% RH; LVP: cathode liquid volume fraction 100% (flooded), anode 38% RH; LLP: cathode and anode flooded, ∆ 1 [bar] note: only a few data points were available to develop each function, which leads to some numerical inaccuracies.

Hydraulic Permeation
Beyond that, hydraulic permeation is only considered if liquid water is present on both sides of the membrane, since the pressure difference has a minor impact on water vapor transport [27].The hydraulic flow  is estimated by a linear correlation with the pressure gradient ∇ -here: difference between cathode and anode-affected by the cell area , dynamic viscosity of water  as well as water concentration inside the membrane  , its thickness  and hydraulic permeability  [16]: For the hydraulic permeability  , a direct dependency on the membrane water content  is assumed [16]: Furthermore, the dynamic viscosity of water is approximated by a function of temperature where  denotes a reference value while  ,  ,  and  are constants.The pressure  has been neglected, since it has a minor impact on  at typical PEMFC operating pressures.Experimental data suggest a nonlinear relation between the hydraulic flow and membrane thickness [24], hence an adjustment-function is applied for the hydraulic permeability  at a reference pressure difference of 0.025  :   24)-( 32), adjusted with data from Adachi et al. [24]); VVP: cathode 96%, anode 38% RH; LVP: cathode liquid volume fraction 100% (flooded), anode 38% RH; LLP: cathode and anode flooded, ∆p 1 [bar] note: only a few data points were available to develop each function, which leads to some numerical inaccuracies.

Osmosis
Water transport due to osmosis j osmo is expressed by a linear correlation composed of the Faraday constant F, the current I [11] and the osmotic coefficient n osmo , which is dependent on the (average) water content of the membrane λ [25]:

Diffusion
In terms of diffusion j diff , the abovementioned functions to differentiate between the liquid and gaseous phase are utilized.The basis for these calculations is given by [17] which considers the cell area A, (average) diffusive coefficient of water D λ , water concentration gradient ∇C w -here: difference between cathode and anode-and membrane thickness δ mem .The latter is treated as a constant (=201 µm) in the above equation and variations are instead considered by the adjustment-functions (Equations ( 26)-( 32)).
To estimate the diffusion coefficient for water through the membrane D λ , the following expression dependent on membrane temperature T mem and (average) water content λ is applied [9,26]: Energies 2019, 12, 3478 7 of 12 Subsequently, the diffusive flow j diff is adjusted for the thickness of the membrane.VVP-correction is applied in the complete absence of liquid water, whilst LVP-correction is used if at least one side is flooded with liquid water.For mixtures of gaseous and liquid phases, a linear scaling is applied:

Hydraulic Permeation
Beyond that, hydraulic permeation is only considered if liquid water is present on both sides of the membrane, since the pressure difference has a minor impact on water vapor transport [27].The hydraulic flow j hyd is estimated by a linear correlation with the pressure gradient ∇p-here: difference between cathode and anode-affected by the cell area A, dynamic viscosity of water µ H 2 O as well as water concentration inside the membrane C w , its thickness δ mem and hydraulic permeability K λ [16]: For the hydraulic permeability K λ , a direct dependency on the membrane water content λ is assumed [16]: Furthermore, the dynamic viscosity of water is approximated by a function of temperature T [28] where µ 0 denotes a reference value while a µ , b µ , c µ and d µ are constants.The pressure p has been neglected, since it has a minor impact on µ w at typical PEMFC operating pressures.Experimental data suggest a nonlinear relation between the hydraulic flow and membrane thickness [24], hence an adjustment-function is applied for the hydraulic permeability K λ at a reference pressure difference of 0.025 [bar]: Subsequently, the hydraulic flow is corrected for the actual pressure difference:

Application
Two modes of operation are available for using this model, as displayed in Figure 3. First, the cell performance can be calculated solely based on physical parameters.Second, the cell performance can be estimated based on the voltage deviation from a supplied polarization table.In the latter case, a distinction can be made between using a single or multiple polarization curves as reference.
For this purpose, the inputs of the MATLAB-function are divided into two categories: one for the state variables mandatory to estimate the cell performance and another for the polarization table references.When running in mode 1, the inputs from the second category are ignored.In mode 2a, the model expects reference values for the experimental conditions of the polarization curve recording.Lastly, in the case of 2b, internal calculations for voltage deviation can be disabled by supplying state Energies 2019, 12, 3478 8 of 12 variables as table references.e.g., when using polarization curves for different temperatures, supplying the current cell temperature will lead to no additional adjustments for this factor, since current and reference values are identical.The operation mode has no impact on the calculations for concentration overpotential and water transport, however, since these heavily depend on stack composition and media supply.

Application
Two modes of operation are available for using this model, as displayed in Figure 3. First, the cell performance can be calculated solely based on physical parameters.Second, the cell performance can be estimated based on the voltage deviation from a supplied polarization table.In the latter case, a distinction can be made between using a single or multiple polarization curves as reference.For this purpose, the inputs of the MATLAB-function are divided into two categories: one for the state variables mandatory to estimate the cell performance and another for the polarization table references.When running in mode 1, the inputs from the second category are ignored.In mode 2a, the model expects reference values for the experimental conditions of the polarization curve recording.Lastly, in the case of 2b, internal calculations for voltage deviation can be disabled by supplying state variables as table references.e.g., when using polarization curves for different temperatures, supplying the current cell temperature will lead to no additional adjustments for this factor, since current and reference values are identical.The operation mode has no impact on the calculations for concentration overpotential and water transport, however, since these heavily depend on stack composition and media supply.
It is also possible to use this model outside of the MATLAB environment via Co-Simulation.As an example, the program code can be run inside of a MATLAB-Function-block within a Simulink model, which can be compiled as C/C++ code.Depending on the target software, the exact Simulink model composition-regarding inputs, outputs and parameters-and compilation procedure may be wary and has to be looked up in the corresponding documentation.Following this practice allows the fuel cell model to be used in any simulation platform with support for Simulink coupling.

Simulation Results
Figure 4 shows a set of polarization curves calculated by the presented model.A variation of the chosen parameters affects the overpotentials and therefore the overall cell voltage (Section 2).In terms of computational time, creating a polarization curve with hundreds of data points only takes a few s on a modern CPU.This shows the suitability of the model at hand for large system simulations or real-time applications.The inputs and most important parameters for the reference case are represented by Table 1, and a complete list of all model parameters is supplied as supplementary material.It is also possible to use this model outside of the MATLAB environment via Co-Simulation.As an example, the program code can be run inside of a MATLAB-Function-block within a Simulink model, which can be compiled as C/C++ code.Depending on the target software, the exact Simulink model composition-regarding inputs, outputs and parameters-and compilation procedure may be wary and has to be looked up in the corresponding documentation.Following this practice allows the fuel cell model to be used in any simulation platform with support for Simulink coupling.

Simulation Results
Figure 4 shows a set of polarization curves calculated by the presented model.A variation of the chosen parameters affects the overpotentials and therefore the overall cell voltage (Section 2).In terms of computational time, creating a polarization curve with hundreds of data points only takes a few s on a modern CPU.This shows the suitability of the model at hand for large system simulations or real-time applications.The inputs and most important parameters for the reference case are represented by Table 1, and a complete list of all model parameters is supplied as Supplementary Material.

Conclusions
The model at hand represents a one-dimensional, dynamic proton exchange membrane fuel cell.To estimate the cell voltage, activation, ohmic and concentration overpotentials are calculated and subtracted from the ideal cell voltage in every time step.Furthermore, the water transport through the membrane by means of osmosis, diffusion and hydraulic permeation is considered.In order to reduce the complexity and computational load, simplified correlations are used to estimate the cell performance.This approach allows the model to be used in complex systems such as complete vehicle models or real-time applications.Two modes of operation allow for flexible use of the fuel cell model by supplying either polarization tables or physical cell parameters, which enables a quick model setup.The program code is written in MATLAB and designed to be used in a MATLAB-functionblock inside of a Simulink model.By compiling the Simulink model as C/C++, this PEMFC model can also be used within any software tool that supports Simulink coupling.It is supplied under an openaccess license to make it available to anyone for free.
However, the use of simplified equations for cell performance estimation also reduces the accuracy of the results.In particular, the consideration of the MEA can be further expanded because water and gas transport as well as concentration, can be significantly affected by the MEA composition.Additionally, cell geometry and local differences in the distribution of temperature, current density and reactants can also affect the overall performance.For future expansions, the

Conclusions
The model at hand represents a one-dimensional, dynamic proton exchange membrane fuel cell.To estimate the cell voltage, activation, ohmic and concentration overpotentials are calculated and subtracted from the ideal cell voltage in every time step.Furthermore, the water transport through the membrane by means of osmosis, diffusion and hydraulic permeation is considered.In order to reduce the complexity and computational load, simplified correlations are used to estimate the cell performance.This approach allows the model to be used in complex systems such as complete vehicle models or real-time applications.Two modes of operation allow for flexible use of the fuel cell model by supplying either polarization tables or physical cell parameters, which enables a quick model setup.The program code is written in MATLAB and designed to be used in a MATLAB-function-block inside of a Simulink model.By compiling the Simulink model as C/C++, this PEMFC model can also be used within any software tool that supports Simulink coupling.It is supplied under an open-access license to make it available to anyone for free.
However, the use of simplified equations for cell performance estimation also reduces the accuracy of the results.In particular, the consideration of the MEA can be further expanded because water and gas transport as well as concentration, can be significantly affected by the MEA composition.Additionally, cell geometry and local differences in the distribution of temperature, current density and reactants can also affect the overall performance.For future expansions, the computational load should be considered, in order to not increase the model's calculation time too much.

Figure 1 .
Figure 1.Model structure and calculation order.Italic functions are implemented as a subfunction.

Table 1 .
Inputs and parameters for the reference case.Modes of operation.1: cell mode; 2a: table mode, single curve polarization table; 2b: table mode, multi curve polarization table; S: state variable; T: polarization table; R: polarization table reference values.

Table 1 .
Inputs and parameters for the reference case.
catalyst-specific area, cm 2 /mg a µ constant in the calculation of µ w , bar −1 b µ constant in the calculation of µ w , J/mol bar c µ constant in the calculation of µ w , K/mol bar C i molar concentration, mol/cm 3 C w membrane water concentration, mol/cm 3 D i diffusion coefficient, cm 2 /s D λ diffusion coefficient of water through the membrane, cm 2 /s d µ constant in the calculation of µ w , J/mol