Investigation of performance heterogeneity of PEMFC stack based on 1+1D and flow distribution models

Energy Conversion and Management 207 (2020) 112502

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Investigation of performance heterogeneity of PEMFC stack based on 1+1D and flow distribution models

T

Zirong Yang, Kui Jiao, Zhi Liu, Yan Yin, Qing Du

⁎

State Key Laboratory of Engines, Tianjin University, 135 Yaguan Road, Tianjin 300350, China

ARTICLE INFO

ABSTRACT

Keywords: PEMFC stack Flow distribution Performance heterogeneity Voltage overestimation Manifold geometric parameter

Flow distributions among proton exchange membrane fuel cell (PEMFC) stacks remain an important issue owing to the great influences on stack performances. To simultaneously consider flow distributions as well as reactions, phase changes, and transport processes inside every fuel cell, a comprehensive stack model is developed based on the integration of a 1+1 dimensional multiphase stack sub-model and a flow distribution sub-model. Frictional and local pressure drops due to diverging, converging, and bending configurations are calculated. After rigorous model validation, differences between the uniform flow assumption and the authentically nonuniform distribution are quantitatively investigated, including the distribution of voltage, mass flow rate, temperature, reactant concentration, and other parameters. Results show that the uniform assumption not only overestimates the stack output performance but also underestimates the cell voltage variations. Besides, the uniform assumption may lead to higher predictions of the overall stack temperature and lower predictions of temperature variations among different fuel cells. Even though the total amount of air seems abundant, it is still possible for some fuel cells to suffer from the local reactant starvation. Therefore, a higher cathode stoichiometry is preferable since it increases the inlet mass flow rate for middle cells. Increasing the inlet pressure contributes negligibly to the uniformity of reactant distribution, but it improves the stack performance. A larger manifold cross-sectional area leads to more uniform reactant distributions among the stack and less local current density variations inside single cells.

1. Introduction During the past decades, there is growing concerns about energy and environment around the world. In addition to the limited storage, fossil fuels may cause environmental problems in actual applications such as automobile vehicles. Therefore, researchers are paying more and more attention to the environmentally friendly vehicles nowadays (e.g. hybrid electric vehicles, battery electric vehicles, and fuel cell vehicles) [1–4]. Among these vehicles, proton exchange membrane fuel cell (PEMFC) has the advantages of long cruising distance, low operating temperature, and zero emissions [5–7]. Several motor corporations have already issued their fuel cell vehicles, such as Toyota Mirai [8], Honda Clarity [9], and Hyundai Nexo [10]. To enhance the output power, numerous PEMFCs are typically stacked together. For Mirai, 370 single fuel cells are stacked in single-line, which achieves the maximum output at 114 kW [8]. During the scaling up process, the uneven flow distribution remains an important issue because it not only affects the overall stack output power, but also significantly influences the stack reliability and durability [11–13]. Therefore, it is necessary to

⁎

comprehensively investigate the mechanisms and influences of flow distributions within the whole stack. After being supplied into the stack inlet manifold, reactants are subsequently distributed into each individual fuel cell where the hydrogen–oxygen reaction occurs [14]. Due to the pressure losses along with reactants movement, the mass flow rates into each fuel cell are typically different, leading to the performance heterogeneity such as output voltage, temperature, and other parameters. To investigate the uneven flow distribution phenomenon and its influences, many experimental and numerical studies have been conducted. Kim and Hong [15] evaluated the effects of operating conditions on the performance of a PEMFC stack with 10 cm2 active area and 10 single cells. The output voltage of each individual fuel cell was measured. It was found that the voltage variations among different cells could reach 5% under the tested conditions. Yi et al. [16] used the same amorphous carboncoated 304 stainless steel as bipolar plates in a single cell and a 100 Wclass short stack. The results showed that the peak power density was 1150.6 mW cm−2 for the single cell and 1040.1 mW cm−2 for the stack, indicating that the maximum power density deteriorated after scaling

Corresponding author. E-mail address: [email protected] (Q. Du).

https://doi.org/10.1016/j.enconman.2020.112502 Received 4 October 2019; Received in revised form 11 December 2019; Accepted 12 January 2020 0196-8904/ © 2020 Published by Elsevier Ltd.

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

Nomenclature a A ASR c Cp D Dh EW F f G h I l k K m M p p R Re s S ST t t T u V X Y

κ

water activity effective geometric area (m2) area specific resistance (Ω m2) mole concentration of gas species (mol m−3) specific heat capacity (J kg−1 K−1) diffusivity (m2 s−1) hydraulic diameter (m) equivalent weight of membrane (kg kmol−1) Faraday’s constant (C mol−1) function Gibbs free energy (J mol−1) heat transfer coefficient (W m−2 K−1) current density (A m−2) length (m) thermal conductivity (W m−1 K−1) permeability (m2) mass flow rate (kg s−1) molecular weight (g mol−1) pressure (Pa) pressure drop (Pa) universal gas constant (J mol−1 K−1) Reynolds number volume fraction source terms (kmol m−3 s−1, kg m−3 s−1, W m−3), entropy (J mol−1 K−1) stoichiometry time (s) time step size (s) temperature (K) velocity (m s−1) voltage (V) mole fraction of gas species mass fraction of gas species

Subscripts and superscripts a act atm bend bot BP c CH cir CL con conc cool div eff EOD eq FC fr g GDL H2 l lq MEM mw MPL N2 O2 ohmic out pc per react rect ref sat surr vp m-l m-v v-l

Greek letters ε ζ λ ξ μ ρ ω δ θ

conductivity (S m−1) friction factor

porosity water transfer rate (s−1), pressure drop coefficient membrane water content stoichiometry ratio dynamic viscosity (kg m−1 s−1) density (kg m−3) volume fraction of ionomer thickness (m) contact angle (°)

up. Pei et al. [17] measured the temperature distribution characteristic in a 200 cm2 stack with 46 cells by embedding micro-thermocouples into the cathode plate of four marked cells (cell 1, 11, 23, and 46). Nine thermocouples were placed at different positions for the aforementioned cells. It was found that the highest temperature difference among the marked cells could reach approximately 11 °C. Besides, the temperature variation inside the single cell could reach around 8 °C. To simultaneously measure the voltage and temperature distributions, Devrim et al. [18] conducted a research with an air-cooled 500 W stack consisting of 24 fuel cells. The highest temperature difference among the individual cells was around 9 °C when the stack was operated at 65 °C. It should be noted that the performance heterogeneity mainly results from the non-uniform reactant distributions within the stack. However, the mass flow rate of each fuel cell is hard to experimentally measure since the mass flow sensor could not be placed in the limited

anode activation atmosphere bending bottom bipolar plate cathode, capillary flow channel circular catalyst layer converging concentration coolant diverging effective electro-osmotic drag equilibrium fuel cell frictional gas phase gas diffusion layer hydrogen liquid phase liquid water membrane membrane water micro-porous layer nitrogen oxygen ohmic output phase change permeation reaction rectangular reference state saturation surroundings water vapor membrane water to liquid water membrane water to water vapor water vapor to liquid water

space of stack manifolds or fuel cell channel inlet regions. Therefore, numerical models have been developed to quantitatively depict the reactant and coolant distributions. Baschuk and Li [19] investigated the distributions of pressure and reactants based on a hydraulic network approach. The frictional pressure drops in stack manifolds and fuel cell channels were calculated. Mass conservation analysis for hydrogen, air, and water vapor was conducted for the zero-dimensional steady-state stack model. It was found that the power density of a stack is lower than that of a single cell, which was resulted from unequal reactant distributions. However, local pressure drops owing to the bending, converging, and diverging configurations were not considered. To supplement the model, the local pressure drops and heat transport processes were added in the work of Park and Li [20]. The number of heat transfer unit (NTU) method was adopted to calculate the temperature of every fuel cell. To investigate the thermal and electrical behaviors of a water-cooled PEMFC stack, Cozzolino et al. [21] developed a one2

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

dimensional electrochemical model, which was integrated with thermal analysis based on mass and energy balance. The experimental tests were also conducted to validate and modify the stack model. However, the pressure drops and non-uniform reactant distribution were not considered in the study. Similar researches were improved by Salva et al. [22], which calculated the frictional pressure drops through the mass and pressure balance. Amirfazli et al. [23,24] established an analytical model to study the temperature distributions and coolant channel designs with the U-type and Z-type stacks consisting of 65 fuel cells. The concept of index of uniform temperature (IUT) was used to evaluate the degree of temperature uniformity. However, the mass transport processes and electrochemical reactions inside fuel cells were not involved. There were several three-dimensional (3D) stack models developed by commercial computational fluid dynamics (CFD) software [25–28]. Macedo-Valencia et al. [26] established a single-phase 3D stack model with five 31.2 cm2 fuel cells. The reactants concentration, local current density, and temperature distributions among the stack were presented. Luo et al. [27] investigated the cold start performance of a 3-cell stack with the transient multiphase model. However, the inlet mass flow rates of all fuel cells were assumed identical, and the stack manifolds were not taken into consideration. To improve the modeling studies, Lim et al. [28] developed a steady-state single-phase stack model, and the manifold structures were further included in the computational domains. The performances of single-cell, dual-cell, quad-cell, and hexa-cell were comparatively studied. Despite aforementioned studies, the researches which simultaneously consider nonuniform flow distributions, electrochemical reactions, and phase changes are rarely presented. Note that the phase changes among water vapor, membrane water, and liquid water have great influences on the components of mixture gases, which affect the pressure drops inside single cells and correspondingly influence the reactant distributions within the whole stack. Besides, comparison between the uniform flow assumption and the non-uniform flow distribution has not been quantitatively investigated, which remains an important issue for stack modeling studies. In the study, a flow distribution sub-model is integrated with a 1+1D multiphase stack sub-model to develop the comprehensive stack model, which is capable of simultaneously demonstrating non-uniform reactant and coolant distributions among the stack as well as electrochemical reactions, phase changes, and transport processes inside every fuel cell. Detailed correlations for frictional and local pressure drop coefficients are presented. The established stack model is rigorously validated against experimental data. Differences between the uniform flow assumption and the authentically non-uniform distribution are quantitatively investigated, including the distribution of cell voltage, mass flow rate, temperature, reactant concentration, and other parameters, which are rarely presented in literature but of great significance for stack researches. To alleviate non-uniform flow distributions and thereby enhance stack performances, the effects of stack operating conditions (e.g. stoichiometry, pressure) and manifold geometric

parameters are further studied, which gives valuable suggestions for the stack performance improvement. 2. Model development The schematic diagram of a proton exchange membrane fuel cell (PEMFC) stack with the Z-type configuration is shown in Fig. 1. The stack configuration is a U-type if the flow direction in the outlet manifold is designed contrary to the inlet manifold. To investigate the uneven flow distributions among different fuel cells, a flow distribution sub-model is developed based on the hydraulic network approach. As regards heat and mass transfer processes inside every fuel cell, a 1+1 dimensional multiphase stack sub-model is developed. The two submodels are subsequently integrated to develop the comprehensive stack model. Detailed methodologies about the two sub-models are explained in Subsection 2.1 and 2.2. The calculation procedure, initial and boundary conditions are explained in Subsection 2.3. 2.1. Flow distribution sub-model During the motion of fluid, two types of pressure losses should be considered, including the frictional pressure loss and the local pressure loss. The frictional losses result from the friction between fluids and walls. The local pressure losses are usually caused by the converging, diverging and bending of the flow. It should be noted that the two types of pressure losses are arbitrarily divided for the sake of calculation while the total pressure losses are inseparable [29].

ptotal = pfr + plocal 1 u2 2 fr 1 p local = 2 local

pfr =

u2

(1)

where ptotal (Pa) is the total pressure loss, pfr (Pa) is the frictional pressure loss, plocal (Pa) is the local pressure loss, fr , local are the pressure loss coefficients, (kg m−3) is the density, and u (m s−1) is the velocity. To determine the total pressure losses, it is necessary to obtain the pressure loss coefficients which are related to flow regimes, mass flow rates, and structural designs. The frictional loss coefficient is defined as [29]: fr

=

l Dh

(2)

where is the friction factor which depends on the Reynolds number and the wall roughness, l (m) is the length, and D h (m) is the hydraulic diameter. For circular tube with smooth walls under stabilized flow conditions, the friction factor is calculated as [29]:

Fig. 1. Schematic diagram of a Z-type PEMFC stack. 3

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

uD h µ

Re =

(Re < 2000)

f1 (Re )

=

f1,div2 = p00 + p10 x + p01y + p20 x 2 + p11xy + p02 y 2 + p30 x 3 + p21x 2y + p12 xy 2 + p03 y 3 + p40 x 4

64 Re cir

x = A2 /A1 , y = m2 /m1

(3)

Re < 4000)

(2000

0.3164 Re 0.25

(4000

(1.8 lg Re

1.64)

+ p31x 3y + p22 x 2y 2 + p13 xy 3 + p04 y 4 + p50 x 5 + p41x 4y + p32 x 3y 2

2

Re <

+ p23 x 2y 3 + p14 xy 4 + p05 y 5 f1,div3 = 0.969(1

105)

g

=

RT

Yi Mi

i

1

=

RT

MiX i

µg = i

j

ij

=

1.0 = 0.9(1 0.55

(6)

The function f1 in Eq. (4) is obtained by polynomial curve fitting methods based on the original reference [29].

f1 (Re ) =

12Re 3

1.49 × 10

1.561 × 10

5Re

+ 0.0331

rect = f cali

)

y2

(11)

(x < 0.35) y ) (x > 0.35, y < 0.4) (x > 0.35, y > 0.4)

(12)

(13)

+ 1.036 × 10 8Re 2 (2000

Re < 4000)

(7)

Table 1 Schematic diagram and calculation of local pressure losses [29]. Type

Correlations

p2div = div 1,2

1 div 2 1,2

div = f1,2

(

u12 ,

A2 m2 , A1 m1

p3div =

),

div 1,3

1 div 2 1,3

div = f1,3

u12

( ) m3 m1

(8)

cir

( ) 0.0781( ) 0.707

y)2

x = A2 A1 , z = a b

rect is the calibration coefficient which typically depends on the where f cali Reynolds number and cross-sectional designs. The calibration coefficient is also obtained by curve fitting methods based on the presented figures and tables in the reference [29].

rect f cali =

2(1

f1bend = p00 + p10 z + p01x + p20 z 2 + p11zx + p02 x 2 + p30 z 3 + p21z 2x + p12 zx 2 + p03 x 3

The sum of squares due to error (SSE) is 7.94 × 10−-8, and the coefficient of determination (R-square) is 0.9994, indicating good fitting accuracy. If the cross-section of stack manifold differs from circular tubes (e.g. rectangular, elliptical, trapezoid), the structural calibration is necessary for the corresponding friction factor. For rectangular conduit, the actual friction factor is calculated as [29]: rect

y 2 x

For the first fuel cell near the inlet manifold or the last fuel cell near the outlet manifold, the local pressure drop owing to the bending configuration should also be taken into consideration. For rectangular section elbows with 90° sharp corners, the friction factor f1bend is calculated as

[1 + (µi µ j )0.5 (Mj Mi )0.25] 2 [8(1 + Mi Mj)] 0.5

(1 + ( )

f1,con 3 = 1.55y

where pg (Pa) is the gas pressure, R (J mol−1 K−1) is the universal gas constant, T (K) is the temperature, Yi is the mass fraction of gas species i, Xi is the mole fraction, and Mi (g mol−1) is the molecular weight. The mixture dynamic viscosity is calculated as [37]:

Xi µi , Xj ij

−1

x = A2 A1 , y = m2 m1

(5)

i

y ) + 0.728

0.282(1

where A (m ) is the cross-sectional area, and m (kg s ) is the mass flow rate. The multinomial coefficients of curving fitting results are listed in Table 2. The R-square of the two functions are 0.9988 and 0.9683. For the converging configuration, the corresponding functions are calculated as [29]:

f1,con 2 =

pg

y) 2

2

where µ (kg m−1 s−1) is the dynamic viscosity. It should be noted that different types of gas are involved in stack manifolds and fuel cell channels, including hydrogen, oxygen, nitrogen, and water vapor. Therefore, the density and dynamic viscosity refer to the corresponding mixed gases. Based on the ideal gas law, the mixture density is calculated as

pg

0.709(1

(10)

(4)

(Re > 4000)

y) 3

b0 a0

3

b0 a0

3

( ) + 0.203( ) + 1.96

b0 a0

b0 a0

2

con 1,2

( ) + 1.50 (Re < 2000) 0.224( ) + 1.10 (Re > 2000)

1.87 2

p2con =

1 con 2 1,2

con = f1,2

(

u12 ,

A2 m2 , A1 m1

p3con =

),

con 1,3

1 con 2 1,3

con = f1,3

u12

( ) m3 m1

b0 a0

b0 a0

(9) where b0 (m) is the cross-sectional width, and a0 (m) is the cross-sectional length. The R-square of the two functions are 0.9998 and 0.9983, respectively. The local pressure loss coefficients due to diverging, converging, and bending configurations are related to the cross-sectional area and the mass flow rate. The schematic diagram and calculation of local pressure losses are summarized in Table 1. The local pressure loss coefficients are determined based on the original tables and figures [29]. It should be noted that the cross-sectional area of the straight passage is identical, namely, A1 is equal to A3 as shown in Table 1. For the diverging configuration, the functions f1,div2 , are calculated as

p1bend = bend 1

1 bend 2 1

= f1bend

(

u2

a A2 , b A1

) (a, b refers to the

length and width of the inlet cross-section A1 )

4

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

Table 2 Multinomial coefficients in the diverging configuration. Parameter

Value

Parameter

Value

Parameter

Value

Parameter

Value

p00 p10 p01 p20 p11 p02

8.747 −139.7 27.15 825.5 −172.1 5.251

p30 p21 p12 p03

−2132 388.8 77.77 −41.55

p40 p31 p22 p13 p04

2485 −475.6 17.53 −119.7 72.27

p50 p41 p32 p23 p14 p05

−1046 225.1 −58.86 71.95 8.671 −25

The multinomial coefficients in the bending configuration are listed in Table 3. The SSE is 0.02034 and the R-square is 0.993. It should be noted that the fitting results apply to the conditions with smooth walls and the Reynolds number larger than 2 × 105. However, appropriate coefficients are not available if the Reynolds number is out of the above range. In stack manifolds, the gas velocity can be obtained by the corresponding mass flow rate, density, and cross-sectional area.

liquid pressure because the distinct hydrophobicity may cause liquid water jump phenomenon at the interface of neighboring layers [36]. The water saturation in porous layers are subsequently calculated by Leverett-J function once the capillary pressure is obtained. Detailed transport properties can be found in our previous studies [35].

m u= A

where s lq is the liquid water saturation, pl (Pa) is the liquid pressure, and Slq (kg m−3 s−1) is the liquid water source term. The hydraulic permeation caused by pressure differences between the anode CL and the cathode CL is also considered [37].

(

(14)

EW

t

=

MEM

EW

·(

1.5D

mw

mw )

+ Smw

Sv

pc = pg

Sreact Sm

Sm v

v

Sm

l

Sm l + SEOD (in cathode CL) SEOD (in anode CL)

pl + Slq

µlq

Sper (in cathode CL)

l l

(17)

+ Sper

(in anode CL) (18)

(in MPL, GDL)

l

pl

(19)

( ) [1.42(1 cos( ) ( ) [1.42s 0.5

cos( )

pc =

s lq )

K0

0.5

K0

slq ) 2 + 1.26(1

2.12(1

s lq )3]( < 900)

2.12s lq 2 + 1.26s lq 3]

lq

( > 900)

(20) For gas species, the ideal gas law is adopted. The governing equations of gas species are solved in porous layers and flow channels. t t

( (1

s lq ) g Yi ) =

( g Yi ) +

·( g Dieff Yi ) + Si

·( g ug Yi ) =

·(

eff g Di

(in CL, MPL, GDL)

Yi ) + Si

(in CH)

(21)

where Yi is the mass fraction of gas species i (H2, O2, N2, vp), g (kg m−3) is the mixture gas density, and Si (kg m−3 s−1) is the source term. The energy equation is calculated in the whole fuel cell.

t

(( cp )eff fl,sl T ) = I 2ASRCL 3 CL

ST =

eff ·(k fl,sl T ) + ST

+ IVact,a + Spc

I 2ASRCL 3 CL

T S I 2F

(22)

(in CLa)

+ IVact,c + Spc

I 2ASR I 2ASR

(in CLc)

(in MEM, BP)

+ Spc

(in MPL, GDL, CH)

(23)

−3

where T (K) is the temperature, and ST (W m ) is the heat source term, including activational heat, reversible heat, ohmic heat, and latent heat.

(15)

where mw is the membrane water content, Dmw (m2 s−1) is the membrane water diffusivity, and Smw (kmol m−3 s−1) is the source term, which consists of electrochemical reactions, electro-osmotic drag effect, and phase changes (listed in Table 5).

Smw =

lq

Sm l MH2O + Sv

A typical PEMFC stack consists of many single cells to increase the output voltage and power. A single PEMFC is usually composed of 11 layers, including proton exchange membrane, catalyst layer (CL), micro-porous layer (MPL), gas diffusion layer (GDL), flow channel (CH), and bipolar plate (BP). The coolant is designed to flow through cooling channels inside every bipolar plate to improve the temperature distribution uniformity [31]. To depict the parameter variations (e.g. reactant concentration, temperature, local current density, and velocity) from the channel inlet region towards the outlet region, the fuel cell is divided into several nodes in the flow direction. The schematic diagram and detailed explanation about the stack sub-model can be found in our previous studies [32,33]. Table 4 lists the structural properties and operating conditions. The governing equations include membrane water, liquid pressure, gas species, and temperature. The generated water is assumed to be membrane water since it is stated that the liquid water and vapor production assumptions should be cautiously adopted [34]. mw )

Klq

·

Slq = Sm l MH2O + S v

2.2. 1+1 dimensional multiphase stack sub-model

(

=

t

The pressure drops in fuel cell channels are calculated in the similar way as shown in Eq. (1). It should be noted that the pressure drops in fuel cell channels are also influenced by the possible presence of liquid water [30]. However, specific correlations between pressure drops and the gas–water two-phase flow seem unavailable in literature. Therefore, only single-phase pressure drops are calculated in the study. Note that the gas velocity varies in the flow direction owing to reactants consumption and water vapor generation. Besides, the mixture gas density and viscosity are also changing. Therefore, the density, viscosity, and velocity for the pressure drop calculation are treated as the average value of the inlet status and the outlet status.

MEM

lq s lq )

Table 3 Multinomial coefficients in the bending configuration.

(16)

In porous layers, the transport process is assumed to be dominated by diffusion, and the previous studies also use the assumption [27,35]. The liquid water volume fraction is solved based on the continuity of 5

Parameter

Value

Parameter

Value

Parameter

Value

p00 p10 p01

0.8257 −0.1395 −0.1586

p20 p11 p02

0.07811 −0.01717 0.1637

p30 p21 p12 p03

−0.04289 −0.004508 0.01079 −0.06378

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

Table 4 Stack structural properties and operating conditions. Parameter

Value

Number of fuel cells Gas manifold cross-sectional length, width Coolant manifold cross-sectional length, width Effective area of single cells Cooling area of single cells Channel length, width, depth, rib width Thickness of MEM, CL, MPL, GDL, BP Density of MEM, CL, MPL, GDL, BP

5 10, 6 mm 10, 5 mm

Specific heat capacity of MEM, CL, MPL, GDL, BP Thermal conductivity of MEM, CL, MPL, GDL, BP Electric conductivity of CL, MPL, GDL, BP Ionomer volume fraction in CL Contact angle of CL, MPL, GDL Porosity of CL, MPL, GDL Expected stack operating temperature Stack manifold inlet pressure Stoichiometry ratio Inlet relative humidity Inlet gas temperature and surrounding temperature Inlet mass flow rate and temperature of coolant Heat transfer coefficient between endplate and environment Heat transfer coefficient between fuel cell and coolant

120 cm2 60 cm2 100, 1.0, 1.0, 1.0 mm 0.025, 0.01, 0.03, 0.2, 2 mm 1980, 1000, 1000, 1000, 1000 kg m−3 833, 3300, 2000, 568, 1580 J kg−1 K−1 0.95, 1.0, 1.0, 1.0, 20 W m−1 K−1 5000, 5000, 5000, 20,000 S m−1 0.4 100°, 120°, 120° 0.3, 0.4, 0.6 60 °C pa = 1.5 atm, pc = 1.5 atm ST a = 1.5, ST c = 2.0 RHa = 1, RHc = 1 60 °C, 25 °C 0.2 kg s−1, 25 °C 100 W m−2 K−1 200 W m−2 K−1

The output voltage of a single fuel cell is calculated based on Tafel equations for simplification [37,38]. The stack output voltage can be obtained by adding all single cell voltages.

Vout = VNernst VNernst

Vact

Vohmic

G S = + (T 2F 2F

Fig. 2. Calculation procedure of the integrated steady-state PEMFC stack model.

(24)

Vconc

The inlet mass flow rate of individual fuel cells is initially guessed which satisfies the mass conservation at each diverging, converging, and bending point. The pressure drops in stack manifolds and fuel cell channels are calculated, including frictional pressure drops and local pressure drops as shown in Fig. 1. It should be noted that the total pressure drop in each closed loop should be equal to zero if the total mass flow rate is correctly distributed into single fuel cells. Hardy cross method is adopted to update the mass flow rate of gases and coolant until the convergence standard is reached [20].

c H2,CLa c O2,CLc RT 1 Tref ) + ln + ln 2F c H2,ref 2 c O2,ref (25)

where Vact , Vohmic , Vconc (V) are the activation loss, ohmic loss, and concentration loss. Detailed calculation equations can be found in our previous studies [35]. 2.3. Calculation procedures

Max { ploop,1 , ..., ploop,N 1 } < 10

As stated previously, the two sub-models are strongly related. The calculation procedure of the integrated stack model is shown in Fig. 2. The total mass flow rate of reactant is obtained based on the required stoichiometry and operating current density.

mHtotal = 2

Nstack IAcell STa MH2 N IA ST M total , mair = stack cell c air 2F × 1000 4F × 0.21 × 1000

3

Pa

(27)

After the flow distribution sub-model is solved, the inlet conditions for each fuel cell are set (e.g. mass flow rate and pressure). The 1+1 dimensional stack sub-model is subsequently calculated. For single cells, the explicit formulation calculation method is adopted. More detailed information is presented in our previous study [14]. In the channel direction inside single cells, the voltage is assumed to be

(26)

Table 5 Phase change functions. Phase change

Correlation

Membrane water and water vapor

Sm

v

=

Membrane water and liquid water

Sm

l

=

Water vapor and liquid water

Sv

l

=

MEM (1 EW MEM (1 m l EW

m v,

s lq )(

m v

v l l v

Phase change rates [5]

Unit

m l

(1 slq

s lq )(

s lq )

= 1.0 ,

m

(pg Xvp

(pg Xvp

eq )

m

eq )

psat )

RT psat )

RT

l v,

6

v l

= 1000

(pg Xvp > psat )

kmol m−3 s−1 kmol m−3 s−1 kg m−3 s−1

(pg Xvp < psat ) s−1

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

identical for all nodes since bipolar plate has good electrical conductivity [14,33]. Therefore, the local current density distribution along the channel direction is obtained through the following nonlinear equations.

data, including the stack polarization curve and the cell voltage distribution [16]. In the experiment, the 100 W-class short stack consists of three fuel cells with 120 cm2 active area each. 30% humidified hydrogen and air are supplied at 80 °C. The anode stoichiometry and cathode stoichiometry are 1.2, 2.0, respectively. The stack temperature is maintained at 60 °C through the water cooling. The overall polarization curve is shown in Fig. 3(c), and it shows good agreement. Meanwhile, the cell voltage variations at 2.1 V, 1.8 V, and 1.5 V are also correspondingly compared as shown in Fig. 3(d). The voltage distribution matches reasonably, and the sum of squares due to error is 0.007. Besides, rigorous validation of the 1+1 dimensional stack submodel has been presented in our previous studies [32,35].

n

(28)

where n is the number of nodes in the channel direction. For different fuel cells among the stack, the total current is identical since the fuel cells are stacked in a single line [8]. The heat loss to the surroundings happens in the end bipolar plates.

Qsurr = hsurr Acell (T

(29)

Tsurr )

The heat exchange between fuel cells and coolant is calculated by the heat convection.

Qcool = hcool A cool (T

3.2. Comparison between uniform assumption and non-uniform distribution As stated previously, reactant and coolant mass flow rate of individual fuel cells are assumed to be equal in numerous stack modeling studies [26,27], which may lead to inaccurate predictions of stack performances. Therefore, differences between the uniform flow assumption and the authentically non-uniform distribution are quantitatively investigated, which are rarely presented in literature. Fig. 4(a) shows the stack voltage under various operating current densities. It is seen that the uniform flow assumption leads to higher performance predictions. The stack voltage difference becomes more significant with the stack current increasing. For current lower than 108 A, the relative difference is calculated to be less than 0.5%. The relative difference reaches 1.2% if the current increases to 156 A, which refers to the current density of 1.3 A cm−2. In other words, assuming the reactants and coolant to be uniformly distributed into individual fuel cells overestimates the stack output performance. For modeling studies, if the correlation of stack voltage difference and current density is obtained, the performance predictions of an authentic stack model can be calculated based on the modification of a simplified stack model with the uniform flow assumption, which largely enhances the calculation efficiency and improves the prediction accuracy. Although the correlation is likely to change with different setting of structural designs and operating conditions, the same method is theoretically applicable. The

(30)

ave Tcool )

ave where Tcool (K) is the average of coolant inlet temperature and outlet temperature. To guarantee the convergence of the explicit calculation method, the time step size for the stack sub-model is set as 10−6 s. To modify the corresponding pressure drops in the flow distribution sub-model, the mixture gas and coolant properties are updated, such as temperature, density, viscosity, and velocity. The flow distribution sub-model is recalculated, and the new inlet mass flow rates are obtained. Iterations continues until the convergence of the integrated stack model is reached.

V , T,

, s lq, Ygas

in in mHin2 , mair , mcool

< 2 × 10

6

(31)

where is the relative error. For the flow distribution sub-model, is calculated as the differences between two iterations divided by the present value. For the stack sub-model, a transient mode rather than a steady-state one is developed. Note that the time step size is set as 10−6 s, and the differences between a time step could be extremely small. Therefore, the relative error is calculated based on the interval of 10−3 s, and it is selected as the maximum value of all nodes in the stack submodel. The influences of convergence criteria are studied, and the results are presented in the model validation subsection. All sub-models are developed with self-written codes in Matlab 2018a.

3.5

8

Relative error (%)

3.0

3. Results and discussion 3.1. Model validation

2.5 2.0 1.5 1.0 0.5

For the integrated proton exchange membrane fuel cell (PEMFC) stack model, the influences of node number in the fuel cell flow direction and convergence criteria are firstly investigated. As shown in Fig. 3(a), the output voltage, inlet mass flow rate, and outlet gas velocity are compared when the node number increases from 1 to 10. Note that the relative error is calculated as the maximum value of all nodes in the stack, which consists of 5 fuel cells. Compared with numerical results in the case of 10 nodes, the relative error falls below 0.5% if the node number reaches 4, and it drops below 0.25% if the node number is increased to 5. Further increasing the node number from 5 to 10 leads to negligible changes. Therefore, the node number in fuel cell flow direction is set as 5 in the study. The influences of convergence criteria are shown in Fig. 3(b). The output voltage, inlet mass flow rate, oxygen concentration, local current density, temperature, and water volume fraction are compared. It is found that the relative error drops below 0.5% if the criteria is set at 2 × 10−6. Decreasing the criteria to 1 × 10−6 requires around 200-hour calculation time, which is almost twice the time than that of the criteria 2 × 10−6. Therefore, the convergence criteria are set as 2 × 10−6. To validate the integrated stack model, simulation results are compared with the experimental

Output voltage Anode inlet flow rate Oxygen concentration Current density Temperature Water volume fraction

7

Output voltage Anode inlet mass flow rate Cathode inlet mass flow rate Anode outlet gas velocity Cathode outlet gas velocity

Relative error (%)

Ii = n·I i

6 5 4 3 2 1

0.0

0 1

10

2 3 4 5 6 7 8 9 10 Node number in fuel cell flow direction

9

8

0.9

Stack voltage (V)

6

5

4

3

2

1

(b)

3.0

0.8 Experimental data Simulation results

2.5

7

Convergency criteria ( x10-6 )

(a)

Cell voltage (V)

f (Ii ) = Vout,

2.0

1.5

Exp. 2.1 V Exp. 1.8 V Exp. 1.5 V

Sim. 2.1 V Sim. 1.8 V Sim. 1.5 V

The sum of squares due to error: 0.007

0.7 0.6 0.5 0.4 0.3

1.0

0

20

40

60

Current (A)

(c)

80

100

0.2

1

2

Cell number

3

(d)

Fig. 3. Comparison between simulation results and experimental data. (a) Effects of node number in the fuel cell channel direction. (b) Effects of convergence criteria. (c) Stack polarization curve [16]. (d) Cell voltage variations at different stack voltages [16]. 7

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al. 4.3

0.84

4.2

Fuel cell voltage (V)

Stack voltage (V)

0.82

Uniform flow assumption Non-uniform distribution

4.1 4.0 3.9 3.8 3.7 3.6

0.78

0.6 A cm-2 1.0 A cm-2 1.3 A cm-2

Uniform Uniform Uniform

Distributed Distributed Distributed

0.76 0.74 0.72 0.70

3.5 3.4

0.80

48

60

72

84

0.68

96 108 120 132 144 156 168 Stack current (A)

1

2

3

4

5

Fuel cell number

(a)

(b)

Fig. 4. Comparison between the uniform flow assumption and the non-uniform distribution. (a) Stack voltage. (b) Fuel cell voltage distribution at 0.6, 1.0, and 1.3 A cm−2.

cell voltage distributions at 0.6, 1.0, and 1.3 A cm−2 are further presented in Fig. 4(b). With the uniform flow assumption, the cell voltages are almost identical under the same current density. For the distributed stack model, the output voltages of middle cells are found to be lower than end fuel cells, which are similar to parabola curves and consistent with previous studies [19,20]. The highest voltage locates in the fuel cell 1, and the lowest voltage happens in the cell 2. Therefore, the largest voltage difference is calculated to be 1.1% at 0.6 A cm−2, and it rises to 4.5% at 1.3 A cm−2. Higher current densities mean higher inlet mass flow rates of the stack, which results in larger pressure drops during the motion of reactants, leading to more uneven flow distribution among individual fuel cells. In general, neglecting the non-uniform flow distribution not only overestimates the overall stack performance but also underestimates the single cell voltage variations. It is well known that cell voltage variations among the stack mainly result from the uneven reactant and temperature distribution. Fig. 5 shows the anode and cathode inlet mass flow rate. It is found that variations of mass flow rate increase with the current density increasing. Besides, the degrees of cathode variations are more significant than that of the anode. The inlet mass flow rates of the end cells (cell 1, 5) are larger than the middle cells (cell 2, 3, 4), which partly leads to higher output voltages, and it is consistent with previous studies [20,23]. At 1.3 A cm−2, the real cathodic stoichiometry for cell 1 is calculated to be around 3.31, and it is approximately 1.25 for cell 2. The non-uniformity is supposed to aggravate if the current density further increases. Therefore, it is inferred that the cathodic reactant distributed to middle fuel cells might be insufficient under extremely high current densities even though the overall cathode stoichiometry is

2 in the study. In other words, supplying abundant air is an important issue for avoiding the possible oxygen starvation caused by the uneven reactant distribution [39,40]. The distribution of coolant follows similar trends as the reactants, and negligible change is seen under different current densities since the total inlet mass flow rate is set as constant. Apart from reactant distributions, stack performances are also influenced by the temperature distribution. Meanwhile, the temperature uniformity is of vital importance for the stack durability because temperature gradients may lead to thermal stresses and consequently result in the deformation of components. Fig. 6 shows the temperature differences between the uniform flow assumption and the non-uniform distribution, which is plotted as the temperature of the latter subtracted by the former. Note that the layer thickness in the figure is not to scale as the real thickness. At 0.6 A cm−2, the temperature of distributed stack model is observed to be around 0.7 °C higher than that of the uniform stack model. Besides, temperature differences do not change noticeably among the five cells. At 1.0 A cm−2, the highest temperature difference reaches around 1.7 °C which happens in the cell 3 and cell 4. The temperature difference in the anode channel of cell 1 and the cathode channel of cell 5 are found to be −0.4 °C. At the current density of 1.3 A cm−2, temperature differences in the cell 1 are approximately −1.0 °C while temperature differences in the cell 3 and cell 4 reach about 3.6 °C. It is obvious that the distributed stack model leads to more significant temperature variations at relatively high current densities. The reason is that reactants and coolant mass flow rate of end cells are much larger than middle cells, which takes away more heat after exiting the fuel cells. Besides, the output voltage of middle fuel

Fig. 5. Inlet mass flow rate of reactants among different fuel cells. (a) Anode. (b) Cathode. 8

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

BPc

0.6 A cm-2

CL MPL GDL CH

2.17

1.0 A cm-2

1.46 0.74 0.03

1.3 A cm-2

-0.69

Distributed, STc = 1.5 9.35

Flow

2.89

Temperature difference (K)

Flow

3.60

8.23 7.10

Distributed, STc = 1.8

5.98 4.85 3.73

Distributed, STc = 2.0

2.60

-1.41

1.48

-2.12

Cell 1

Cell 2

Cell 3

Cell 4

Oxygen concentration (mol m-3)

BPa

0.35

Cell 5

Cell 1

Fig. 6. Temperature differences between the uniform flow assumption and the non-uniform distribution.

Cell 2

Cell 3

Cell 4

Cell 5

Fig. 8. Effects of cathode stoichiometry on the oxygen distribution, 1.2 A cm−2.

cells is lower than end cells, which implies that more heat is generated owing to larger polarization voltage losses. Hence, neglecting the nonuniform flow distribution may lead to higher predictions of the overall stack temperature and lower predictions of the temperature variations among individual fuel cells.

changes. The corresponding cathode inlet mass flow rate is shown in Fig. 7(b). With the stoichiometry increasing, the inlet mass flow rate for all fuel cells rises although the non-uniformity becomes more obvious. At the stoichiometry of 1.5, the real stoichiometry for cell 1 is around 2.23 while it is approximately 1.07 for cell 2. At the stoichiometry of 1.8, the real stoichiometry for cell 1 is about 2.82, and it is about 1.20 for cell 2. The reason that a higher cathode stoichiometry is preferable because it increases the amount of reactant distributed into middle fuel cells, which consequently results in the major performance improvement. Generally speaking, the key to stack performance enhancement is to guarantee sufficient reactants for every fuel cell. The oxygen concentration distribution inside each single fuel cell is further presented in Fig. 8. Note that the layer thickness in the figure is not to scale as the real thickness of porous layers and gas channel. The oxygen concentration keeps decreasing towards the channel outlet owing to electrochemical consumptions. Besides, the oxygen concentration in the cell 1 and cell 5 are much higher than other cells, which results from larger inlet mass flow rates as shown in Fig. 7(b). At the stoichiometry of 1.5, the oxygen concentration of cell 2 and cell 3 in the outlet regions are approximately 0.5 mol m−3, indicating that oxygen is almost completely consumed. Suppose that the stack stoichiometry is less than 1.5, the cell 2 and cell 3 is possible to suffer from the local reactant starvation. The outlet oxygen concentration of cell 2 rises to 1.5 mol m−3 at the stoichiometry of 1.8, and it increases to 2.6 mol m−3 at the stoichiometry of 2.0. Therefore, a larger cathode stoichiometry also helps reduce the possibility of air starvation, which

3.3. Effects of stack operating conditions It is widely acknowledged that the operating conditions have great influences on stack performances. Based on Fig. 5, it is found that the cathode inlet mass flow rate varies significantly among different fuel cells while variations in the anode are much less noticeable. Therefore, effects of cathode stoichiometry on the flow distribution and stack performances are investigated while the current density is set as 1.2 A cm−2. Fig. 7(a) shows the fuel cell voltage distribution at the cathode stoichiometry of 1.5, 1.8, and 2.0 under the uniform flow assumption and the non-uniform distribution. It is seen that increasing the stoichiometry is beneficial to the stack output voltage. The performance improvement is more noteworthy when the stoichiometry increases from 1.5 to 1.8 in comparison with the cases when the stoichiometry rises from 1.8 to 2.0. It implies that further increasing the cathode stoichiometry may no longer lead to performance enhancement once it has reached the specific value. For the distributed stack model, the voltage of middle fuel cells rises remarkably with the increment of stoichiometry while the voltage of end cells only increases slightly, indicating that middle cells are more sensitive to stoichiometry

Fig. 7. Effects of cathode stoichiometry on stack performances, 1.2 A cm−2. (a) Fuel cell voltage. (b) Cathode inlet mass flow rate. 9

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

Fig. 9. Effects of inlet pressure on stack performances, 1.2 A cm−2. (a) Fuel cell voltage. (b) Cathode inlet mass flow rate.

Fig. 10. Effects of manifold geometric parameter on stack performances, 1.2 A cm−2. (a) Fuel cell voltage. (b) Cathode inlet mass flow rate.

is important for the stack durability. Since the reactant distribution is associated with pressure drops, the effects of stack inlet pressure are also studied. The anode and cathode inlet pressure are simultaneously changed from 1.4 atm to 1.6 atm. The cell voltage distribution is presented in Fig. 9(a), and the cathode inlet mass flow rate is shown in Fig. 9(b). It is seen that a higher inlet

pressure contributes to the stack performance improvement, which results from the increased reactant concentration. However, the increase of inlet pressure does not lead to better uniformity of cell voltage distribution. For different inlet pressures, the cathode inlet mass flow rates are almost identical. The reason is that the reactant distribution is strongly coupled with pressure drops rather than the absolute pressure. 3.4. Effects of stack manifold geometric parameters Since the reactant distribution is strongly dependent on the stack manifold structural designs, the effects of manifold geometric parameters are investigated to alleviate the non-uniform distribution phenomena. The operating current density is selected as 1.2 A cm−2. The cross-sectional sizes of reactant manifold are changed, and the corresponding results are presented in Fig. 10. In the case of 10 × 6 mm2, the largest voltage difference among the five fuel cells is calculated to be 3.5% while it decreases to about 0.8% in the case of 10 × 12 mm2. The cell voltage variations significantly decrease with larger reactant manifold areas. This is because increasing the sectional dimensions reduces the gas velocities, which consequently decreases the frictional pressure drops in stack manifolds. Besides, the local pressure drops caused by the diverging and converging configurations are also influenced. Note that the maximum voltage difference between the uniform assumption and the case of 10 × 12 mm2 is reduced to 0.4%. The air mass flow rates of individual fuel cells are shown in Fig. 10(b). With the increment of manifold sizes, the air distribution is found to be more uniform, and it is the same with the anode reactant distribution. The variations of reactant flow rate among the stack decreases, remarkably enhancing the output voltage of middle fuel cells.

Fig. 11. Local current density distribution along the channel direction of the fuel cell 2, 1.2 A cm−2. 10

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

The local current density distribution of the fuel cell 2 in the channel direction is shown in Fig. 11. Note that the node number 1 represents the channel inlet region while the number 5 represents the outlet region. The local current density keeps decreasing from the upstream channel towards the downstream channel, and it is caused by the reactants consumption. It is obvious that the distributed stack model leads to more severe non-uniformity of the local current density distribution compared with the uniform assumption. However, the current density variations are decreased with the cross-sectional area increasing, which actually results from the increased inlet reactant mass flow rate. In other words, supplying enough reactants for every fuel cell contributes to the uniformity of local current density inside single cells.

References [1] Nilsson M, Nykvist B. Governing the electric vehicle transition–near term interventions to support a green energy economy. Appl Energy 2016;179:1360–71. https://doi.org/10.1016/j.apenergy.2016.03.056. [2] Sabri MFM, Danapalasingam KA, Rahmat MF. A review on hybrid electric vehicles architecture and energy management strategies. Renew Sustain Energy Rev 2016;53:1433–42. https://doi.org/10.1016/j.rser.2015.09.036. [3] Wilberforce T, El-Hassan Z, Khatib FN, et al. Developments of electric cars and fuel cell hydrogen electric cars. Int J Hydrogen Energy 2017;42(40):25695–734. https:// doi.org/10.1016/j.ijhydene.2017.07.054. [4] Fathabadi H. Fuel cell hybrid electric vehicle (FCHEV): novel fuel cell/SC hybrid power generation system. Energy Convers Manag 2018;156:192–201. https://doi. org/10.1016/j.enconman.2017.11.001. [5] Jiao K, Li X. Water transport in polymer electrolyte membrane fuel cells. Prog Energy Combust Sci 2011;37(3):221–91. https://doi.org/10.1016/j.pecs.2010.06. 002. [6] Lü X, Qu Y, Wang Y, et al. A comprehensive review on hybrid power system for PEMFC-HEV: issues and strategies. Energy Convers Manag 2018;171:1273–91. https://doi.org/10.1016/j.enconman.2018.06.065. [7] Ehsani M, Gao Y, Longo S, et al. Modern electric, hybrid electric, and fuel cell vehicles. CRC Press 2018. https://doi.org/10.1201/9780429504884. [8] https://ssl.toyota.com/mirai/fcv.html (Accessed on June 21. 2019). [9] https://automobiles.honda.com/clarity-fuel-cell (Accessed on June 21. 2019). [10] https://www.hyundai.co.uk/new-cars/nexo (Accessed on June 21. 2019). [11] Miller M, Bazylak A. A review of polymer electrolyte membrane fuel cell stack testing. J Power Sources 2011;196(2):601–13. https://doi.org/10.1016/j.jpowsour. 2010.07.072. [12] Pei P, Chen H. Main factors affecting the lifetime of Proton Exchange Membrane fuel cells in vehicle applications: a review. Appl Energy 2014;125:60–75. https:// doi.org/10.1016/j.apenergy.2014.03.048. [13] Hu Z, Xu L, Li J, et al. A novel diagnostic methodology for fuel cell stack health: performance, consistency and uniformity. Energy Convers Manag 2019;185:611–21. https://doi.org/10.1016/j.enconman.2019.02.031. [14] Yang Z, Du Q, Jia Z, et al. Effects of operating conditions on water and heat management by a transient multi-dimensional PEMFC system model. Energy 2019. https://doi.org/10.1016/j.energy.2019.06.148. [15] Kim S, Hong I. Effects of humidity and temperature on a proton exchange membrane fuel cell (PEMFC) stack. J Ind Eng Chem 2008;14(3):357–64. https://doi.org/ 10.1016/j.jiec.2008.01.007. [16] Yi P, Peng L, Feng L, et al. Performance of a proton exchange membrane fuel cell stack using conductive amorphous carbon-coated 304 stainless steel bipolar plates. J Power Sources 2010;195(20):7061–6. https://doi.org/10.1016/j.jpowsour.2010. 05.019. [17] Pei H, Liu Z, Zhang H, et al. In situ measurement of temperature distribution in proton exchange membrane fuel cell I a hydrogen–air stack. J Power Sources 2013;227:72–9. https://doi.org/10.1016/j.jpowsour.2012.11.027. [18] Devrim Y, Devrim H, Eroglu I. Development of 500 W PEM fuel cell stack for portable power generators. Int J Hydrogen Energy 2015;40(24):7707–19. https:// doi.org/10.1016/j.ijhydene.2015.02.005. [19] Baschuk JJ, Li X. Modelling of polymer electrolyte membrane fuel cell stacks based on a hydraulic network approach. Int J Energy Research 2004;28(8):697–724. https://doi.org/10.1002/er.993. [20] Park J, Li X. Effect of flow and temperature distribution on the performance of a PEM fuel cell stack. J Power Sources 2006;162(1):444–59. https://doi.org/10. 1016/j.jpowsour.2006.07.030. [21] Cozzolino R, Cicconardi SP, Galloni E, et al. Theoretical and experimental investigations on thermal management of a PEMFC stack. Int J Hydrogen Energy 2011;36(13):8030–7. https://doi.org/10.1016/j.ijhydene.2011.01.052. [22] Salva JA, Iranzo A, Rosa F, et al. Experimental validation of the polarization curve and the temperature distribution in a PEMFC stack using a one dimensional analytical model. Int J Hydrogen Energy 2016;41(45):20615–32. https://doi.org/10. 1016/j.ijhydene.2016.09.152. [23] Amirfazli A, Asghari S, Koosha M. Mathematical modeling and simulation of thermal management in polymer electrolyte membrane fuel cell stacks. J Power Sources 2014;268:533–45. https://doi.org/10.1016/j.jpowsour.2014.06.073. [24] Amirfazli A, Asghari S, Sarraf M. An investigation into the effect of manifold geometry on uniformity of temperature distribution in a PEMFC stack. Energy 2018;145:141–51. https://doi.org/10.1016/j.energy.2017.12.124. [25] Chen D, Ding K, Chen Z, et al. Physics field distributions within fuel cell stacks with manifolds penetrating through the plane zone and open outlet features. Energy Convers Manag 2018;178:190–9. https://doi.org/10.1016/j.enconman.2018.10. 034. [26] Macedo-Valencia J, Sierra JM, Figueroa-Ramírez SJ, et al. 3D CFD modeling of a PEM fuel cell stack. Int J Hydrogen Energy 2016;41(48):23425–33. https://doi.org/ 10.1016/j.ijhydene.2016.10.065. [27] Luo Y, Guo Q, Du Q, et al. Analysis of cold start processes in proton exchange membrane fuel cell stacks. J Power Sources 2013;224:99–114. https://doi.org/10. 1016/j.jpowsour.2012.09.089. [28] Lim BH, Majlan EH, Daud WRW, et al. Three-dimensional study of stack on the performance of the proton exchange membrane fuel cell. Energy 2019;169:338–43. https://doi.org/10.1016/j.energy.2018.12.021. [29] Idelchik IE. Handbook of hydraulic resistance. Washington, DC: Hemisphere Publishing Corp.; 1986, 662 p. Translation, 1986. [30] Pei P, Li Y, Xu H, et al. A review on water fault diagnosis of PEMFC associated with

4. Conclusion In the study, the flow distribution and performance heterogeneity are investigated based on a 1+1 dimensional multiphase proton exchange membrane fuel cell (PEMFC) stack sub-model integrated with a flow distribution sub-model. Frictional pressure drops as well as local pressure drops due to the diverging, converging, and bending configurations are calculated. For every fuel cell, transport processes are considered in the through-plane direction and the gas channel direction. Phase changes among the membrane water, liquid water, and water vapor are taken into account. The integrated stack model has been validated against experimental data, including the polarization curve and the cell voltage distributions. Differences between the uniform flow assumption and the nonuniform distribution are quantitatively compared, including the distribution of cell voltage, mass flow rate, temperature, reactant concentration, and local current density. Results show that the uniform flow assumption not only overestimates the stack output performance but also underestimates the fuel cell voltage variations. Besides, neglecting the nonuniform flow distribution may lead to higher predictions of the overall stack temperature and lower predictions of the temperature variations among different fuel cells. The inlet reactant mass flow rate of middle fuel cells is found to be lower than end fuel cells, and the cell voltage distribution follows the similar trend. Even though the total amount of air supplied to the stack is abundant, it is still possible for some fuel cells to suffer from the local reactant starvation. Therefore, a higher cathode stoichiometry is preferable because it increases the inlet mass flow rate for middle fuel cells, which significantly contributes to the stack performance improvement. Increasing the inlet pressure contributes negligibly to the uniformity of reactant distribution, but it is favorable for the stack performance owing to enhanced electrochemical reactions. A larger stack manifold cross-sectional area is suggested because it not only contributes to the more uniform reactant distribution among the stack but also helps decreasing the local current density variations inside individual fuel cells. CRediT authorship contribution statement Zirong Yang: Methodology, Software, Validation, Writing - original draft. Kui Jiao: Conceptualization, Validation, Resources, Writing review & editing. Zhi Liu: Data curation, Writing - review & editing. Yan Yin: Writing - review & editing, Funding acquisition. Qing Du: Writing - review & editing, Supervision, Project administration. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work is supported by the National Key Research and Development Program of China (2018YFB0105601) and the National Natural Science Foundation of China (grant No. 51976138). 11

Energy Conversion and Management 207 (2020) 112502

Z. Yang, et al.

[31] [32]

[33] [34] [35]

the pressure drop. Appl Energy 2016;173:366–85. https://doi.org/10.1016/j. apenergy.2016.04.064. Baek SM, Yu SH, Nam JH, et al. A numerical study on uniform cooling of large-scale PEMFCs with different coolant flow field designs. Appl Therm Eng 2011;31(8–9):1427–34. https://doi.org/10.1016/j.applthermaleng.2011.01.009. Yang Z, Liu Z, Fan L et al. Modeling of proton exchange membrane fuel cell system considering various auxiliary subsystems. In: The energy and sustainability 2018 symposium. Springer, Cham; 2018. p. 18–33. https://doi.org/10.1007/978-3-03000105-6_2. Wang B, Wu K, Yang Z, et al. A quasi-2D transient model of proton exchange membrane fuel cell with anode recirculation. Energy Convers Manag 2018;171:1463–75. https://doi.org/10.1016/j.enconman.2018.06.091. Huo S, Jiao K, Park JW. On the water transport behavior and phase transition mechanisms in cold start operation of PEM fuel cell. Appl Energy 2019;233:776–88. https://doi.org/10.1016/j.apenergy.2018.10.068. Yang Z, Du Q, Huo S, et al. Effect of membrane electrode assembly design on the

[36] [37] [38] [39] [40]

12

cold start process of proton exchange membrane fuel cells. Int J Hydrogen Energy 2017;42(40):25372–87. https://doi.org/10.1016/j.ijhydene.2017.08.106. Zhang G, Fan L, Sun J, et al. A 3D model of PEMFC considering detailed multiphase flow and anisotropic transport properties. Int J Heat Mass Trans 2017;115:714–24. https://doi.org/10.1016/j.ijheatmasstransfer.2017.07.102. Li X. Principles of fuel cells. CRC Press; 2005. Luo Y, Jiao K, Jia B. Elucidating the constant power, current and voltage cold start modes of proton exchange membrane fuel cell. Int J Heat Mass Trans 2014;77:489–500. https://doi.org/10.1016/j.ijheatmasstransfer.2014.05.050. Hu Z, Xu L, Li J, et al. A cell interaction phenomenon in a multi-cell stack under one cell suffering fuel starvation. Energy Convers Manag 2018;174:465–74. https://doi. org/10.1016/j.enconman.2018.08.062. Chen H, Zhao X, Zhang T, et al. The reactant starvation of the proton exchange membrane fuel cells for vehicular applications: a review. Energy Convers Manage 2019;182:282–98. https://doi.org/10.1016/j.enconman.2018.12.049.