Numerical Study of the Dynamic Response of Heat and Mass Transfer to Operation Mode Switching of a Unitized Regenerative Fuel Cell

Knowledge concerning the complicated changes of mass and heat transfer is desired to improve the performance and durability of unitized regenerative fuel cells (URFCs). In this study, a transient, non-isothermal, single-phase, and multi-physics mathematical model for a URFC based on the proton exchange membrane is generated to investigate transient responses in the process of operation mode switching from fuel cell (FC) to electrolysis cell (EC). Various heat generation mechanisms, including Joule heat, reaction heat, and the heat attributed to activation polarizations, have been considered in the transient model coupled with electrochemical reaction and mass transfer in porous electrodes. The polarization curves of the steady-state models are validated by experimental data in the literatures. Numerical results reveal that current density, gas mass fractions, and temperature suddenly change with the sudden change of operating voltage in the mode switching process. The response time of temperature is longer than that of current density and gas mass fractions. In both FC and EC modes, the cell temperature and gradient of gas mass fraction in the oxygen side are larger than that in the hydrogen side. The temperature difference of the entire cell is less than 1.5 K. The highest temperature appears at oxygen-side catalyst layer under the FC mode and at membrane under a more stable EC mode. The cell is exothermic all the time. These dynamic responses and phenomena have important implications for heat analysis and provide proven guidelines for the improvement of URFCs mode switching.


Introduction
Science and technology are developing by leaps and bounds, leading to higher energy storage demands.The unitized regenerative fuel cell (URFC) is perceived as one of the cleanest and most effective energy storage and conversion device, whose special energy is several times higher than that of the lightest secondary battery [1,2].Compared with secondary batteries, URFCs have the advantages of no self-discharge or cell capacity limitation.The high special energy of URFCs indicates limitless applications for some weight-critical and time-consuming portable applications such as space energy systems (>40 Wh•kg −1 ) [3].In addition, URFCs can be developed for and applied in high-altitude long-endurance solar aircraft, hybrid energy storage spacecraft propulsion systems, remote area energy storage systems without relying on the grid, power systems for power grid peak adjustment, and portable power source systems [4][5][6].
Hamour et al. [35] adopted a transient model to investigate the heat conducting characteristics of the GDL for a PEM FC and a PEM EC.However, no article analyzing the heat response for URFCs has been found.Heat and mass transfer significantly affect the performance of FC, water electrolysis, and URFCs.Numerical simulation is a cost-effective and time-effective approach in analyzing the intricate change inside the URFCs.
A transient, non-isothermal, single-phase, multi-physics, and full-cell mathematical model for PEM URFC is generated in the present study to analyze the dynamic responses of a PEM URFC under the operation mode switching and enhance the fundamental understandings of the complex interactions of heat and mass transfer during the cell transient operations further.The model is coupled with electrochemical reactions.Heat generation of URFC without phase change includes irreversible heat attributed to activation polarizations and ohmic heat and reversible entropic heat of reactions.The influences of temperature on the parameters of materials and equilibrium potential are considered in this model.The multi-component species in hydrogen and oxygen electrodes transmit in porous medium.The URFC is switched from the FC mode to the EC mode at 2 s.This study aims to explore the intricate changes of heat and mass transfer as well as current density under operation mode switching.

Model Description
Figure 1 shows the computational domain, which is optioned along with the gas flow channels of a PEM URFC.The geometric model consists of oxygen-side and hydrogen-side gas flow channels (OGFC and HGFC, respectively), oxygen-side and hydrogen-side gas diffuse layers (OGDL and HGDL, respectively), oxygen-side and hydrogen-side catalyst layers (OCL and HCL, respectively), and a PEM-based electrolyte.
Energies 2016, 9, 1015 3 of 19 characteristics of the GDL for a PEM FC and a PEM EC.However, no article analyzing the heat response for URFCs has been found.Heat and mass transfer significantly affect the performance of FC, water electrolysis, and URFCs.Numerical simulation is a cost-effective and time-effective approach in analyzing the intricate change inside the URFCs.A transient, non-isothermal, single-phase, multi-physics, and full-cell mathematical model for PEM URFC is generated in the present study to analyze the dynamic responses of a PEM URFC under the operation mode switching and enhance the fundamental understandings of the complex interactions of heat and mass transfer during the cell transient operations further.The model is coupled with electrochemical reactions.Heat generation of URFC without phase change includes irreversible heat attributed to activation polarizations and ohmic heat and reversible entropic heat of reactions.The influences of temperature on the parameters of materials and equilibrium potential are considered in this model.The multi-component species in hydrogen and oxygen electrodes transmit in porous medium.The URFC is switched from the FC mode to the EC mode at 2 s.This study aims to explore the intricate changes of heat and mass transfer as well as current density under operation mode switching.

Model Description
Figure 1 shows the computational domain, which is optioned along with the gas flow channels of a PEM URFC.The geometric model consists of oxygen-side and hydrogen-side gas flow channels (OGFC and HGFC, respectively), oxygen-side and hydrogen-side gas diffuse layers (OGDL and HGDL, respectively), oxygen-side and hydrogen-side catalyst layers (OCL and HCL, respectively), and a PEM-based electrolyte.

Model Assumptions
Suitable assumptions are considered based on the purpose of investigating the characteristics of heat and mass transfer to make the numerical simulation simple and tractable.The main assumptions are listed as follows:

•
Single phase.No phase change occurs in this model, and water, including the reaction product water, is maintained in the gaseous state to simplify the model [23,36,37].

•
Heat transfer.The MEA (including GDLs, CLs and a PEM) is treated as heat transfer in a solid, and the effect of gas flow in a porous medium on heat transfer is ignored.

•
Species transport.The flow state is laminar because of the low Reynolds number.All gas mixtures are treated as incompressible ideal gasses [23,37,38].

•
Material characteristics.The PEM is assumed to be impermeable to gas species [23].In addition, GDLs and CLs are homogeneous and isotropic, and the contact resistance between different layers is ignored [39]. .OGDL: oxygen-side gas diffuse layer; HGDL: hydrogen-side gas diffuse layer; OGFC: oxygen-side gas flow channels; and HGFC: hydrogen-side gas flow channels.

Model Assumptions
Suitable assumptions are considered based on the purpose of investigating the characteristics of heat and mass transfer to make the numerical simulation simple and tractable.The main assumptions are listed as follows:

•
Single phase.No phase change occurs in this model, and water, including the reaction product water, is maintained in the gaseous state to simplify the model [23,36,37].

•
Heat transfer.The MEA (including GDLs, CLs and a PEM) is treated as heat transfer in a solid, and the effect of gas flow in a porous medium on heat transfer is ignored.

•
Species transport.The flow state is laminar because of the low Reynolds number.All gas mixtures are treated as incompressible ideal gasses [23,37,38].

•
Material characteristics.The PEM is assumed to be impermeable to gas species [23].In addition, GDLs and CLs are homogeneous and isotropic, and the contact resistance between different layers is ignored [39].

Governing Equations
In this study, a two-dimensional, non-isothermal, single-phase, multi-physics, transient, and full cell model of URFC is generated numerically, which accounts for the transfer of charge, species, and energy.Equations are presented below to describe these processes.

Charge Transfer and Current Distribution
From Ohm's law, the ionic and electronic currents of electrodes and membrane can be written as: where φ M and φ S are the ionic and electronic potentials, respectively.σ eff M and σ eff S are the effective ionic and electronic conductivities of the electrodes and membrane, respectively.For isotropic GDLs and CLs, σ eff S = σ S ; the value of σ S is listed in Table 1.For the membrane, σ eff M = σ M ; for the CLs, , where ε l is the porosity of catalyst layers.The value of σ M is determined by using the following empirical correlation: where λ is the water content, which is defined as a function of the water activity a. a is assumed to be equal to 1 in this model [27].
In the HCL and OCL, the conservation charge can be described as follows: where S φ S and S φ M are the source terms, which are related to the volume current density in CLs.In the FC mode, S φ S = i V in the HCL, S φ S = −i V in the OCL, S φ M = −i V in the HCL, and S φ M = i V in OCL.In the EC mode, S φ S = −i V in the HCL, S φ S = i V in the OCL, S φ M = i V in the HCL, and S φ M = −i V in the OCL.i V in the hydrogen and oxygen sides is expressed by the Butler-Volmer equation, as follows: where i 0 is the exchange current density, which is dependent on reactant concentrations and reference exchange current density i ref 0 .a υ and η refer to the active specific surface area and over-potential, respectively.They can be described as follows [27]: where E eq is the equilibrium potential, which can be defined as E eq = 0 in the hydrogen side and E eq = 1.23 − 9 × 10 −4 (T − 298.15) in the oxygen side [27].

Gas Flow Equations
In the gas flow channel the Navier-Stokes equation is used to govern momentum transfer: There is an equation of mass conservation: For the porous medium of GDLs and CLs, the following Brinkman equations can describe the flow: where ε and κ are the porosity and permeability of the porous electrodes, respectively.S m is the source term, which is expressed by Equation ( 21) in CLs.For GDLs, S m = 0.

Heat Transfer
In thermodynamics, the total heat generation of a URFC can be classified into reversible heat and irreversible heat.Reversible heat is the chemical reaction heat, whose value is positive in the FC mode and negative in the EC mode.Reversible heat can be calculated as follows: where ∆S is the entropy change of the overall reaction ( 2 in the EC mode), U 0 is the thermodynamic equilibrium potential, and I is the current of the cell.Meanwhile, irreversible heat released from a cell is caused by ohmic polarization and activation polarization.This part of the heat generation rate can be written as follows: where V cell represents the cell operating voltage.The energy conservation equation can be expressed as follows: where k eff is the effective or equivalent thermal conductivity.In porous media, k eff is calculated as follows: where k g and k s are the thermal conductivities of the porous solid frame and the gases, respectively.The materials of the HGDL, OGDL, and CLs are carbon, titanium, and platinum, respectively.Their thermal conductivities are listed in Table 1 Heat conduction through the ribs of the flow field plate makes the most significant contribution to heat transfer.The equivalent thermal conductivity k eq in the y-direction is used for the gas flow channels to simulate the temperature of MEA accurately, and in the flow direction, k eff = k g .As shown in Figure 2, the derivation process of k eq is as follows: where the d rib and d ch are the width of a rib and the gas flow channel.Both the values of parameter d rib and d ch are 2 mm.The heat source S T , the last term of Equation ( 24), considers three kinds of heat sources in a single-phase model, namely, reversible chemical reaction heat, heat attributed to activation polarizations, and Joule heat.In CLs, the three types of heat sources are presented as follows: Energies 2016, 9, 1015 7 of 19 where is the effective or equivalent thermal conductivity.In porous media, is calculated as follows: where and are the thermal conductivities of the porous solid frame and the gases, respectively.The materials of the HGDL, OGDL, and CLs are carbon, titanium, and platinum, respectively.Their thermal conductivities are listed in Table 1 Heat conduction through the ribs of the flow field plate makes the most significant contribution to heat transfer.The equivalent thermal conductivity in the y-direction is used for the gas flow channels to simulate the temperature of MEA accurately, and in the flow direction, = .As shown in Figure 2, the derivation process of is as follows: where the and are the width of a rib and the gas flow channel.Both the values of parameter and are 2 mm.The heat source T S , the last term of Equation ( 24), considers three kinds of heat sources in a single-phase model, namely, reversible chemical reaction heat, heat attributed to activation polarizations, and Joule heat.In CLs, the three types of heat sources are presented as follows: In the GDLs and membrane, only Joule heat caused by electronic and protonic current flows exists.Equations ( 28) and ( 29) are the expressions of heat source in the GDLs and membrane, respectively.For the gas flow channel, the heat sources are equal to zero: In the GDLs and membrane, only Joule heat caused by electronic and protonic current flows exists.Equations ( 28) and ( 29) are the expressions of heat source in the GDLs and membrane, respectively.For the gas flow channel, the heat sources are equal to zero:

Boundary Conditions and Initial Conditions
Boundary conditions.The boundary conditions are listed in Table 2.At all of the outside borders, zero heat or mass flux condition is applied, except for the specific description in Table 2.In addition, the model considers laminar flow in the inlet of the gas flow channel and no slip wall in the boundary of y = 0 and y = 5.134 mm.Initial conditions.The initial values of velocity, pressure, and temperature for the entire cell are expressed as follows: u = 0 m/s, p = 1 atm, and T = 353.15K. Potentials in the GDL, CL at the hydrogen side, and PEM are φ M = 0 and φ S = 0. Potentials in the GDL, and CL at the oxygen side are φ M = 0 and φ S = V cell .The mass fractions of H 2 in the hydrogen side and O 2 in the oxygen side are 0.9 and 0.8, respectively.The mass fraction of H 2 O are 0.1 in the hydrogen side and 0.2 in the oxygen side.

Element Independence Test
The finite volume method is applied to the discrete model equations presented previously.For calculating the errors, a convergent solution is considered when the relative differences of all variables between two consecutive iterations are lower than 10 −5 .The quality of the grid and the number of elements have an effect on the numerical results.A local mesh refinement exists in MEA because of the complex electrochemical reactions and intricate interactions of heat and mass transfer.A large number of elements will lead to an increase in the calculation amount, and small number of elements will lead to a large error.Thus the steady-state models of the FC and EC modes are calculated seven times with different numbers of elements.Figure 3a,b shows the temperature distribution of the central line of the membrane with an operating voltage of 0.6 V and the temperature distribution of the central line of the OCL with an operating voltage of 1.8 V, respectively, under the seven different numbers of elements.As shown in Figure 3, the temperature curves almost coincided with each other when the number of elements is larger than 5300.Thus, 5300 is selected to be the most suitable number of elements in the model to meet the requirements of reducing the computation time and ensuring the accuracy of the numerical results simultaneously.
distribution of the central line of the membrane with an operating voltage of 0.6 V and the temperature distribution of the central line of the OCL with an operating voltage of 1.8 V, respectively, under the seven different numbers of elements.As shown in Figure 3, the temperature curves almost coincided with each other when the number of elements is larger than 5300.Thus, 5300 is selected to be the most suitable number of elements in the model to meet the requirements of reducing the computation time and ensuring the accuracy of the numerical results simultaneously.

Model Validation
Measuring the concentration and temperature fields in the GDLs, CLs and the membrane in the experiment is difficult because of its small thickness.The results computed at the steady state of this model in the FC and EC modes are applied to compare with the experimental data [2,3] in the

Model Validation
Measuring the concentration and temperature fields in the GDLs, CLs and the membrane in the experiment is difficult because of its small thickness.The results computed at the steady state of this model in the FC and EC modes are applied to compare with the experimental data [2,3] in the literatures.The basic conditions of the numerical models match with those of the experimental data to validate the numerical models.The operating temperatures in Figure 4a,b are 353.15 and 365.15 K, respectively.As shown in Figure 4a, the polarization curves (voltage vs. current density) of the simulation results, particularly under the FC mode, are consistent with experimental data.The performance of the experiment in the EC mode is better than that of the present model.This difference is likely because of gaseous state of the water in the present model.As shown in Figure 4b, the simulation results at 365.15 K are relatively consistent with the experimental data without any adjustment on the parameters.The slight difference between the two may be attributed to the constant value of water activity in the membrane.Parameters, such as the exact cell geometries, are unknown in [2,3].

Results and Discussion
This study deals with the understanding of the dynamic response of an operating PEM-based URFC during the mode switching process.The cell operating mode in FC or EC is controlled by cell voltage, as shown in Figure 5. Cell voltage in the FC and EC modes are set as 0.6 V in the first 2 s and 1.5 V during the subsequent 3 s.In the course of the switching process, the key transient responses of the current density, mass transfer, and heat transmission would occur at approximately 2 s, and the results will be shown and discussed in this section.

Results and Discussion
This study deals with the understanding of the dynamic response of an operating PEM-based URFC during the mode switching process.The cell operating mode in FC or EC is controlled by cell voltage, as shown in Figure 5. Cell voltage in the FC and EC modes are set as 0.6 V in the first 2 s and 1.5 V during the subsequent 3 s.In the course of the switching process, the key transient responses of the current density, mass transfer, and heat transmission would occur at approximately 2 s, and the results will be shown and discussed in this section.

Results and Discussion
This study deals with the understanding of the dynamic response of an operating PEM-based URFC during the mode switching process.The cell operating mode in FC or EC is controlled by cell voltage, as shown in Figure 5. Cell voltage in the FC and EC modes are set as 0.6 V in the first 2 s and 1.5 V during the subsequent 3 s.In the course of the switching process, the key transient responses of the current density, mass transfer, and heat transmission would occur at approximately 2 s, and the results will be shown and discussed in this section.The water mass fractions along y-axis with varying times are displayed in Figure 7.The solid lines represent the temperature distribution along A-A line and the dotted lines represent the temperature distribution along the C-C line.Under the FC mode, the water mass fraction along the C-C line is evidently higher than that along the A-A line.The main reasons are the consumption of reactants and generation of H 2 O in oxygen side.Therefore, the H 2 O mass fraction increases gradually along the gas flow direction.Moreover, the consumption and the generation of gases occur in the CLs.H 2 O mass fraction increases gradually from gas flow channels to CLs because of gas diffusion.However, the distribution of H 2 O mass fraction under the EC mode is different from that under the FC mode because of the opposite electrochemical reactions.The mass fraction gradient in the hydrogen side is smaller compared with that in the oxygen side because of the generation of water and the small molecular weight of hydrogen.The water mass fraction at 0.5 and 1.99 s did not change at the A-A line.Meanwhile, at the C-C line, the water mass fraction at 1.99 s is slightly larger than that at 0.05 s, which illustrates that the mass fraction near the inlet of the gas flow channel reaches a steady state faster.The H 2 O mass fraction decreases in different degrees within 0.02 s from t = 1.99 s to t = 2.01 s.Along the y-direction from gas flow channel to CLs, the H 2 O mass fraction significantly changes.At 2.01 s, the maximum value of H 2 O mass fraction appears on the GDLs.The H 2 O mass fraction dramatically declines in the CLs because of mode switching from the FC mode to the EC mode at 2.0 s and the opposite electro-chemical reactions of the two modes.Eventually, the H 2 O mass fraction declines in the GDLs and gas flow channels under the action of diffusion.
mA•cm at 2 s because of mode switching.Then, the current density rapidly increases to approximately 15 mA•cm −2 corresponding to the EC mode.At the same time, O2 mass fraction and temperature at that point change dramatically.These changes can be explained as follows: The accumulated water, which is produced in the FC mode, is the reactant in the EC mode.The concentration of H2O reduces rapidly at the beginning of the EC mode, which leads to a slight decrease in the current density.The O2 mass fraction increases sharply because of the consumption of H2O and the production of O2 when the mode is switched from FC mode to EC mode.Temperature decreases rapidly at 2 s mainly because the exothermic reaction in the FC mode converts into the endothermic reaction in the EC mode.During the subsequent 3 s, the temperature gradient diminishes gradually.Temperature does not reach a steady state during the process.Current density, O2 mass fraction and temperature in the OCL suddenly change under operation mode switching.The response time of temperature is longer than that of current density and O2 mass fraction.The water mass fractions along y-axis with varying times are displayed in Figure 7.The solid lines represent the temperature distribution along A-A line and the dotted lines represent the temperature distribution along the C-C line.Under the FC mode, the water mass fraction along the C-C line is evidently higher than that along the A-A line.The main reasons are the consumption of reactants and generation of H2O in oxygen side.Therefore, the H2O mass fraction increases gradually along the gas flow direction.Moreover, the consumption and the generation of gases occur in the CLs.H2O mass fraction increases gradually from gas flow channels to CLs because of gas diffusion.However, the distribution of H2O mass fraction under the EC mode is different from that under the FC mode because of the opposite electrochemical reactions.The mass fraction gradient in the hydrogen side is smaller compared with that in the oxygen side because of the generation of water and the small molecular weight of hydrogen.The water mass fraction at 0.5 and 1.99 s did not change at the A-A line.Meanwhile, at the C-C line, the water mass fraction at 1.99 s is slightly larger than that at 0.05 s, which illustrates that the mass fraction near the inlet of the gas flow channel reaches a steady state faster.The H2O mass fraction decreases in different degrees within 0.02 s from t = 1.99 s to t = 2.01 s.Along the y-direction from gas flow channel to CLs, the H2O mass fraction significantly changes.At 2.01 s, the maximum value of H2O mass fraction appears on the GDLs.The H2O mass fraction dramatically declines in the CLs because of mode switching from the FC mode to the EC mode at 2.0 s and the opposite electro-chemical reactions of the two modes.Eventually, the H2O mass fraction declines in the GDLs and gas flow channels under the action of diffusion.Figure 8 displays the distribution of H2 and O2 mass fractions along the central lines of the CLs.The results indicate that both O2 and H2 mass fractions are depressed distinctly along the flow direction at 1.99 s.The O2 and H2 mass fractions increase when the mode is switched to the EC mode.Eventually, the O2 and H2 mass fractions increase along the flow direction under a more stable EC mode at 2.5 s, and the magnitudes of O2 and H2 mass fractions are greater than those at the inlets of gas flow channels because of the production of O2 and H2 in the CLs.As with the H2O mass fraction, the gradient of the O2 mass fraction is significantly larger than that of the H2 mass fraction.
Figure 9 shows the distribution of temperature and current density along the central line of the OCL.At the end of the FC mode (1.99 s), the general trends of temperature and current density decline gradually along the x-axis.The distribution of current density probably results from the concentrations of the reactants.As shown in Figure 8, the O2 mass fraction declines along the x-axis.The temperature is determined by the heat source and heat flow.For the homogenous and isotropic CL, a relatively high temperature is observed where the heat source is large.The heat source in the OCL consists of heat attributed to activation polarizations, Joule heat, and reaction heat, all of which are positively correlated with current density.Therefore temperature is depressed along the x-axis.Temperature and current density show remarkable reduction during 0.02 s from 1.99 s to 2.01 s.The  Eventually, the O 2 and H 2 mass fractions increase along the flow direction under a more stable EC mode at 2.5 s, and the magnitudes of O 2 and H 2 mass fractions are greater than those at the inlets of gas flow channels because of the production of O 2 and H 2 in the CLs.As with the H 2 O mass fraction, the gradient of the O 2 mass fraction is significantly larger than that of the H 2 mass fraction.
Figure 9 shows the distribution of temperature and current density along the central line of the OCL.At the end of the FC mode (1.99 s), the general trends of temperature and current density decline gradually along the x-axis.The distribution of current density probably results from the concentrations of the reactants.As shown in Figure 8, the O 2 mass fraction declines along the x-axis.The temperature Energies 2016, 9, 1015 is determined by the heat source and heat flow.For the homogenous and isotropic CL, a relatively high temperature is observed where the heat source is large.The heat source in the OCL consists of heat attributed to activation polarizations, Joule heat, and reaction heat, all of which are positively correlated with current density.Therefore temperature is depressed along the x-axis.Temperature and current density show remarkable reduction during 0.02 s from 1.99 s to 2.01 s.The causes of temperature reduction are varied.The main reasons are the endothermic electrochemical reaction heat source and low current density under the EC mode.However, the current density distribution increases slightly along the x-axis at 2.01 s and decreases at 2.5 s.The main reason for these trends is the change of H 2 O concentration, which has an opposite distribution trend compared with O 2 .From the O 2 mass fraction shown in Figure 8, the concentration of H 2 O, which is reactant for the reaction of the EC mode, has a distribution trend similar to current density along the x-axis.The temperature at the OCL exhibits an evident decrease in the process of operation mode switching.The temperature difference along the x-axis is less than 0.1 K, except at the left end.The current density, temperature, and mass fraction of the reactant influence each other to some extent.Figure 10 displays the temperature profiles of the cell along line B-B at different times.Figure 11 shows the detailed temperature distribution of the MEA.As shown in Figure 10, temperature in the hydrogen side is lower than that in the oxygen side, which can be explained by the following three factors: (1) The heat source in OCL accounts for more than 90% of all heat sources according to Table 3; (2) the thickness of OGDL is double that of HGDL, which causes double Joule heat; and (3) the heat conductivity of oxygen is smaller than that of hydrogen.In spite of the endothermic reaction in the EC mode, the temperature in the hydrogen side is lower than that in the oxygen side during the limited time of this mode because of the heat attributed to activation polarizations, Joule heat, and heat accumulation in the FC mode.The maximum value of temperature appears at the OCL and    Figure 10 displays the temperature profiles of the cell along line B-B at different times.Figure 11 shows the detailed temperature distribution of the MEA.As shown in Figure 10, temperature in the hydrogen side is lower than that in the oxygen side, which can be explained by the following three factors: (1) The heat source in OCL accounts for more than 90% of all heat sources according to Table 3; (2) the thickness of OGDL is double that of HGDL, which causes double Joule heat; and (3) the heat conductivity of oxygen is smaller than that of hydrogen.In spite of the endothermic reaction in the EC mode, the temperature in the hydrogen side is lower than that in the oxygen side during the limited time of this mode because of the heat attributed to activation polarizations, Joule heat,    11 shows the detailed temperature distribution of the MEA.As shown in Figure 10, temperature in the hydrogen side is lower than that in the oxygen side, which can be explained by the following three factors: (1) The heat source in OCL accounts for more than 90% of all heat sources according to Table 3; (2) the thickness of OGDL is double that of HGDL, which causes double Joule heat; and (3) the heat conductivity of oxygen is smaller than that of hydrogen.In spite of the endothermic reaction in the EC mode, the temperature in the hydrogen side is lower than that in the oxygen side during the limited time of this mode because of the heat attributed to activation polarizations, Joule heat, and heat accumulation in the FC mode.The maximum value of temperature appears at the OCL and increases gradually under the FC mode.During the 0.02 s from 1.99 s to 2.01 s, the temperature of the oxygen electrode decreases sharply, but the temperature of other areas almost remains constant.The main reason for these results is that the OCL is the place where heat source exhibits the most significant changes.Furthermore, the low conductivity of the membrane hinders heat transfer to the hydrogen side.After the EC process lasting 0.5 s, at 2.5 s, the maximum of temperature appears at the membrane on account of the largest value of heat source according to Table 3.The temperature of the entire B-B line decreases to different degrees.Similar phenomena and conclusion can be observed and drawn from Figure 12, which shows transient variations of the entire cell temperature.Under the FC mode, the temperature in the OCL is significantly higher than that in other areas.The temperature difference decreases when the mode is switched to the EC mode.The maximum temperature difference is less than 1.5 K.By the time the EC reaches a more stable state, the cell temperature is slightly higher than the operating temperature, and the highest temperature appears at the membrane.In addition, the cell temperature distribution trends along the x-axis are similar at different times.The temperature gradient near the inlet is relatively large and exhibit a small temperature change along the x-axis behind x = 0.05 m. Figure 13 shows the transient responses of average temperature in different layers.The OCL, which generates the most heat under the FC mode but absorbs heat under the EC mode, has the strongest temperature variation.In the FC mode, the temperature of each layer increases with time because of heat accumulation.Then the temperature variation becomes small and decreases rapidly when the cell mode is switched from FC mode to EC mode at 2 s.The temperatures of the CLs are close to that of the membrane under the EC mode.Eventually, the temperature difference of entire Figure 13 shows the transient responses of average temperature in different layers.The OCL, which generates the most heat under the FC mode but absorbs heat under the EC mode, has the strongest temperature variation.In the FC mode, the temperature of each layer increases with time because of heat accumulation.Then the temperature variation becomes small and decreases rapidly when the cell mode is switched from FC mode to EC mode at 2 s.The temperatures of the CLs are close to that of the membrane under the EC mode.Eventually, the temperature difference of entire cell is small because of the endothermic reaction and the small current density in the EC mode.
Table 3 lists the heat sources and its percentages in each layer before and after mode switching.Under the FC mode, heat attributed to activation polarizations and reaction heat are the main sources of cell total heat production, and Joule heat contributes less than 10% of all heat.When the mode is switched to EC, cell total heat production comes mainly from Joule heat in the membrane, and the value of the heat source in the OCL is negative because of the endothermic reaction.However, the cell is still exothermic, although the total heat is small under the EC mode.

Conclusions
In this study, a transient, non-isothermal, and multi-physics model for a PEM URFC is developed, in which the electrochemical reaction is also coupled.The model is validated, and results are consistent with the previously published experimental results.Numerical simulations are conducted to investigate the transient responses of heat and mass transfer as well as current density in a URFC when the operating mode is switched from FC mode to EC mode.The results, which considered the temperature, gas mass fractions, current density, and heat source, are shown and discussed.The following conclusions can be obtained: (1) Mode switching has a significant effect on heat and mass transfer in PEM URFCs.The dynamic responses of temperature, mass fraction and current density are different.The response time of temperature is longer than that of current density and gas mass fractions.(2) The mass fraction of H2O along the direction from gas flow channel to CL and the gas flow direction changes into a decreasing trend in the EC mode from an increasing trend in the FC mode.H2 and O2 reveal an opposite mass fraction distribution trend compared with H2O.The gradient of gas mass fractions in the oxygen side is larger than that in the hydrogen side along x-and y-axes.(3) Cell temperature in the hydrogen side is lower than that in the oxygen side under the FC and EC modes.The temperature difference of the entire cell is less than 1.5 K.The highest temperature appears at the OCL under the FC mode.Under the EC mode, the highest temperature initially exists at the OCL and then at the membrane.(4) The OCL, which generates the most heat under the FC mode but absorbs heat under the EC mode, has the strongest temperature variation when the cell operates from the FC mode to the EC mode.(5) The cell is exothermic under the FC and EC modes.The heat attributed to activation polarizations and reaction heat in the OCL occupies a dominant part of heat production under

Conclusions
In this study, a transient, non-isothermal, and multi-physics model for a PEM URFC is developed, in which the electrochemical reaction is also coupled.The model is validated, and results are consistent with the previously published experimental results.Numerical simulations are conducted to investigate the transient responses of heat and mass transfer as well as current density in a URFC when the operating mode is switched from FC mode to EC mode.The results, which considered the temperature, gas mass fractions, current density, and heat source, are shown and discussed.The following conclusions can be obtained: (1) Mode switching has a significant effect on heat and mass transfer in PEM URFCs.The dynamic responses of temperature, mass fraction and current density are different.The response time of temperature is longer than that of current density and gas mass fractions.(2) The mass fraction of H 2 O along the direction from gas flow channel to CL and the gas flow direction changes into a decreasing trend in the EC mode from an increasing trend in the FC mode.H 2 and O 2 reveal an opposite mass fraction distribution trend compared with H 2 O.The gradient of gas mass fractions in the oxygen side is larger than that in the hydrogen side along xand y-axes.(3) Cell temperature in the hydrogen side is lower than that in the oxygen side under the FC and EC modes.The temperature difference of the entire cell is less than 1.5 K.The highest temperature appears at the OCL under the FC mode.Under the EC mode, the highest temperature initially exists at the OCL and then at the membrane.
(4) The OCL, which generates the most heat under the FC mode but absorbs heat under the EC mode, has the strongest temperature variation when the cell operates from the FC mode to the EC mode.(5) The cell is exothermic under the FC and EC modes.The heat attributed to activation polarizations and reaction heat in the OCL occupies a dominant part of heat production under the FC mode.Meanwhile, in the EC mode, Joule heat in the membrane mainly contributes to heat production.

Figure 2 .
Figure 2. Schematic of the rib and the gas flow channel.

Figure 2 .
Figure 2. Schematic of the rib and the gas flow channel.

Figure 3 .
Figure 3. Element independence test for the present model: (a) temperature distribution of the central line of the membrane with an operating voltage of 0.6 V; and (b) temperature distribution of the central line of the OCL with an operating voltage of 1.8 V.

Figure 3 .
Figure 3. Element independence test for the present model: (a) temperature distribution of the central line of the membrane with an operating voltage of 0.6 V; and (b) temperature distribution of the central line of the OCL with an operating voltage of 1.8 V.

Figure 5 .
Figure 5. Operating voltage setting under mode switching.

Figure 6
Figure 6 displays the transient responses of current density, O2 mass fraction and temperature at the central point of the OCL.At the beginning of the FC mode, the O2 mass fraction declines sharply,

Figure 5 .
Figure 5. Operating voltage setting under mode switching.

Figure 6
Figure 6 displays the transient responses of current density, O 2 mass fraction and temperature at the central point of the OCL.At the beginning of the FC mode, the O 2 mass fraction declines sharply, whereas the current density and temperature increase gradually because of the consumption of O 2 and the generation of heat.The O 2 mass fraction and current density reach a steady state at approximately 1 s.At 2 s the temperature continues to increase.The current density decreases to 2.77 mA•cm −2 at 2 s because of mode switching.Then, the current density rapidly increases to approximately 15 mA•cm −2 corresponding to the EC mode.At the same time, O 2 mass fraction and temperature at that point change dramatically.These changes can be explained as follows: The accumulated water, which is produced in the FC mode, is the reactant in the EC mode.The concentration of H 2 O reduces rapidly at the beginning of the EC mode, which leads to a slight decrease in the current density.The O 2 mass fraction increases sharply because of the consumption of H 2 O and the production of O 2 when the mode is switched from FC mode to EC mode.Temperature decreases rapidly at 2 s mainly because the exothermic reaction in the FC mode converts into the endothermic reaction in the EC mode.During the subsequent 3 s, the temperature gradient diminishes gradually.Temperature does not reach a steady state during the process.Current density, O 2 mass fraction and temperature in the OCL suddenly change under operation mode switching.The response time of temperature is longer than that of current density and O 2 mass fraction.The water mass fractions along y-axis with varying times are displayed in Figure7.The solid lines represent the temperature distribution along A-A line and the dotted lines represent the temperature distribution along the C-C line.Under the FC mode, the water mass fraction along the C-C line is evidently higher than that along the A-A line.The main reasons are the consumption of reactants and generation of H 2 O in oxygen side.Therefore, the H 2 O mass fraction increases gradually along the gas flow direction.Moreover, the consumption and the generation of gases occur in the CLs.H 2 O mass fraction increases gradually from gas flow channels to CLs because of gas diffusion.However, the distribution of H 2 O mass fraction under the EC mode is different from that under the FC mode because of the opposite electrochemical reactions.The mass fraction gradient in the hydrogen side is smaller compared with that in the oxygen side because of the generation of water and the small molecular weight of hydrogen.The water mass fraction at 0.5 and 1.99 s did not change at the A-A line.Meanwhile, at the C-C line, the water mass fraction at 1.99 s is slightly larger than that at 0.05 s, which illustrates that the mass fraction near the inlet of the gas flow channel reaches a steady state faster.The H 2 O mass fraction decreases in different degrees within 0.02 s from t = 1.99 s to

Figure 6 .
Figure 6.Time-dependent evolution of relevant parameters in the central point of the OCL.

Figure 6 .Figure 7 .
Figure 6.Time-dependent evolution of relevant parameters in the central point of the OCL.Energies 2016, 9, 1015 11 of 19

Figure 7 .
Figure 7. Water mass fraction distribution along lines A-A and C-C: (a) oxygen side; and (b) hydrogen side.

Figure 8
Figure 8 displays the distribution of H 2 and O 2 mass fractions along the central lines of the CLs.The results indicate that both O 2 and H 2 mass fractions are depressed distinctly along the flow direction at 1.99 s.The O 2 and H 2 mass fractions increase when the mode is switched to the EC mode.Eventually, the O 2 and H 2 mass fractions increase along the flow direction under a more stable EC mode at 2.5 s, and the magnitudes of O 2 and H 2 mass fractions are greater than those at the inlets of gas flow channels because of the production of O 2 and H 2 in the CLs.As with the H 2 O mass fraction, the gradient of the O 2 mass fraction is significantly larger than that of the H 2 mass fraction.Figure9shows the distribution of temperature and current density along the central line of the OCL.At the end of the FC mode (1.99 s), the general trends of temperature and current density decline gradually along the x-axis.The distribution of current density probably results from the concentrations of the reactants.As shown in Figure8, the O 2 mass fraction declines along the x-axis.The temperature

Figure 8 .
Figure 8. H2 and O2 mass fraction distribution at different time.

Figure 9 .
Figure 9.The distribution of temperature and current density.

2 )Figure 8 .
Figure 8. H 2 and O 2 mass fraction distribution at different time.

Figure 8 .
Figure 8. H2 and O2 mass fraction distribution at different time.

Figure 9 .
Figure 9.The distribution of temperature and current density.

2 )Figure 9 .
Figure 9.The distribution of temperature and current density.

Figure 10
Figure10displays the temperature profiles of the cell along line B-B at different times.Figure11shows the detailed temperature distribution of the MEA.As shown in Figure10, temperature in the hydrogen side is lower than that in the oxygen side, which can be explained by the following three Figure10displays the temperature profiles of the cell along line B-B at different times.Figure11shows the detailed temperature distribution of the MEA.As shown in Figure10, temperature in the hydrogen side is lower than that in the oxygen side, which can be explained by the following three

Figure 10 .
Figure 10.Temperature profiles along the B-B plane of the entire cell.

Figure 11 .
Figure 11.Temperature profiles along the B-B plane of the membrane electrode assembly (MEA).

Figure 10 .Figure 10 .
Figure 10.Temperature profiles along the B-B plane of the entire cell.MEA: membrane electrode assembly.

Figure 11 .
Figure 11.Temperature profiles along the B-B plane of the membrane electrode assembly (MEA).

Figure 11 .
Figure 11.Temperature profiles along the B-B plane of the MEA.

Figure 13 .
Figure 13.Transient responses of average temperature in different layers.

Figure 13 .
Figure 13.Transient responses of average temperature in different layers.

Table 1 .
Physical parameters and basic conditions.CL: Catalyst layer; and GDL: gas diffuse layer. 2

Table 3 .
Heat sources of electrodes and percentages in different times.

Table 3 .
Heat sources of electrodes and percentages in different times.

Table 3 .
Heat sources of electrodes and percentages in different times.