Regression analysis of PEM fuel cell transient response

To develop operating strategies in polymer electrolyte membrane (PEM) fuel cell-powered applications, precise computationally efficient models of the fuel cell stack voltage are required. Models are needed for all operating conditions, including transients. In this work, transient evolutions of voltage, in response to load changes, are modeled with a sum of three exponential decay functions. Amplitude factors are correlated to steady-state operating data (temperature, humidity, average current, resistance, and voltage). The obtained time constants reflect known processes of the membrane heat/water transport. These model parameters can form the basis for the prediction of voltage overshoot/undershoot used in computational-based control systems, used in real-time simulation. Furthermore, the results provide an empirical basis for the estimation of the magnitude of temporary voltage loss to be expected with sudden load changes, as well as a systematic method for the analysis of experimental data. Its applicability is currently limited to thin membranes with low to moderate humidity gases, and with adequately high reactant-gas stoichiometry.


Introduction 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 [1,2]. 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 [3]. 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. [4] 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 [5] incorporated HFR (high-frequency resistance) measurements to infer changes in MEA hydration. Hou [6] 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. [7] modeled temperature changes occurring during this transient period. Yan et al. [8] 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 [9] 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 [10]. Later, they showed the effects of a degraded gas diffusion layer (GDL) on transient response, repeating many of their earlier conclusions [11].
The multiphase model of Loo [12] reproduced several key findings of these experiments. Takaichi et al. [13] 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. [14] 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 [15,16] 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 [3,17] 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 leastsquares fit to experimental data. They rejected theory-based ''physical'' models of dynamic response as being too complex and requiring far too much computing time [3]. Hussaini and Wang [18] 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 [19]. Later works incorporated changing MEA water content [20], 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 [21]. 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 [22]: (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.

Experimental study Test setup and instrumentation
The testing equipment consisted of a Scribner Associates model 850e Fuel Cell Test System. Figure 1 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.

Fuel cell
Experiments were conducted on a fuel cell with an effective area of 25 cm 2 . 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 lm thick with 15 lm catalyst layers on either side. Gas diffusion layers were the commercially available SGL Carbon 10BC from the same supplier. A break-in procedure [23] 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 9 0.75 mm 9 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. 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 1. 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 [8,15,24]. 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 3 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 V OC ðtÞ, as well as activation losses g act , masstransport losses g mt , and total ohmic resistance iR, where i is the current density (A cm -2 ) and R the area-specific resistance (mX cm 2 ).
VðtÞ ¼ V OC ðtÞ À g act ðtÞ À g mt ðtÞ À iRðtÞ: ð1Þ Table 2 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 2 .

Steady-state results
Steady-state polarization curves were measured at all of the temperature/humidity combinations in Table 1. Ohmic losses significantly impact the voltage output of an operating PEMFC, and they can be measured in real time as a diagnostic tool. Equation (1) describes the total ohmic resistance in a PEMFC. It results from the sum of electronic contact resistances and current flow resistances (R eÀ ) in addition to protonic (R H þ ) resistance losses. The term (R eÀ ) 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 (R H þ m ). The ionic resistances of the protonically conducting phase of the catalyst layers R H þ acl and R H þ ccl are not accounted for by the measurement [25 -27]. All of the protonic resistances are hydration dependent [28], 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 4 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 lowhumidity measurements show steeper drops in membrane resistance with increasing current density. Protonic conduction resistance R H þ m 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.

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 1 is the steady-state value which is reached at the end of the recorded data, and DR u represents the drop in resistance (current step-up) or increase (current stepdown). s R represents the time constant of decay, in s.  Figure 7 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. T A/C = 27/27 [A3] T A/C = 32/32 [A4] T A/C = 37/37 T A/C = 32/32 [B2] T A/C = 43/43 [B3] T A/C = 51/51 [B4] T A/C = 56/56 T A/C = 48/48 [C2] T A/C = 60/60 [C3] T A/C = 70/70 [C4] T A/C = 76/76 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): ð4Þ V(t) represents the experimentally measured voltage during recovery. V ? is the steady-state value. V i , V ? , and s i are parameter fits. The V 1 , V 2 and V 3 terms are amplitudes (positive is overshoot and negative is undershoot). Time constants s i are ordered as s 1 \s 2 \s 3 . Figure 8 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 V 1 [ 0 and V 2 \0 should have opposite signs. The    Table 3. 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 mX cm 2 shift in ohmic resistance (i.e., 0.033 X cm 2 9 0.4 A/cm 2 = 13.2 mV) is approximately half of the 25 mV amplitude associated with the second component.
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 2 , is used, though its application to nonlinear regression is admittedly problematic [29][30][31]. Table 4 shows uniformly high r 2 values. The second test is Akaike's information criterion [32]-corrected, abbreviated AIC c, calculated from Eq. (5): 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 yvalues. 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 [33].
Discussion of voltage transient model The parameters in the three-term model can be related to changes in operating conditions, such as current step size and reactant gas humidity. Figure 9 shows the first parameters, V 1 and s 1 , plotted against the current step-up (I 2 À I 1 ).
Data for the largest current step at 85 % RH are excluded due to suspected flooding. The amplitude increases with both the size of the current step and the dryness of the reactant gases. V 1 is positive for the lower RH, serving to lessen the undershoot. It becomes small and negative, contributing to undershoot, with higher RH. The time constant s 1 is essentially scattered around 1 s. Many of the outliers are points where V 1 indicates small amplitude; hence the time constant cannot be estimated accurately.  The second parameter pair, V 2 and s 2 , exhibited a similar trend. Figure 10 shows a linear relationship between V 2 and I 2 DR u . The coefficient of determination is 0.91; in this case. V 2 increases with both the size of the current step and the dryness of the reactant gases. It is negative, representing undershoot in the voltage recovery. The low-humidity data (24 % RH) has time constant s 2 that follows the time constant of cell ohmic resistance s R . The high-humidity data did not show a resistance change; a measurement of s R was not possible. Amplitude V 2 was correspondingly very small. Low-humidity voltage response has typically been associated with resistance changes, and component V 2 has a time constant s 2 similar to the resistance measurement. Consider Eq. (2) which includes both membrane and catalyst layer protonic resistances.
The third parameter pair, V 3 and s 3 , could be correlated to changes of the rate of heat generation, per unit area, in the MEA. The amplitude V 3 is significantly greater with dry reactant inputs and varies monotonously with the final current step. The time constant s 3 is 30-150 s and rises with reactant gas humidity. Wu [34] suggested that the effective heat capacity of the membrane contains contributions from both the dry membrane and its water content. Figure 11 shows that V 3 increases with changing heat generation, D _ Q,between the initial and final states. The rate of MEA heat generation in either state is [35]: Problems of mathematical ill-conditioning are not thought to be responsible for the data scatter, because the third time constant was significantly removed from the first two in all measurements. The greatest outliers from the trend line are the exceedingly dry measurements of group C1, which have the largest V 3 amplitudes. These most violate the implicit assumption in Eq. (6) that there is a single sudden change in levels of heat generation pre-step and post-step. The automatic action of the unit's temperature control system, attempting to compensate, may have also added scatter to all readings. The best parabolic fits to    The typical resistance transient displays timescales s m;d and s m;h in its measurement of ðR eÀ þ R H þ m Þ. s m;d indicates the amount of time post-current step for water imbalances through the membrane thickness to equilibrate: back diffusion occurs, anode side water content levels normalize, and membrane conductivity improves. Wang et al. [21] proposed the equation:

Discussion of timescales
The hydration time constant s m;h reflects the time required for the MEA to gain or lose water content [20]. It was given at *5-20 s, in line with s R measured here.
Voltage transients can be examined with a 0-dimensional model equation: The cell voltage (V) is expressed in terms of the opencircuit voltage (V oc ), and various losses. The activation overpotential (g act ) drives both the anode and cathode reactions; both are commonly assumed to follow Butler-Vollmer kinetics. This term has been assumed to remain constant with current throughout the transient event of voltage recovery [18]. It will change with catalyst layer temperature, however [7]. Fuel cell inner-component temperatures have been observed to vary with first-order response characteristic due to varied step inputs, for example, in SOFC devices [38]. The third component V 3 , s 3 matches s T estimates at 40-120 s. Presumably, activation overpotential (g act ) changes due to temperature variations occurring at this time constant. Ionomeric-phase conductivity also has temperature dependence, potentially imparting a change in the resistance terms. Direct in situ temperature measurements have shown a linear temperature rise of the cathode catalyst layer (CCL), above the end-plate or current collector, of 0°C to 5-7°C, as current density increased from 0 to 1 A/cm 2 [39,40]. This 3°C temperature rise from a cell block temperature between (50-80°C) can be combined with the long-established conductivity relationship [41] to yield 3 % conductivity change. At a post-step current of 0.5 A/cm 2 with total ionomeric phase resistance of 0.090 X cm 2 , the change in voltage caused by this temperature-induced resistance change is only 1.3 mV. Hence, no significant s 3 component is seen in the ohmic resistance measurement, and s 3 is not indicated as affecting the catalyst layer resistances in Eq. (8).
The first component V 1 , s 1 also matches the reported value of s R H þ CL at *1 s [21]. The catalyst layer resistances ðR H þ acl þ R H þ ccl Þ vary with time and contribute to cell dynamic response, particularly in the low-humidity case. Temporary flooding of the cathode catalyst layer has been suggested as an explanation of the 1 s delay observed in voltage recovery. If flooding was causing this delay, they would become more severe, or of longer duration, with increasing  Figure 12 shows measurements at 80°C and four incrementally varying humidity levels for a 100-400 mA/cm 2 current step. The higher humidity measurements have been vertically shifted only for purposes of clarity. Delay is greatest where humidity is lowest, and delay is effectively suppressed by increasing humidity levels. Table 6 shows the first fit parameters and the observed delay in voltage response. The delay is greatest at low humidity, where amplitude V 1 is positive. At 85 % RH, amplitude V 1 becomes negative.

Conclusions
Dynamic responses of a typical single-cell PEMFC with a thin MEA to step changes in current load have been studied experimentally. Operating conditions were simplified to use excess stoichiometry and avoid water accumulation. Resistance and voltage transient responses are examined. These transient responses were well fitted by mono-exponential functions for the resistance and tri-exponential functions for the voltage. Previous studies have mostly used non-physical ''black box'' models to estimate dynamic response to step changes. The dynamic model's optimal fitted parameters would change somewhat with operating conditions. With the tri-exponential model, a strategy which incorporates variable-fitted parameters might offer improved results. This work showed how a very simplified understanding of the physics of the MEA can explain some of the variations in amplitudes and timescales. The first component varies with temperature, humidity and current step size, at the fixed stoichiometry employed here. The second component showed amplitude variations, which correlate to ohmic resistance changes, from membrane hydration occurring during the current step. The third component is consistent with variable MEA heat generation, impacting the temperature-dependent activation losses in the cathode.