Regression analysis of PEM fuel cell transient response

Problem considered

Polymer electrolyte membrane (PEM) fuel cells are electrochemical devices that produce electrical energy from the chemical energy present in hydrogen fuel, with water and heat being the byproducts. They have high efficiency at relatively low operating temperatures, and so have been investigated for transportation and electronic applications. Durability and cost issues are seen as the major barriers to commercialization. The major focus of PEM fuel cell cost reduction and performance improvement strategies is on issues of (1) heat and water management and (2) new materials development [ ¹ , ² ]. Heat and water management processes determine transient response characteristics of PEM fuel cell power systems, which need to produce power with varying load conditions. In practice, sophisticated models have been used, which combine a nonlinear static model for steady-state voltage with a linear dynamic subsystem to describe the voltage evolution, with time, between steady points [ ³ ]. Steady voltage might be estimated by a lookup table that references operating parameters. Transient voltage behavior, as the result of a current step between steady-state operating points, might be estimated with linear dynamic transfer functions.

Previous works

Transient response has been the subject of experimental and modeling efforts, which are both qualitative and quantitative in nature. Hamelin et al. [ ⁴ ] suggested that the hysteresis seen in swept load commutations was due to changing membrane ionic conductivity, which resulted from changes in the membrane’s hydration level. Yu and Zeigler [ ⁵ ] incorporated HFR (high-frequency resistance) measurements to infer changes in MEA hydration. Hou [ ⁶ ] follows a similar approach. He concluded that water redistribution in the PEM membrane was the slowest transient process in the operating PEMFC, but that heat transfer also plays an important role.

Several experimental investigations examined the response to step load changes in the PEMFC. Pathapati et al. [ ⁷ ] modeled temperature changes occurring during this transient period. Yan et al. [ ⁸ ] investigated voltage responses under many operating conditions. The work was very broad and mainly reported general trends without illuminating underlying causes of dynamic behaviors. Kim and Min [ ⁹ ] measured the PEMFC voltage response under fully humidified conditions at low temperatures (30–50 °C) and high current density. Liquid water accumulation (flooding) was observed, leading to restrictions in reactant gas transport and subsequent influence of the voltage response of the PEMFC. PEMFC transient response was surmised to be controlled by the diffusion processes associated with restricted gas transport. Subsequent experiments investigated the effects of stoichiometric ratios, humidity levels, and flooding intensity [ ¹⁰ ]. Later, they showed the effects of a degraded gas diffusion layer (GDL) on transient response, repeating many of their earlier conclusions [ ¹¹ ]. The multiphase model of Loo [ ¹² ] reproduced several key findings of these experiments. Takaichi et al. [ ¹³ ] measured transient redistribution of water content through the thickness of the membrane. They reported that the cathode side experiences a drop in resistivity/increase in water content when current is increased. The anode side experiences the converse effect. Reactant gas humidity suppressed these effects. Moçotéguy et al. [ ¹⁴ ] investigated changes in transient response of a PEMFC stack with different aging times. Cell resistivities and response times were not systematically influenced by aging. Kim and Shimpalee [ ¹⁵ , ¹⁶ ] investigated single-cell current transients with fully humidified reactants while operating at varying fuel stoichiometry. The measured current fluctuations reflect rapid gas-phase transients that produce anode reactant starvation. MEA hydration and thermal transients are not considered. Longer transients were seen when air penetrated the anode flow channels, deactivating the affected anode areas until hydrogen flow could increase.

Experimentally inspired empirical approaches have been developed that do not consider the underlying physics of the problem. Meiler and co-workers [ ³ , ¹⁷ ] modeled dynamic voltage response with “black box” models constructed from linear or nonlinear transfer functions from experimental data. The model structure was chosen only by suitability and the model parameters were found by least-squares fit to experimental data. They rejected theory-based “physical” models of dynamic response as being too complex and requiring far too much computing time [ ³ ]. Hussaini and Wang [ ¹⁸ ] empirically correlated maximum voltage undershoot to the size of change in HFR.

Some transient modeling approaches involved numerical solution of transport equations within the MEA, starting with models devised for the fully humidified case [ ¹⁹ ]. Later works incorporated changing MEA water content [ ²⁰ ], necessary to model low-humidity transients. The model was refined to incorporate an electron transport equation and could then simulate transient step changes in load current [ ²¹ ]. They showed how the transient response of the PEMFC, in low-humidity operation, was impacted by water storage in the Nafion membrane, and described various time scales of PEMFC response.

Present experimental study

This work investigates a decomposition of the transient response of a PEMFC to a step change in operating current. It follows a previous strategy, but here fit parameters are directly linked to underlying physical phenomena, which should allow improved parameter estimation (in the operation control application) and better understanding. Mechanisms and time constants of transient response can be correctly assessed without excessive speculation.

This work fits load step responses to a nearly equivalent model. It is intended for thin MEAs in low- to mid-humidity operation with adequate reactant stoichiometry. The “thin” membrane presents two key operational characteristics [ ²² ]: (1) membrane water concentration profiles become approximately linear (anode dryout due to electroosmotic drag does not occur) and (2) cell resistance drops with increasing load (from water production). This work is not intended for low-temperature and high-humidity response where liquid water flooding develops a post-current step, changing the characteristic shape of the voltage recovery.

Test setup and instrumentation

The testing equipment consisted of a Scribner Associates model 850e Fuel Cell Test System. Figure  ¹ shows a schematic diagram. This system controls load, reactant flow rate, humidification, and cell temperature. There is a PC-driven data acquisition system for recording cell temperature, anode and cathode flow rates, humidity levels, voltage, current, and resistance measurements conducted with the current-interrupt (CI) technique. All measurements are acquired at a sampling frequency of 10 Hz. This approach has better transient response than the alternative high-frequency resistance (HFR) method.

Fig. 1

Fuel cell

Experiments were conducted on a fuel cell with an effective area of 25 cm ² . The MEA used in this study consisted of a commercially available Nafion NR211 (PO#2430, Lot#1619) supplied by Ion Power, Inc. It is about 30 μm thick with 15 μm catalyst layers on either side. Gas diffusion layers were the commercially available SGL Carbon 10BC from the same supplier. A break-in procedure [ ²³ ] was utilized. Humidified hydrogen and air flowed in a 6-turn triple-pass serpentine flow field etched in graphite bipolar plates. The flow channels have dimensions of 1 mm × 0.75 mm × 30 mm. The backpressure was 34 kPa. The fuel cell hardware was supplied by Fuel Cell Technologies. Components of the fuel cell are shown in Fig.  ² .

Fig. 2

Operating conditions and test procedure

Tests were conducted to measure both the steady-state polarization curves and a number of transient load (current density) steps. Cell temperatures and relative humidity levels were varied. The operating conditions are summarized in Table ¹ . Note that a capital letter is used to describe the cell temperature ( A  = 40 °C, B  = 60 °C, C  = 80 °C) and a subsequent number is used to describe the humidity levels of the reactant gases (1 = 24 %, 2 = 48 %, 3 = 64 %, 4 = 85 %). For example, at 60 °C cell temperature and 24 % reactant gas humidity, measurement [B1] has anode and cathode saturator temperatures T _A/C , each maintained at 32 °C. Gas flow rates were controlled in terms of a predefined stoichiometric ratio.

Steady-state polarization curves were obtained at all of these operating conditions. Reactant flow rates were set to a constant stoichiometry of 1.25 at the anode and 2.5 at the cathode in all tests. Steady-state conditions were achieved by operating the cell at a given condition for a period of 15 min and observing constant cell voltage, temperature, current, and ohmic resistance.

Transient measurements were performed in which the fuel cell’s response to a step change in load current was observed. Flow stoichiometry values of 1.25 anode/2.5 cathode were set according to the maximum current. These had previously been found sufficient to avoid retarding the voltage recovery [ ⁸ , ¹⁵ , ²⁴ ]. First, a steady-state is reached by maintaining a constant current load for several minutes. Then, the control system imposed a step change in current on the cell, which was not a perfect step, but could be completed in about 0.2–0.3 s. The time response of the cell’s voltage and resistance was recorded for 190 s in each case.

Figure  ³ shows typical step events and the voltage responses. When the current increases, there is a voltage undershoot, which then recovers and reaches steady-state value during the recovery time. When the current decreases, an overshoot occurs, with subsequent recovery to an equilibrium value. Voltage transients can be expressed in terms of overall cell voltage and mean losses. The observed voltage V(t) depends on the time-varying open-circuit voltage VOC(t) , as well as activation losses ηact , mass-transport losses ηmt , and total ohmic resistance iR , where i is the current density (A cm ⁻² ) and R the area-specific resistance (mΩ cm ² ).

Fig. 3

Table  ² summarizes the five different current step increases and decreases utilized here. These varied in the size of the step and in the initial average current density. Lowercase letters indicate the step sequences. Step increases are shown in the left column end as ‘+’, and step decreases in the right end as ‘−’. The codes are combined to refer to specific measurements. For example, A3[c+] refers to the measurement at 40 °C, 64 % RH reactant gas feeds, with a current step from 0.1 to 0.4 A/cm ² .

Steady-state results

Steady-state polarization curves were measured at all of the temperature/humidity combinations in Table  ¹ . Ohmic losses significantly impact the voltage output of an operating PEMFC, and they can be measured in real time as a diagnostic tool. Equation ( ¹ ) describes the total ohmic resistance in a PEMFC. It results from the sum of electronic contact resistances and current flow resistances ( Re- ) in addition to protonic ( RH+ ) resistance losses. The term ( Re- ) is a given constant for a cell assembly and does not contribute to cell dynamic response. It is typically assumed that, in a well-built PEMFC, the measured ohmic resistance will be dominated by membrane ionic resistance ( RmH+ ). The ionic resistances of the protonically conducting phase of the catalyst layers RaclH+ and RcclH+ are not accounted for by the measurement [ ²⁵ – ²⁷ ]. All of the protonic resistances are hydration dependent [ ²⁸ ], decreasing with rising water content. Thus, membrane ohmic resistance has been used to infer the overall water content of the membrane during transients,

The ohmic resistance measurements explain the observed hydration effects in transient response and will be discussed in a later sub-section.

Figure  ⁴ shows the polarization curves and ohmic resistance measurements at 60 °C with varying humidity levels. Both are typical of low-humidity operation with a thin MEA, the intended application for this work. The low-humidity measurements show steeper drops in membrane resistance with increasing current density. Protonic conduction resistance RmH+ drops as current (water production) increases. Increases in membrane resistance due to anode dryout (by electroosmotic drag) are not observed. In high-humidity measurements, the moisture to hydrate the MEA can come from humidified reactant gases. Therefore, resistance values were observed to be less sensitive to current density. Observations at 80 °C are not shown here, but are similar.

Fig. 4

Modeling of resistance response

Measured resistance showed a brief spike ( t  ~0.3 s) followed by a first-order exponential decay characteristic, before reaching a new steady-state value. The resistance data could be modeled with a single amplitude and time constant, after omitting the first 0.3 s. High-humidity measurements showed small resistance changes, which took place immediately and were then constant.

The experimentally measured resistance vs. time curves could be fitted with a single-term exponential regression model:

The term R ( t ) represents the experimentally measured resistance. R∞ is the steady-state value which is reached at the end of the recorded data, and ΔRu represents the drop in resistance (current step-up) or increase (current step-down). τ _R represents the time constant of decay, in s.

Figure  ⁶ a, b shows the model fits for cases C1[a+] and C1[a−], respectively. The trust-region-reflective algorithm as implemented in the “lsqcurvefit” function of MATLAB was utilized.

Fig. 6

Figure  ⁷ shows predicted time constants as functions of the higher current density in the step, for three of the low-humidity operating conditions. For the step-up conditions, time constants were higher at lower temperatures and higher humidity levels. They also show an inverse logarithmic relationship to the final current value. For step-down conditions, the same trend was observed at 80 °C and 48 % relative humidity. However, other cases showed no discernible trends.

Fig. 7

Modeling the voltage response

The observed voltage showed a drop to minimum value, followed by an asymptotic recovery. The voltage response could be modeled with multiple exponential terms. Nonlinear regression techniques were employed to fit experimentally measured voltage data with n  = 1, 2, or 3 exponential terms as in Eq. (4): V ( t ) represents the experimentally measured voltage during recovery. V _∞ is the steady-state value. V _i , V _∞ , and τ _i are parameter fits. The V ₁ , V ₂ and V ₃ terms are amplitudes (positive is overshoot and negative is undershoot). Time constants τi are ordered as τ1<τ2<τ3 .

Figure  ⁸ shows one-, two-, and three-term regression fits for a representative low-humidity measurement: C2[c+]. Fits with only one or two exponential terms systematically deviate from the measured data. The three-term fit is needed to match the short time constant behavior at 2 s as well as the long time constant behavior between 40 and 80 s. The ‘delay’ behavior is seen: the voltage reaches a minimum of nearly 2 s after the current step is complete. The delay behavior necessitates that the first two amplitude terms V1>0 0$$\end{document}]]> and V2<0 should have opposite signs. The fitted parameters are presented in Table  ³ . Note that the second time constant have a similar value to that of the resistance, indicative of a connection in origin. The change in voltage resulting from a 33 mΩ cm ² shift in ohmic resistance (i.e., 0.033 Ω cm ² × 0.4 A/cm ²  = 13.2 mV) is approximately half of the 25 mV amplitude associated with the second component.

Fig. 8

Table 3

C2[c+] fit values:	Resistance	Voltage
		V1 = +23.5 mV	τ1 = 0.82 s
ΔRu = −33 mΩ cm2	τR = 5.6 s	V2 = −25.4 mV	τ2 = 7.0 s
		V3 = −13.2 mV	τ3 = 51.7 s

In addition to the qualitative picture, statistical tests can be used to determine the most correct model fit of the data. The coefficient of determination, commonly abbreviated as r ² , is used, though its application to nonlinear regression is admittedly problematic [ ²⁹ – ³¹ ]. Table  ⁴ shows uniformly high r ² values. The second test is Akaike’s information criterion [ ³² ]-corrected, abbreviated AIC _c, calculated from Eq. ( ⁵ ): where N is the number of time history points in the fitted portion of the curve, K is the number of model parameters plus 1, and SS represents the sum of the squares of the distances of the fit y -values from the measured y -values. Models are compared by looking at the difference in AIC _c between two models applied to the same dataset. The smallest value is most likely to be correct, and a difference of 6 indicates that the model with the lowest score has a 95 % chance of being correct. The difference of 14 is greater than this threshold and so the three term fits are most correct. This approach was used in the field of pharmacokinetics to determine the best number of exponential terms [ ³³ ].

AICC=NlnSSN+2K+2KK+1N-K-1.

Table 4

Goodness of fit	C2[c]		C2[f]
	r2	ΔAICc	r2	ΔAICc
1 Term	0.9747	2165	0.9831	691
2 Terms	0.9863	1436	0.9904	14
3 Terms	0.9959	0	0.9906	0