Semi-Empirical Models for Stack and Balance of Plant in Closed-Cathode Fuel Cell Systems for Aviation

Abstract

In recent years, there has been a growing interest in utilizing hydrogen as an energy carrier across various transportation sectors, including aerospace applications. This interest stems from its unique capability to yield energy without generating direct carbon dioxide emissions. The conversion process is particularly efficient when performed in a fuel cell system. In aerospace applications, two crucial factors come into play: power-to-weight ratio and the simplicity of the powerplant. In fact, the transient behavior and control of the fuel cell are complicated by the continuously changing values of load and altitude during the flight. To meet these criteria, air-cooled open-cathode Proton Exchange Membrane (PEM) fuel cells should be the preferred choice. However, they have limitations regarding the amount of thermal power they can dissipate. Moreover, the performances of fuel cell systems are significantly worsened at high altitude operating conditions because of the lower air density. Consequently, they find suitability primarily in applications such as Unmanned Aerial Vehicles (UAVs) and Urban Air Mobility (UAM). In the case of ultralight and light aviation, liquid-cooled solutions with a separate circuit for compressed air supply are adopted. The goal of this investigation is to identify the correct simulation approach to predict the behavior of such systems under dynamic conditions, typical of their application in aerial vehicles. To this aim, a detailed review of the scientific literature has been performed, with specific reference to semi-empirical and control-oriented models of the whole fuel cell systems including not only the stack but also the complete balance of plant.

1. Introduction

The hydrogen economy represents a highly promising opportunity for acquiring cleaner energy and propulsion technologies in today’s context. Within the realm of transportation, there is a strong fascination with fully electric propulsion systems across various domains, including aircraft, unmanned aerial vehicles (UAVs), ground vehicles, trains, ferries, and more. Fuel cell systems present an attractive solution for harnessing hydrogen energy due to their significantly higher efficiency compared to internal combustion engines and the superior energy density of hydrogen in comparison to batteries. In the field of aerospace engineering, this becomes particularly critical in enhancing flight endurance and range. Initially, fuel cells were introduced in aerospace as components of auxiliary power units. However, there is a growing interest within the European Union in adopting hybrid electric power systems that combine fuel cells, batteries, and possibly supercapacitors to power electric propulsion systems in small aircraft [1]. As pointed out in [2], the climate impact reduction potential of a hydrogen fuel cell system is in the range of 75–90%. This range was obtained by considering direct emissions of CO₂ and NO_x, water vapor, contrails, and cirrus. Another advantageous aspect of utilizing fuel cells in aerospace applications is that their by-products, namely water, heat, and oxygen-depleted air, can be harnessed to conserve resources and reduce overall take-off weight [3]. The disadvantages of these systems are mainly the high cost and the lack of a hydrogen refilling infrastructure [4].

There are three main types of fuel cells: Proton exchange membrane fuel cells (PEMFC), direct methanol fuel cells (DMFC), and solid oxide fuel cells (SOFC) [5]. PEMFCs have the highest power density of all fuel cell types, have good start-stop capabilities, and operate at low temperatures. On the other hand, they are more expensive because of precious catalysts (platinum), polymer membranes, and high-cost ancillaries. Moreover, they require very high purity for hydrogen, and high-temperature technologies must be adopted to improve tolerance to CO. Since PEMFCs are the most used fuel cell in the aerospace field, this review concentrates on this technology.

The structure of a PEM-FC is shown in Figure 1. The core of the cell is the membrane electrolyte assembly (MEA), which includes the membrane, the catalyst layers, and the gas diffusion layers. To obtain the cell, the MEA is inserted between bipolar plates that are used as electrods and to provide flow channes for the reactants. Depending on the voltage requirement, a certain number of cells are connected in series to form the stack. Current collectors are used to connect the cells and gaskets are added to obtain a gas-tight seal. Cooling plates can be introduced in the stack in case of liquid-cooled systems. The cell converts the chemical energy of hydrogen into electricity through the following overall reaction:

$$H_2+1/2O_2\to H_2{O}_{}+286 \mathrm{k}\mathrm{J}{\mathrm{m}\mathrm{o}\mathrm{l}}^{-1}$$

(1)

Since the average voltage of a single cell is quite low (0.6–0.7 V), to reach the desired voltage, a certain number of cells are connected in series to form the stack [6].

To obtain a complete fuel-cell system, the stack must be integrated with various other components, namely, the hydrogen supply system, the air management system, the cooling system, and the electric power control system. These subsystems are commonly referred to as the Balance of Plant (BOP) [6]. The overall efficiency of a fuel cell system, and consequently, the range and endurance of a hydrogen-powered aerial vehicle, are influenced not only by the specifications of the fuel-cell stack but also by the management and power consumption of the auxiliary systems necessary for its operation [4].

Three operation strategies are usually adopted for fuel supply: flow-through, dead-end anode, and recirculation [7]. In the flow-through mode, usually used for testing fuel cells, excess hydrogen is supplied to the anode to remove water diffused from the cathode, resulting in lower hydrogen utilization with a waste of fuel and H₂ emissions. In the recirculation mode, non-reacted residual hydrogen is recirculated back to the supply line by a pump or an ejector. The excess fuel and impurities pass through a separator to isolate the fuel from the impurities. Then, the hydrogen is re-humidified and re-fueled to the stack. This system is quite complex and causes additional parasitic loads (power absorbed by the pump) [8,9,10]. In the dead-end mode, the anode is sealed off to achieve 100% fuel utilization. A pressure regulator is installed at the inlet of the anode to reduce the pressure of hydrogen from a gas tank, while a solenoid valve is placed at the outlet to seal the anode. In the dead-ended anode operating condition, the cell voltage gradually decreases with time because of the accumulation of impurities (inert gases, liquid water) in the anode, until the purge valve is used. The frequency of activation of the purge valve according to the current density is another critical issue in this system [11].

Different configurations can be used for the cathode as well. Commonly, a fuel cell is called an open-cathode PEM fuel cell (OCPEMFC) when the same airstream is used as a reactant and as a cooling medium. The air can be freely breathed or forced into cathode flow channels through one or more fans [12]. Free-breathing open cathode fuel cells have the simplest configuration and no parasite load, but the stack temperature cannot be controlled, and an external heater is required for a cold start. Moreover, there is a high risk of membrane dehydration and flooding. Forced-air open cathode systems are used in applications where the required power is below than a certain threshold (5–6 kW). The performances of the OCPEMFC are strongly affected by environmental conditions. For higher levels of power, the heat released from the electrochemical reaction may not be sufficiently removed without a separate cooling circuit. Moreover, the reactant supplied to the cell may not be sufficient to sustain the optimum electrochemical reaction. These problems are solved by adopting a liquid-cooled closed-cathode configuration (CCPEMFC) [13,14]. Closed cathode configurations require a more complex balance of plants that determines higher parasitic power, but their behavior is less affected by surrounding conditions, and heat dissipation can be better controlled. Moreover, they are less prone to membrane dehydration and flooding. Research on CCPEMFCs has mainly been carried out on ground applications. Most automotive fuel cell stacks (50–90 kW) are liquid-cooled using deionized water or a water-glycol mixture. The cooling fluid is continuously recirculated and only periodically replenished. The most relevant recent research topic for CCPEMFC is the control of the oxygen excess ratio (OER). If the OER is less than one, oxygen starvation occurs, which in turn produces cell voltage reduction, stack floodings, and membrane damage. On the other hand, if OER is too high, excessive power is consumed by the air compressor so that the net power of the fuel cell decreases [15]. A third configuration that solves the problem of air starvation at very high altitudes is represented by the so-called Closed Cathode and Anode fuel cell (CCAPEMFC) [16] or hydrogen-oxygen fuel cell [17] where both the cathode and anode are sealed, and the oxidizer used in the reaction is not ambient air, but pure oxygen stored in a vessel. However, this solution requires a further storing vessel for oxygen that adds to the weight of the system and is rarely adopted in aerospace applications.

In the aerospace field, open cathode fuel cells are generally considered unsuitable for power higher than 6 kW. In fact, liquid-cooled PEMS are used in research investigations and projects like ENFICA-FC [18], and the Sigma-4 ultralight aircraft [19]. However, there is some research on the possibility of extending the application range of air-cooled fuel cells [20] since using air as a coolant requires less components and increases the overall system size by less than 50%, while water-cooled systems can increase the system size up to 200% [21]. If the fuel cell stack is located in the front part of the fuselage, the intake air can be used to feed the cathodes and cool down the entire system [22,23]. An example of an air-cooled closed-cathode fuel cell is the Nexa Ballard [24]. Passive cooling with a heat spreader and/or heat pipes is also explored on the medium scale (<10 kW) [14]. Edge-cooling and phase change cooling have been proposed for the control of temperature [25,26,27]. Edge-cooling makes use of fins on the stack and can improve the maximum stack power by more than 15% [13]. This increases the volume of the stack and requires a high flow of air but can be effective for UAVs where the ambient flow can be utilized for heat removal [28].

The performances of the fuel cells are affected by the drastic changes in elevation and load during the flight [29]. The increase in flight altitude determines a decrease in air pressure, air temperature, air density, etc., that changes the flow and pressure working point of the air compressor. When the load increases suddenly, the chemical reactions inside the fuel cell accelerate, increasing hydrogen consumption. However, because of the time lag of the compressor and its motor, protons migrating from the anode to the cathode cannot participate in the oxygen reduction reaction due to a lack of oxygen. Moreover, a reduction reaction of hydrolyzed hydrogen occurs in the air-deficient region, causing a rapid decrease in voltage. The hydrogen precipitated from the cathode can also react directly with oxygen, releasing more heat that can irreversibly damage the membrane [30]. The change in aircraft flight power and elevation also causes fluctuations in the cathode pressure. An increase in cathode pressure is positive since it increases the gas diffusion rate and improves the air distribution on the surface of the membrane electrode, resulting in higher stack performance and faster dynamic response. On the other hand, if the pressure is too high, the proton exchange membrane can be damaged, and the service life of the PEMFC is shortened. As pointed out in [29], the change of cathode pressure will increase or decrease the output voltage of PEMFC, and then result in the fluctuation of DC bus voltage. Moreover, the change in bus voltage will affect the speed regulation of the air compressor, and then affect the supply of air mass flow. Finally, ambient air can damage the fuel cell membranes due to the presence of dust and very dry and cold ambient air [28]. For this reason, a protective dome with controlled ventilation is adopted with the additional advantage of the potential recovery of water vapor produced by the fuel cell.

From the above considerations, it can be easily understood that modeling a fuel cell for aerospace applications is a very challenging task, not only because of the interaction between different physical phenomena but also because of the nonlinear behavior of the whole system [31].

The models for Proton Exchange Membrane Fuel Cells (PEMFCs) can be categorized into three types: black box, grey box, and white box models, as described in reference [32]. Black box models associate experimental input and output data without considering the underlying physical processes within the fuel cell system. Examples of this model type include artificial neural networks and fuzzy logic. They are known for their computational efficiency and are suitable for real-time identification [32]. In [33], a data-driven model is employed for temperature control, utilizing active disturbance rejection control, which requires the tuning of only two parameters. White box models, on the other hand, are built upon algebraic and differential equations based on the principles of thermodynamics, electrochemistry, and fluid mechanics. These models are valuable for studying the impact of geometric and operational parameters on fuel cell performance. Grey box models represent a middle ground, incorporating physical relationships, such as the polarization curve, supported by experimental data. They strike a balance between complexity and simplicity, and are particularly well-suited for energy management and control applications [32]. Nòbrega [13] recently reviewed physics-based models for PEM fuel cells, focusing on stationary models applied to liquid-cooled fuel cell systems and their behavior under specific operating conditions.

Quasi-static and dynamic can be used for fuel cell systems. Most of the models proposed in the scientific literature aim to predict the current-voltage and current-power curves of the PEMFC [8] under very slow variation of the current (quasi-static operation). A group of resistors and capacitors is used to describe the behavior of electrochemical reaction kinetics, ohmic conduction processes, and mass transport. These models are often tuned at a single operating condition, and a variety of tuning procedures have been proposed. Dynamic models are generally implemented in simulation studies, using equivalent electric circuit models implemented in simulation environments like MATLAB Simulink [34] and PSPICE [35]. The Aspen simulation platform has been also proposed for the implementation of dynamic models of CCPEMFC BOPs in [36] and Rabbani [37]. As pointed out in [38], dynamic models are very useful in emulators for testing the system control, evaluating responses to fast load changes, and estimating maximum current capabilities and other critical operational conditions. However, the correct evaluation of hydrogen consumption and parasitic power is important for evaluating the actual performance of the fuel cell in terms of power and overall efficiency. From our analysis of the scientific literature, only a few papers address this issue.

Finally, fuel cell models can be classified according to the level of analysis: single cell, stack, or system. Single-cell and stack models primarily focus on the voltage-current characteristics of the fuel cell, while system-level models encompass the overall plant balance.

The scientific literature on the modeling and control of the BOP of a CCPEMFC is quite rich, but most works come from the automotive field [4,36]. An overview of hydrogen fuel cell architectures, including the DC/DC converter, and intelligent controls has recently been proposed by Guo et al. [39], while a review of physical-based and data-driven models for real-time control of PEMS is proposed in [33]. Understanding and modeling the dynamic behavior of a PEMFC fuel cell system are particularly critical in applications like small drones or Unmanned Aerial Vehicles (UAVs), where nonlinear factors, such as the variable space environments and different load conditions, must be taken into account [29].

The present study explores quasi-static and dynamic system-level grey-box models. This choice is motivated by the fact that data-driven models, although highly accurate for predicting specific fuel cell system performance within a certain range of operational parameters, lack generality and cannot be applied beyond the range of variation used for the inputs in the training phase.

The goal of this investigation is to review the semi-empirical approaches used in the literature, in particular, the possibility of accounting for the effect of ambient pressure, temperature, and relative humidity on each subsystem, in order to define a complete simulation protocol for the whole CCPEMFC fuel cell system for aerial applications.

The innovative contribution of this review, compared to review papers recently presented in the scientific literature, can be found in these aspects:

• The review includes models for the whole fuel cell system and not only the stack.
• Both quasi-static and dynamic models are analyzed in detail.
• The interaction between the complex phenomena is illustrated, including describing current trends in control strategies.
• The effect of variable altitude conditions is specifically addressed.

The paper is organized as follows. Section 2 describes the typical configuration of a CCPEMFC power plant for road and air transportation. The main performance indexes and the experimental facilities adopted in the scientific literature are also discussed in this section. Section 3 addresses all the transient phenomena taking place in a CCPEMFC system and mentions some control techniques. Section 4 and Section 5 describe, respectively, the state-of-the-art quasi-stationary and dynamic sub-models. Section 6 addresses the effect of variable ambient pressure and humidity caused by a variable flight altitude and how this can be modeled.

2. Balance of Plant and Performance Indexes

A fuel cell system essentially converts the chemical content of hydrogen directly into electric power. At the stack level, the output power can be expressed as:

$$P_{el}=V_{st}{I}_{st}$$

(2)

where $V_{st}$ and $I_{st}$ are the voltage and the current of the stack, which consists of $N$ elementary cells connected in series. As the current increases, the voltage drops due to a series of phenomena described later in the paper. Actually, the main goal of a fuel cell model is to predict the dependence of $V_{st}$ on $I_{st}$, the so-called voltage-current curve or polarization curve. Since $V_{st}$ depends on $I_{st}$, the electric power shows a non-linear monotonical dependence on $I_{st}$. Abu-Hawa et al. [40] performed a numerical sensitivity analysis to evaluate the effect of the operational variables (temperature, pressure, and water content) on the maximum power output of a fuel cell. If the stack temperature varies from 80 °C to 64 °C, the output power is reduced by about 9%. If pressure is decreased from 1 bar to 0.8 bar, the gross power reduction is small (less than 1%).

At the system level, it is necessary to account for the parasitic power of the auxiliaries needed for the management of the mass flow, humidity, and thermal control. Figure 2 reports the typical BOP of a fuel cell system for automotive applications [41].

Note that the power plant of Figure 2 includes several operating fluid machines: the air compressor, the hydrogen circulation pump, the humidifier water circulation pump, the coolant pump, and the cooling fan of the heat exchanger. Simpler configurations are often used in the aerospace field where self-humidified stacks are adopted, and the cooling fan can be omitted [42]. In some cases, a two-level cooling system is adopted consisting of high-temperature and low-temperature circuits [43].

The configuration of Figure 3 refers to an aerospace application and shows the air circuit in red, the hydrogen circuit in magenta, and the electric flows in green. The cooling system is represented in blue. The control volume identified in the figure will be used in this work to model the stack voltage and temperature. The main nomenclature used in the description of the subsystems is also reported in the figure.

More complex configurations are sometimes adopted in aerospace applications. In [44], a heater is used to preheat the hydrogen stored at cryogenic temperatures, while a cooler is used to reduce the temperature at the compressor outlet.

All operating fluid machines included in the BOP are moved by an electric motor that absorb a part of the gross electric power generated by the stack, known as “parasitic power”. Brushless DC motors (BLDCMs) are generally adopted to move these machines [29]. An additional bias power is considered by Guzzella et al. [41] to account for the minimum flow of reactants necessary to keep the fuel cell from shutting down.

To account for these effects, the gross electrical power expressed by Equation (2) must be corrected to obtain the net power of a fuel cell system, which can be written as:

$$P_{net}=V_{st}{I}_{st}{\eta }_{DC-DC}-P_{aux}$$

(3)

where ${\eta }_{DC-DC}$ is the efficiency of the DC-DC converter and $P_{aux}$ is the parasitic power of the whole BOP.

In a CCPEMFC like that shown in Figure 3, the only contribution to the parasitic losses are the compressor and the pump of the thermal management module. In [37], the compressor consumption was estimated to be 18.5% of the gross power. In [4], for a fuel cell stack with a gross power of 100 kW, the largest compressor consumption is 12 kW, while the consumption of the thermal management module is estimated to be 1.9 kW.

Like power, efficiency can be defined at the stack level or the system level. At the stack level:

$${\eta }_{stack}={\mu }_F\frac{V_{st}}{N·E_0}$$

(4)

where ${\mu }_F$ is the fuel utilization rate [40], defined as the ratio between the flow rate of hydrogen that reacts in the stack, divided by the inlet fuel flow rate. This value is assumed to be equal to 0.95 in [45] and 0.9 in [46], depending on the anode configuration. In a dead-mode anode with a purging valve, the fuel utilization depends on the managing of the purging process. $E_0$ is the electric equivalent of the hydrogen heating value that corresponds to 1.462 V or 1.254 V if calculated from the higher heating value (HHV) [47] or the lower heating value [40], respectively. The voltage of the fuel cell, $V_{cell}$, is significantly lower than $E_0$, mainly because of the three typical loss phenomena of electrochemical devices, namely activation, ohmic, and concentration losses [6].

At the system level, the net power of a fuel cell system can be expressed as [1]:

$${\eta }_{net}=\frac{P_{net}}{{\dot{m}}_{H_2}HHV}$$

(5)

where ${\dot{m}}_{H_2}$ is the hydrogen mass flow rate consumed in the process. This expression of efficiency is self-evident since the fuel cell system is a converter from chemical energy to electric energy. The net efficiency can be quite lower than the fuel stack in a CCPEMFC at low load, because of the higher effect of parasitic losses (see Figure 4).

To reduce the contribution of the parasite power, the compressor is coupled with a turbine to form a turbocharger like that used in piston engines [48]. However, differently from piston engines, the low enthalpy of exhaust gases makes the turbine incapable of producing all the power needed by the compressor. To solve this problem, a combustor can be included in the cathode sub-system. In [44], the combustion of hydrogen is proposed to increase the enthalpy of the exhaust gas before entering the turbine. In this case, the net efficiency is evaluated as:

$${\eta }_{net}=\frac{P_{net}}{{\dot{(m}}_{H_2}+{\dot{m}}_{fuel})HHV}$$

(6)

where ${\dot{m}}_{fuel}$ is the flow rate of hydrogen burned in the combustor.

For a given application, FC system designers can select the operating range of the fuel cell voltage according to [47]:

• Maximum power point (lower stack voltage)
• Maximum efficiency (higher stack voltage)
• Optimal fuel utilization

Other parameters can be introduced to account for the performance of the cooling system. In particular, the stack cooling effectiveness is defined as the ratio between the heat removal rate and the generated thermal power by the stack [49]. The effectiveness of a particular cooling device can be also evaluated by considering the amount of heat removal it accomplishes compared to the electrical power it consumes. An effectiveness in the range of 20–40 is generally attainable for well-designed cooling systems.

Experimental Facilities

The tuning of the grey models needs a comparison to experimental data that is usually obtained from electrochemical in situ characterization techniques. In these techniques, variables like voltage, current, temperature, etc., are used to characterize the performance of the fuel cell under typical operating conditions. The most useful methods for the characterization of fuel cell behavior and the fitting of the quasi-static and dynamic models are [1]:

• Current-voltage measurements.
• Current Interrupt Tests (able to separate the losses into ohmic and non-ohmic).
• Electrochemical Impedance Spectroscopy (it can distingue between ohmic, activation, and concentration voltage drops).
• Cyclic Voltammetry (provides insight into fuel cell reaction kinetics).

Current-voltage measurements are usually considered sufficient for the validation of stack voltage models. However, they cannot allow the separate analysis of the activation, ohmic, and concentration losses.

The test bench usually consists of sensors, data acquisition systems for the relevant measurements, controllers for the electric motors and valves, and a programmable electronic load to simulate a variable current load. In transport applications, the load can be applied by connecting the fuel cell to the traction motor [50]. In this case, the power input is directly controlled by the throttle. Hydrogen supply systems and water and heat management systems complete the setup [51]. The values of ambient pressure, temperature, and humidity, and the relaxation times at which the voltage-current tests are performed are not reported in most of the papers analyzed in this review.

Altitude chambers have been proposed to study the effects of reduced pressure and elevation when the fuel cell is operated at high (<11,000 m) and very high altitudes (>11,000 m). Some of these chambers control pressure, while temperature and humidity are only monitored [52]. In [53], a climate chamber with control of temperature, pressure and humidity is proposed to simulate atmospheric conditions up to 5500 m altitude. In [28], a wind tunnel is used to verify the effect of using edge-cooling in a fixed-wing drone. A similar approach is used in [54], where the wind tunnel is used to determine the convective heat transfer coefficient of a high-temperature PEM stack. In [31], the experimental setup includes a high-accuracy air flow sensor to measure the inlet air velocity at the cathode. The stack, moreover, is installed in an environmental chamber where it is possible to control humidity, temperature, and oxygen concentration. An environmental test chamber was created in a wind tunnel by the authors of [55] to investigate the effect of elevation on a single PEM cell. The chamber was equipped with a test system, including a frequency response analyzer to collect the impedance spectrum. Electrochemical Impendence Spectroscopy (EIS) has been applied by [8] to consider the diffusion and charge-transfer phenomena.

3. Transient Behavior of a Fuel Cell System and Control Actions

Analyzing the dynamic behavior of a fuel cell is of paramount importance in applications like automotive and UAVs [56] that are characterized by rapid load variations. One limitation of fuel cell systems, compared to battery-powered configurations, is the slower response to dynamically variable load profiles. According to [46], PEMFCs have a no-load-to-rated-power transient of 20–30 s. As pointed out in [57], a supercapacitor can be charged in 1–10 s, while a few minutes are needed to charge the battery to 70%. In [58], the dynamic response to transient loads was evaluated for a 1.2 kW Ballard fuel cell and compared to that of a battery. The step-down (20 A–0 A) and step-up (0 A–20 A) transient response times were found to be 1.2 s and 0.6 s, respectively. In [59], the change rate of the load demand for a 62 kW fuel cell is set between −10 kW/s and 7 kW/s. A response rate of up to 8 kW/s was found in a 42 kW system for automotive applications [60]. For a 200 kW PEMFC designed for marine propulsion, a typical value of 20 kW/s is assumed in [61]. This explains the need for a battery or a supercapacitor to provide sufficient responsiveness to load changes.

The dynamic response of a fuel cell system to load variation is slow because of the necessity to monitor and control the flow rate and pressure of the reactants (OER and HER), the temperature of the stack, and the humidity of the membrane. In particular, the slow dynamic of the air compressor affects the response of the whole system to fast-changing load demands [62]. This makes the optimization of the performance of the whole system very complex [40,63], particularly at high altitudes, also due to the trade-off between fast response and minimum reactant consumption [64].

The control of temperature is a critical issue because of the strong effect that this operating parameter has on the performance of fuel cells. When the temperature increases, the electrochemical reaction becomes faster, improving efficiency. However, if the temperature increases too much, the internal resistance increases because of damage to the membrane [65], and both performance and efficiency are degraded. Excessively high temperatures produce drying of the catalyst layer and a decrease in voltage [31]. The temperature of the stack, according to [40], should keep at 80 °C, while in [65], a range of 65–85 °C is considered. This issue, however, is more critical in open-cathode fuel cells, where the control of the temperature is coupled with the management of the air flow rate through the fan [66]. In liquid-cooled closed-cathode systems, on the other hand, temperature can be easily controlled by regulating the flow rate of cold water.

The transient behavior of a PEMFC is also relevant for start-up (to achieve a rapid increase in stack temperature up to the optimum value) and shut-down operations (to avoid degradation of the fuel cell, as reported in [63]). At a cold start, the efficiency of the fuel cell is quite low, and heaters can be used to rapidly heat the stack. Moreover, fluctuations in temperature must be avoided to ensure a long cell life [67]. As already explained, the cooling circuit is bypassed to achieve a fast reaching of the desired stack temperature.

During transients, it is paramount to avoid the starvation of fuel and oxygen, which reduces cell voltage [65,68]. A common method is to use more fuel and air than needed for the specified power. However, this can lead to a waste of hydrogen, which is an expensive and valuable energy carrier. An excess of air (OER), on the other hand, does not cause a waste of precious fluid. However, since air is fed through a compressor in a CCPEMFC, it causes an increase in the power of these auxiliaries and reduces the net power of the fuel cell. In [65], stack orientation is also proposed as a potential variable to be controlled to solve this problem. The excess water accumulated at the anode is removed through a purge valve that also causes a loss of hydrogen and a further subsystem to be controlled.

The pressure of the hydrogen must also be controlled so that the difference of pressure across the membrane does not exceed a certain value; otherwise, reactants may cross over, and membrane damage can potentially occur [40]. This is particularly critical in high-altitude applications because of the low pressure at the cathode.

Humidity must also be kept within optimal bounds because a dried membrane presents higher resistance, but excessive humidification can cause flooding in the stack. In both cases, the performance of the fuel cell is reduced. Control of water content can be achieved by humidifying the gas stream of the reactants with an external circuit. However, most fuel cells used in the transportation sector are “self-humidifying” [31] to avoid the burden of an external humidifier. In these fuel cells, internal humidification is accomplished by retaining the water produced by the reaction, thanks to the use of additives that can be integrated either into the catalyst layer or the membrane [66]. Water content is not uniform in the stack, being higher at the cathode where protons drag some water molecules. This effect is counterbalanced by back-diffusion. During transients, the anode dries out until back-diffusion starts to take effect. The estimated time constant for this phenomenon depends on the rapidity of current changes and the final value of the current. Six seconds are needed for a step change in current density from 0.1 to 0.8 A/cm² [40]. As shown in [56], the dynamic polarization curve shows hysteresis: a sudden increase in the load reduces the water content in the membrane and generates a lower voltage at a given current than in the case of reduced load. The faster the load variation, the higher the magnitude of the hysteresis.

At the cell level, a further dynamic effect is the so-called double-layer capacitor effect. During the normal operation of the cell, the anode-electrolyte layer and the electrolyte-cathode layer can be charged by the polarization effect and behave as a supercapacitor. The double-layer discharge has a very small time constant, whereas the thermal management system has a time constant of a few minutes [66,69]. The time constants for electrochemical double-layer and water membrane hydration were obtained from [70] with a numerical approach, but the order of magnitude was also confirmed experimentally [56].

Control Methods

Up to five interconnected controls can be included in a complete CCPEMFC system (Figure 2) since the whole balance of plan needs to be coordinated to meet variable load requests [41,71]:

• The control of the hydrogen regulation valve;
• The control of the compressor air flow through the motor voltage;
• The control of the cathode pressure;
• The control of the temperature through the cooling pump speed;
• The control of humidity.

Moreover, power management control is needed to obtain the load power since the voltage of the stack is a function of the current and the parasitic power of the BOP. Coordination between the DC/DC converter and the air compressor driver circuit is suggested to avoid starvation [39].

As already explained, the main control action to protect the cell and optimize its performance is the Oxygen Excess Ratio (OER) [36]. The typical curves of the net power as a function of the OER for each value of stack current are shown in Figure 5. At constant current (green lines), the net power initially increases with OER because of the improved reaction rate until reaching a maximum (blue circle). For higher values of the OER, the net power declines due to the increased power absorbed by the compressor. The net power increases with current, but the range of feasible OER becomes narrower because of the compressor limits. This leads to identifying an optimal value of OER for each value of the current that is used as a setpoint for the control of the CCPEMFC [72,73]. The optimization of OER was found to increase the net power by about 3.5% in [74].

By interpolating the optimal values of OER for each value of stack current, an empirical equation is derived and implemented in [75] to obtain the optimal value of oxygen rate at any time:

$$OE{R}_{opt}\left(t\right)=c_1{I}_{st}^3+c_2{I}_{st}^2+c_3{I}_{st}+c_4{I}_{st}^{-1}+c_5$$

(7)

where $c_1, \dots c_5$ are obtained by fitting the typical curves of performance shown in Figure 5.

The setpoints guarantee the optimal net power in steady-state operation, but the power consumption during the transient OER also needs to be considered, as discussed in [76]. During a rapid increase in load current, the acceleration of the chemical reaction can cause the OER to deviate from the optimal value and cause oxygen starvation [76] or excessive parasitic power. Moreover, there is a coupling between the control of the OER and the control of the hydrogen flow. In fact, the difference between the anode and the cathode pressure must be controlled to avoid membrane damage. This is performed in [36] by tracking the anode pressure with hydrogen flow manipulation using the cathode pressure signal as a reference value. The HER is also considered in [77,78] as a control parameter to keep the pressure difference between the cathode and the anode within a reasonable range and to minimize the fuel consumption of the fuel cell. The authors highlight that the influence of air and hydrogen excess ratio on power is more relevant at high current levels. The contemporary control of OER and anode pressure is particularly critical in UAV applications because of the changeable load power and variable atmospheric environments [29]. The trade-off choice between fast response and reactants consumption at high altitude is considered in [64].

The temperature of the reactants is regulated in the heat exchanger by adjusting the cooling water flow rate. The regulation of the temperature is quite simple for hydrogen since its temperature is almost constant during the operation, while in the case of air cooling, the temperature of the air is strictly correlated with the operating point of the compressor [36]. Gas humidity control is obtained by varying the flow rate of the water in the humidifier to achieve saturation.

A detailed review of control methods applied to CCPEMFCs is beyond the goal of this investigation. Nevertheless, it is noteworthy to report some examples for two reasons. On the one hand, the dynamic behavior of the fuel cell is strictly related to the control actions performed on the system. On the other hand, the number of recent papers on advanced control techniques shows the importance of developing accurate and comprehensive control-oriented models. To deal with the non-linear behavior of a CCPEMFC, fuzzy-PID techniques, adaptive methods, and sliding mode controllers have been proposed in the scientific literature [36,51,75,76]. A cascade controller based on disturbance and uncertainty estimation is proposed by Liu et al. [51]. As pointed out in [76], conventional PID controllers are not able to deal with the complexity and non-linearity of the air supply system. A model reference adaptive control is proposed in [36] to regulate the air temperature through the flow of cold water. In [4], two different approaches are compared to define the optimal setpoints. In the first case, the setpoints are obtained by optimizing the system offline under steady-state conditions. In the Model Predictive Control case, the optimal setpoints are continuously computed using an incoming time horizon. The dynamic effects of the systems are classified as fast and slow dynamics. The potential for energy savings during transient operating conditions by controlling the overshoot of OER is addressed in [76].

4. Quasi-Static Models

The goal of a quasi-static model of a CCPEMFC is to calculate hydrogen consumption by using the net power of the system as the main input variable, which must meet the load power at any time. Modeling the behavior of a fuel cell is a challenging task because of the multivariate, non-linear, and complex nature of such a system [32]. According to [79], five interacting sub-models (stack voltage, manifold flow, anode flow, cathode flow, and membrane hydration) should be considered. The sub-models for the voltage-current curve are the same for open-cathode and closed-cathode configurations [30], while the control of the reactant flows, the analysis of pressure dynamics, and the thermal balance are quite different in the two cases. In all models retrieved in the scientific literature, the gases are treated as ideal gases. Moreover, for small fuel cells, the temperatures of the gases in the stack are considered uniform and equal to the stack temperature [79].

As mentioned in the introduction, different kinds of models have been proposed for modeling the stack: white, gray, and black models [32]. White models address the thermodynamic, electrochemical, and fluid dynamical processes that take place in the stack using complex algebraic and differential equations, like the Navier-stokes equations for fluid dynamics. Fluid dynamic processes can be solved with one-dimensional, bi-dimensional, or even three-dimensional approaches [80]. Computational Fluid Dynamic (CFD) models are valuable for studying the impact of geometric and technological parameters on fuel cell performance, like the optimization of cooling channels to reduce pumping power [81] or the analysis of the distribution of temperature, pressure, and concentration of reactants in the stack [82]. On the other hand, they require too much computational effort and are not suitable for the analysis and control of the whole fuel cell system, which is the focus of this review. Therefore, CFD models will not be analyzed here.

4.1. Stack Voltage

The goal of stack-level models is to estimate the quasi-stationary voltage-current or polarization curve of the fuel cell, which is a function of current, partial pressure of the reactants, temperature, water content of the membrane, and oxygen excess ratio [83].

The Nernst voltage is the theoretical voltage that the fuel cell can develop at a certain pressure of hydrogen and oxygen. When the fuel cell is connected to the load and a current is drawn, the voltage decreases because of reaction activation, crossover, internal currents, ohmic losses, and mass concentration limitations [1,6]. These phenomena take place at all values of the current, but as the current increases, activation, ohmic, and mass concentration limitations become more relevant in this order, as shown in Figure 6. Note that the open-circuit voltage is slightly lower than the Nernst voltage.

The operational voltage curves can be expressed in terms of current density/cell voltage or current/stack voltage. The first approach is used when the goal of the investigation is to evaluate or optimize the single cell, where the current is directly proportional to the active area or to highlight the technological limits of a fuel cell [1]. For a given stack, it is indifferent to express the voltage in terms of current ($I_{st}$) or current density (${i=I}_{st}/S_{act}$), where $S_{act}$ is the active area of the cell.

Equivalent resistors are used in equivalent electric circuit models for activation, ohmic voltage, and concentration voltage drops. The most common equivalent electric circuit for a PEMFC stack is shown in Figure 7. The model is usually referred to as a single cell, but by multiplying it by the number of cells in series (N), it is possible to obtain the voltage of the whole stack [84].

The scheme in Figure 7 is used in several works in literature [6,12,49,66,83,85,86] to evaluate the stack voltage as follows:

$$V_{st}=N\left(E_{Nernst}-∆V_{ohm}-∆V_{act}-∆V_{conc}\right).$$

(8)

In this equation, the voltage corrections are expressed in terms of voltage drop $∆V$ due to ohmic $(ohm$), activation $\left(act\right),$ and concentration $\left(conc\right)$ losses. In the scientific literature however, these are also called overvoltage, polarization, irreversibility, or losses. Sometimes they are also expressed as equivalent resistances.

The voltage of the stack is assumed to be the same as that of a single cell multiplied by the cell number $N$. However, since the temperature of the stack is not uniform [87], this assumption cannot be correct. To account for this effect, a voltage fluctuation rate is expressed and measured in terms of fluctuation rate $C_v$:

$$C_v=\sqrt{\frac{{{\sum }_{i=1}^N\left(\frac{V_i-V}{V}\right)}^2 }{N}·100\%}$$

(9)

where $V$ is the average cell voltage and $V_i$ is the voltage of cell $i$.

The authors of [87] found that the cell voltage tends to be more uniform at low power than at high power and when working at high temperatures because of membrane dehydration.

4.1.1. Nernst Voltage

The Nernst voltage is commonly expressed as:

$$E_{Nernst}=E_0+\frac{R{T}_{st}}{2F}\left[\ln \left(p_{H_2}\right)+\frac{1}{2}\ln \left(p_{O_2}\right)\right]$$

(10)

where $R/2F$ can be approximated as $4.308·{10}^{-5}\frac{J}{C·K}$.

The partial pressures of the reactants in Equation (10) can be expressed as $p_{H_2}=p_a{x}_{H_2,a}$ and $p_{O_2}=p_c{x}_{O_2,c}$, where the under scripts $a$ and $c$ for pressure $p$ stand for anode and cathode, respectively, while $x_{H_2,a}$ and $x_{O_2,c}$ are the molar fractions of hydrogen at the anode and of oxygen at the cathode. The partial pressure of the cathode air is usually considered constant in open-cathode configurations, while in a CCPEMFC, the pressure at the anode and at the cathode are controlled to optimize the fuel cell performance. In CCAPEMFCs, the pressure of the oxygen can be easily regulated.

Equation (10) is valid for $T_{st}\le 100 °\mathrm{C}$. For higher temperatures, the partial pressure of water vapor must be considered too:

$$E_{Nernst}=E_0+\frac{R{T}_{st}}{2F}\left[\ln \left(p_{H_2}\right)+\frac{1}{2}\ln \left(p_{O_2}\right)-\ln \left(p_{H_{2}O}\right)\right].$$

(11)

$E_0$ is the reversible single cell open circuit voltage at standard pressure, which can be assumed constant (equal to 1.2 V according to Mahjoubi et al. 2019 [66]) or expressed as a function of the stack temperature [87,88]:

$$E_0=1.229-8.46·{10}^{-4}\left(T_{st}-298.15\right)$$

(12)

where the first term (1.229 V) is the voltage obtained by the electrochemical conversion of the free reaction enthalpy at 298 K, while the second term is associated with the reaction entropy (see [12,49] for details). However, the effect of temperature on the Nernst voltage was found to be negligible when the stack temperature is between 30 and 60 °C [31].

4.1.2. Activation Losses

Different approaches are used in the literature to estimate the activation loss $∆V_{act}$. The first one is based on the Tafel equation [6]:

$$∆V_{act}=\frac{RT}{2\alpha F}ln\left(\frac{I_{st}}{S_{act}{ I}_0}\right)$$

(13)

or its modified form:

$$∆V_{act}=\frac{RT}{2\alpha F}ln\left(\frac{I_{st}+I_n}{S_{act}{ I}_0}\right)$$

(14)

where ${ I}_0$ is the exchange current density, i.e., the minimum current at which the overvoltage begins to move from zero. In fact, this equation is valid only for $I_{st}>I_0$ [6]. The parameter $\alpha$ is the transfer coefficient that accounts for the part of electrical energy harnessed in changing the rate of the electrochemical reaction. $I_n$ neglected in [66], accounts for fuel crossover and internal currents, i.e., fuels and electrons that diffuse through the electrolyte. The crossover term can be determined empirically or calculated experimentally.

Note that activation polarization loss has two dependencies on temperature. In fact, these losses increase linearly with temperature, as shown in Equation (13), but are indirectly affected by the temperature through the exchange current density that increases exponentially with temperature [31]. An increase in temperature of 2 °C, corresponding to an increase in ${ I}_0$ by 17%, caused an improvement of the cell voltage of 86 mV.

As shown in [88], the exchange current can be expressed as a function of the temperature:

$$\frac{{ I}_0}{S_{act}}=1.256-8\exp [0.10988(T_{st}-323.15)].$$

(15)

Note that the improvement in the activation losses due to increased temperature is largely guaranteed by the higher exchange current density.

In [50], the temperature dependence of the exchange current density is approximated with a quadratic fitting curve:

$$\frac{{ I}_0}{S_{act}}=5.52\times {10}^{-7}{T}_{st}^2-3.21\times {10}^{-4}{T}_{st}+0.04674.$$

(16)

The exchange temperature is also dependent on the square root of the oxygen partial pressure, as pointed out in [31]. In [44], the exchange current is calculated as a function not only of the temperature but also of the catalyst specifications and the electrode roughness. The transfer coefficient $\alpha$ depends on temperature too and can be calculated as:

$$\alpha ={\alpha }_0+2.3\times {10}^{-3}(T_{st}-298.15)$$

(17)

with ${\alpha }_0=0.175$.

Using pure O₂ instead of air (CCAPEMFC) is a way to reduce the activation voltage drop, but this is not evident in the Tafel equation. Moreover, the Tafel approximation (Equation (14)) is valid only when $I_{st}>I_0$. When $I_{st}\le I_0$, a detailed modeling the low-current region is proposed in [1] with the Butler-Volmer equation. However, Dickinson et al. [89] discussed this equation and showed that it is unsuitable for PEMFCs and offers no advantages over empirical approaches that express the activation losses at the anode (A) and cathode (C) as:

$$∆V_{act}=a_A+b_A\ln \left(\frac{I_{st}}{S_{act}}\right)-a_C-b_C\ln \left(\frac{I_{st}}{S_{act}}\right)$$

(18)

with $a_A$, $b_A$, $a_C$, $b_C$ empirical variables to be tuned.

In fact, activation losses take place at both the cathode and anode even if cathode activation is usually neglected because the anode kinetics are fast compared to cathode activation kinetics [1,38].

In [35], the authors use the following empirical equation for the activation loss:

$$∆V_{act}=∆V_{act,1}+∆V_{act,2}={\eta }_0+a\left(T_{st}-298\right)+b{T}_{ST}\mathrm{l}\mathrm{n}\left(I_{st}\right)$$

(19)

where incremental $∆V_{act,1}$ is a term of loss affected only by the fuel cell internal temperature, while $∆V_{act,2}$ includes the effect of both current and temperature; ${\eta }_0$ is an empirical parameter.

A dependence of activation loss on oxygen partial pressure at the catalyst/gas interface is put into evidence in [88] by proposing a slightly different expression for $∆V_{act}$.

The effect of oxygen concentration is also accounted for in the so-called Amphett’s model [90], and the following empirical correlation is used [49,91,92,93] for the activation loss:

$$∆V_{act}={\left(R_{act1}+R_{act2}\right)I_{st}=\xi }_1+{\xi }_2{T}_{st}+{\xi }_3{T}_{st}\ln \left(c_{O_2}\right)+{\xi }_4{T}_{st}\ln \left(I_{st}\right)$$

(20)

where ${\xi }_i$ with $i=1, 4$ are empirical parameters obtained by fitting experimental data. $c_{O2}$ denotes the concentration of oxygen:

$$c_{O_2}=\frac{P_{O2}}{5.08·{10}^6\exp \left(-\frac{498}{T_{st}}\right)}$$

(21)

Even if the terms ${\xi }_i$ originally had a physical interpretation in the so-called Amphlett’s model [90], these parameters are treated as empirical parameters to be tuned. From the analysis of the scientific literature [32,40,49,79,83,91,92,94,95,96,97,98,99], the lower and upper bounds of

Table 1. Upper and lower bounds of the empirical parameters of Equation (20).

Parameter	${\mathit{\xi }}_1$	${\mathit{\xi }}_2$	${\mathit{\xi }}_3$	${\mathit{\xi }}_4$
Upper bound	−0.801	0.00386	$9.3·{10}^{-5}$	$-0.954·{10}^{-4}$
Lower bound	−1.1997	0.002283	$3.4·{10}^{-5}$	$-2.6·{10}^{-4}$

are found for the parameters of the activation loss Equation (20).

Note that values of this parameter outside the bounds of Table 1 were obtained by [38]. However, in that case, the parameters were tuned using dynamic data together with many variables, including the dynamics of reactants and the action of the cooling system.

In [92,98], the second coefficient is calculated as proposed by [95] in an attempt to generalize the model:

$${\xi }_2=0.00286+0.002\ln \left(S_{act}\right)+4.3·{10}^{-5}+\ln c_{H2}$$

(22)

where $c_{H2}$ is the hydrogen concentration calculated as:

$$c_{H_2}=9{.17·{10}^{-7}P}_{H_2}\exp \left(-\frac{77}{T_{st}}\right).$$

(23)

After tuning the model over a Ballard fuel cell stack, Lee et al. [100] tried to generalize the model to other stacks but found it did not accurately predict the effect of temperature on the voltage-current curve of fuel cells developed by other manufacturers.

Recently, a different approach has been proposed in [85] where the activation loss is expressed as:

$${∆V}_0+∆V_1[1-\exp (-c_1{I}_{st})]$$

(24)

where $c_1$, ${∆V}_0$, and $∆V_1$ are empirical parameters that depend on the stack temperature, oxygen partial pressure, cathode pressure, and saturated vapor pressure (see [85] for details). Despite the complexity of the model, the comparison of the modeled and experimental voltage curves of the 500 W fuel cell tested in [85] reveals a large error in the activation region.

As pointed out in [6], it is possible to assume that when the pressure is changed from $p_1$ to $p_2$, the variation in the cell voltage is equal to:

$$∆V_{gain}=c·\ln \left(\frac{p_2}{p_1}\right)$$

(25)

where $c$ is a constant that depends on how the exchange current density is affected by pressure, fuel cell design, and humidification. A value of 0.06 V is a realistic value according to [6].

4.1.3. Ohmic Losses

Ohmic losses are universally expressed through an equivalent resistance $R_{ohm}$:

$$∆V_{ohm}=R_{ohm}{I}_{st}.$$

(26)

The ohmic resistance includes two effects: the resistance met by the ions through the membrane, $R_{mem}$, and the contact resistance $R_c$ between different PEM fuel cell layers [66]:

$$R_{ohm}=\left(R_{mem}+R_c\right)$$

(27)

where $R_{mem}$ depends on some specification of the membrane, like the thickness $t_m$ and operational parameters like ${ \lambda }_m$ (the membrane water content), the stack temperature ${ T}_{st}$, and the current density [49,86,98]:

$$i=\frac{{ I}_{St}}{S_{act}}$$

(28)

$$R_{mem}=\frac{181.6t_m·\left[1-0.03i+0.062{\left(\frac{{ T}_{st}}{303}\right)}^2{i}^{2.5} \right]}{\left[{\lambda }_m-0.634-3\frac{{ I}_{St}}{S_{act}}\right]·S_{act}·\exp \left(4.18\frac{T_{st}-303}{T_s}\right)}$$

(29)

This equation is sometimes expressed in terms of membrane-specific resistivity $R_{mem}/{t_m·S}_{Act}$ [95].

From the literature data collected in [32], ${\lambda }_m$ varies between 10 and 24, while $R_c$ is in the interval [$8·{10}^{-5},8·{10}^{-4}]$. According to some authors (see [40,70], the membrane water content can be calculated as:

$${\lambda }_m=\left\{\begin{array}{ll}0.043+17.81a_{H_{2}O}-39.85a_{H_{2}O}^2-36.0a_{H_{2}O}^3& for 0<a_{H_{2}O}<1\\ 14+1.4\left(a_{H_{2}O}-1\right)& for 1\le a_{H_{2}O}\le 3.\end{array}\right.$$

(30)

For $a_{H_{2}O}>3$, a value of 16.8 is suggested in [101]. $a_{H_{2}O}$ is the water activity in the membrane:

$$a_{H_{2}O}=\frac{C_{w}R{T}_{st}}{p_{sat}}$$

(31)

where $C_w$ is the molar concentration of water in the electrolyte membrane. In other words [84], the water activity is the ratio of water partial pressure to saturation pressure $a_{H_{2}O}=p_{H_{2}O}/p_{sat}$. The saturation pressure dependence on the stack temperature can be empirically expressed as:

$${\log }_{10}{p}_{sat}=-2.1794+0.02953·T_{St}-9.1837·{10}^{-5}·T_{st}^2+1.4454·{10}^3·T_{st}^3.$$

(32)

The membrane water content is also dependent on the temperature as shown in [88], which correlated the membrane water content to the water activity. There is a tradeoff between increasing temperature and water content to minimize ohmic losses. In fact, the membrane water content depends on the water saturation pressure, which in turn is a function of the stack temperature.

In [83] the ohmic loss is expressed as:

$$∆V_{ohm}=\frac{t_m}{{\sigma }_m}$$

(33)

where $t_m$ is the membrane thickness, and ${\sigma }_m$ is the ion conductivity, which is a function of the water content of the membrane and of the stack temperature:

$${\sigma }_m=\left(b_1{\lambda }_m-b_2\right)\exp \left[b_3\left(\frac{1}{303}-\frac{1}{T_{st}} \right)\right].$$

(34)

The values of the empirical parameters $b_i$ are not reported in [83], while in [85], the following values are used: $b_1=0.005139$, $b_2=0.00326$, $b_3=350$. The same correlation is also used in [44].

In [32,91], the internal resistance is expressed as:

$$R_{ohm}=-\frac{I_{st}}{S_{act}}({\zeta }_1+{\zeta }_2{T}_{st}+{\zeta }_3\frac{I_{st}}{S_{act}}).$$

(35)

The values suggested for the empirical coefficients ${\zeta }_1,{\zeta }_2$, and ${\zeta }_3$ are $0.0375$, $1·{10}^{-6}$, and $-0.00027$ [91]. Even with a different notation, a similar approach is used by Wang et al. [35].

As pointed out in [6], the internal resistance is also dependent on the ambient pressure, but this effect is not accounted for in the models available in the scientific literature.

4.1.4. Concentration Losses

The loss of voltage due to mass transport is often expressed as [49,66,83,86,91,92]:

$$∆V_{conc}=\mathrm{B} \ln \left(1-\frac{I_{st}}{I_{max}}\right)$$

(36)

where B is assumed to be equal to $\frac{R{T}_{st}}{2F}$ in [35,66,91].

$I_{max}$, also called the limiting current [83], is the current at which the reactant concentration at the surface of the cell falls to zero. As pointed out in [32,102], $I_{max}$ is treated as a constant in the scientific literature; however, it is highly dependent on degradation phenomena and operating conditions.

According to [98,103], $B=\frac{R{T}_{st}}{2F}\left(1+\frac{1}{\alpha }\right),$ where $\alpha$ is the charge transfer coefficient that describes the amount of electrical energy applied to change the rate of the electrochemical reaction. In [56], the term $\mathrm{B}$ is expressed as a function of two more tunable parameters ${\mathsf{\alpha }}_{\mathrm{B}}$ and $G$:

$$\mathrm{B}={\mathsf{\alpha }}_{\mathrm{B}}{I}_{st}^G.$$

(37)

Since fuel cells exhibit a voltage drop larger than that predicted by Equation (36), $B$ is often treated as an empirical variable [49,83,91,92,93].

The maximum current density can be expressed as [84]:

$$\frac{I_{max}}{S_{act}}=2F{h}_m\left(\frac{C_{in}-C_{out}}{\mathrm{l}\mathrm{n}\frac{{\mathrm{C}}_{\mathrm{i}\mathrm{n}}}{C_{out}}}\right)$$

(38)

where $h_m$ is the mass transfer coefficient, and $C_{in}$ and $C_{out}$ are the concentrations of the reactants at the inlet and outlet, respectively.

The maximum current is evaluated in [38] as:

$$\frac{I_{max}}{S_{act}}=2F·D·\frac{C_B}{\delta }$$

(39)

where $D$ is the effective diffusion coefficient of the reacting species, $C_B$ is its bulk concentration, and $\delta$ the thickness of the diffusion layer.

From the literature data collected in [32], $B$ varies between 0.01 and 0.355 V, and $I_{max}/S_{act}$ ranges between 0.50 and 1.5 $\mathrm{A}/\mathrm{c}{\mathrm{m}}^2$.

An alternative way to express the concentration losses is given in [12]:

$$∆V_{conc}=m·\exp (n{I}_{st}/S_{act})$$

(40)

with $m=0.095741V$, and $n=0.0076301\mathrm{c}{\mathrm{m}}^2/\mathrm{A}$.

In some studies ([68,85,104]), the approach proposed by Guzzella [48] is adopted:

$$∆V_{conc}=-\frac{Ist}{S_{act}}{\left(\frac{c_2{I}_{st}}{I_{max}}\right)}^{c_3}.$$

(41)

The concentration losses are evaluated in [45] as:

$$∆V_{conc}=m·\exp \left(N{I}_{st}\right)+b·\ln \left(\frac{P_{O_2}}{a}\right)$$

(42)

where $m$, $b,$ and $a$ are empirical coefficients.

It is necessary to point out that concentration losses can be neglected if the fuel cell is not allowed to work in its high-current region where this kind of loss is most relevant [105].

4.1.5. Parametric Models

From the analysis above, it is evident that the performance of a fuel cell stack is affected by a large number of geometrical and operational parameters. There is uncertainty about some geometrical parameters, like the effective area. In fact, its value is used as a parameter to model the degradation of the fuel cell in [36]. In the absence of data about the geometrical features, the voltage curve can be estimated with parametric models where only the main operational variables are considered. The operational variables include pressures, humidity, control of reactants, stack temperature, and water content [40].

One of the most commonly used simulation approaches is the so-called Amphlett’s model, which consists of Equations (11) and (12) for the Nernst voltage, Equations (20) and (21) for the activation polarization, Equation (35) for the ohmic losses, and Equation (36) for the concentration voltage drop. In Larminie and Dicks [6], the following parametric model is proposed:

$$V=E_0-A ·ln\left(i\right)-R_i·i-m·exp\left(ni\right)$$

(43)

where all the five parameters $E_0$, $A$, $R_i$, $m$, and $n$ can be treated as tunable variables even if they are a function of the fuel cell operating conditions (temperature, humidity, pressure, oxygen concentration, etc.) as evident from the more detailed models presented above. A slightly different version can be found in Lee [100] where the effect of pressure is added:

$$V=E_0-A ·ln\left(i\right)-R_i·i-m·exp\left(ni\right)-C·ln\left(\frac{P}{P_{O_2}}\right).$$

(44)

The constant $C$ is assumed to be equal to A in [100], while in [46] a value between 0.03 and 0.06 V at the cell level is suggested (see also [6]).

As pointed out in [32], the number of parameters in the semi-empirical models ranges from 3 [41] to 22 [38]. In particular, nine parameters are needed for Amphlett’s model [95], four for Squadrito’s [54], and six for Lee’s [95]. A sensitivity analysis of the tuning parameters of Equation (43) can be found in [46]. This investigation shows how the results are strongly affected by each parameter. In particular, the Nernst voltage $E_0$ and the internal resistance present the highest impact, but also the parameter $A$ that quantifies the activation losses has a relevant effect on the whole voltage-current curve. Other parameters like $m$ and $n$ affect only the high-current part of the polarization curve where the concentration losses play a significant role. A similar analysis is performed in [47] for parameters ${\xi }_i$ of Equation (20) where the authors pointed out how numerical and approximation errors are introduced if the tuning of the parameters is performed improperly. In fact, due to this large sensitivity of the model, it can be concluded that a polarization curve obtained experimentally can be fitted with a different combination of the parameters of the same model if the tuning is performed on a single curve, all the more so if obtained under unspecified operating conditions. This can lead to a perfectly tuned model but to unrealistic values of the losses estimated by the model. This means that the model is valid only for the specific operating condition used in the test and has no prediction capability.

To solve this problem, online identification, based on the processing of data in real-time, has been suggested in assessing fuel cell degradation and aging. In [32], an online identification method is proposed for a fuel cell system to counteract undermines due to fuel cell aging for integration with energy management in hybrid electric power systems. The paper considers three inputs (current, temperature, and pressure) and eight tuning parameters: ${\xi }_i$ with $i=1, 4$ for Equation (20), ${\zeta }_1,{\zeta }_2,{\zeta }_3$ for the internal resistance, Equation (35), and ${\mathsf{\alpha }}_{\mathrm{B}}$ in the concentration losses, Equation (37). The output is the fuel cell voltage. An adaptive identification method is used in [102] to correct the maximum current $I_{max}$ which changes over time owing to degradation and operating conditions variation. Note that $I_{max}$ affects the concentration losses (see Equation (36)), and, consequently, the maximum gross power that the fuel cell can produce.

4.2. Cooling System

Under stationary conditions, the temperature of the stack is established by the flow rate of cool water. With reference to the control volume of Figure 3, the required mass flow rate of the pump is computed as:

$${\dot{m}}_{cool}=\frac{{\dot{\mathrm{Q}}}_{cool}}{c_{cool}(T_{cool,out}-T_{cool,in}) } \cong \frac{{\dot{\mathrm{Q}}}_{cool}}{c_{cool}(T_{st}-T_{cool,in}) }$$

(45)

where ${\dot{\mathrm{Q}}}_{cool}$ is the heat to be removed from the stack, and $c_{cool}$ is the specific heat of the cooling fluid. $T_{cool,in}$ is the inlet temperature of the cooling fluid, while $T_{cool,out}$ is the outlet temperature, often considered equal to the stack temperature.

The power absorbed by the cooling pump motor is computed as [41]:

$$P_{aux, pump}=\frac{{\dot{m}}_{cool}\mathsf{\Delta }{p}_{pump}}{{\eta }_{pump}}$$

(46)

where $\mathsf{\Delta }{p}_{pump}$ is the head required to move the cooling water in the cooling system because of major and minor head losses, ${\dot{m}}_{cool}$ is the mass flow rate of the coolant fluid (generally water), and ${\eta }_{pump}$ is the pump efficiency inclusive of its motor.

In the most general form, ${\dot{\mathrm{Q}}}_{cool}$ can be calculated as:

$${\dot{\mathrm{Q}}}_{cool}={\dot{\mathrm{Q}}}_{gen}-{\dot{\mathrm{Q}}}_{nat}-{\dot{\mathrm{Q}}}_{rad}-{\dot{\mathrm{Q}}}_{sens+lat}+\sum _i{\dot{m}}_{i,in}{c}_{p,i,in}{T}_{i,in}-{\dot{m}}_{i,out}{c}_{p_i,out}{T}_{i,out}$$

(47)

where ${\dot{m}}_{i,in}{c}_{p,i,in}{T}_{i,in}$ and ${\dot{m}}_{i,out}{c}_{p_i,out}{T}_{i,out}$ represent the mass flow entering and leaving the stack. The subscript $i$ is used to represent the flows at the anode and the cathode. The different heat contributions are explained below.

The heat generated from the chemical reactions taking place within the fuel cell stacks can be expressed as ${\dot{\mathrm{Q}}}_{gen}=\left(E_{nernst}·N-V_{st}\right)·I_{st}$ [31,66], or simply as ${\dot{\mathrm{Q}}}_{gen}=\left(1.25·N-V_{st}\right)I_{st}$ [40]. In [35,83], the generated heat is calculated as ${\dot{\mathrm{Q}}}_{gen}={\dot{\mathrm{Q}}}_{reac}-P_{st}$ where ${\dot{\mathrm{Q}}}_{reac}=\frac{N{I}_{st}}{2F}∆G$ is the total energy released from the electrochemical reaction, while $P_{st}=V_{st}{I}_{st}$ is the gross electric power generated from the stack. In [83,84], the reaction heat is calculated from the enthalpy variation as ${\dot{\mathrm{Q}}}_{reac}=\frac{N{I}_{st}}{2F}∆H$.

The calculation of the reaction heat flux is detailed in [85]:

$${\dot{\mathrm{Q}}}_{reac}={\dot{m}}_{H_2,reac}\left(H_{f,H_2}^0+\mathsf{\Delta }{H}_{H_2}\right)+{\dot{m}}_{O_2,reac}\left(H_{f,O_2}^0+\mathsf{\Delta }{H}_{O_2}\right)-{\dot{m}}_{H_{2}O,gen}\left(H_{f,H_{2}O}^0+\mathsf{\Delta }{H}_{H_{2}O}\right)$$

where $H_f^0$ is the standard enthalpy of formation (zero for elements like $H_2$ and $O_2$).

In [86], the reaction heat is computed with the following equations where $M_j$ is the molar mass of species $j$:

$${\dot{\mathrm{Q}}}_{reac}={\dot{m}}_{H_2,reac}[H_{f,H_2}^0+c_{H_2}\left({\mathrm{T}}_{\mathrm{a}\mathrm{m}\mathrm{b}}-T_0\right)]+{\dot{m}}_{O_2,reac}[H_{f,O_2}^0+c_{O_2}\left({\mathrm{T}}_{\mathrm{a}\mathrm{m}\mathrm{b}}-T_0\right)]-{\dot{m}}_{H_{2}O,gen}[H_{f,H_{2}O}^0+c_{H_{2}O}\left({\mathrm{T}}_{\mathrm{a}\mathrm{m}\mathrm{b}}-T_0\right)]$$

(48)

$${\dot{m}}_{H_2,reac}=\frac{M_{H_2}N{I}_{st}}{2F}$$

(49)

$${\dot{m}}_{O_2,reac}=\frac{M_{O_2}N{I}_{st}}{4F}$$

(50)

$${\dot{m}}_{H_{2}O,gen}=\frac{M_{H_{2}O}N{I}_{st}}{2F}$$

(51)

${\dot{\mathrm{Q}}}_{nat}$ is the natural convection heat flux through the fuel cell surface exposed to the external environment and can be expressed as ${\dot{\mathrm{Q}}}_{nat}=h_{nat}·S_{nat}·\left(T_{st}-T_{amb}\right)$ [49], while the radiant flow rate ${\dot{\mathrm{Q}}}_{rad}$ can be calculated by [86]:

$${\dot{\mathrm{Q}}}_{rad}=c_1{T}_{st}^2-c_2{T}_{ca,in}^2$$

(52)

where $T_{ca,in}^{}$ is the temperature of the air entering the cathode. In Wang et al. [85], the following values are used for the parameters: $c_1=2.179·{10}^{-3}$, $c_2=2.436·{10}^{-3}$.

According to Musio et al. [84], the contribution of natural convection and radiation is an order of magnitude lower than that of forced convection. In Strahl et al. [31], heat loss due to natural convection and radiation is found to be around 10% of the total waste heat. Therefore, the radiant flow rate and the natural convection contributions are often considered negligible. ${\dot{\mathrm{Q}}}_{sens+lat}$ is the sensible and laten heat absorbed during the process. It includes ${\dot{\mathrm{Q}}}_{l/g}$, the heat removal by phase change from liquid to gas of the water produced inside the PEM fuel cell [66]. This contribution is also considered negligible.

In automotive applications (Figure 2), the thermal power is removed by water with the help of cooling fans operated with on/off control, as in internal combustion engines [41]. This introduces a further contribution to parasitic power when the fan is turned on, which can be calculated with an expression like Equation (46). The parasitic power absorbed by the fan can be assumed to be linearly proportional to the coolant mass flow rate and proportional to the square of the stack current, according to Guzzella et al. [41].

In aerospace applications (Figure 3), the fan is not required because the required air flow rate is obtained from an air inlet positioned in the frontal part of the nacelle (see [106]). The airflow through the heat exchanger (HEX) is subject to a pressure loss that can be expressed empirically as [42]:

$$∆p_{HEX}=k_1{\dot{m}}_{air,intake}^2+k_2{\dot{m}}_{air,intake}+k_3$$

(53)

where ${\dot{m}}_{air,intake}$ is the intake flow rate of atmospheric air, and $k_3$ is an empirical parameter.

During steady-state fuel cell operation, the pressure loss in the HEX is equal to the pressure gain provided by the intake air:

$$∆p_{intake}=\frac{1}{2}{\rho }_{amb}{V}_{\infty }^2+\frac{T}{\pi R^2}-\frac{1}{2}{\rho }_{amb}{\left(\frac{{\dot{m}}_{air,intake}}{S_{HEX}}\right)}^2$$

(54)

where $V_{\infty }$ is the aircraft speed, $T$ is the thrust provided by the propeller (with radius $R$), $\frac{{\dot{m}}_{air,intake}}{S_{HEX}}$ is the speed of flow that actually enters the heat exchanger with the inlet surface $S_{HEX}$.

4.3. Anode Sub-System

The anode subsystem (magenta circuit in Figure 3) provides the flux of hydrogen from the tank to the fuel cell. The pressure at the anode inlet needs to be accurately regulated because it must not exceed the cathode pressure by a certain quantity (0.2 bar in [4]) to avoid damages.

In equilibrium conditions, the flow rate of hydrogen, stoichiometric oxygen, and water are directly related to the stack current [6,40]. According to the basic reaction (1), the theoretical consumption of hydrogen is expressed by Equation (49). To obtain a uniform feeding of hydrogen, a larger mass flow is sent to the anode. The hydrogen excess ratio (HER) [78] is used to account for this:

$${\dot{m}}_{H_2, in}=HER \frac{M_{H_2}}{2F}·I_{st}.$$

(55)

The excess of hydrogen is often named “fuel utilization” and represented with the symbol ${\mu }_F$, which is the ratio of the fuel used by the cell to generate electric current versus the total fuel provided to the fuel cell according to [6].

Four possible hydrogen delivery schemes can be adopted for the cathode: flow through, recirculation, dead-end, and recirculation with a dead-end. In the first case, the excess hydrogen is discharged into the atmosphere, leading to a waste of fuel. In the recirculation scheme, the anode outlet gases pass through a separator to eliminate water and impurities and then are recirculated to the anode inlet. This solution increases the complexity and the parasitic power due to the need for a circulation pump [107]. In dead-end systems like those considered in this investigation, the anode outlet is sealed, and the accumulated water and nitrogen in the anode are eliminated by the periodical opening of a purge valve. In fact, nitrogen and water are generated at the cathode but diffuse back to the anode because of the concentration gradient and accumulate in the dead-end anode, causing channel clogging and voltage reduction (because hydrogen delivery is prevented) and stack degradation [12].

The volumetric flow rate of gas through the purge valve, when open, $\dot{V}$, depends on the specification of the purge valve. In [108] it is calculated as:

$$\dot{V}=\frac{5140}{3600}{k}_V\sqrt{\frac{{(p}_{an}-p_{amb}) p_{amb} }{\rho T_{st} }}$$

(56)

where $k_V$ is the valve constant, while the gas density $\rho$ is calculated from the molar fraction of nitrogen, hydrogen, and water vapor. The nitrogen molar fraction was found to reach 0.14 when the cell voltage decreases by 0.1 V during the dead anode mode. When the valve is open, the gaseous mixture is purged together with liquid vapor. Therefore, the fraction of water is calculated by the saturated water vapor.

The hydrogen pressure drop is expressed as a function of the squared hydrogen flow rate:

$$\mathsf{\Delta }{p}_{H_2}=k_f{\dot{V}}_{H_2}^2.$$

(57)

A value $k_f=2.22·{10}^{-3}\mathrm{a}\mathrm{t}\mathrm{m} \mathrm{m}{\mathrm{i}\mathrm{n}}^2/\mathrm{n}{\mathrm{L}}^2$ is used in [99].

The purging process affects the fuel utilization ${\mu }_F$ particularly at low load. To filter and smooth the peaks in the fuel-flow rate due to purging, a moving average over a sample range of 10% of the total number of measurements is performed in [56] for each operating condition. In this way, a fuel utilization factor of 90% is obtained in a wide range of loads (see Figure 8). Optimization of the purging process is proposed in [109,110].

To recover the loss of hydrogen from the purging valve, this valve can be combined with a recirculation circuit. The flow rate of hydrogen is in part lost because of the purging process, and in part recirculated to the control valve at the anode inlet in the most complex configuration of Figure 2. An optimization of the purging valve opening frequency must be performed to minimize the losses of hydrogen. The subsystem adopted in [4] also presents a water trap that collects liquid and a drain valve that discharges the water collected in the water trap.

In the cases of recirculation, the inlet mass flow rate is controlled through the hydrogen supply valve and is the sum of a quantity recirculated from the anode outlet, ${\dot{m}}_{H_2,rec}$, and a quantity extracted from the hydrogen tank. Therefore, the outlet flow rate of hydrogen from the anode can be expressed as the difference between ${\dot{m}}_{H_2, in}$ and ${\dot{m}}_{H_2, reac}$.

A circulation pump is needed to compensate for head losses in the anode circuit that determine a reduction of the pressure at the anode outlet ($p_{an,out}$) with respect to the inlet value $p_{an,in}$. The power absorbed by the recirculation pump, $P_{hp}$, can be approximated as [41]:

$$P_{hp}={\dot{m}}_{H_2,rec}\left(\frac{p_{an,in}-p_{an,out}}{{\rho }_{H_2}{\eta }_{hp}}\right)$$

(58)

where ${\rho }_{H_2}$ is the average hydrogen density, and ${\eta }_{hp}$ is the overall efficiency of the recirculation pump and its drive motor. Since the pressure loss $p_{an,in}-p_{an,out}$ in a laminar flow is proportional (through the constant $k_{hp}$) to the flow rate ${\dot{m}}_{H_2,rec}$, which, in turn, depends on the stack current, it is possible to write:

$$P_{hp}=k_{hp}{I}_{st}^2\left(t\right).$$

(59)

4.4. Cathode Subsystem

The air supply system of Figure 3 includes a supply manifold, a compressor, an electric motor driving the compressor, a return manifold, and a control valve. Such a backpressure valve is used to control the cathode pressure [4]. Under stationary conditions, the cathode pressure is constant, and the flow rate in all these component is the same (${\dot{m}}_{air,in}={\dot{m}}_{cp}\cong {\dot{m}}_{rm,out}$).

The consumption of oxygen is given by Equation (50). However, since oxygen comes from the air, it is necessary to refer to the stoichiometric air [6]:

$${\dot{m}}_{air,st}=3.57·{10}^{-7}{P}_e/V_{cell}$$

(60)

or:

$${\dot{m}}_{air,st}=\frac{5N}{4F}·I_{st}.$$

(61)

To optimize the stack performance, an excess of oxygen is adopted. The OER is defined as the actual air flow into the cathode divided by ${\dot{m}}_{air,st}$:

$${\dot{m}}_{air,in}=OER·{\dot{m}}_{air,st}.$$

(62)

Values of OER below one determines starvation. An OER of 2 is usually considered suitable for transport applications; however, recent studies have highlighted the need to optimize OER according to operating conditions to optimize net power and, consequently, the overall efficiency of the system.

Under quasi-static operation, the OER and the cathode pressure are kept at the desired values by selecting appropriate values for the speed of the compressor motor and the opening of the cathode pressure-control valve. The pressure at the cathode outlet can be calculated as:

$$p_{rm}=p_{amb}+{\xi }_a{\dot{m}}_{rm,out}^2$$

(63)

where ${\xi }_a$ is the opening of the pressure-control valve, and ${\dot{m}}_{air,out}^{}$ is equal ${\dot{m}}_{air,in}$ minus the flow rate of oxygen consumed by the reaction (Equation (50)).

The compressor working characteristic is described by a performance map, which provides the relationships of the speed, mass flow, pressure ratio $\frac{p_{cp}}{p_{amb} }$, and efficiency. A generic performance map is shown in Figure 9. The blue continuous line represents the surge limit. The pressure ratio curves corresponding to four different values of shaft speeds are shown as red continuous lines. The dotted iso-lines represent the values of compressor efficiency.

For a given speed, the mass flow rate of air, ${\dot{m}}_{air,in}$, and the supply manifold pressure are obtained by intersecting the compressor map with the pressure drops in the entire air circuit, defining the so-called “system characteristic curve”.

The mass flow rate must satisfy the OER requirement and be accurately controlled to avoid the compressor operating beyond the surge limit and becoming unstable [46]. Consequently, the compressor capacities limit the power delivery range of the whole system [29,111].

Compressor maps are often considered unaffected by ambient conditions. However, a more correct approach is to use the so-called universal compressor map that reports the pressure rate over the air mass flow for defined standard conditions. Under different operating conditions, the formula for mass flow and compressor speed correction can be applied [52].

The compressor needs the following electrical power to increase the air pressure from $p_{amb}$ to $p_{sm}$:

$$P_{comp}={\dot{m}}_{air,in}\frac{c_{p,air}{T}_{amb}}{{\eta }_{cp}}\left[{\left(\frac{p_{sm}}{p_{amb} }\right)}^{\frac{\gamma -1}{\gamma }}-1 \right]$$

(64)

where $p_{sm}$ is the pressure in the supply manifold, $T_{amb}$ and $p_{amb}$ are the ambient values of temperature and pressure at the working altitude, and ${\eta }_{cp}$ is the efficiency of the compressor, inclusive of its motor. The flow rate of the air from the compressor, ${\dot{m}}_{cp}$, and its overall efficiency ${\eta }_{cp}$ depend on the compressor speed and pressure ratio [36].

For a given value of the cathode pressure, the power absorbed by the compressor motor can be considered directly proportional to the stack current (through Equations (64) and (62)).

The outlet temperature of the compressor can be calculated as:

$$T_{cp}=T_{amb}+\frac{1}{{\eta }_{cp}}{T}_{amb}\left[{\left(\frac{p_{sm}}{p_{amb} }\right)}^{\frac{\gamma -1}{\gamma }}-1 \right].$$

(65)

Since this temperature can be quite higher than the stack temperature, a cooler is sometimes introduced in the supply manifold.

In some applications, an expander is used to recuperate part of the enthalpy from the exhaust air [41]. In this case, $p_{rm}$ coincides with the backpressure generated by the expander. However, this solution is seldom applied in the aerospace field because of the reduced enthalpy of the exhaust gas in a PEM FC because of its low working temperature. Moreover, such expanders have low efficiency and increase complexity and costs. To increase the enthalpy of exhaust gases, a hydrogen combustor is adopted in [44] before sending the exhaust gases to the turbine, while a heat exchanger is adopted to reduce the pressure of the compressed air, as shown in Equation (65). The temperature at the combustor outlet is calculated as:

$$T_{comb, out}=T_{comb,in}+\frac{{{\eta }_{comb}·\dot{m}}_{fuel}·LHV }{{\dot{m}}_{exh}·c_{p,exh}}$$

(66)

where $T_{comb, out}$ and $T_{comb,in}$ are the temperatures at the combustor outlet and inlet, respectively; ${\dot{m}}_{fuel}$ is the mass flowrate of hydrogen burnt in the combustor, $LHV$ is the lower heating value of hydrogen, ${\dot{m}}_{exh} \mathrm{a}\mathrm{n}\mathrm{d} c_{p,exh}$ are the mass flow rate and the specific heat of the exhaust gas. The combustion efficiency ${\eta }_{comb}$ is assumed to be equal to 0.98 in [44].

The power produced by the turbine is calculated as [44]:

$$P_{turb}={\eta }_{turb}{·\dot{m}}_{exh}·c_{p,exh}\left[1-{\left(\frac{p_{amb}}{p_{exh} }\right)}^{\frac{{\gamma }_{exh}-1}{{\gamma }_{exh}}}\right]$$

(67)

where $p_{exh}$ is the pressure of the exhaust gases at the turbine inlet, and ${\gamma }_{exh}$ is their heat capacity ratio. As in the case of the compressor, the isentropic efficiency of the turbine ${\eta }_{turb}$ can be calculated from the turbine maps as a function of the turbine speed and pressure ratio.

4.5. Water Circuit

The water circuit (Figure 2) includes the recirculation of the water separated from the anode outlet and the humidification of the air flow needed to control the humidity of the membrane. The humidifier also serves the additional function of reducing the temperature of the cathode flow rate after the increase in temperature caused by compression. In some cases, the hydrogen flow is also humidified. To improve the stability of humidification during dynamic load variations, different solutions can be adopted [53].

The amount of water produced by the chemical reaction at the cathode is given by Equation (51). The water vapor produced at the cathode is recirculated in the wet circuit of the humidified (${\dot{m}}_{H_{2}O,h}$). The power of the air-humified water pump, $P_{ahp},$ can be calculated as:

$$P_{ahp}=\frac{{\dot{m}}_{H_{2}O,inj}·\mathsf{\Delta }{p}_{ahp}}{{\rho }_{H_{2}O·}{\eta }_{ahp}}$$

(68)

where $\mathsf{\Delta }{p}_{ahp}$ is the pressure drop in the circuit, and ${\eta }_{ahp}$ is the efficiency of the water pump. ${\dot{m}}_{H_{2}O,inj}$ is the injected flow rate of water that is regulated to achieve the desired value of RH in the intake air and depends on the humidity of ambient air.

As already explained, such a circuit is absent in a typical aerospace application, thanks to the usage of self-humidified stacks (Figure 3).

4.6. Overall Parasitic Power

Bringing together the contributions of power absorbed by the compressor, the pumps, and the fan of a BOP like those shown in Figure 2, the total auxiliary power can be expressed as a function of the fuel cell current with the following empirical formula:

$$P_{aux}=P_0+k_{aux1}{I}_{st}+k_{aux1}{I}_{st}^2$$

(69)

where $P_0$ is an additional bias power considered by Guzzella et al., 2013 [41], to account for the minimum flow of reactants necessary to keep the fuel cell from shutting down.

As already explained, the goal of a quasi-static model is to obtain the stack current and, consequently, the hydrogen consumption from the required load power that must coincide with the fuel cell net power.

$$I_{st}{V}_{st}-P_{aux}=P_{load}$$

(70)

By integrating Equations (69) and (43), the stack current can be obtained by solving the following nonlinear equation:

$$P_{load}=N\left(E_0-A ·ln\left(i\right)-R_i·i-m·exp\left(ni\right)\right)I_{st}-k_{aux1}{I}_{st}-k_{aux1}{I}_{st}^2.$$

(71)

4.7. DC/DC Converter

A DC/DC converter is needed to adjust the variable voltage of the fuel cell to the load voltage. A detailed procedure for the design of the DC/DC converter is reported in Lee et al. [112], while an efficiency map is used to model the converter in [72] and to calculate the net power in Equation (3).

5. Dynamic Models

A comprehensive dynamic model simulating all transient phenomena taking place in the entire fuel cell system is not available in the literature. This is justified by the presence of different time constants which make it possible to neglect or not some dynamic effects depending on the focus of the control action.

5.1. Double-Layer Charging Effect

Under dynamic operation, the voltage-current curve is not constant but shows hysteresis due to the double-layer effect, i.e., the formation of two charged layers of opposite polarity across the boundary between the porous cathode and the membrane that acts as a capacitor [56].

To model the double-layer effect, the scheme of Figure 10 is applied in Larminee et al. [6], Donateo et al. [105], Kandidayeni et al. 2019 [91], Qingshan et al. [47], and Wang et al. [35], with small differences in the position of the resistors associated with activation, ohmic, and concentration losses. All these schemes include an equivalent capacitor $V_{st}=N(E_{Nernst}-∆V_{Ohm}-∆V_{act,1}-∆V_C)$, where $∆V_C=\left(I_{st}-C\frac{dV}{dT}\right)\left(R_{conc}+R_{act2}\right)$.

A similar scheme, but with $∆V_{act,1}=0,$ is proposed in [92]:

$$∆V_C=\left(I_{st}-C\frac{d{V}_{st}}{dt}\right)\left(R_{act2}+R_{con}\right)$$

(72)

$$\frac{d{V}_C}{dt}=\zeta (V_{act}-V_c).$$

(73)

In other schemes, the capacitor is put in parallel with activation and concentration losses [98].

A value of $C=6.107$ F is adopted in [91], but figures between 0.1 F and 15 F are suggested in the scientific literature [38]. In [99], a specific capacitance of $0.01974 \mathrm{F}/\mathrm{c}{\mathrm{m}}^2$ is considered. In [98]:

$$\frac{d{V}_C}{dt}=\frac{I_{st}}{C}-\frac{∆V_{act}+∆V_{conc}}{\tau }$$

(74)

where the time constant $\tau$ can be expressed as a function of the activation overvoltage and the concentration overvoltage, $\tau =C\frac{∆V_{act}+∆V_{conc}}{I_{st}}$.

Less common is the adoption of the Randles circuit [8] that includes a frequency-dependent resistor and capacitor in parallel to a conventional capacitor. A further resistance $R_s$ is considered in series. This model is used to interpret the impedance spectra of a fuel cell [8].

5.2. Dynamics of Compressor Speed

The dynamic of angular speed ${\omega }_{cp}$ of the compressor is expressed as [51]:

$$J_{cp}\frac{d{\omega }_{cp}}{dt}={\tau }_{cp}-{\tau }_{cm}$$

(75)

where $J_{cp}$ is the combined rotational inertia of the compressor and the motor [76], and ${\tau }_{cp}$ and ${\tau }_{cm}$ are the loading torque and the driving torque, respectively.

The compressor is usually driven by a BLDC motor that can be modeled through three main variables: the constant $k_a$ that correlates the back electromotive force to the motor speed, and the constant $k_t$ that links the motor electric torque to the current, and the internal resistance $R_{cm}$ [76,113].

In a DC motor, the relation between the armature voltage $V_{cm}$, the armature current $i_a$, and the shaft angular speed, ${\omega }_{cp},$ is expressed by the following two equations:

$$V_{cm}\left(t\right)=L\frac{d{i}_a\left(t\right)}{dt}+R_{cm}{i}_a\left(t\right)+k_a{\omega }_{cp}(t)$$

(76)

$${\tau }_{id}=k_t{i}_a\left(r\right)$$

(77)

where ${\tau }_{id}$ is the electrical ideal torque. A loss of torque is determined by friction and pneumatic losses in the electric machine, ${\tau }_l$, that can be expressed as a function of the compressor speed ${\omega }_{cp}$ and output pressure $p_{cp}$:

$${\tau }_l=A_{00}+A_{10}{\omega }_{cp}+A_{20}{\omega }_{cp}^2+A_{01}{p}_{cp}+A_{11}{\omega }_{cp}{p}_{cp}+A_{02}{p}_{cp}^2$$

(78)

where $A_{00}$, $A_{10}$, $A_{20}$, $A_{01}$, $A_{11}$, and $A_{02}$ are empirical variables to be tuned.

The losses can also be quantified by introducing the efficiency of the electric machine, ${\eta }_{cm}=\frac{{\tau }_{cm}}{{\tau }_e}.$ Therefore, the torque of the motor can be written as [76]:

$${\tau }_{cm}={\eta }_{cm}\frac{k_t}{R_{cm}}\left(V_{cm}-k_a{\omega }_{cp}\right)$$

(79)

where $V_{cm}$ is the control voltage of the compressor.

The loading torque can be expressed as:

$${\tau }_{cp}=\frac{c_{p,air}{T}_{amb}}{{\omega }_{cp}{\eta }_{cp}}\left[{\left(\frac{p_{sm}}{P_{amb} }\right)}^{\frac{\gamma -1}{\gamma }}-1 \right]{\dot{m}}_{cp}$$

(80)

where $p_{sm}$ is the pressure in the supply manifold that depends on the mass balance of the supply manifold.

The flow rate of the air from the compressor, ${\dot{m}}_{cp}$, is calculated from the balance of the cathode sub-system, as explained in the next section. The compressor mass flow rate and outlet temperatures are calculated by using steady-state performance maps like that of Figure 9 in the form of lookup table. As pointed out in [114], there is a deviation between the dynamic and steady-state maps of a dynamic compressor during rapid transients. This aspect is generally neglected in the scientific literature.

Instead of supplying the compressor flow maps in the form of a lookup table that is not well-suited for the simulation of dynamic systems, standard interpolation techniques can be adopted [71]. In [51], the following correlation is considered:

$${\dot{m}}_{cp}=-5.395+0.116{\omega }_{cp}-639.833\frac{p_{sm}}{p_{amb}}+174.71{\left(\frac{p_{sm}}{p_{amb}}\right)}^2+0.173\frac{p_{sm}}{p_{amb}}{\omega }_{cp}.$$

(81)

By controlling the voltage command to the compressor motor, $V_{cm}$, it is possible to control the compressor speed and, consequently, the OER [76].

5.3. Dynamics of Pressure in the Cathode Sub-System

The assumption made in the scientific literature is that all cathodes (anodes) in the stack are lumped as one stack cathode (anode) volume [79]. The models available in the literature are all based on the mass conservation equation but differ in the number of lumped volumes considered in the analysis and, consequently, on the number of pressure level introduced. The most common adoption is to consider a single pressure level for the cathode system ${ p}_{ca}$ that does not change from inlet to outlet but is a function of the mass flow rate ${\dot{m}}_{air,in}$ [41]. In other more recent works, up to three levels are considered; the pressure in the supply manifold, $p_{sm}$, the pressure at the cathode, $p_{ca},$ and the pressure in the return manifold, $p_{rm}$ [51,76].

The supply manifold volume includes the volume of the pipes between the compressor and the fuel cell stack, including the cooler and the humidifier, if present [71]. The pressure dynamics in the manifold are expressed by the following mass conservation equation:

$${\dot{p }}_{cp}=\frac{R{T}_{sm}}{M_{air}{V}_{sm}}({\dot{m}}_{cp}-{\dot{m}}_{air,in})$$

(82)

where $T_{sm}$ is the temperature in the supply manifold, and $V_{sm}$ is the volume of the supply manifold. Note that the temperature in the supply manifold can be considered equal to the stack temperature [41].

Since the supply manifold is directly connected to the stack, the pressure drop between them is very small, and the mass flow rate can be assumed to have a linear dependence on the pressure drop:

$${\dot{m}}_{air,in}=k_{ca,in}\left(p_{sm}-p_{ca}\right)$$

(83)

where $k_{ca,in}$ is the flow constant, and $p_{ca}$ is the cathode pressure. This equation is very important in the control of the system because it relates the OER to the state variable $p_{sm}$ and $p_{ca}$.

The pressure dynamics in the cathode are dependent on the air flow into and out of the cathode and the oxygen consumed in the reaction (see Equation (90)).

The return manifold represents the pipeline at the fuel cell stack exhaust. The mass flow out of the cathode is calculated using the nozzle flow equations for compressible flows and depends on the pressure in the return manifold $p_{rm}$:

$$\begin{gathered} {\dot{m}}_{rm,out}=C_D{A}_T{p}_{rm}\sqrt{\frac{\gamma }{R{T}_{st}}}{\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma +1}{2\left(\gamma -1\right)}}, for \frac{p_{amb}}{p_{rm}} \le {\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma }{\gamma -1}} (chocked flow) \\ {\dot{m}}_{rm,out}=C_D{A}_T{p}_{rm}\sqrt{\frac{\gamma }{R{T}_{st}}}\sqrt{\frac{2}{\gamma -1}\left[{\left(\frac{p_{amb}}{p_{rm}}\right)}^{\frac{2}{\gamma }}-{\left(\frac{p_{amb}}{p_{rm}}\right)}^{\frac{\gamma +1}{\gamma }}\right]}, for \frac{p_{amb}}{p_{rm}}>{\left(\frac{2}{\gamma +1}\right)}^{\frac{\gamma }{\gamma -1}} \end{gathered}$$

(84)

where $C_D$ is the discharge coefficient, and $A_T$ is the opening area of the nozzle. In [51], only the chocked condition is considered, while [29] presents only the formula for non-choked flows.

The partial pressure of the reactants during transients is not constant [98]. During transient operation, the mass flow rate of reactants is usually greater than that consumed in the reaction to avoid starvation or short circuit. Therefore, a reactant flow model is needed to determine the effective partial pressure of hydrogen and oxygen at the anode and the cathode [98].

The oxygen partial pressure at the cathode side can be estimated by summing the partial pressures of oxygen, nitrogen, and water vapor [76]. The dynamics of the partial pressure of oxygen are given by [68,96]:

$$\frac{d}{dt}{m}_{O_2}=\frac{V_c}{RT}\frac{d{p_O}_2}{dt}={\dot{m}}_{O_2,in}-{\dot{m}}_{O_2,react}-{\dot{m}}_{O_2,out}$$

(85)

where ${\dot{m}}_{O_2,in}$ is the mass flow rate provided by the compressor, ${\dot{m}}_{O_2,react}$, is the oxygen consumed by the reaction (see Equation (50)), and $V_c$ is the cathode volume.

The oxygen flow rate leaving the stack is ${\dot{m}}_{O_2,out}=k_c{m}_{O_2}$ according to [68], and detailed as ${\dot{m}}_{O_2,out}=k_c\frac{m_{O_2}}{m_{O_2}+m_{N_2}}\left(p_{ca}-p_{amb}\right)M_{O_2}$, according to [96]. In [76], it is expressed as ${\dot{m}}_{O_2,out}=k_c\left(p_{ca}-p_{rm}\right)$, where $p_{ca}$ is the cathode pressure, $M_{O_2}$ is the oxygen molar mass, and $k_c$ represents the effects of the outlet conditions.

A similar mass balance equation can be written for nitrogen that does not take part in the reactions in the fuel cell:

$$\frac{d}{dt}{m}_{N_2}={\dot{m}}_{N_2,in}-{\dot{m}}_{N_2,out}={{\dot{m}}_{N_2,in}-k}_c\frac{m_{N_2}}{m_{O_2}+m_{N_2}}\left(p_{ca}-p_{amb}\right)M_{N_2}.$$

(86)

The dynamics of the water formed at the cathode and passed through the layers to reach the anode is expressed by [68]:

$$\frac{d}{dt}{m}_{H_{2}O}=\frac{m_{H_{2}O}}{p_{ca}{k}_{W1}}{m}_c+k_{w2}{I}_{st}+\left(k_{w3}{I}_{st}-k_{w4}\right)k_{w5}-\frac{m_{H_{2}O}}{I_{st}{k}_{W1}{k}_{W6}}{m}_{c,out}{k}_{w7}$$

(87)

where $k_{w1}$, …, $k_{w7}$ are empirical variables, $m_c$ is the mass at the cathode inlet, and $m_{c,out}$ is the mass at the cathode outlet.

By summing up the contributions of oxygen, nitrogen, and humidity at the cathode, it is possible to write [79]:

$$\frac{d{P}_{ca}}{dt}=\frac{R{T}_s}{V_c}\left(\frac{{\dot{m}}_{air,in}}{M_{air,in}}-\frac{{\dot{m}}_{air,out}}{M_{air,out}}-\frac{{\dot{m}}_{O_2,react}}{M_{O_2}}+\frac{{\dot{m}}_{H_{2O}gen}-{\dot{m}}_{H_{2}O,mem}}{M_{H_{2}O}} \right)$$

(88)

where ${\dot{m}}_{H_{2}O,mem}$ is the mass flow rate of water vapor across the electrolyte membrane:

$${\dot{m}}_{H_{2}O,mem}=2M_{H_{2}O}{S}_{act}\left(n_d\frac{i}{F}-D_w\frac{{\phi }_c-{\phi }_a}{l}\right)$$

(89)

with $n_d$ as the electro-osmotic drag coefficient (number of water molecules dragged with each proton of hydrogen transported from the anode to the cathode), $D_w$ as the diffusion coefficient, ${\phi }_c$ and ${\phi }_a$ are the relative humidity at the cathode and the anode, respectively. $D_w$ and $n_d$ are functions of the membrane water content (Equation (30)) that varies between 0 and 14:

$$n_d=0.0029{\lambda }_m^2+0.05{\lambda }_m-3.4\times {10}^{-19}$$

(90)

5.4. Dynamics of Pressure at the Anode

The mass balance at the supply manifold of the anode is given by:

$$\frac{d}{dt}{m}_{H_2}=\frac{V_{an}}{RT}\dot{\frac{d{p}_{H_2}}{dt}}=m_{H_2,in}-{\dot{m}}_{H_2,react}-{\dot{m}}_{H_2,out}$$

(91)

where ${\dot{m}}_{H_2,react}$ can be obtained from Equation (61), while $V_{an}$ is the anode volume. The inlet hydrogen flow can be controlled by means of an electromagnetic valve [78]. The flow rate of hydrogen is controlled by the feedback pressure difference between the cathode and anode in [115]. Using a control valve with constant $k_v$, it is possible to write [68]:

$${\dot{m}}_{H_2,in}=\left(1-\frac{p_{v,a}}{p_a}\right)k_v{I}_{st}$$

(92)

where $p_{v,a}$ and $p_a$ are the partial pressure of the water vapor and the pressure at the anode, respectively. In [96,98], the outflow of hydrogen is expressed as:

$${\dot{m}}_{H_2,out}=k_a\left(P_{H_2}-P_{amb}\right)M_{H_2}$$

(93)

while in [68], it is simply ${\dot{m}}_{H_2,out}=k_a{m}_{H_2}$.

The dynamic of water content in a self-humidified stack is expressed in [31] as a function of the liquid water saturation, $s_{st}$, defined as the ratio of the liquid volume to the total volume of void space in the porous structure:

$$K_s\frac{d{s}_{st}}{dt}={\dot{m}}_{H_2,react}-{\dot{m}}_{H_2,evap}-{\dot{m}}_{H_2,diff}$$

(94)

where ${\dot{m}}_{H_2,evap}$ and ${\dot{m}}_{H_2,diff}$ are the loss of water due to evaporation and diffusion, respectively. The evaporation rate and the water production rate can be of the same order of magnitude:

$${\dot{m}}_{H_2,evap}=\frac{K_{evap}{M}_{H_{2}O}}{R{·T}_{stack{D}_{pore}}}\left(p_{sat}-p_v\right) if p_v<p_{sat}$$

(95)

where $p_v$ is the vapor pressure, and $p_{sat}$ is the saturated vapor pressure, which is a function of the temperature:

$$p_{sat}=p_0\exp \left(-\frac{E_a}{K_B{T}_{stact}}\right)$$

(96)

In this equation, $p_0$ is the pre-exponential factor for a temperature range of 0–100 K, and $E_a$ is the activation energy of evaporation. Since vapor pressure does not usually exceed the saturated vapor pressure, there is no condensation in the stack. Moreover, the changing of local vapor pressure because of liquid water evaporation and diffusion, is negligible [31].

The diffusive liquid water flux is calculated by Strahl et al. [31] by combining capillary theory with Darcy’s law. The overall model includes three empirical variables to be tuned in comparison with experimental data. The authors of [31] pointed out the uncertainties related to the initial state of humidification that is always unknown after shutdown and restart.

In [35,38], the dynamics of reactant flows are accounted for by using a delay in the Nernst equation:

$$E_{Nernst}=E_0+\frac{R{T}_{st}}{2F}\left[\ln \left(p_{H_2}\right)+\frac{1}{2}\ln \left(p_{O_2}\right)\right]-k_E\left(T-298.15\right)-{\lambda }_e{I}_{st}(s)\frac{{\tau }_{e}s}{{\tau }_{e}s+1}$$

(97)

where $k_E$ is an empirical variable, ${\lambda }_e$ is a constant (expressed in Ohm), and ${\tau }_e$ is s the constant time delay due to fuel and oxidant delays during load transients. For the PEM SR-12 FC, the values suggested in [35] for ${\tau }_e$ and ${\lambda }_e$ are 80 s and $3.33 \mathrm{m}\mathsf{\Omega }$, respectively, while their values were calibrated in [38] for a 1.2 kW Nexa PEMFC using different fitting functions. The following values were suggested: ${\tau }_e$ = 115.3 s and ${\lambda }_e=4.48 \mathrm{m}\mathsf{\Omega }$.

5.5. Humidification System

The degree of humidification is regulated by adjusting the water temperature within the humidifier in schemes like that of Figure 2. To model this process, two steps are considered in [113]. In step one, the dynamics of this subsystem are assumed to be dominated by the pressure change in the air humidifier, without considering the effect of the vapor injected into the gas. In the second step, the effect of injected vapor is taken into consideration. This subsystem is considered not relevant in this investigation because, as already explained, self-humidified stacks (Figure 3) are preferred in aerospace applications [116].

5.6. Thermal Dynamics and Coolant Pump

Unlike an open-cathode fuel cell, in CCPEMFC, the control of the temperature is decoupled from the control of the OER. Therefore, it is often neglected, and the stack temperature is assumed constant in most investigations on the control of CCPEMFC (see, for example, [76]).

The thermal dynamic process can be modeled by using the following equation [50]:

$$m_{FC}{C}_{FC}\frac{d{T}_{st}}{dt}=\left({\dot{\mathrm{Q}}}_{gen}-{\dot{\mathrm{Q}}}_{nat}-{\dot{\mathrm{Q}}}_{rad}-{\dot{\mathrm{Q}}}_{cool}-{\dot{\mathrm{Q}}}_{sens+lat}\right)+\sum _i{\dot{m}}_{i,in}{c}_{p,i,in}{T}_{i,in}-{\dot{m}}_{i,out}{c}_{p_i,out}{T}_{i,out}$$

(98)

where $m_{FC}$ and $C_{FC}$ are the stack mass and thermal capacity, respectively.

The forced convection from the fuel cell body to the liquid cooling system can be calculated as [50]:

$${\dot{\mathrm{Q}}}_{cool}=h_{forc}{S}_{forc}(T_{st}-T_{cool,0})$$

(99)

where $T_{cool,0}$ is the inlet temperature of the cooling fluid.

The change in the coolant temperature in the coolant channel, $i$, is represented by:

$$m_{cool}{c}_{cool}\frac{d{T}_{cool}}{dt}={\dot{\mathrm{Q}}}_{cool}+\sum {\left({\dot{m}}_i{h}_i\right)}_{cool}$$

(100)

The inlet temperature $T_{cool,0}$ can be calculated from:

$$m_{cool}{c}_{cool}\frac{d{T}_{cool,0}}{dt}=h_{forc}{S}_{forc}\left(T_{st}-T_{cool,0}\right)+{\dot{m}}_{cool}{c}_{cool}(T_{cool,in}-T_{cool,out})$$

(101)

where $c_{cool}$ is the specific heat of the coolant (generally water), and ${\dot{m}}_{cool}$ is its flow rate. The inlet temperature $T_{cool,in }$ is obtained from the outlet of the reservoir tank, while the outlet temperature $T_{cool,out }$ can be calculated using the logarithmic mean temperature difference [50]:

$$∆T=\frac{\left(T_{cool,out}-T_{cool,0}\right)-(T_{cool,in}-T_{cool,0})}{\ln \frac{T_{cool,0}-T_{cool,0}}{T_{cool,0}-T_{cool,0}}}$$

(102)

The heat exchanger, the coolant tank, and the pump are also modeled in [50]. The required mass flow rate of the pump is computed from Equation (45). The flow rate ${\dot{m}}_{cool}$ can be expressed as a linear function of the coolant pump speed:

$${\dot{m}}_{cool}=k_{bl}{\omega }_{pump}.$$

(103)

The dynamics of the pump can be written as:

$$J_{pump}\frac{d{\omega }_{pump}}{dt}=k_m-k_p{\omega }_{pump}^2.$$

(104)

The values of $k_{bl}$, $k_m$, $k_p$, and ${\eta }_{pump}$ can be obtained from the pump manufacturer [50].

When the stack is operated at an ambient temperature below 0 °C, the water generated by the reaction could freeze and destroy the structure of the MEA because of the increase in water volume. To avoid such a problem, the heating rate of the stack must be faster than the freezing rate of water. This can be achieved by different methods like shutdown purging, external auxiliary heating, hydrogen pre-heating, etc. [49,64,117,118]. The start-up behavior of PEM fuel cells below 0 °C requires very accurate control of the stoichiometric flow ratio, and possibly an additional electric heater can be included in the bypass circuit of the cooling system [119]. Repeated cold-start operations degrade the stack. as shown in [120], where the rated point performance of the CCPEMFC was found to be reduced by about 1% after 30 cold start-up tests performed in an environmental chamber.

6. Modeling Effects of Altitude

It is well known that ambient temperature, pressure, and relative humidity vary with elevation. According to the International Standard Atmosphere model (ISA), the following equations can be used to assess the ambient temperature and pressure at altitude $h_L$ in meters from sea level to 11,000 m:

$$T_{amb}\left(h_L\right)=T_0-6.5 h_L$$

(105)

$$p_{amb}\left(h_L\right)=p_0{\left(1-6.5\frac{h_L}{T_0}\right)}^{5.2561}$$

(106)

where $T_0=288.15 \mathrm{K}$ and $p_0=1.013 \mathrm{b}\mathrm{a}\mathrm{r}$.

From 110,000 m to 250,000 m, the pressure at altitude $h_L$ is expressed as [88]:

$$p_{amb}\left(h_L\right)=p_0\exp \left(-g\frac{h_L-\mathrm{11,000}}{RT}\right).$$

(107)

The humidity of the air from 0 to 15 km, expressed in ppmv (parts per million by volume), can be fitted by [88]:

$$H_{2}O\left(ppmv\right)=\frac{\mathrm{16,176}}{1+\exp \left(\frac{h_L}{2299}\right)}+20.51\frac{h_L}{1000}-320$$

(108)

In [111], empirical equations in $h_L$ are adopted to fit the variation of atmospheric pressure (in bar) and density with altitude:

$$p_{amb}\left(h_L\right)=7.782·{10}^{-9}{h}_L^2-0.0001508h_L+1.1$$

(109)

$${\rho }_{amb}\left(h_L\right)=6.332·{10}^{-9}{h}_L^2-1.46·{10}^{-4}{h}_L+1.289$$

(110)

A review of PEMS fuel cells utilized in high-altitude environments is reported in [64], where the influence of vibration, low pressure, cathode air starvation, and cold start is analyzed. In a closed-cathode PEMFC with a rated power of 1.2 kW tested at an altitude of 1524 m, a drop in the power of 25% compared to sea level was found. This is coherent with the information reported in the fuel cell manual [24], which estimates a power drop of 15 W for every 100 m increase in altitude and 10 watts for every degree Celsius above 30 °C, up to a maximum of 40 °C [58]. Other works on the performance drop of closed-cathode liquid-cooled fuel cells are described in [64].

A high cell temperature and low atmospheric values of relative humidity worsen the performance of the stack because of the decreased proton conductivity and increased ohmic losses. By elaborating on the voltage-current curves of [88], the author of this paper estimated a loss of 20% and 30% of the single cell power at 3000 m and 8000 m, respectively, when the working temperature of the cell was 353.15 K (80 °C). The loss of power was significantly reduced when the fuel cell works at the temperature of 323.15 K (50 °C). As pointed out above, the investigation of Gonzàlez-Espasadìn et al. [88] refers to the single cell, so the model is not suitable to analyze the behavior of the whole fuel cell system. A summary of the information available in the scientific literature about the reduction of the maximum power of the fuel cell at altitudes up to 8000 m is reported in Table 2.

The effect of the operating conditions in terms of temperature, pressure, and water content on the stack voltage ca be expressed as [87]:

$$V_{cell}=E_{Nernst}\left(p_{amb}, T_{stack}\right)-\mathsf{\Delta }{V}_{act}\left(p_{amb}, T_{stack}\right)-\mathsf{\Delta }{V}_{Ohm}\left(T_{stack}, {\lambda }_m\right)-\mathsf{\Delta }{V}_{conc}(p_{amb}, T_{stack}).$$

(111)

However, according to the equations reported in Section 4, the correlation between the stack behavior and the ambient conditions is more complex. As pointed out in [6], the internal resistance is also dependent not only on the temperature and water content but also on the ambient pressure, but this effect is not accounted for in the models available in the scientific literature. The temperature of the $T_{stack}$ under stationary conditions is set by the control unit of the cooling circuit while in transient operation is affected by the ambient temperature.

To account for the effect of ambient conditions on a PEM fuel cell, Barelli et al. [121] developed a semi-empirical model by using existing functions of AspenPlus and a user-defined Fortran subroutine. The effect of changes in pressure, RH, temperature, and humidity is modeled as a deviation from an experimentally fitted polarization curve using literature data for calibration.

Regarding the balance of plant, the variation of elevation can lead to the degradation of the stack performances and create challenges for thermal and OER management [98]. As pointed out in [111], under low-pressure operation, the performance and efficiency of the stack are reduced also because of the effect of altitude on the air compressor map [76]. The optimal values of OER are significantly affected by altitude, as pointed out in [111], where the authors performed off-line quasi-static experiments under altitudes 0–4000 m. They found that the maximum OER that can be achieved is also reduced, and the optimal OER curves move towards lower values. In fact, the increase in altitude, for the same current, results in a reduction of the net power of the CCPEMFC. The operating range of the compressor is also reduced because the slope of the system characteristic curve increases [122], and the operating point of the compressor at low-flow conditions can reach the surge zone [111]. This shows how critical the control of OER at high altitudes is. The effect of altitude on the compressor map is shown in Figure 11.

Table 2. Review of studies on the loss of performance at high altitude.

Fuel Cell	Type of Cooling	Max Altitude	Power Drop	Ref.	Type of Analysis
Nexa Ballard 1.2 kW	Closed-cathode air cooled	1525 m	22%	[58,64]	Experimental (altitude chamber)
Nexa Ballard 1.2 kW	Closed-cathode air cooled	2256 m	24%	[123]	Experimental
Open-cathode single cell	No cooling considered (cell temperature experimentally set)	3239 m	17% (from 2.21 to 1.83)	[55]	Experimental (altitude chamber)
Generic fuel cell stack with 50 cells	Open-cathode air cooled	4000 m	8%	[98]	Numerical
Single-cell	No cooling considered (cell temperature arbitrarily set)	8000 m	30%	[88]	Numerical

Table 2. Review of studies on the loss of performance at high altitude.

Fuel Cell	Type of Cooling	Max Altitude	Power Drop	Ref.	Type of Analysis
Nexa Ballard 1.2 kW	Closed-cathode air cooled	1525 m	22%	[58,64]	Experimental (altitude chamber)
Nexa Ballard 1.2 kW	Closed-cathode air cooled	2256 m	24%	[123]	Experimental
Open-cathode single cell	No cooling considered (cell temperature experimentally set)	3239 m	17% (from 2.21 to 1.83)	[55]	Experimental (altitude chamber)
Generic fuel cell stack with 50 cells	Open-cathode air cooled	4000 m	8%	[98]	Numerical
Single-cell	No cooling considered (cell temperature arbitrarily set)	8000 m	30%	[88]	Numerical

In [29,111], some modeling approaches are proposed to account for the effect of altitude on the air supply system of a CCPEMFC. In these investigations, the temperature of the stack and the gas flowing into the anode are considered constant, uniformly distributed, and equal to 80 °C. However, the effect of different ambient conditions can be accounted for by introducing the corrected mass flow and the corrected revolution speed, and expressing the compressor performance maps as a function of these variables [52]:

$${\dot{m}}_{cp}\left(h_L\right)={\dot{m}}_{cp}\left(s.l.\right)·\frac{p_{amb}\left(h_L\right)}{p_0}\sqrt{\frac{T_0}{T_{amb}\left(h_L\right)}}$$

(112)

$${\omega }_{cp}\left(h_L\right)=\omega \left(s.l.\right)·\sqrt{\frac{T_{amb}\left(h_L\right)}{T_0}}.$$

(113)

The authors of [52] verified the validity of this procedure experimentally and obtained a map of pressure rate vs. corrected air mass flow for altitudes up to 4000 m. The experiments were performed in a custom-made climatic chamber. The same authors also pointed out that in some applications, the operation map of the compressor is limited by the maximum power that the motor can provide. Therefore, the performance map of the electric machine is also relevant when analyzing the behavior of a fuel cell system at high altitude.

In [111], Li et al. detailed the friction losses and the shock losses on the impeller and diffuser. The losses due to reverse flow, volute, and diffuser are also considered in the calculation of the compressor efficiency. The model is validated at zero level and at 3000 m. The effect of ambient conditions is also to be considered in the dynamics of the cathode subsystem by replacing sea-level ambient conditions with pressure and temperature calculated at altitude $h_L$ in the relevant Equations (82), (86), and (88).

Li et al. [111] also model the minimum flow to avoid surges at different altitudes by fitting experimental data:

$${\dot{m}}_{air,surge}=4.594·{10}^{-10}{h}_L^2+1.25·{10}^{-6}{h}_L+0.01494.$$

(114)

This value is used to correct the value of the OER by ensuring that the system is far from the surge zone:

$$OER>\frac{0.92\left({\dot{m}}_{air,surge}+∆\right){\dot{m}}_{air}}{N I_{st}{M}_{O_2}}$$

(115)

where $∆$ is the safety margin.

7. Conclusions and Future Directions

The state-of-the-art in semi-empirical models for closed-cathode fuel cell systems was analyzed to identify the best simulation approach under steady-state and dynamic operating conditions typical of low-power aerial vehicles. The analysis was performed at the system level, considering not only the stack (with its voltage-current characteristic) but also the balance of plant, particularly the air-supply system.

The analysis of the quasi-static stack models at the stack level revealed a consolidated approach in modeling the different polarization losses that affect the fuel cell voltage, at least for applications at sea-level operation. These models are nominally suitable to account for the effect of variable ambient conditions characteristic of aerospace applications. However, some inconsistencies were found in the tuning and validation procedures. The models are usually fitted over a single voltage-current curve and therefore, are expected to have limited prediction capability.

The modeling of the auxiliary systems, while relevant for the determination of actual fuel consumption and net power of the systems, is a topic not sufficiently addressed. In fact, the evaluation of the net power and efficiency of the fuel cell system is seldom performed at a numerical level, and a standardized quasi-static simulation framework for the balance of the plant at variable-altitude applications has not been developed yet. Often, the models are validated at sea level, and then extended to simulate different altitudes, while altitude chambers are rarely used to obtain an experimental evaluation of the fuel cell behavior at high altitudes.

About dynamic modeling, a comprehensive dynamic model simulating all transient phenomena taking place in the whole fuel cell system is not available in the literature. Works concentrate, in fact, on a single subsystem, particularly on the cathode circuit. The limitations found in the quasi-static and dynamic models prove the necessity of developing a standard validation procedure under different operating conditions using altitude chambers or on-board measurements for fuel cell systems adopted in aerospace applications.

As far as control issues are concerned, recent research interest is in online identification methods and advanced control strategies to overcome the limits of standard PID controllers. Moreover, the oxygen excess ratio of the fuel cell is no longer considered a constant value but a parameter to be optimized according to operating conditions. Also, in this case, research on the effect of altitude is still quite limited, and numerical and experimental investigations are still needed to satisfactorily address this important issue for the aerospace applications of fuel cells.

According to the author, one of the difficulties of developing an overall simulation framework for the stack and its balance of plant could be the multi-physics nature of such a system that includes electrochemical devices, compressible and incompressible fluid machinery, electric machines, heat exchangers, electronic devices, etc., and therefore, requires multi-disciplinary knowledge. Therefore, as future directions, the author suggests the development of a comprehensive simulation approach for variable-altitude applications, the elaboration of accurate tuning procedures, and the acquisition of vast experimental data using altitude chambers or on-board acquisition systems.