Simulation of the Dynamic Characteristics of a PEMFC System in Fluctuating Operating Conditions

: A greater understanding of the dynamic processes inside the stack is urgently needed to optimize the PEMFC (proton exchange membrane fuel cell). In this study, we examined the gas, water and electrochemical processes inside the stack, studied the physical dynamics of system accessories such as gas supplement, ﬂow and pressure-regulating devices, then used Simulink to build a mathematical model of a complete PEMFC system; a segmented testing platform was built to test the spatial distribution of RH (relative humidity) and pressure, which was used to verify the simulation model; based on this model, the complicated phenomena occurring inside the stack during ﬂuctuating operating states were calculated. Our ﬁndings showed that the pressure in the gas channel and exhaust manifolds decreased when the external load increased, changing sharply at the moment of load change. The transient pressure di ﬀ erence between the cathode and anode sides (several kPa) had a huge impact on the MEA (membrane electrode assembly); when the load current increased, RH in cathode and cathode channel increased gradually, and the increasing rate of anode side was bigger than that in cathode side. The inﬂuence of variance magnitude and change interval of external load were also studied based on the model. stack’s steady state was disrupted with load fluctuation. Based on these findings, we used the model to study how the load change magnitude and the time interval of load variation affected the stack.


Introduction
It is widely known that the development of new-energy vehicles (NEV), especially those powered by lithium-ion batteries (LIBs), has made great strides in recent years [1]. However, large-scale commercialization of LIB vehicles is hindered by their inadequate driving mileage and long charging times. In this context, the proton exchange membrane fuel cell (PEMFC) is considered a more promising solution for NEV applications [2]. The PEMFC uses hydrogen fuel (143 M J kg −1 ), can be driven further before recharging, and can be refueled in a few minutes. In addition, it is environmentally benign, as water is its only by-product. After the release of the Toyota Mirai and Honda Clarity, PEMFC research gained huge support from governments and powerful syndicates. A 2019 report by the Chinese government [3] states that the country aims to produce 5.5 million PEMFC systems per year, and that the production of PEMFC vehicles will reach 5.2 million per year by 2050.
Currently, PEMFCs still face many challenges before widespread industrial commercialization will be possible. A reliable system that experiences fewer failures and has a long lifetime is vital. Such a The physical phenomena taking place in the cathode and anode channels are relatively easy to understand and can be described using a modified ideal gas state equation [4]. The mass transfer process, which is related to the distribution of moisture, pressure, and temperature inside the stack, is complicated, and numerous papers have been published on the subject. To calculate the mass transfer of hydrogen and air in a porous GDL, the Stefan-Maxwell equation is used [5]. However, for the catalyst layers, it consists of a complex, multi-phase, nanostructured porous material that is difficult to characterize, and it hosts an array of coupled transport phenomena including the flow of gases, liquid water, heat and charged occurring in conjunction with electrochemical reactions. Therefore, the multi-physics dynamics inside the catalyst layers still need to be studied comprehensively. In some research, to simplify the analysis, the catalyst layer (CL) and membrane are usually simplified, using only a permeance parameter related to inner pressure [6].
With respect to water content, different methods have been adopted to monitor water movement in multilayer structures. Steven Holdcroft's group studied liquid water and vapor absorption by a Nafion 211 membrane by measuring weight change, then quantifying the relationship between proton conductivity and membrane humidity [7]. In Japan, Phengxay Deevanhxay's team used high-resolution soft X-ray radiography to monitor the dynamic process of water diffusion in the GDL, CL, and membrane and found that more water gathered at the interface of the GDL and CL than at the interface of the CL and membrane [8]. In reference [9], soft X-ray radiography was used to compare the liquid water distribution in carbon paper and carbon cloth GDL, and the effect of microporous layer (MPL) in liquid water transportation between GDL and CL was also observed by the high resolution method. Other methods, such as magnetic resonance imaging, have also been used to study water distribution [10]. Actually, water sorption on the membrane [11] and transfer across the porous GDL, CL are very complex process. Jay Benziger's team in Princeton University observed the water drops, slugs and flooding process in the GDL and flow channel by experimental methods, and also revealed their influence on PEMFC performance such as current ignition and remaining aspects [12]. On the aspect of simulation method, a mathematical description of changes in water between the liquid and gaseous states under different pressure and temperature conditions can be found in ref. [13]. Studies of electrochemical characteristics-especially in terms of cell voltage-have been published in many journals and books [14], examining Nernst voltage, activation overvoltage, concentration overvoltage, and ohmic resistance loss. The physical phenomena taking place in the cathode and anode channels are relatively easy to understand and can be described using a modified ideal gas state equation [4]. The mass transfer process, which is related to the distribution of moisture, pressure, and temperature inside the stack, is complicated, and numerous papers have been published on the subject. To calculate the mass transfer of hydrogen and air in a porous GDL, the Stefan-Maxwell equation is used [5]. However, for the catalyst layers, it consists of a complex, multi-phase, nanostructured porous material that is difficult to characterize, and it hosts an array of coupled transport phenomena including the flow of gases, liquid water, heat and charged occurring in conjunction with electrochemical reactions. Therefore, the multi-physics dynamics inside the catalyst layers still need to be studied comprehensively. In some research, to simplify the analysis, the catalyst layer (CL) and membrane are usually simplified, using only a permeance parameter related to inner pressure [6].
With respect to water content, different methods have been adopted to monitor water movement in multilayer structures. Steven Holdcroft's group studied liquid water and vapor absorption by a Nafion 211 membrane by measuring weight change, then quantifying the relationship between proton conductivity and membrane humidity [7]. In Japan, Phengxay Deevanhxay's team used high-resolution soft X-ray radiography to monitor the dynamic process of water diffusion in the GDL, CL, and membrane and found that more water gathered at the interface of the GDL and CL than at the interface of the CL and membrane [8]. In reference [9], soft X-ray radiography was used to compare the liquid water distribution in carbon paper and carbon cloth GDL, and the effect of microporous layer (MPL) in liquid water transportation between GDL and CL was also observed by the high resolution method. Other methods, such as magnetic resonance imaging, have also been used to study water distribution [10]. Actually, water sorption on the membrane [11] and transfer across the porous GDL, CL are very complex process. Jay Benziger's team in Princeton University observed the water drops, slugs and flooding process in the GDL and flow channel by experimental methods, and also revealed their influence on PEMFC performance such as current ignition and remaining aspects [12]. On the aspect of simulation method, a mathematical description of changes in water between the liquid and gaseous states under different pressure and temperature conditions can be found in ref. [13]. Studies of electrochemical characteristics-especially in terms of cell voltage-have been published in many journals and books [14], examining Nernst voltage, activation overvoltage, concentration overvoltage, and ohmic resistance loss. Drawing upon such earlier research on the physical characteristics and processes inside the stack, various groups have tried to build models of the stack, combining different components and physical phenomena to calculate the stack state, explain certain special experimental results, and improve the PEMFC system design [15][16][17]. For example, reference [18] reviewed the macroscopic and microscopic models of catalyst layers. This also examined modeling of liquid water transport in the catalyst layer and its implications on the overall transport properties. Review [19] introduced recent progress in terms of cathodic and anodic reactions with a focus on rational design. Xu Liangfei built a model of a PEMFC to explore the reliability of a self-humidification system with cathodic exhaust gas recirculation [20]. A pseudo-2D model of a cathode CL together with a 1D model of a membrane have been proposed to analyze all of the various kinds of voltage loss inside the stack [21]. Similarly, an anode fuel supply system was built [22] to study the influence of a purging valve on system stability. Compared to a lumped parameter calculation, field spatial distribution is a more straightforward approach to reveal the spatial distribution of different fields inside the stack. For example, the temperature distribution of a PEMFC cooling system was simulated using a 1D and 3D combined simulation method to study an optimized thermal management strategy [23]. In reference [24], a non-isothermal, two-phase model was developed to investigate simultaneous heat and mass transfer in the cathode gas diffusion layer GDL. The anisotropy in the GDL properties was taken into account and found to be an important factor controlling the temperature distribution in the GDL. A review [25] in 2009 summarized the development of multi-phase mass transport models using CFD method, also pointed out the problems it faced. These are very useful references when building our own models. Although CFD simulation usually reveals more spatial details of field distribution, it is too complicated for industry use, especially when it is used in system simulation.
Instead of simulation methods, many researchers have tried to reveal the inner state of a PEMFC through experimental testing. Traditional testing can only determine the temperature, humidity, and pressure outside the stack, so a segmented cell method has been developed to monitor the stack's inner state. In segmented cell testing, the stack is divided into dozens of parts on the bipolar plates, and sensors for humidity, current density, pressure, and resistance are inserted in each part of the plates, yielding direct measurements of the stack's inner state [26,27]. Segmented cell testing can provide detailed information about the stack's operating state, which can be used to study dynamic stack characteristics [28], start-stop strategies [29], and aging mechanisms [30]. However, current segmented cell testing is limited by the sensitivity of the sensors, and the sensors inevitably influence the field distribution. Hence, better experimental methods are needed to achieve accurate results, especially when so many physical phenomena are influencing one another. Previous research has made significant progress in understanding physical mechanisms, establishing models, and conducting experimental testing on PEMFCs, but how the dynamic characteristics inside the stack are altered by sharp changes in load has seldom been reported which has great significance for optimizing PEMFC control.
In this paper, we describe studying pressure, gas flow, water content, and electrochemical processes using a mathematical model of a PEMFC system that we built using MATLAB/Simulink. With this model we were able to simulate the instantaneous state inside the stack by making a series of changes in the operating conditions. Applying this dynamic model, we studied water transport and phase transition, the impact of transient pressure on the membrane, and electrochemical changes during the stack operating process. By combining our findings with experimental results, we revealed the effects of sharp changes in load. However, some critical issues, such as the mass transport and phase transition at different interface, water droplets and slugs formation were not considered in this paper because we have difficulty building the complex model of those accurately and quantitatively currently. We aimed to build a simplified model which can be used in industry, so this model has its own limits upon revealing more subtle details.

Mathematical Model and Experimental Platform
We claim the following assumptions made behind the model:

Operational Voltage
A single cell is calculated in this paper, which is consistent with the experimental method. The stack voltage, V cell , can be described as Equation (1): where I cell is the stack current, V Nernst is the Nernst voltage, V act is the activation overvoltage loss, V conc is the concentration overvoltage loss, and R ohm is the ohmic resistance. The value of V Nernst can be calculated by the Nernst equation, as shown in Equation (2): where V 0 is the reversible cell potential at standard pressure, R = 8.31 J mol −1 K −1 and is the ideal gas constant, F = 96,483 C mol −1 and is the Faraday constant, T fc is the cell temperature, the average values of temperature from experimental results under different current were calculated and applied as known parameters in the simulation. P H 2 ,an and P O 2 ,ca are the partial pressures of hydrogen and oxygen gases in the anode and cathode channels, respectively, and P 0 = 100 kPa is the normalization coefficient. The activation voltage loss, V act , can be calculated using a Tafel equation: where α is the charge transfer coefficient, I n is the internal current due to fuel crossover, and I 0 is the exchange current. V conc is the concentration voltage loss and can be calculated using Equation (4): where I 1 is the limiting current at which the fuel is used at the maximum supply rate. The cell ohmic resistance (R ohm ) consists of the resistance of the membrane, R cell , and the conducting resistance, R 0 (see Equation (5)). R cell is a function of water content and cell temperature, and it can be fitted using the experimental data from Peron et al. [7]: where L m is membrane thickness, σ m is membrane conductivity, A cell is the active area of the cell, and λ is the water content in the PEM, the value of which can be obtained using Dutta et al.'s study [31].

Fluid Physics
The dynamic model for the cathode side includes air compressor, humidifier, supply and exhaust manifolds, cathode channel, and back valve; the anode side has similar structures except that hydrogen Energies 2020, 13, 3596 5 of 17 is discharged to the atmosphere directly. The mathematical models of different parts of the anode and cathode will be introduced below.

Gas Flow into the Stack
The hydrogen comes out from a high-pressure tank on the anode side and a flowmeter is used to regulate the fuel flow. On cathode side, air flow is controlled by the air compressor. To reduce calculations, the gas flow is simplified as Equation (6): where W in,tg is the target value for the mass flow of hydrogen and air, and τ is the time delay to reach the target value. W in,tg in anode and cathode sides can be calculated by as follows: W in,antg = 0.00697Nε an Iρ H 2 /60 (8) where N is the cell number, ε ca = 2.0 and ε an = 1.3 are the stoichiometric ratios in cathode and anode, ρ air and ρ H 2 are density of air and hydrogen.

Humidifier
The pressure difference between the humidifier and the supply manifold is negligible. The mole and mass fractions for the humidifier can be calculated as follows: where x vap,hum is the mole fraction of vapor at the output port of the humidifier, P hum is the gas pressure in the humidifier, y gas,hum is the mass fraction of air or hydrogen at the output port of the humidifier, T dew,hum is the dew point of the humidified air, P sat is the saturated vapor pressure, RH is the relative humidity in %, and M gas and M H 2 O are the molecular weights of reactant gases and water, respectively. We used Equation (11) to calculate the flow rate of mixed air and vapor in the humidifier: W hum = W cp y gas,atm /y gas,hum (11) where W cp is the rate of air flowing out from the air compressor, and y gas,atm can be calculated as shown in Equation (10) using the atmospheric parameters.

Supply Manifolds
Assuming the temperature and pressure in the supply manifolds are uniform, the mass flow at the output port of the supply manifolds, W sm,ca , can be expressed as follows (taking cathode side as example): W sm,ca = k sm,ca (P sm,ca − P ca ) (12) where k sm,ca is the nozzle orifice coefficient of the supply manifolds, P sm,ca and P ca are the gas pressures in the supply manifolds and the cathode volume, respectively, and P sm,ca is the supply manifold pressure, which can be expressed as follows: Energies 2020, 13, 3596 6 of 17 In Equation (13), v sm,ca is the volume of the supply manifolds on the cathode side, M sm,ca is the molar mass of the air mixture in the cathode manifolds, and W hum is the mass flow rate of air through the humidifier. The molar mass of the mixed gas in the cathode manifolds M sm,ca can be expressed as: 2.2.4. Cathode Channel The cathode channel contains oxygen, nitrogen, water vapor, and liquid water. The mass balance equations are as follows: The vapor change rate (Equation (17)) in cathode channel equals the algebraic sum of humid gas supplement, electro-osmotic drag water from anode to cathode, water diffusion from cathode to anode, flowing-out vapor, water condensation and membrane sorption (from left to right). Electrophoretic drag from cathode to anode due to concentration difference was neglected in this model because it is very small compared to other factors. Equation (18) shows the liquid water dynamics in cathode. The liquid water change rate is the algebraic sum of reaction product, water condensation and flowing-out liquid water. Actually, the influence of membrane water sorption usually can be neglected due to its slow change.
In these equations, m O 2 , m N 2 , m vap,ca , and m liq,ca are the respective masses of oxygen, nitrogen, water vapor, and liquid water in the cathode volume; y O 2 ,sm , y N 2 ,sm , y vap,sm , y O 2 ,ca , y N 2 ,ca , and y vap,ca are the respective mass fractions of oxygen, nitrogen, and water vapor in the supply manifolds and the cathode volume; W ca is the mass flow rate of the mixed gas at the output port of the cathode volume; W liq,ca is the liquid water mass flow at the output port of the cathode volume; m phase is the water condensation in the cathode volume; M O 2 and M H 2 O are the respective molar masses of oxygen and water vapor; N H 2 O,diff is the water diffusion parameter which can be calculated according to ref [14] and n d describes the number of water molecules dragged per H + ion moved by electric field through the membrane, according to Ref. [5], n d = 2.5λ/22.
The water content in the membrane λ is calculated by Equation (19). a is the water activity calculated by the mean value of relative humidity in the anode and cathode channels (Equation (20)) [14]. The water sorption equilibrium process in membrane usually costs more than 10 5 s, so a delay module (1/(1 + 10 5 s)) was applied in the model to simulate the slow process of water sorption. λ = 0.043 + 17.18a − 39.85a 2 + 36a 3 (19) a = (P ca x vap,ca /P sat,ca + P an x vap,an /P sat,an )/2 (20) In addition, the water droplet and slug formation process would affect the mass transport across GDL and flow channel, which in turn affects the PEMFC performance [12,32]. However, the model built in this paper is a lumped parameter model and it is limited to describe the local clogging phenomenon during mass transport, so it has to be assumed that the the gases and fluid water transport successfully inside the stack. The pressure in the cathode volume, P ca , can be calculated as follows: The water condensation rate, m phase , can be described as RT ca x vap,ca (P ca x vap,ca − P sat )M H 2 O P ca x vap,ca ≥ P sat k e m liq,ca (P ca x vap,ca − P sat ) P ca x vap,ca < P sat (22) where k c is the steam condensation rate, k e is the evaporation rate of liquid water, x vap,ca is the molar fraction of water vapor in the cathode, and P sat is the saturated water vapor pressure. The liquid water outflow from the cathode volume, W liq,ca , can be expressed as follows: where k liq,ca is a coefficient.

Anode Channel
The anode volume is comprised of hydrogen and water. The mass balance equations are as follows: The pressure in the cathode volume, P an , can be calculated as follows:

Cathode Exhaust Manifolds
The gas pressure in the exhaust manifolds of the cathode side, P em,ca , can be expressed as follows: where v em,ca is the volume of the cathode exhaust manifolds, M em,ca is the molar mass of the mixed gas, and W em,ca is the mass flow rate at the output port of the manifolds. W ca can be expressed as follows: where k ca is the nozzle orifice coefficient at the interface of the cathode volume and exhaust manifolds.

Anode Exhaust Manifolds
The gas pressure in the exhaust manifolds of the cathode side, P em,an , can be expressed as follows: dP em,an dt = RT cell v em,an M em,an (W an − W em,an ) Energies 2020, 13, 3596 8 of 17 The two mass flows, W em,an and W an , can be expressed as follows: W an = k an (P an − P em,an ) W em,an = k em,an (P em,an − P atm )

Back-Pressure Valve
The back-pressure valve is used to regulate the stack's cathode inlet pressure. Laval nozzle flow theory is used to describe its mathematical features: where s is the opening area of the valve and is tested by experimental platform; ρ em,ca is the gas density before the valve; and P em,ca and P atm are the pressures before and after the valve, respectively. All the parameters used in the model are listed in Table 1.

Segmented Cell Testing Experiment to Verify the Model
A segmented cell testing experiment to test the pressure and humidity distribution on bipolar plates was conducted and is shown in Figure 2. In this study, each plate (graphite plate) was divided into 36 zones, and in each zone, humidity and pressure sensors were installed to measure the humidity and pressure distribution for different current densities. The relative humidity (RH) measurement is based on sensitive change of polymer material capacitance under different humidity, while the pressure testing is based on the resistance change of pressure sensitive resistor. The test errors of RH and pressure are 1.5% and ±4.0 kPa respectively. PCB (print circuit board) is used to collect signals from sensors and communicate with the upper computer with Labview. To reduce the influence of so many sensors, only one field was tested in each experiment. into 36 zones, and in each zone, humidity and pressure sensors were installed to measure the humidity and pressure distribution for different current densities. The relative humidity (RH) measurement is based on sensitive change of polymer material capacitance under different humidity, while the pressure testing is based on the resistance change of pressure sensitive resistor. The test errors of RH and pressure are 1.5% and ±4.0 kPa respectively. PCB (print circuit board) is used to collect signals from sensors and communicate with the upper computer with Labview. To reduce the influence of so many sensors, only one field was tested in each experiment.

Calculations and Experimental Conditions
In the simulation and experimental processes, the inlet pressures were set at 203 kPa and 201 kPa at the anode and cathode sides, respectively, with relative humidity (RH) values of 100% and 50%. The hydrogen and air were in co-flow mode. The requisite fuel and oxidant flows were calculated using stack current and stoichiometric ratios. The stoichiometric ratios were 1.3 and 2.0, respectively, consistent with the experimental platform.
In the simulation process, a current of 55.8 A was selected as the initial condition for the calculations, then the load current was increased with time, as shown in Table 2, and the dynamic process inside the stack was calculated.

Pressure Dynamics
The pressure dynamics in the supply manifolds, channel and exhaust manifolds can be obtained, as can be seen in Figure 3. The pressure fluctuated when load current changed with time, and it shows obvious declines in the pressure from the supply manifolds, to the channel, to the exhaust manifolds. The pressure inside the channel and exhaust manifold declined gradually as the load current rose, because of the increasing oxygen consumption and the action of the back-pressure valve. On the anode side, the pressure in the channel and exhaust manifolds also decreased with time.
50%. The hydrogen and air were in co-flow mode. The requisite fuel and oxidant flows were calculated using stack current and stoichiometric ratios. The stoichiometric ratios were 1.3 and 2.0, respectively, consistent with the experimental platform.
In the simulation process, a current of 55.8 A was selected as the initial condition for the calculations, then the load current was increased with time, as shown in Table 2, and the dynamic process inside the stack was calculated.

Pressure Dynamics
The pressure dynamics in the supply manifolds, channel and exhaust manifolds can be obtained, as can be seen in Figure 3. The pressure fluctuated when load current changed with time, and it shows obvious declines in the pressure from the supply manifolds, to the channel, to the exhaust manifolds. The pressure inside the channel and exhaust manifold declined gradually as the load current rose, because of the increasing oxygen consumption and the action of the back-pressure valve. On the anode side, the pressure in the channel and exhaust manifolds also decreased with time.
(a) (b) The pressure distribution from the segmented cell platform is shown in Figure 4 (in the interest of space, only the cathode channel pressure is shown as an example). It can be seen that the pressure declined gradually from stack inlet (right) to stack outlet (left) which consistent with the simulation results.  The pressure distribution from the segmented cell platform is shown in Figure 4 (in the interest of space, only the cathode channel pressure is shown as an example). It can be seen that the pressure declined gradually from stack inlet (right) to stack outlet (left) which consistent with the simulation results.  The average pressure from the experiment inside the channel was calculated and compared with the simulation results (the current was increased every 3 minutes considering the operating condition of the stack), as shown in Figure 5. It can be seen that the experimental pressure and calculated pressure showed the same trend with current. In addition, the average relative error of calculated results is less than 5%, which verifies the simulation model built in this study to some extent. The average pressure from the experiment inside the channel was calculated and compared with the simulation results (the current was increased every 3 min considering the operating condition of Energies 2020, 13, 3596 11 of 17 the stack), as shown in Figure 5. It can be seen that the experimental pressure and calculated pressure showed the same trend with current. In addition, the average relative error of calculated results is less than 5%, which verifies the simulation model built in this study to some extent.  The average pressure from the experiment inside the channel was calculated and compared with the simulation results (the current was increased every 3 minutes considering the operating condition of the stack), as shown in Figure 5. It can be seen that the experimental pressure and calculated pressure showed the same trend with current. In addition, the average relative error of calculated results is less than 5%, which verifies the simulation model built in this study to some extent.

Water Dynamics
Water content and its transition between states in the channel are important indicators for monitoring states and diagnosing faults in a PEMFC. Values for vapor, liquid water, and the water condensation are shown in Figure 6. On the cathode side, when the load current increased, the amount of vapor and liquid water increased gradually with current. Water condensation rate fluctuated only slightly. On the anode side, vapor increased with time slightly under the multi-process of water diffusion and electro-osmotic drag, while water condensation and liquid water seldom existed.

Water Dynamics
Water content and its transition between states in the channel are important indicators for monitoring states and diagnosing faults in a PEMFC. Values for vapor, liquid water, and the water condensation are shown in Figure 6. On the cathode side, when the load current increased, the amount of vapor and liquid water increased gradually with current. Water condensation rate fluctuated only slightly. On the anode side, vapor increased with time slightly under the multi-process of water diffusion and electro-osmotic drag, while water condensation and liquid water seldom existed.  The RH distribution under different load currents was obtained by the experimental method, the results of which are shown in Figure 7. RH clearly increased from stack inlet to stack outlet, based on the field distribution data.  The RH distribution under different load currents was obtained by the experimental method, the results of which are shown in Figure 7. RH clearly increased from stack inlet to stack outlet, based on the field distribution data.
The average RH value was calculated from the experimental data, and Figure 8 compares the experimental and simulation results (the RH results from the simulation were calculated using Equation (9)). It can be seen that RH in cathode channel and anode channel both increased with current under the conditions set in this paper. The RH distribution under different load currents was obtained by the experimental method, the results of which are shown in Figure 7. RH clearly increased from stack inlet to stack outlet, based on the field distribution data. The average RH value was calculated from the experimental data, and Figure 8 compares the experimental and simulation results (the RH results from the simulation were calculated using Equation (9)). It can be seen that RH in cathode channel and anode channel both increased with current under the conditions set in this paper. The RH distribution under different load currents was obtained by the experimental method, the results of which are shown in Figure 7. RH clearly increased from stack inlet to stack outlet, based on the field distribution data. The average RH value was calculated from the experimental data, and Figure 8 compares the experimental and simulation results (the RH results from the simulation were calculated using Equation (9)). It can be seen that RH in cathode channel and anode channel both increased with current under the conditions set in this paper. With the simulation model above, we were able to obtain the membrane water and pressure impact on the MEA (i.e., the difference between the cathode and anode pressure), which we could not measure experimentally. Water content in the membrane increased very slowly with load current according to the calculation results shown in Figure 9, because the water sorption process in the membrane usually cost more than 10 5 s. With load variation, the pressure difference between the cathode and anode channel pressure, and hence the pressure impact on the MEA, also fluctuated. The impact was particularly large, in the order of kPa.
Based on the above results, it can be concluded that the PEMFC underwent obvious state change inside the stack when the load current varied over time, potentially leading to unpredictable faults and even adversely influencing the stack's lifespan. In the aspect of pressure dynamics, the anode and cathode pressure decreased with time because more oxygen and hydrogen were consumed when a larger current was needed. On the cathode side, the back-pressure valve connected to the atmosphere usually had a bigger opening when a larger load occurs, causing the channel pressure to drop again. Changes in anode and cathode pressure meant the MEA was subjected to pressure fluctuations, and steep pressure increases threatened the MEA's structure and materials.
Water content is also a very important indicator for evaluating a stack's operating state. The water dynamics in the flow channel comprise a complicated balancing process involving reactive synthesis, bidirectional diffusion, phase transition, and so on. Under the conditions described in this paper, the amount of water vapor and the water condensation rate in the cathode channel changed slightly with the load current, while the amount of liquid water gradually increased. On the anode side, water vapor increased with load, while there were nearly no liquid water and water condensation existing. With the simulation model above, we were able to obtain the membrane water and pressure impact on the MEA (i.e., the difference between the cathode and anode pressure), which we could not measure experimentally. Water content in the membrane increased very slowly with load current according to the calculation results shown in Figure 9, because the water sorption process in the membrane usually cost more than 10 5 s. With load variation, the pressure difference between the cathode and anode channel pressure, and hence the pressure impact on the MEA, also fluctuated. The impact was particularly large, in the order of kPa. Based on the above results, it can be concluded that the PEMFC underwent obvious state change inside the stack when the load current varied over time, potentially leading to unpredictable faults and even adversely influencing the stack's lifespan. In the aspect of pressure dynamics, the anode and cathode pressure decreased with time because more oxygen and hydrogen were consumed when a larger current was needed. On the cathode side, the back-pressure valve connected to the atmosphere usually had a bigger opening when a larger load occurs, causing the channel pressure to drop again. Changes in anode and cathode pressure meant the MEA was subjected to pressure fluctuations, and steep pressure increases threatened the MEA's structure and materials.
Water content is also a very important indicator for evaluating a stack's operating state. The water dynamics in the flow channel comprise a complicated balancing process involving reactive synthesis, bidirectional diffusion, phase transition, and so on. Under the conditions described in this paper, the amount of water vapor and the water condensation rate in the cathode channel changed slightly with the load current, while the amount of liquid water gradually increased. On the anode side, water vapor increased with load, while there were nearly no liquid water and water condensation existing.

Discussion
As described above, we were able to verify the model's reliability to some extent, and we found that the stack's steady state was disrupted with load fluctuation. Based on these findings, we used the model to study how the load change magnitude and the time interval of load variation affected the stack.

Discussion
As described above, we were able to verify the model's reliability to some extent, and we found that the stack's steady state was disrupted with load fluctuation. Based on these findings, we used the model to study how the load change magnitude and the time interval of load variation affected the stack.

Role of Load Change Magnitude
A current of 55.8 A was chosen as the initial value. It was then increased to 111.6 A, 167.4 A, 223.2 A, and 279 A over a period of 10 s. Figure 10 presents a curve showing the relationship between state change and current change magnitude. It can be seen from Figure 10a that as the current was increased from 55.8 A, the magnitude of the pressure impact and output voltage increased almost linearly, while membrane water change only increased slightly. Figure 10b shows that the RH change in the cathode and anode channels also increased with current change, and the increasing rate of anode was bigger in this paper.

Role of Load Change Magnitude
A current of 55.8 A was chosen as the initial value. It was then increased to 111.6 A, 167.4 A, 223.2 A, and 279 A over a period of 10 s. Figure 10 presents a curve showing the relationship between state change and current change magnitude. It can be seen from Figure 10a that as the current was increased from 55.8A, the magnitude of the pressure impact and output voltage increased almost linearly, while membrane water change only increased slightly. Figure 10b shows that the RH change in the cathode and anode channels also increased with current change, and the increasing rate of anode was bigger in this paper.

Role of Time Interval of Load Change
In actual operation, a PEMFC's operational mode changes over time, usually at very short intervals, which leads to sharp state changes in the PEMFC. According to our calculated results, the trends in pressure, RH in cathode and anode at different time intervals were similar when the time was greater than 5 s. However, when the interval was reduced to 1 s, the PEMFC struggled to maintain a steady state (Figure 11b). It can also be concluded from Figure 11 that a shorter time interval altered the normal trend of the state parameters, resulting in even greater stress on the MEA.

Role of Time Interval of Load Change
In actual operation, a PEMFC's operational mode changes over time, usually at very short intervals, which leads to sharp state changes in the PEMFC. According to our calculated results, the trends in pressure, RH in cathode and anode at different time intervals were similar when the time was greater than 5 s. However, when the interval was reduced to 1 s, the PEMFC struggled to maintain a steady state (Figure 11b). It can also be concluded from Figure 11 that a shorter time interval altered the normal trend of the state parameters, resulting in even greater stress on the MEA.