A Control-Oriented Model for Predicting Variations in Membrane Water Content of an Open-Cathode Proton Exchange Membrane Fuel Cell

1. Introduction

Due to the negative environmental impacts of fossil fuels, alternatives are being sought and developed particularly in the transportation sector where reliance on petroleum-based fuels is high. While electrification of light-duty vehicles is growing, medium and heavy-duty vehicles can be harder to electrify due to their high energy demands. For medium to heavy-duty vehicles, proton exchange membrane (PEM) fuel cells have received considerable attention. Fuel cells provide an electrification solution for vehicles that alleviates range anxiety due to their refuelling capability. They can also function with a better efficiency than internal combustion engines and produce little to no pollutants on board. Among various fuel cell types, the PEM fuel cell stands out as a popular choice, primarily due to its high-power density and low operating temperature.

There are two types of PEM fuel cells: the closed-cathode and open-cathode varieties. The closed-cathode design is usually equipped with a compressor, a manifold, a dedicated cooling circuit for temperature management and an external humidifier for water management. These auxiliary components in the closed-cathode design induce parasitic loads that impact the stack performance. In contrast, the open-cathode design reduces the system’s complexity by integrating the stack’s air supply and temperature control and consequently water management.

Reliable and sustainable fuel cell operation requires effective water management as the membrane humidity has been shown experimentally to influence the cell performance. (Zeng et al., 2019) investigated the performance of an open cathode PEM fuel cell under variable speed control. They established that strong air flow rate leads to a decrease in cell voltage due to excessive loss in membrane water content. Further (Morner & Klein, 2001) evaluated the dynamic behavior of an air-breathing fuel cell stack experimentally and concluded that the performance of the stack is more dependent on membrane humidity than the stack temperature.

The membrane water content needs to be in the right range to promote ionic conductivity but must not interfere with reactant transport. Due to its dependence on temperature and reactant inlet humidity, control of membrane water content is inherently challenging. Internal cell variations that occur during transients can also lead to losses in fuel cell performance. To guarantee the stack’s performance and continued operation, these variations must be carefully monitored and controlled. As the fuel cell’s internal humidity status cannot be measured, control strategies depend on models to predict the internal humidity dynamics.

There are a number of existing models that aim to predict humidity distribution in a fuel cell stack. For instance, (Zhang & Jiao, 2018) developed a three-dimensional (3D) mathematical model of a PEM fuel cell. The established model can describe the relative humidity distribution of the gases in the anode and cathode channels as well as the membrane water content. (Dutta et al., 2001) established a 3D model using Navier Stoke’s equations with a multi-species mixture. They found that species mass flow direction significantly depends on mass consumption on the membrane electrolyte assembly. To study the water and thermal behaviour of an open-cathode PEM fuel cell, (Sagar et al., 2020) developed a 3D CFD model to predict the spatial distribution of relative humidity, temperature and membrane water content. They concluded that high water content variations in the fuel cell limits the cell performance. The studies provide insight into the fact there are distributions in the fuel cell stack that impact performance and hence need to be controlled. However, these computationally fluid dynamics (CFD)-based models are not suitable for control purposes as they are far too computationally expensive.

A number of control-oriented models have been created over the years. For instance, (Pukrushpan et al., 2004) developed a control-oriented lumped-parameter model for automotive fuel cell systems. The model was proven to be useful for analysis of transient effect of step inputs and analysis of system observability. (Meyer & Yao, 2006) developed a control-oriented model for a self-humidifying fuel cell stack and the model was capable of characterizing the transient response of the fuel cell stack. (X. Chen et al., 2022) developed a physics-based model to control the temperature and humidity of a PEM fuel cell. To study and control water in a fuel cell stack, (X. Chen et al., 2020) developed a dynamic water management model based on the equations developed by (Pukrushpan et al., 2004). (F. Chen et al., 2021) developed a lumped parameter control-oriented model for predicting humidity dynamics for an open-cathode PEM fuel cell. Based on their model, the output performance of an open-cathode PEMFC, taking into account the internal humidity state of the stack, is predicted. All these models adopt the lumped-parameter approach which assumes minimal spatial variations. However, the CFD models reviewed demonstrate that variations exist; therefore, control-oriented models must account for these variations to ensure effective control of the membrane water content in the PEM fuel cell stack.

For predicting humidity spatial variations, (Headley et al., 2016) developed and experimentally validated a control-oriented model that captures variation of relative humidity along the cathode channel but did not explore humidity variation in the anode channel, which can occur due to back diffusion. It is not possible to accurately capture variations in membrane water content with this approach as the membrane water content depends on humidity dynamics in both the cathode and anode channels. In fact, a more conservative approach of determining the membrane water content is to model it as a function of the water activity in the anode alone. (Nguyen & White, 1993) adopted this approach in their work where the membrane water content was taken to be the anode water content to reflect the importance of the anode side water activity as it is often the limiting electrode due to drying conditions. Therefore, for more accurate predictions of variations in membrane water content, a model that predicts both variations of humidity along the cathode and anode channels is needed.

The main contribution of this paper is that it overcomes the limitation of existing control-oriented models for water management by developing a model that predicts humidity variations along both the cathode and anode channels, with the goal of providing more accurate insights into the humidity dynamics in the PEM fuel cell membrane. In this study, a reduced order model that integrates thermal effects to predict spatial variations in membrane water content of a 5kW open-cathode PEM fuel cell stack is presented. First, a lumped parameter model for the stack was developed. This model was then further discretized into smaller control volumes. The model was validated by comparing its predictions to simulation results from GT-Suite and achieved good agreement. The developed model will be useful for real-time control of membrane water content variation in PEM fuel cells.

2. Experimental Setup

The model is based on the experimental set up shown in Figure 1. The fuel cell stack is a commercially available 5kW open-cathode stack from Horizon. It has 120 cells connected in series. The active area of each cell is 150cm². Air is supplied to the stack by four blowers connected to the casing of the stack. Hydrogen of 99.99% purity is supplied to the stack from a pressurized tank. The stack is equipped with an in-built controller from the manufacturer that monitors and records the current, voltage and stack temperature. A mass flow controller from Aalborg is used to measure the hydrogen inlet pressure and control the supply of hydrogen to the anode. There are also four thermocouples which are attached on the other side of the stack casing (opposite the fans) with the aim of monitoring temperature variation along the stack (Figure 2). Data taken on this experimental setup was used for model validation.

3. Model Development

3.1. Physics-Based Model

The model is made up of four interconnected sub models. An air flow model computes the mass flow rate of air (${\dot{m}}_{air}$) from the fan’s PWM signal. The air flow rate is then fed into both the thermal and water management models. The thermal model generates the stack temperature $\left(T_{st}\right)$ ,which is used as an input to the water management and stack voltage models. The water management model’s prediction of membrane water content $\left({\lambda }_{mem}\right)$ then serves as an input to the stack voltage model. Figure 3 depicts the model’s structure.

3.1.1. Stack Voltage Model

The fuel cell reaction is split into the cathodic and anodic half reactions. To facilitate these reactions, air is supplied to the cathode and hydrogen is supplied to the anode. The cathode and the anode are separated by the polymer electrolyte membrane, which is permeable to protons but not electrons, so electrons are transferred via an electrical circuit generating electric current. The electrons then recombine with the protons in the cathode side to generate water. Due to the transfer of electrons, the anodic reaction is referred to as the hydrogen oxidation reaction (HOR) and the cathodic reaction is called the oxygen reduction reaction (ORR) due to the addition of electrons. The electrochemical reactions involved are:

$${\mathrm{H}}_2\to 2{\mathrm{H}}^{+}+2\mathrm{e}\left(\mathrm{H}\mathrm{O}\mathrm{R},\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}\right)$$

$$2{\mathrm{H}}^{+}+\frac{1}{2}{\mathrm{O}}_2+2\mathrm{e}\to {\mathrm{H}}_2\mathrm{O}(\mathrm{O}\mathrm{R}\mathrm{R},\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{h}\mathrm{o}\mathrm{d}\mathrm{e})$$

(1)

$${\mathrm{H}}_2+\frac{1}{2}{\mathrm{O}}_2\to {\mathrm{H}}_2\mathrm{O}(\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n})$$

The voltage of a single cell is expressed by

$$V_{cell} \left(V\right)=V_{oc}-V_{act}-V_{ohmic}-V_{conc}$$

(2)

where $V_{oc}$ is the open circuit voltage, $V_{act }$ is the activation overvoltage, $V_{ohmic}$ is the ohmic overvoltage and $V_{conc}$ is the concentration overvoltage.

$V_{oc}$ is the thermodynamic reversible potential or the fuel cell voltage when no load is connected. The reversible fuel cell voltage can be expressed as a function of temperature and pressure as follows (O’hayre et al, 2014):

$$V_{oc}=E^o+\frac{∆s_{rxn}\left(T_{st}\right)}{nF}\left(T_{st}-T_{amb}\right)-\frac{R{T}_{st}}{nF}In\left(\frac{1_{ }}{P{H}_2\ast P{O}_2^{0.5}}\right)$$

(3)

in which $E^o$ is the non-standard reversible voltage, $∆s_{rxn}$ is the entropy change of the fuel cell reaction, $T_{st}$ is the stack temperature, $T_{amb}$ is the ambient temperature, $R$ is the ideal gas constant, $n$ is the number of moles of electrons transferred from the anodic reaction, $F$is the Faraday constant, $P{H}_2$ is the hydrogen inlet pressure and $P{O}_2$ is the oxygen partial pressure.

The activation overvoltage includes losses in voltage that arise due to the kinetic reaction. The activation overvoltage can be determined by

$$V_{act}=\begin{array}{c}\frac{R{T}_{st}}{anF}In \left(\frac{i}{i_o}\right)\end{array}$$

(4)

where $a$ is the charge transfer coefficient, $i$ is the cell current density and $i_o$ is the exchange current density.

The exchange current density can be computed by (F. Chen et al., 2021)

$$i_o=i_o^{ref}{\left(\frac{P{O}_2}{P_{ref}}\right)}^{y}exp\left[\frac{i}{i_o}\left(-\frac{E_c}{R}-\frac{1}{T_{ref}}\right)\right]$$

(5)

where $i_o^{ref}$ is the reference exchange current density, $P_{ref}$ is the reference pressure of 1atm, $y$ is the pressure coefficient, $T_{ref}$ is the reference temperature of 298K and $E_c$ is the activation energy for the cathodic reaction. Due to the slower reaction kinetics of the ORR as compared to the HOR, the activation loss in this model is only attributed to the cathodic reaction.

Ohmic overvoltage captures the losses due to ionic and electronic conduction and can be calculated by

$$V_{ohmic}=i\left(R_{ionic}+R_{elec}\right)$$

(6)

where $R_{ionic}$ and $R_{elec}$ are the ionic resistance and electronic resistance respectively. $R_{ionic}$ can be computed by (Wang et al., 2014) (Huo & Hall, 2023)

$$R_{ionic}=\frac{181.6[1+0.03i+0.062{\left(\frac{T_{st}}{303}\right)}^2\left(i^{2.5}\right)]t_{mem}}{\left[{\lambda }_{mem}-0.634-3i\right]\exp \left[4.18\left(T_{st}-\frac{303}{T_{st}}\right)\right]}$$

(7)

in which $t_{mem}$ is the membrane thickness and ${\lambda }_{mem}$is the membrane water content and is determined from the water management model as discussed in Section 3.1.4.

The concentration overvoltage includes losses attributed to mass transport in the fuel cell stack. The concentration overvoltage can be determined by

$$V_{conc}=mexp\left(ni\right)$$

(8)

where $m$ and $n$are empirical coefficients that can be tuned according to experimental polarization data.

The total stack voltage is then computed by

$$V_{st}=V_{cell}\ast N_{st}$$

(9)

in which $V_{cell}$ is the voltage of a single cell as computed from Eq. (2) and $N_{st}$ is the total number of cells in the stack.

3.1.2. Air Flow Model

Air is delivered to the cathode channels via fans mounted on the stack casing. The flow to the cathode channel is assumed to be laminar in this model. The cathode air mass flow rate is modeled as

$${\dot{m}}_{air}\left(t\right)={\rho }_{air}\ast Q_{air}\left(t\right)$$

(10)

where ${\rho }_{air}$ is the density of air in kg/m³ at Standard Temperature and Pressure (STP) and $Q_{air}\left(t\right)$ is the volumetric air flow rate in m³/s.

For an arbitrary fan rotational speed, $Q_{air}\left(t\right)$ is determined by (Huo & Hall, 2023)

$$Q_{air}\left(t\right)=w\left(t\right)\frac{Q_{nom}}{w_{nom}}$$

(11)

in which $w\left(t\right)$ denotes the arbitrary fan speed in revolutions per minute (rpm), $Q_{nom}$ is the nominal volumetric flow rate of the fan in m³/s and $w_{nom}$ is the nominal fan speed in rpm. $Q_{nom}$ and $w_{nom}$ are determined in a manner as described by (Ishaku et al., 2014). The pressure drop in the cathode channel is calculated as a function of the Reynolds number and friction factor. $Q_{nom}$ is then determined from the intersection of the determined pressure drop on the fan performance curve. $w_{nom}$ is the speed at which the performance curve is expressed in the datasheet. For this study $w_{nom}$ is 6500rpm.

3.1.3. Thermal Model

The thermal dynamics of the stack can be described as follows (Huo & Hall, 2023) (Ishaku et al., 2014)

$$C_t\frac{d{T}_{st}}{dt}=P_{total}\left(t\right)-P_{st}\left(t\right)-{\dot{Q}}_{coolant}(t)$$

(12)

where $C_t$ represents the thermal capacitance of the stack, $P_{total}\left(t\right)$ is the power from the electrochemical reaction and is computed by

$$P_{total}\left(t\right)=\frac{N_{st}\ast I_{st}\left(t\right)\ast ∆h}{2F}$$

(13)

in which $∆h$ is the molar enthalpy change of the electrochemical reaction. $P_{st}\left(t\right)$ is the total output power from the stack and is determined by

$$P_{st}\left(t\right)=V_{st}\left(t\right)\ast I_{st}\left(t\right)$$

(14)

${\dot{Q}}_{coolant}\left(t\right)$ is the heat dissipated from the air to the environment and it is calculated by

$${\dot{Q}}_{coolant}\left(t\right)={\eta }_{fan}\ast {\dot{m}}_{air}\left(t\right)\ast C_p(T_{st}\left(t\right)-T_{amb})$$

(15)

where ${\eta }_{fan}$ is the efficiency of the blower, ${\dot{m}}_{air}\left(t\right)$ is the air mass flow rate as determined from Section 3.1.2, $C_p$ is the specific heat capacity of air and $T_{amb}$ is the ambient temperature.

3.1.4. Water Management Model

The water management model is divided into three submodels: the cathode channel, the anode channel and the polymer electrolyte membrane. A block diagram of the model is shown in Figure 4.

Cathode Channel

Conservation of mass was utilized to track the mass of species in the cathode channel over time. The conservation of mass of water, oxygen and nitrogen are given by Eqs. (16), (17) and (18), respectively (Pukrushpan et al., 2004). It is assumed that water enters and leaves the channel in vapor form only.

$$\frac{d{m}_{H_{2}Oca}}{dt }={\dot{m}}_{H_2{O_{in},ca}_{ }}+{\dot{m}}_{H_2{O}_{gen}}-{\dot{m}}_{H_2{O_{out},ca}_{ }}+{\dot{m}}_{mem}$$

(16)

$$\frac{d{m}_{O_2}}{dt }={\dot{m}}_{{O_{2,in}}_{ }}-{\dot{m}}_{{O_{2,reacted}}_{ }}-+{\dot{m}}_{{O_{2,out}}_{ }}$$

(17)

$$\frac{d{m}_{N_2}}{dt }={\dot{m}}_{{N_{2,in}}_{ }}-{\dot{m}}_{{N_{2,out}}_{ }}$$

(18)

In these equations, ${\dot{m}}_{H_2{O_{in},ca}_{ }}$ ${\dot{m}}_{{O_{2,in}}_{ }}$and ${\dot{m}}_{{N_{2,in}}_{ }}$ are the mass flow rates of water, oxygen, and nitrogen into the cathode channel, computed by Eqs. (24), (30) and (31), respectively. The mass flow rate of water vapor, oxygen and nitrogen leaving the cathode channel are depicted by ${\dot{m}}_{H_2{O_{out},ca}_{ }}$, ${\dot{m}}_{{O_2,out}_{ }}$and ${\dot{m}}_{{N_{2,out}}_{ }}$ and are calculated by Eqs. (33), (43) and (44), respectively. The net flow of water across the membrane, ${\dot{m}}_{mem}$ is determined from the membrane model in Section 3.1.4.3.

The water generated due to the electrochemical reaction is given by

$${\dot{m}}_{H2O_{gen}}=\frac{N_{st}\ast I_{st}\ast {M_{H_{2}O}}_{ } }{2F}$$

(19)

where $M_{H_{2}O}$ is the molar mass of water in kg/mol. The mass flow rate of oxygen that partakes in the electrochemical reaction computed by

$${\dot{m}}_{{O_{2,reacted}}_{ }}=\frac{N_{st}\ast I_{st}\ast {M_O}_2 }{2F}$$

(20)

in which ${M_O}_2$ is the molar mass of oxygen in kg/mol.

Anode Channel

Similar to the cathode channel, conservation of mass was utilized to track the mass of water and hydrogen in the anode channel. The anode of the self-humidifying open-cathode PEMFC under consideration is fed with dry hydrogen and it is typically dead-ended except for the occasional purging of water and unused hydrogen by the purge valve. It is assumed in this model that the anode runs exclusively in the dead-ended mode.

The mass balance of water is determined by

$$\frac{d{m}_{H2Oan}}{dt }=-{\dot{m}}_{mem}$$

(21)

The mass balance of hydrogen is given by

$$\frac{d{m}_{H_2}}{dt }={\dot{m}}_{{{H_2}_{in,an}}_{ }}+{\dot{m}}_{{{H_2}_{out,an}}_{ }}-{\dot{m}}_{{H_2}_{reacted}}$$

(22)

The mass rate of hydrogen reacted is determined by

$${\dot{m}}_{{H_2}_{reacted}}=\frac{N_{st}\ast I_{st}\ast {M_H}_2 }{4F}$$

(23)

The mass flow rate of water entering the channels is given by

$${\dot{m}}_{H_2{O_{in}\left(.\right)=y_{H_2{Ov}_{in}(.)}\ast {\dot{m}}_{in(.)}}_{ }}$$

(24)

where $(.)$ represents either the cathode (ca) or anode (an) channel. The term ${\mathrm{y}}_{H_2{Ov}_{in}}$is the mass fraction of vapor in the inlet flow and ${\dot{m}}_{in}$ is the inlet flow rate. For the cathode channel, ${\dot{m}}_{in,ca}$ is determined from the air flow rate model in Section 3.1.2 and ${\dot{m}}_{in,an}$ is a measured input to the model.

The mass fraction of vapor into the channel is computed by (Karnik et al., 2007)

$$\begin{array}{c}{\mathrm{y}}_{H_2{Ov}_{in}\left(.\right)}=\hfill \\ \frac{P_{H_{2}O{v}_{in}(.)}\ast M_{H_{2}Ov}}{\left(P_{H_{2}O{v}_{in}(.)}\ast M_{H_{2}Ov}\right)+\left(\right({P_{(.)}}_{in}-P_{H_{2}O{v}_{in}(.)})\ast M_g)}\hfill \end{array}$$

(25)

The subscript, $g$ represents the gas in the channel. For the cathode, the gas is air, and it is hydrogen for the anode.

The inlet vapor pressure is calculated by

$$P_{H_{2}O{v}_{in}(.)}=R{H_{(.)}}_{in}\ast P_{sat}\left({T_{(.)}}_{in}\right)$$

(26)

where $R{H_{(.)}}_{in}$ is the channel inlet relative humidity and $P_{sat}\left({T_{(.)}}_{in}\right)$ is the saturation pressure evaluated at the channel inlet temperature.

${P_{(.)}}_{in}$ is the channel inlet pressure and $M_g$ is the molar mass of gas in the channel. Once the mass flow rate of water in the inlet stream has been determined, the inlet flow rate of the gas is determined by subtracting the water mass flow rate from the total inlet flow rate.

$${\dot{m}}_{g_{in(.)}}={\dot{m}}_{in(.)-}{\dot{m}}_{H_2{O}_{in,(.)}}$$

(27)

To determine the mass flow rate of each of the nitrogen and oxygen species in the cathode, the mass fraction of each of the gases needs to be computed and this is given by (Pukrushpan et al., 2004)

$${x_O}_{2}c{a}_{in}=\frac{{y_O}_{2}c{a}_{in}\ast {M_O}_2}{({y_O}_{2}c{a}_{in}\ast {M_O}_2)+(\left(1-{y_O}_{2}c{a}_{in}\right)\ast {M_N}_2)}$$

(28)

where ${y_O}_{2}c{a}_{in}$ is the mole fraction of the oxygen in the inlet air which is 0.21.

The mass fraction of nitrogen is then given by

$${x_N}_{2}c{a}_{in}=1-{x_O}_{2}c{a}_{in}$$

(29)

The mass flow rate of oxygen into the cathode is

$${\dot{m}}_{{O_2}_{in.ca}}={x_O}_{2}c{a}_{in}\ast {\dot{m}}_{{air}_{in,ca}}$$

(30)

and the mass flow rate of nitrogen is

$${\dot{m}}_{{N_2}_{in,ca}}={x_N}_{2}c{a}_{in}\ast {\dot{m}}_{{air}_{in,ca}}$$

(31)

The outlet flow rate is determined from the linearized nozzle equation (Karnik et al., 2007)

$${\dot{m}}_{ou{t}_{(.)}}={k_{\left(.\right)}}_{out}(P_{\left(.\right)}-P_{ds})$$

(32)

in which ${k_{\left(.\right)}}_{out}$ is the outlet anode or cathode nozzle constant coefficient, $P_{\left(.\right)}$ is the anode or cathode channel pressure and $P_{ds}$ is the pressure downstream of the stack and for the lumped model this is equal to standard atmospheric pressure.

Knowing the outlet flow rate, the mass flow rate of vapor out can be computed from

$${\dot{m}}_{H_2{O}_{out}(.)}={\mathrm{y}}_{H_{2}Ov(.)}\ast {\dot{m}}_{out(.)}$$

(33)

The mass fraction of water vapor is given by

$${\mathrm{y}}_{H2Ov(.)}=\frac{m_{H_{2}Ov(.)}}{m_{H_{2}Ov(.)}+m_g\left(.\right)}$$

(34)

where $m_{H_{2}Ov(.)}$ and $m_g\left(.\right)$are the masses of vapor and gas respectively in the cathode or anode channel. The total mass of gas in the cathode channel is obtained by summing the masses of nitrogen and oxygen.

Based on the calculated quantity of water exiting, the amount of gas out of the channels can be computed by

$${\dot{m}}_{{g_{out,(.)}}_{ }}={\dot{m}}_{ou{t}_{(.)}}-{\dot{m}}_{H_2{O}_{out}(.)}$$

(35)

To determine the flow rate of the individual species of nitrogen and oxygen out of the cathode, the mass fractions have to be recomputed as oxygen is consumed in the reaction.

The mass fraction of oxygen leaving the cathode is (Pukrushpan et al., 2004)

$${x_O}_2=\frac{{y_O}_2\ast {M_O}_2}{({y_O}_2\ast {M_O}_2)+(\left(1-{y_O}_2\right)\ast {M_N}_2)}$$

(36)

where ${y_O}_2$ is the mole fraction of oxygen in the exit flow rate and it is determined by

$${y_O}_2=\frac{{P_O}_2}{Pai{r}_{ }}$$

(37)

Assuming the gases to be ideal, the partial pressures of oxygen and nitrogen in the cathode channel are determined from the ideal gas law as

$${P_O}_2=\frac{m_{O_2}\ast R_{O_2}\ast T_{ca} }{V_{ca}}$$

(38)

$${P_N}_2=\frac{m_{N_2}\ast R_{N_2}\ast T_{ca} }{V_{ca}}$$

(39)

where $m_{O_{2}ca}$ and $m_{N_{2}ca}$ are the masses of oxygen and nitrogen obtained by solving Eqs. (18) and (19) respectively. $R_{O_2}$ is the gas constant of oxygen and $R_{N_2}$ is the gas constant of nitrogen. $T_{ca}$ is the cathode channel temperature and it is taken to be the same as the stack temperature $\left(T_{st}\right)$ in this study. $V_{ca}$ is the lumped volume of the cathode channel.

Similarly, the partial pressure of hydrogen in the anode channel is calculated by

$${P_H}_2=\frac{m_{H_2}\ast R_{H_2}\ast T_{an} }{V_{an}}$$

(40)

where $m_{H_2}$ is the mass of hydrogen found in Eq. (22), $R_{H_2}$ is hydrogen gas constant, $T_{an}$ is the anode temperature, which is equal to $T_{st}$ in this model and $V_{an}$ is the lumped anode channel volume.

The partial pressure of air in the channel is then given by the summation of Eqs. (38) and (39) as

$$P_{air}={P_O}_2+{P_N}_2$$

(41)

The mass fraction of nitrogen in the exit flow rate is

$${x_N}_2{=1-}_{ }{x_O}_2$$

(42)

Knowing the mass fraction of oxygen and nitrogen in the exit flow rate, the exit mass flow rate of oxygen and nitrogen can then be determined from Eqs. (43) and (44), respectively.

$${\dot{m}}_{{O2\_outca}_{ }}={x_O}_2\ast {\dot{m}}_{{air,out}_{ }}$$

(43)

$${\dot{m}}_{{N2\_outca}_{ }}={{x_N}_2}_{ }\ast {\dot{m}}_{{air,out}_{ }}$$

(44)

If the mass of water in the channels computed from Eqs. (16) and (21) is greater than the maximum amount of vapor the gas can hold (saturation mass), the extra amount is assumed to condense into liquid form instantaneously.

The saturation mass is computed from

$$m_{sat(.)}=\frac{P_{sat}\left(T_{st}\right)\ast V(.)}{R_{H2Ov}\ast T_{st}}$$

(45)

where $V\left(.\right)$ is the lumped volume of the anode or cathode channel and $P_{sat}\left(T_{st}\right)$ is the saturation pressure and can be determined as a function of $T_{st}$ by

$$log{P}_{sat}=-{5.609\ast 10}^{-10}{T}_{st}^4+9.8172\ast {10}^{-7}{T}_{st}^3-6.7687\ast {10}^{-4}{T}_{st}^2+0.22471T_{st}-28.365$$

(46)

where $T_{st}$ is in K and the calculated $P_{sat}$ is in kPa.

The relative humidity inside the channels is then calculated by

$$RH(.)=\frac{P_{H2O_v\left(.\right)}}{P_{sat}\left(T_{st}\right)}$$

(47)

where $P_{H2O_v\left(.\right)}$ Is the partial pressure of water vapor in the channel and is given by

$$P_{H2O_v\left(.\right)}=\frac{m_{H2O_v(.)}\ast R_{H2O_v}\ast T_{st} }{V(.)}$$

(48)

in which $m_{H2O_v}$ is the minimum of the dynamically computed mass and the saturation mass. For the anode, $m_{H2O_v}$ is the minimum of Eqs. (21) and (45). For the cathode, it is the minimum of Eqs. (16) and (45).

The total anode and cathode channel pressures are then determined by Eqs. (49) and (50) respectively

$$P_{an}={P_H}_2+P_{H2O_{v}an}$$

(49)

$$P_{ca}=P_{air}+P_{H2O_{v}ca}$$

(50)

Proton Exchange Membrane

Water is transported across the PEM fuel cell membrane by two main phenomena: electro-osmotic drag (EOD) and back diffusion. As protons move from the anode, they drag water molecules along and this is termed as EOD. This causes a water gradient to form, causing a back diffusion of water from the cathode to the anode.

The net mass flow rate of water across the membrane is given by (Pukrushpan et al., 2002)

$${\dot{m}}_{mem}=M_{H2O_v}\ast A_{fc}\ast N_{st}\ast (\frac{n_d\ast i}{F}-D_w\frac{ (C{v}_{ca}-C{v}_{an})}{tm})$$

(51)

where $n_d$ is the electroosmotic drag coefficient, $D_w$ is the diffusion coefficient, $C{v}_{ca}$is the concentration of water in the cathode and $C{v}_{an}$ is the concentration of water in the anode.

The first term, $\frac{nd\ast i}{F}$ represents the molar rate of EOD and the second term, $D_w\frac{ (C{v}_{ca}-C{v}_{an})}{tm}$ represents the molar rate of water vapor by back diffusion.

The electroosmotic drag is determined by (Springer et al., 1991)

$$n_d \frac{ 2.5\ast {\lambda }_{mem}}{22}$$

(52)

The water concentrations in the cathode and anode are given by Eqs. (53) and (54), respectively,

$$C{v}_{ca}=\frac{ {\rho }_{mem}}{M_{mem}}\ast {\lambda }_{ca}$$

(53)

$$C{v}_{an}=\frac{ {\rho }_{mem}}{M_{mem}}\ast {\lambda }_{an}$$

(54)

where ${\rho }_{mem}$ and $M_{mem}$ are the membrane dry density and dry equivalent weight, respectively and ${\lambda }_{ca}$ and ${\lambda }_{an }$are the water content in the cathode and anode channels, respectively. The water content is defined as the number of water molecules per sulfonic acid sites and is determined by the following empirical relation (Springer et al., 1991)

$${\lambda }_{\left(.\right)}=\left\{\begin{array}{c}\begin{array}{c}0.043+17.81RH\left(.\right)-39.85RH{\left(.\right)}^2+36RH{\left(.\right)}^3, 0<RH\left(.\right)\le 1\\ 14+1.4\left(RH\left(.\right)-1\right), 1<RH\left(.\right)\le 3\end{array}\end{array}\right.$$

(55)

where the subscript (.) represents either the membrane (mem), cathode (ca) or anode (an) and $RH$ is the relative humidity.

The relative humidity of the membrane is modeled as the average of the cathode and anode relative humidities.

$$R{H}_{mem}=\frac{ ({RH}_{ca}+R{H}_{an})}{2}$$

(56)

The diffusion coefficient is expressed by

$${D_w=D}_{\lambda }\exp (2416(\frac{1}{303}-\frac{ 1}{T_{st}}))$$

(57)

in which

$$D_{\lambda }=\left\{ \begin{array}{c}{10}^{-6}, {\lambda }_{mem}<2 \\ {10}^{-6}\left(1+2\left({\lambda }_{mem}-2\right)\right),2 \le {\lambda }_{mem}\le 3\\ {10}^{-6}\left(3-1.67\left({\lambda }_{mem}-3\right)\right), 3 \le {\lambda }_{mem}\le 4.5 \\ 1.25x{10}^{-6},{\lambda }_{mem}\ge 4.5\end{array} \right.$$

(58)

Together these submodels allow the water content in the stack to be monitored.

3.2. GT-Suite Model

The fuel cell stack was modeled in GT Suite v2023 using the stack’s specification. The model was calibrated using experimental voltage and temperature data. Throughout the GT modeling process, a steady inlet hydrogen tank pressure was assumed. In addition, the cathode’s input was a scaled version of the air mass flow rate that was determined in Section 3.1.2 which was optimized to provide accurate voltage predictions. The model, like the experimental PEM fuel cell stack, is devoid of a separate cooling channel. A heat transfer multiplier (HTM) is utilized in its place to account for cooling.

4. Model Discretization

4.1. Physics-Based Model

To predict localized effects, the lumped model was divided into four control volumes (CVs), each with 30 cells, and the output of each CV fed into the next in order to capture variations in membrane water content along the fuel cell stack. The choice of the number of CVs was contingent upon experimental data available for model validation. More details on validation of the model are covered in Section 5.2. The governing equations remain unchanged from the lumped model, with slight modifications made to specific equations to account for the subdivision.

In the discretization, a co-flow technique in which the fuel and air both flow in the same direction was employed. The inlet flow to CV1 is the same as the inlet flow to the lumped model. Subsequent CVs however receive their inlet flows from the outlet of the CV immediately preceding it. The outlet flow is determined by using the linearized nozzle equation in (32) which is adjusted as

$${\dot{m}}_{{ou{t}_{\left(cvi\right)}}_{ }}={k_{ca}}_{out}(P_{cvi}-P_{cv(i+1)})$$

(58)

where $P_{cvi}$ is the total pressure of the current CV cathode channel and $P_{cv(i+1)})$ is the total pressure of the subsequent CV, which is the downstream pressure. For the last CV, $P_{cv(i+1)}$ is equal to the standard atmospheric pressure.

For the anode, it is assumed that the mass flow controller is able to regulate the flow such that pressure differences in the channel are kept to a minimum. As a result, the flow rate to each CV remains unchanged. To guarantee that there is sufficient flow of oxygen for the fuel cell reaction to occur, the anode channel is assumed to have a low flow resistance. Due to this low resistance, the mass flow controller can regulate the flow without significant pressure variation.

Regarding the stack voltage model, the number of cells was modified to reflect the number of cells in a CV rather than the whole stack. The voltage output for each CV was then determined by

$${\mathrm{V}}_{{st,c{v}_i}_{ }}=V_{cell}\ast {N_{cells}}_{c{v}_i}$$

(59)

where $V_{cell}$ was computed by Eq. (2) as described in Section 3.1.1 and ${N_{cells}}_{c{v}_i}$ is the number of cells in the control volume.

The total stack voltage is then computed by adding the voltages in each CV.

$${\mathrm{V}}_{{st}_{ }}=V_{st,cv1}+V_{st,cv2}+V_{st,cv3}+V_{st,cv4}$$

(60)

With regards to discretization of the temperature dynamics model, the efficiency of the fan, ${\eta }_{fan}$ in Eq. (15) was optimized at each current input for each CV using Simulink’s parameter estimation feature in the design optimization toolbox to align the model’s prediction with experimental data.

For the water management model, each CV had a distinct thermal subsystem to account for the temperature variations. Also, the anode and cathode channel volumes were divided according to the number of cells in each CV. Another modification to the discretized water management model is that while the anode is fed with dry hydrogen, this was solely the case for CV1. Back diffusion from the cathode channel results in humid hydrogen entering subsequent CVs. This inlet relative humidity is the output relative humidity from the previous CV which is computed by Eq. (47).

4.2. GT-Suite Model

The developed GT-Suite model was discretized by applying the same method used to discretize the lumped physics-based model. The GT-Suite model was discretized using four distinctive fuel cell stacks, one for each control volume. The input for the first stack was the known measured input, much like the physics-based model, and the input for successive stacks was the output of the stack preceding it. Each stack’s thermal and electrical domains were also appropriately linked. The thermocouple readings were used to calibrate each stack of the discretized model. Similar to the lumped model, the HTM was applied to account for cooling and the HTM used for each stack was iteratively optimized to match the temperature to the experimental thermocouple reading.

5. Results and Discussion

5.1. Lumped Model Validation

The stack voltage and thermal models were validated using data from the experimental setup as described in Section 2. Figure 5 highlights the steady state calibration result for the physics-based stack voltage model where errors are within 5% of the experimental data.

The calibration results for the lumped thermal model are depicted in Figure 6. The fan efficiency was used as a tuning parameter in the thermal model calibration. In Simulink, an optimization search for best efficiency was carried out, and an optimal efficiency of 35.75% was reached. When compared to the experimental data, the maximum steady state error for the thermal model is only 0.4%, highlighting the model’s accuracy in predicting the stack temperature.

Due to the lack of experimental data to validate the water management model’s prediction, simulation results from GT-Suite model provided a basis for comparison. The GT-Suite model was calibrated using the experimental voltage and temperature data and the results are shown in Figure 7 and Figure 8 respectively, with a maximum relative error of 2.38% for the voltage model and 0.81% for the thermal model.

After validating the thermal and voltage predictions of GT-Suite, its membrane water content was used to validate the physics-based water management model. The comparison results are shown in Figure 9. The root mean square error (RMSE) between the predictions was found to be 0.72 membrane water content. Notably, the GT-Suite model predicted higher water content for all current cases than the physics-based model. One plausible explanation for this could be due to the temperature predictions of the GT-Suite model. For current inputs ranging from 0A to 40A, the GT-Suite model overestimates the stack temperature which could lead to an increase in water content. On the contrary, for these same current values, the physics-based model slightly underestimates the stack temperature. The discrepancy in the predictions could be due to variations in the models’ initialization. Even though the GT-Suite model underestimates the water content from 60A to 85A, the water content remains high. This could be due to water saturation as a result of increasing water generation in these high current zones. As a result, despite the modest shift in trend of stack temperature, the water content may remain elevated.

The membrane water content dynamics from the physics-based model is as shown in Figure 10. As can be observed, an increase in current input leads to a decrease in the membrane water content. The dynamics of the membrane water content is influenced by the relative humidity dynamics in the anode and cathode channels as depicted in Figure 11. The cathode relative humidity remains constant at 100% at all current cases due to water generation at the cathode electrode. The relative humidity of the anode channel however decreases with increase in current. This can be due to an increase of electroosmotic drag from the anode to the cathode channel. The parameters used for simulation of the models are shown in Table 1.

5.2. Discretized Model

While the lumped model captures the overall stack behavior, variations are present internally. The steady state experimental results of the temperature recorded by the four thermocouples is illustrated in Figure 12. As depicted in the figure, for low current regions (0A-20A), the temperature recorded by the thermocouples is comparable with a variation within 2K. As the current increases, the temperature variation becomes more pronounced due to higher electrochemical reaction kinetics. Particularly at 85A current, there is a maximum temperature variation of 10K recorded. Notably, the left thermocouple (which is farthest from the direction of airflow) recorded the highest temperature for the medium to high current regions. This may be due to reduced cooling efficiency as a result of the location of the thermocouple. Interestingly, the mid-left thermocouple (just before the left thermocouple) recorded the lowest temperature at the higher current region. This suggests that the airflow may be directed in such a manner that it promotes effective cooling in this region.

The discretized model was validated using the experimental thermocouple measurements. For the discretized model, the temperature in CV1 corresponds to the right thermocouple reading, in CV2 it is the mid-right, in CV3 it is the mid-left and in CV4 it is the left thermocouple. Figure 13, Figure 14, Figure 15 and Figure 16 display the thermocouple calibration results for each of the four thermocouples. The maximum relative error for each case is shown in Table 2.

The membrane water content dynamics for each CV is as shown in Figure 17. It is evident that the membrane water content varies along each CV. These variations are most noticeable in the medium to high current regions (40A to 85A) where there is increased water activity. Notably, the trends in the temperature variation results and the variation in membrane water content were the same. The CV with the highest temperature also had the highest water content .

To gain a better understanding of the membrane water activity, the relative humidity dynamics occurring along the cathode and anode channels must be highlighted. The relative humidity dynamics in the anode and cathode channels are shown in Figure 18 and Figure 19 respectively. As can be seen, the cathode channel remained constant at 100% relative humidity for all CVs for all current cases. Conversely, the anode relative humidity varied along each CV. Hence, it is evident then that the observed variation in membrane water content stemmed solely from variations in relative humidity in the anode channel. This explains why some researchers opt to model the membrane water content as a function of simply the anode water content since it is a limiting factor in determining the membrane water content.

To validate the discretized water management model, the model’s steady state predictions were compared to simulation results from the discretized GT-Suite model. Firstly, the GT-Suite model was calibrated and validated using the thermocouple readings and the calibration results are illustrated in Figure 20, Figure 21, Figure 22 and Figure 23. The maximum relative error for each case is shown in Table 2. Following validation of the GT-Suite temperature variations prediction, its membrane water content predictions were used to validate the membrane water content variation of the physics-based model.

The steady state variation of membrane water content in the GT-Suite model is shown in Figure 24. As can be seen, the trend for the physics-based and GT-Suite models are comparable. The variations in membrane water content are much more evident in the medium to high current regions. The highest water content was recorded in the CV with the highest temperature. In the GT-Suite model however, CV3 has a higher membrane water content than CV1 despite having the lowest temperature. Nonetheless the difference in membrane water content between CV1 and CV3 for both the physics-based and GT-Suite model is only about 3%.

Comparison of the steady state result of the physics-based model to the GT-Suite simulation results for each CV is shown in Figure 25, Figure 26, Figure 27 and Figure 28. For each CV, the trend in membrane water content prediction is the same. Table 3 shows the RMSE values between the models’ predictions. The physics-based model’s ability to predict the distribution of membrane water content in the stack is demonstrated by the RMSE, which was within 1.3 membrane water content for all cases.

As shown in Figure 29, temperature variation and its ensuing variation in membrane water content resulted in variation in voltage output along the stack. As can be observed, at low current density, the voltage output for each CV is identical but as the current density increases, the voltage output varies more noticeably. Particularly, CV4 recorded the lowest voltage output, while CV3 and CV1 recorded the highest voltage, with CV3 slightly exceeding CV1 in voltage. As determined from the thermal and water management physics-based models, CV3 had the lowest stack temperature and water content followed by CV1 while CV4 recorded the highest stack temperature and the highest water content. One possible explanation for CV3 and CV1’s enhanced performance could be due to the improved cooling system in these regions which resulted in the lower stack temperatures. Better cooling implies higher air flow rate, which increases oxygen availability and in turn, increases performance. It is, however, well documented that lower membrane water content leads to higher ohmic overvoltage in a PEM fuel cell stack. Hence, this improvement in performance only occurs because the increased ohmic resistance is smaller than the increased performance caused by the availability of oxygen.

6. Conclusion and Future Work

A reduced-order control-oriented physics-based model for predicting variation in membrane water content in an open-cathode PEM fuel cell was developed in this study. The water management model was validated by comparing its prediction to GT-Suite simulation results, and RMSE were found to be within 1.3 membrane water content for all cases. The results indicate that the water content of the PEM fuel cell membrane varies along the fuel cell stack. Variations abound, particularly in medium to high current density regions where there is significant water activity. The variations of the membrane water content coupled with the temperature variations leads to variation of the voltage output along the stack. Control of these variations is therefore critical especially during transients. The developed model thereby presents a computationally inexpensive means for prediction and subsequent control of membrane water content variations. Compared to the lumped physics-based model which took ~12 seconds to simulate about 40 minutes of experimental data, the computational time for the 4-CV discretized model was ~80 seconds. Despite the increased computational time, it is still relatively low for control. By contrast, the lumped GT-Suite model took ~240 seconds to compute, which is almost 20 times slower than the lumped physics-based model. The computational time for the discretized GT-Suite model was ~600 seconds, which is about 7.5 times slower than the discretized physics-based model.

Future research will include validating the water management model experimentally to improve the model’s accuracy and reliability by aligning its predictions with real-world observations. In addition, the model will be improved to include liquid water effects to predict flooding conditions particularly in the cathode channel.