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

Abstract

[...]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. [...]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. Concentration Overvoltage Using the Nernst-equation, concentration overvoltage Econi is described as a function of temperature T , current density i and limiting current density iL [11] Econi=RTinFln(iLiiLi−i), 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]. [...]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.

1. 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 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.

2. Mathematical Model

In this fuel cell model, the real cell voltage _Ecell is estimated by subtracting the voltage losses inside the cell from the ideal cell voltage _ET,p . These are summarized as activation _Eacti , ohmic _Eohm and concentration _Econi overpotentials [11]:

_Ecell=_ET,p−_Eacti−_Eohm−_Econi.

2.1. Ideal Voltage

The theoretical maximum voltage of a PEMFC under reference conditions _E25C,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]:

_E25C,1atm=−ΔGnF=2373402×96485=1.23 V.

Variations in temperature T and partial pressure of reactants _pi can be accounted for by using the Nernst-Equation [11]:

_ET,P=−(ΔHnF−TΔSnF)+RTnFln[_{pH2 pO20.5pH2O}].

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 _ET,P can be expressed as follows [11]:

_ET,p=1.482−0.000845 T+0.0000431 T ln(_{pH2 pO20.5}).

2.2. Activation Overvoltage

Based on the Butler–Volmer equation, the activation overpotential _Eact,i can be estimated as a function of current density i , exchange current density _i0 as well as temperature T and charge transfer coefficient _αi [12,13]:

_Eact,i=R_TiαiFarcsinh(i2_i0,i),

where R denotes the universal gas constant and F the Faraday constant.

The exchange current density _i0 of a platinum electrode as a function of partial pressure _pi , temperature T , catalyst loading _Li and specific area _ai can be calculated on the basis of a reference value _i0,ref [11]:

_i0,i=_{i0,ref,i ai Li} ^(_pipref)γexp[−Δ_GiR_Ti(1−_TiTref)],

_Tref=298.15 K; _pref=1.0125 bar.

For simplicity, a constant value was used for the activation energy Δ_Gi , even though it can vary under real operating conditions [14].

2.3. Membrane Water Content

For estimating the water content λ inside of a Nafion-membrane, a function of water activity a is used [15]:

λ={0.043+17.81 a−39.85 ^a2+36 ^a3fora≤114+1.4(a−1)for1<a≤3.

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]

a=RH+2 s,

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.

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]:

_λsat{=4.837 for T<223.15 K=[−1.304+0.01479 T^{−3.594·10−5 T2]−1} for 223.15 K≤T<_Tfrost>λ for T≥_Tfrost.

To estimate the water concentration _Cw , a simple proportional correlation between the membrane density _ρmem , equivalent weight and water content λ is used [11]:

_Cw=_ρmemEWλ.

2.4. 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]:

_Eohm=_δmem i_σmem.

To estimate the membrane conductivity _σmem , the following correlation is used [9,12,19]:

_σmem=1.16 max^{{0,f−0.06}1.5}exp[15000R(1_Tref−1T)],

_Tref=353.15 K; f=λ_VWλ_VW+_Vm,_Vm=EW_ρmem; _VW=18.01528_ρW(T).

It is mostly dependent on membrane water content λ and temperature _Tmem , 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]:

_ρw(T)=999.972−7·^{10−3(T−4)2}·¹⁰⁻³,

T in °C.

2.5. Concentration Overvoltage

Using the Nernst-equation, concentration overvoltage _Econi is described as a function of temperature T , current density i and limiting current density _iL [11]

_Econi=R_TinFln(_iLiiLi−i),

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].

For the limiting current density _iL , a simplified expression including the Faraday F constant, number of electrons n , diffusion coefficient _Di , gas concentration _Ci and diffusion distance (here: electrode thickness) _δe is used [11]:

_iLi=nF_{Di Ciδe}.

To calculate the diffusion coefficient _Di , the model expects an external reference value _Di,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]:

_Di=ϵ^τ2(1−s)3 _Di,ref ^(T_Tref)1.5_prefp,

_Tref=353.15 K; _pref=1.01325 bar.

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 _Di,ref as a model input for the calculation of _Di 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 _Clm . If desired, this can be disabled by supplying the same value for input and output [5]:

_Clm=_Cin−_Coutln_CinCout.

2.6. Cell and Stack Performance

What is labeled as the electrical efficiency _ηelectric in this work relates to the lower heating value of hydrogen _ELHV [11]:

_ηelectric=_EcellELHV,_ELHV=1.254 V,

and, furthermore, is used to calculate the heat flow Q˙ of the cell/stack based on the power delivered _Pstack . As a simplification, it is assumed that the product water fully evaporates before leaving the cell/stack [11]:

_Q˙stack=(_{Pstackμelectric})−_Pstack.

2.7. Water Transport

In this model, the estimation of the flow of water _jw from anode to cathode is divided in three categories: osmotic _josmo and diffusive _jdiff flow as well as hydraulic permeation _jhyd . A positive value means an increase in water concentration:

_jwanode=_jdiff+_jhyd−_josmo,

_jwcathode=_jgen+_josm−_jdiff−_jhyd.

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).

2.7.1. Osmosis

Water transport due to osmosis _josmo is expressed by a linear correlation composed of the Faraday constant F , the current I [11]

_josmo=_nosmoIF

and the osmotic coefficient _nosmo , which is dependent on the (average) water content of the membrane λ [25]:

_nosmo=0.0029 ^λ2+0.05 λ−3.4·¹⁰⁻¹⁹.

2.7.2. Diffusion

In terms of diffusion _jdiff , the abovementioned functions to differentiate between the liquid and gaseous phase are utilized. The basis for these calculations is given by [17]

_jdiff=A_Dλ∇_Cwδmem,

which considers the cell area A , (average) diffusive coefficient of water _Dλ , water concentration gradient ∇_Cw —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 _Tmem and (average) water content λ is applied [9,26]:

_Dλ=3.842 ^λ3−32.03 ^λ2−67.74 λ^λ3−2.115 ^λ2−33.013 λ+103.37·¹⁰⁻⁶exp[20·20000R(1_Tref−1_Tmem)],_Tref=353.15K.

Subsequently, the diffusive flow _jdiff 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:

_jdiffVVP=0.9178·_jdiff(_{δmem |0.0201})(−947.5 _δmem2−6.198 _δmem+1.508),

_jdiffLVP=3.592·_jdiff(_{δmem |0.0201})(−687 _δmem2−21.73 _δmem+1.714).

2.7.3. 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 _jhyd 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 _μH2O as well as water concentration inside the membrane _Cw , its thickness _δmem and hydraulic permeability _Kλ [16]:

_jhyd=A_{Cw KλμH2O δmem}∇p·¹⁰⁵.

For the hydraulic permeability _Kλ , a direct dependency on the membrane water content λ is assumed [16]:

_Kλ=_Kwλ.

Furthermore, the dynamic viscosity of water is approximated by a function of temperature T [28]

_μw=_μ0exp[_aμp+_dμ−_bμpR(T−_θμ−_cμp)],

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] :

_KλLLP{=0.1158 _Kλ(5.749·¹⁰⁻³ _δmemexp[−1.326]) for _δmem≥0.0056cm=0.1158 _Kλ(2.518·¹⁰⁻⁴ _δmemexp[−1.872]) for _δmem<0.0056cm.

Subsequently, the hydraulic flow is corrected for the actual pressure difference:

_jhydLLP=_jhydref(32.41 Δp+0.06016).

3. 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.

4. 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.

5. 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.

[Image omitted. See PDF.]

[Image omitted. See PDF.]); 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."]

[Image omitted. See PDF.]

[Image omitted. See PDF.]; (d) O2 concentration [mol/cm3]."]

Name/type	Value	Unit
Input
Temperature	70	[°C]
Pressure (absolute)	1.013	[bar]
Relative humidity cat/an	96/38	[%]
O2 concentration	1.2 × 10−6	[mol/cm3]
H2 concentration	5 × 10−5	[mol/cm3]
Liquid water volume fraction	0	[%]
O2 diffusive reference	0.36	[mol/cm2 s]
H2 diffusive reference	1.24	[mol/cm2 s]
Parameter
Membrane thickness	201	[µm]
Membrane density	1.97	[g/cm3]
Membrane EW	1020	[g/mol]

Supplementary Materials

The supplementary materials are available online at https://www.mdpi.com/1996-1073/12/18/3478/s1.

Author Contributions

Formal analysis, investigation, methodology, software, and visualization, writing-original draft, A.L.L.; data curation, S.C.K.; supervision, resources, H.R.; validation, A.L.L. and S.C.K.; writing-review and editing, S.C.K. and H.R.; conceptualization, project administration, A.L.L., S.C.K. and H.R.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

a water activity

_ai catalyst-specific area, cm2/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

_Ci molar concentration, mol/cm3

_Cw membrane water concentration, mol/cm3

_Di diffusion coefficient, cm2/s

_Dλ diffusion coefficient of water through the membrane, cm2/s

_dμ constant in the calculation of _μw , J/mol

_Ei voltage, V

EW membrane equivalent weight, g/mol

f membrane liquid water volume fraction

F Faraday constant, A s/mol

ΔG Gibbs free energy, J/mol

Δ_Gi activation energy, J/mol

H enthalpy, J

i current density, A/cm2

_i0 exchange current density, A/cm2

_i0,ref reference exchange current density, A/cm2

_iL limiting current density, A/cm2

_ji water flow, mol/s

_Khyd hydraulic permeability of the membrane (liquid water), cm2

_Kw hydraulic permeability coefficient (liquid water)

_Li catalyst loading, mg/cm2

_Ni molar flux, mol/cm2 s

n number of electrodes involved

_ni water transport coefficient

P power, W

p pressure, bar

_pi partial pressure, bar

_pref reference pressure, bar

Q˙ heat flow, W

R universal gas constant, J/mol K

s liquid water volume fraction

S entropy, J/K

_Ti temperature, K

_Tref reference temperature, K

_Vm acid equivalent volume of the membrane, cm3/mol

_VW molar volume of water, cm3/mol

_xi percentage

Greek Letters

_αi charge transfer coefficient

γ pressure dependency coefficient

_δi thickness, cm

ϵ electrode porosity

_ηi efficiency

λ membrane water content

_μw,0 reference water dynamic viscosity, Pa s

_μw water dynamic viscosity, Pa s

_ρi density, g/cm3

_σi conductivity, A/V cm

τ electrode tortuosity

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

References

1. Harvey, D.B.; Pharaoh, J.G.; Karan, K. FAST-FC: Open source fuel cell model: An Open Source Performance and Degradation Fuel Cell Model for PEMFCs. Available online: https://www.fastsimulations.com/source-code.html (accessed on 3 September 2019).

2. Falcão, D.S.; Pinho, C.; Pinto, A.M.F.R. Water management in PEMFC: 1-D model simulations. Ciênc. Tecnol. dos Mater. 2016, 28, 81-87.

3. Secanell, M.; Putz, A.; Wardlaw, P.; Zingan, V.; Bhaiya, M.; Moore, M.; Zhou, J.; Balen, C.; Domican, K. OpenFCST: An Open-Source Mathematical Modelling Software for Polymer Electrolyte Fuel Cells. ECS Trans. 2014, 64, 655-680.

4. Falcão, D.S.; Gomes, P.J.; Oliveira, V.B.; Pinho, C.; Pinto, A.M.F.R. 1D and 3D numerical simulations in PEM fuel cells. Int. J. Hydrog. Energy 2011, 36, 12486-12498.

5. Spiegel, C. PEM Fuel Cell Modeling and Simulation Using Matlab; Elsevier: Amsterdam, The Netherlands, 2008; ISBN 978-0-12-374259-9.

6. Chakraborty, U.K.; Abbott, T.E.; Das, S.K. PEM fuel cell modeling using differential evolution. Energy 2012, 40, 387-399.

7. Conti, B.; Bosio, B.; McPhail, S.J.; Santoni, F.; Pumiglia, D.; Arato, E. A 2-D model for Intermediate Temperature Solid Oxide Fuel Cells Preliminarily Validated on Local Values. Catalysts 2019, 9, 36.

8. Chakraborty, U.K. A New Model for Constant Fuel Utilization and Constant Fuel Flow in Fuel Cells. Appl. Sci. 2019, 9, 1066.

9. Vetter, R.; Schumacher, J.O. Free open reference implementation of a two-phase PEM fuel cell model. Comput. Phys. Commun. 2019, 234, 223-234.

10. Chakraborty, U.K. Reversible and Irreversible Potentials and an Inaccuracy in Popular Models in the Fuel Cell Literature. Energies 2018, 11, 1851.

11. Barbir, F. PEM Fuel Cells. Theory and Practice, 2nd ed.; Academic Press: Waltham, MA, USA, 2013; ISBN 9780123877109.

12. Yigit, T.; Selamet, O.F. Mathematical modeling and dynamic Simulink simulation of high-pressure PEM electrolyzer system. Int. J. Hydrog. Energy 2016, 41, 13901-13914.

13. Carmo, M.; Fritz, D.L.; Mergel, J.; Stolten, D. A comprehensive review on PEM water electrolysis. Int. J. Hydrog. Energy 2013, 38, 4901-4934.

14. Song, C.; Tang, Y.; Zhang, J.L.; Zhang, J.; Wang, H.; Shen, J.; McDermid, S.; Li, J.; Kozak, P. PEM fuel cell reaction kinetics in the temperature range of 23-120°C. Electrochim. Acta 2007, 52, 2552-2561.

15. Springer, T.E. Polymer Electrolyte Fuel Cell Model. J. Electrochem. Soc. 1991, 138, 2334.

16. Jiao, K.; Li, X. Water transport in polymer electrolyte membrane fuel cells. Prog. Energy Combust. Sci. 2011, 37, 221-291.

17. Abdin, Z.; Webb, C.J.; Gray, E.M. PEM fuel cell model and simulation in Matlab-Simulink based on physical parameters. Energy 2016, 116, 1131-1144.

18. Chemours. Chemours Product-Bulletin Nafion N115-N117-N1110. Available online: https://nafionstore-us.americommerce.com/Shared/Bulletins/N115-N117-N1110-Product-Bulletin-Chemours.pdf (accessed on 13 May 2019).

19. Weber, A.Z.; Newman, J. Transport in Polymer-Electrolyte Membranes. J. Electrochem. Soc. 2003, 150, A1008-A1015.

20. Kersting, A.; Aeschbach-Hertig, W. Vorlesungen Physik aquatischer Systeme I. Vorlesung. Available online: http://www.iup.uni-heidelberg.de/institut/studium/lehre/AquaPhys/docMVEnv3_11/MVEnv3_2011_p1_2_Dichte.pdf (accessed on 13 May 2019).

21. Musio, F.; Tacchi, F.; Omati, L.; Gallo Stampino, P.; Dotelli, G.; Limonta, S.; Brivio, D.; Grassini, P. PEMFC system simulation in MATLAB-Simulink® environment. Int. J. Hydrog. Energy 2011, 36, 8045-8052.

22. Rosén, T.; Eller, J.; Kang, J.; Prasianakis, N.I.; Mantzaras, J.; Büchi, F.N. Saturation Dependent Effective Transport Properties of PEFC Gas Diffusion Layers. J. Electrochem. Soc. 2012, 159, F536-F544.

23. Bird, R.B. Transport phenomena. Appl. Mech. Rev. 2002, 55, R1-R4.

24. Adachi, M.; Navessin, T.; Xie, Z.; Li, F.H.; Tanaka, S.; Holdcroft, S. Thickness dependence of water permeation through proton exchange membranes. J. Membr. Sci. 2010, 364, 183-193.

25. Dutta, S.; Shimpalee, S.; van Zee, J.W. Numerical prediction of mass-exchange between cathode and anode channels in a PEM fuel cell. Int. J. Heat Mass Transf. 2001, 44, 2029-2042.

26. Mittelsteadt, C.K.; Staser, J. Simultaneous Water Uptake, Diffusivity and Permeability Measurement of Perfluorinated Sulfonic Acid Polymer Electrolyte Membranes. ECS Trans. 2011, 41, 101-121.

27. Karpenko-Jereb, L.; Innerwinkler, P.; Kelterer, A.-M.; Sternig, C.; Fink, C.; Prenninger, P.; Tatschl, R. A novel membrane transport model for polymer electrolyte fuel cell simulations. Int. J. Hydrog. Energy 2014, 39, 7077-7088.

28. Likhachev, E.R. Dependence of water viscosity on temperature and pressure. Tech. Phys. 2003, 48, 514-515.