Effect of water vapor condensation on the flow distribution in a PEM fuel cell stack

International Journal of Heat and Mass Transfer 151 (2020) 119471

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/hmt

Effect of water vapor condensation on the flow distribution in a PEM fuel cell stack Ahmad Rezaei Sangtabi a, Ali Kianifar a,∗, Ebrahim Alizadeh b a b

Department of Mechanical Engineering, Ferdowsi University of Mashhad, Mashhad, Iran Malek Ashtar University of Technology, Tehran, Iran

a r t i c l e

i n f o

Article history: Received 25 November 2019 Revised 20 January 2020 Accepted 2 February 2020

Keywords: Maldistribution PEMFC Water management Manifold Flow uniformity

a b s t r a c t In this study, the effect of water vapor condensation on the flow distribution in the manifold of a PEMFC stack is investigated. A solver is developed in OpenFOAM software and is validated against experimental and numerical data. Due to the condensation, the water vapor concentration in the mixture decreases near the manifold walls. Thus, beginning and last cells receive more oxygen than the middle cells. The condensed water forms small droplets on the manifold walls. The small droplets coalesce into larger droplets. Descending large droplets into the last cell of the 26-cell stack leads to the formation of the slug. This process results in blockage of flow through the last cell, causing a significant increase in the flow non-uniformity. The non-uniformity index increases significantly from 0.046 to 0.94 in the water removal process. When the number of cells increases from 26 to 39, the large droplet detaches from the manifold and sinks into seven middle cells. The flow passing through these cells significantly reduces, resulting in a severe flow maldistribution. A comparison between flow distribution parameters indicates that the root mean square value of the deviation from the mean flow rate is a better criterion than the non-uniformity index to determine the severity of maldistribution. © 2020 Elsevier Ltd. All rights reserved.

Abbreviations PIMPLE SLPM

Merged PISO-SIMPLE Algorithms Standard Liter Per Minute

1. Introduction Fuel cells are electrochemical engines that directly convert the chemical energy stored in fuels to electricity. Due to their higher efficiency, durability, reliability and lower emissions, fuel cells are considered as an alternative energy to conventional thermal power systems. Among various types of fuel cells, polymer electrolyte membrane fuel cells (PEMFCs) with a low operating temperature have a high power density, fast start-up and zero carbon emission [1,2]. A single PEMFC has a potential of 0.6–1.0 V depending on the load, therefore several individual fuel cells are serially connected to create a PEM fuel cell stack. Manifold system is used to feed reactant gasses (fuel and oxidant) to each individual cell of a stack. Since all cells in a stack are thermally and electrically connected

∗

Corresponding author. E-mail address: [email protected] (A. Kianifar).

https://doi.org/10.1016/j.ijheatmasstransfer.2020.119471 0017-9310/© 2020 Elsevier Ltd. All rights reserved.

in series, the overall performance of a stack depends on a satisfactory operation of all individual cells. The reactants must be uniformly distributed from the inlet manifold to the individual cells. Gas flow maldistribution from cell to cell could introduce water flooding, membrane drying, localized hot spots in the membrane, material degradation, which has a significant impact on the fuel cell efficiency [1,3]. The flow maldistribution in the PEMFC stack has been studied using numerical and experimental approaches. Baschuk and Li [4] used a two part mathematical model based on the hydraulic network analysis to predict the performance of the PEMFC stack with 50 cells. They used the distribution of mass flow rate and pressure in the stack from the first part calculation as an input parameter for the second part of their mathematical model, which determined the voltage of each cell in the stack. The results showed that the output voltage of cells within a stack being lower when compared to a fuel cell operating individually. They suggested three methods for improving the flow distribution in the stack: increasing the cross section area of the stack manifold, decreasing the number of gas flow channels per bipolar plate and varying the resistance to mass flow in the gas flow channels. Koh et al. [5] used the microscopic mechanical energy balance with the pressure loss by geometrical effects and wall friction to analyze the gas flow in the cathode and anode manifolds of a stack.

2

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

Nomenclature Cp D f F1 F2 g h Hfg k kcond kevap m˙ m˙ i ¯˙ m M n N p R t t T U Ur x y

Specific heat, J.kg−1 .K−1 Diffusion coefficient, m2 .s−1 Body force, kg.m−2 .s−2 Non-uniformity index Flow distribution parameter Gravity, m.s−2 Heat convection coefficient, W.m−2 .K−1 Latent heat of water, J.kg−1 Thermal conductivity, W.m−1 .K−1 Condensation rate constant, s−1 Evaporation rate constant, Pa−1 .s−1 Phase change mass transfer rate, kg.m−3 .s−1 Mass flow rate, kg.m−3 .s−1 Average mass flow rate, kg.m−3 .s−1 molecular weight, kg.mol−1 unit vector Number of Channels Pressure, Pa Gas constant, J.kg−1 .K−1 Time, s unit vector Temperature, K Velocity, m.s−1 Compressive velocity, m.s−1 Molar fraction Mass fraction

Greek symbols ρ Density, kg. m−3 α Volume fraction μ Dynamic viscosity, kg.m−1 .s−1 ϑ Permeability, m2 κ Interface curvature, m−1 σ surface tension coefficient, J.s−2 θ contact angle, ° subscript 1 2 c amb g i j mix oxy sat vapor w

Liquid phase Gas phase Critical Ambient Gas Species i Species j Mixture Oxygen Saturation Water vapor Wall

superscript H2 o Water They showed that the non-uniform distribution is more significant at the cathode side. Karimi et al. [6] determined the molar flow rate and pressure for the air and hydrogen streams in the stack based on a flow network model incorporating the minor losses. They studied the effect of inlet-outlet configuration on the stack performance and showed that a symmetric double inlet-single outlet configuration provides better performance than a single inletsingle outlet configuration. Maharudrayya et al. [7] used a onedimensional analysis for the evaluation of relative flow distribution and pressure drop in the PEMFC stack. They showed that the

geometric factors such as the manifold dimension, the rib width between cells and the channel dimension play an important role in the stack performance. Wang [8] developed an analytical model based on the mass and momentum conversation to solve the pressure and flow distribution of U type PEMFC stacks. They found that momentum and friction effects work in opposite directions and the proper balance of two effects can result into an optimal design of the stack and less flow maldistribution. Two phase flow, turbulence, heat transfer and variation of fluid properties were neglected in the most analytical methods [1,4–10]. Due to a large number of exits to/from the cells, pipe network approximations are inappropriate for use to calculate gas flows in PEMFC stacks [11]. The computational fluid dynamics is used to the simulation of flow distribution in the manifold of the stack to take account of the effects of these factors. The effect of cathode manifold configuration on the overall performance of a 10-cell stack was studied by Kim and Kim [12]. Numerical results showed that flow separation occurs at the inner wall in a manifold with a 90° turn. Elimination of separation zones provides an 8.5% more uniform flow distribution. Chen et al. [13] used computational fluid dynamics to investigate the effects of flow channel resistance, manifold widths and air feed on the flow distribution in the PEMFC stack composed of 72 cells. Their numerical results indicated that lesser air feed, larger flow resistance and manifold widths are beneficial for flow distribution. Lebaek et al. [14,15] experimentally and numerically studied the inlet effect on the flow distribution in a 70 cell PEMFC stack. They found that a very distinct jet forms in the manifold when using the circular inlet configuration and proper inlet conditions are a key factor for an even flow distribution. Maldistribution reduced significantly by implementing a simple diffuser at the inlet of the manifold. Liu et al. [16] simulated the air side manifold of an industrial 5 kW PEMFC stack consisting of 67 cells. Numerical results showed that adding a baffle in the inlet manifold reduces the pressure variation in the channels from 39% to 2%. Dong et al. [17] designed a tree-type modular manifold based on the principle of minimizing entropy generation. CFD simulations showed that better uniformity achieves when the ratio of channel length to equivalent diameter in each stage increases to larger than 7. Jackson et al. [18] modified a simple discrete method to deduced a mathematical relationship for increasing manifold width to reduce flow maldistribution. Although flow distribution improves by their introduced relationship, the manifold width increases significantly. Alizadeh et al. [19] designed a cascade type dead end PEMFC stack. They found that water formation and accumulation at cathode sides is more than anode sides, and the fluctuation of voltage in cathode purge cells is higher than anode’s one. In the reviewed numerical studies, flow in the manifold of the PEMFC stack is considered as a single phase. Two phase flow simulations were used to investigate the performance of a single PEMFC [20–23] and liquid water removal process from the flow channels [24–29]. In this paper, the effect of the water vapor phase change on the flow distribution in the PEMFC stack is investigated at the different operating conditions. For this purpose, a new solver is developed in the OpenFOAM software. First, the flow channels of a bipolar plate are simulated to calculate the pressure drop across the flow channels of the bipolar plate. Then, flow channels are replaced by straight parallel channels filled with porous media in the simulation of the stack, due to reducing the computational time. The effects of the water vapor condensation, operating temperature, water removal process and number of fuel cells in the stack on the flow distribution are investigated. Also, the effect of the number of cells on the water removal process is discussed. The flow uniformity under different conditions is calculated using two parameters.

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

3

Fig. 1. Computational domain. Table 1 Dimensions of the stack and bipolar plate. Stack dimensions Number of cells Width (mm) Depth (mm) Rib (mm) Channel width (mm) Channel height (mm) Inlet manifold height (mm) Outlet manifold height (mm)

Bipolar plate 26 191.36 37.2 2.76 0.8 57.2 21.7 21.7

Number of channel Channels width (mm) Channels depth (mm) Rib (mm) Active area (cm2 )

45 1 0.8 1 500

reduce the computational time. Then, 2D straight parallel channels filled with porous media are used to describe the pressure drop in the stack. Fig. 2 shows a schematic representation of the bipolar plate with a multi-pass serpentine flow channel configuration with 45 channels per plate and 12 U-turns per channel. The performance of a fuel cell improves by multi-pass serpentine flow fields [30–32]. The dimensions of the stack and bipolar plates are listed in Table 1. 3. Governing equations Fig. 2. Pressure drop in the flow channels of bipolar plates.

2. Computation domain A two-dimensional computational domain is displayed in Fig. 1. The domain represents the cathode side of a 26-cell PEM fuel cell stack consisting of the inlet and outlet manifolds, ribs and parallel channels. The dimensions for this PEMFC stack are obtained from a 6-kW fuel cell stack and are based on an active area of 500 cm2 . Both manifolds have a cross-section of 21.7 mm × 37.2 mm. Due to the high computational cost of two-phase 3D simulations, only the two-dimensional simulations are performed to determine flow distribution in the fuel cell stack. The dimensions for the geometry are shown in Fig. 1. As shown in the cut through its mid-plane (Fig. 2), 3D simulation of flow in the gas flow channels of one bipolar plate is performed to evaluate the pressure drop across the flow channels to

In this paper, attention is concerned only on the flow distribution according to the condensation of water vapor. Therefore, the electrochemical reaction is not taken into account. The following assumptions are considered to simplify the PEMFC model in this work: (1) The incompressible flow regime is laminar due to the low mixture flow velocity. (2) The stack inlet flow is humidified oxygen. (3) The gas mixture behaves as an ideal gas. (4) The flow channel is assumed to be homogeneous and isotropic with uniform permeability. (5) The gas and liquid phases move at the same velocity. The proton conductivity of the membrane depends on its water content. Fuel and oxidant gasses typically must be humidified before entering the cell [33]. Based on these assumptions, the gov-

4

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

erning equations of unsteady multiphase flow in the cathode side of PEMFC stack are given as below: Continuity equation:

∂ρ + ∇ . (ρ U ) = 0 ∂t

(1)

Where ρ and U are density and velocity. The density is calculated by volume fraction weighted averaging:

ρ = ρ1 α1 + ρ2 α2

(2)

Where α denote the volume fraction. The Volume of Fluids (VOF) technique is used to calculate the volume fraction of gas and liquid phases. The interface between gas and liquid is tracked using the volume fraction of liquid water in the computational cell volume. Among different surface reconstruction schemes, the VOF method is widely used to describe the phase transition mechanism between liquid and vapor phases [24–27,34–36]. The subscripts 1 and 2 refer to the liquid and gas phase, respectively. The density of liquid water is considered to be constant and the density of the mixture is obtained by the ideal gas law:

ρ2 =

P

(3)

Rmix T

Momentum equation:

(4)

f = σ κ∇ α1

(5)

Where σ is the surface tension coefficient. κ is the interface curvature and is defined as:

κ = ∇ .(nw cos θw + tw sin θw )

(6)

Where nw and tw are the unit vectors normal and tangential to the wall, respectively. In the all simulations, the static contact angle (θ w ) is considered to be constant. Energy equation

∂ (ρ C p T ) + ∇ .(ρC pUT ) = ∇ .(k∇ T ) + m˙ H f g ∂t

(7)

Where Cp is the heat capacity, T the temperature, k the heat conduction coefficient,m˙ the phase change rate, Hfg the latent heat of water. The phase changes rate is expressed as [36]:

m˙ =

MgH2 o



2o 2o P xH − xH g sat ⎩k α ρ PxH2Ro T− xH2 o evap 1 1 g sat



2o 2o xH ≥ xH g sat

(8)

2o 2o xH < xH g sat

Where kcond and kevap are condensation and evaporation rates constants. The saturation pressure of water vapor (Psat ) as a function of local temperature is calculated by [36]:

Psat = −2846.4 + 411.24(T − 273.15 ) − 10.554(T − 273.15 )2 + 0.166636(T − 273.15 )3

(9)

The heat latent of water is evaluated as [38]:

H f g = 307090(647.15 − T )0.35549



xμ i i x j φi j

i

kmix =

 i

xk i i x j φi j

(12)

j

1+

φi j =

(11)

j

0.5 0.25 2 μi μj



8 1+

Mj Mi

Mi Mj

0 . 5

(13)

Species equation:

∂ ( ρ yi ) + ∇ .(ρU yi ) = ∇ .(ρ Di j ∇ yi ) − m˙ i ∂t

(14)

Where yi is the mass fraction of species i and Dij is a binary diffusion coefficient of species i in the species j and is defined as [39]:

a T 1 1 Di j = ( )b (Pc,i Pc, j )1/3 (Tc,i Tc, j )5/12 + P M M i j Tc,i Tc, j

Where ϑ is the permeability of the porous media. The highest Reynolds number in the manifold of the 26-cell stack is lower than 20 0 0, and the average Reynolds number in the parallel channels is less than 200. The range of Reynolds number shows that the regime of the flow inside the manifold is laminar. The numerical result shows that the pressure drop in the flow channels is about 1142 Pa and the corresponding permeability is 1.765E-09 m2 . f is the source term corresponding to the surface tension effect. The surface tension effect is considered using the continuum surface force (CSF) model proposed by Brackbill et al. [37]:

kcond (1 − α1 )

μmix =



∂ (ρ U ) + ∇ .(ρU U ) = −∇ P + ρ g + ∇ .(μ(∇ U ) ∂t μ + ( ∇ U )T ) + f − U ϑ

⎧ ⎨

Viscosity, thermal conductivity and heat capacity of the two phase flow are calculated by volume fraction weighted averaging. The ideal gas mixture relation is used for calculating the viscosity and thermal conductivity of the mixture [36]:

(10)

1/2 (15)

Where Pc,i , Pc,j , Tc,i and Tc,j are the critical pressure and temperature of species of i and j, respectively. The constants are a = 0.0 0 0364 and b = 2.334 for water vapor [39]. Volume fraction equation:

1 ∂ α1 1 1 + ∇ .(U α1 ) + ∇ .(Ur α1 (1 − α1 )) = −m˙ ( − α1 ( − ) ∂t ρ1 ρ1 ρ2 (16) Where Ur is compressive velocity [40] and is calculated in the normal direction to avoid any dispersion. 4. Boundary condition and numerical method A Velocity inlet boundary condition is applied at the stack inlet. The inlet velocity is calculated by using the gas mass flow rate and the inlet surface area. A no-slip boundary condition is applied to the cells and manifolds walls. The inlet temperature is equal to the cell operating temperature. The temperature on the internal surface of cells is specified as the cell operating temperature. Convective heat transfer is applied to the external manifold walls with a heat transfer coefficient of 50 W/m2 k. The partial pressure of water vapor is a function of temperature. The water vapor concentration in the inlet mixture is calculated by using the partial pressure and the relative humidity. The relative humidity of the oxygen at the inlet is set to 100% and it is assumed that there are no liquid water droplets in the inlet flow. Different static contact angles are applied for walls. The static contact angles of the cell and manifold walls are set at 135° and 45°, respectively, based on typical PTFE treated carbon paper GDL materials and carbon plate [41,42]. A pressure outlet boundary condition is specified at the stack outlet and the other variables are extrapolated at the outlet. The governing equations and boundary conditions are implemented and solved in the finite volume software OpenFOAM. The pressure-velocity coupling is accomplished by the PIMPLE algorithm [43]. Central differencing and VanLeer [44] schemes are used for discretization of the diffusion and convection terms of the conversion equations. For the discretization of temporal terms, the first order accurate Euler scheme was adopted. The simulating parameters are listed in Table 2.

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

5

Fig. 3. Mesh independency (N = 26, T = 353.15 K).

Table 2 Simulation parameters [36,41,42]. Operating conditions and simulation parameters P (atm) T (K) Oxygen flow rate (SLPM/Cell) Current density (A/cm2 ) Relative humidity (%) Tamb (K)

2 343.15 - 353.15 1.6 2 100 298.15

σ (J/s2 ) θ manifold (°) θ channel (°) kcond (1/s) kevap (1/Pa.s) hamb (W/m2 .K)

0.07 45 135 100 9.8e-6 50

5. Mesh independency and validation Four different meshes are tested by reporting the gas vertical velocity at the center of fuel cells to ensure that the solution was grid independent. The results of grid independence are shown in Fig. 3. The vertical velocity does not change by more than 1% when using cell number 60,0 0 0 and 80,0 0 0. Hence, a mesh of about 60,0 0 0 elements is selected for further simulations. There are not any experimental and numerical investigations on the effect of water vapor condensation on the flow distribution in a fuel cell stack. For this reason, the developed code is validated by two different problems: flow distribution in a 10-cell stack and interaction between liquid water and air in the anode channel of a fuel cell. The experimental and numerical data for normalized air velocity in the fuel cells of a stack given by Kim and Kim [12] are used to validate the present numerical results. The geometry and operating conditions are adopted from Ref. [12]. The air at ambient pressure and temperature flowed through the 10-cell stack with a velocity of 0.57 m/s. Due to low velocity and relative humidity, the flow is single phase and incompressible. Fig. 4 shows that the present numerical result has a good agreement with the experimental and numerical results of Kim and Kim [12]. Interaction between liquid water with air flow in the anode channel of a fuel cell and water removal process is simulated and compared with numerical results obtained by Soopee et al.

[26] and Zhu et al. [45]. Air flow enters the anode channel at a velocity of 10 m/s and water with a velocity of 1 m/s enters the channel from the pore on the bottom wall. A static contact angle of 45° is used for the upper wall. The bottom wall is set to be hydrophobic with a contact angle of 140°. The boundary, operating and initial conditions are similar to the studies of Soopee et al. [26] and Zhu et al. [45]. The different factors in the present study are given in Table 3. The computational domain, boundary conditions and the water droplet growth profiles at different simulation times are shown in Fig 5. The comparison between numerical results shows a quite small difference and indicates the results of the present study are reasonable. 6. Results and discussion In this section, the influence of vapor condensation on the flow distribution in a 26-cell PEM fuel cell stack is investigated. Then, the effects of the operating temperature and number of fuel cells in a stack on the flow distribution are studied. In this paper, the severity of maldistribution is determined by two parameters proposed by Maharudrayya et al. [7]. The non-uniformity index (F1 ) is based on the maximum and minimum values of the mass flow rates, defined as:

F1 =

max(m˙ 1 . . . .m˙ N ) − min(m˙ 1 . . . .m˙ N ) max(m˙ 1 . . . .m˙ N )

(17)

The mass flow rates are computed via integrating over the cell cross section. The value of F1 varies between 0 and 1 and shows the extent of flow maldistribution among cells. F1 =0 indicates uniform flow distribution and F1 →1 shows the severity of nonuniformity. The second parameter is the root mean square value of the deviation from the mean flow rate (F2 ) and given by:



F2 =

n  i=1

2

(m˙ i − m¯˙ ) ¯˙ Nm

(18)

6

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

Fig. 4. Comparison between present simulation and experimental data. Table 3 List of differences between the present study and the simulations of Soopee et al. [26] and Zhu et al. [45]. Differencing parameters

Soopee et al. [26]

Zhu et al. [45]

Present study

Software Number of cells VOF formulation Spatial discretization Transient formulation

Fluent 92,250 Implicit Compressive Non-Iterative advancement

Fluent 15,810 Explicit Geometric reconstruction Iterative Advancement

OpenFOAM 60,750 Semi-implicit Compressive Non-Iterative advancement

Fig. 5. Comparison of water droplet growth profiles as simulated by the present study and by Zhu et al. [45] and Soopee et al. [26].

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

7

Fig. 6. Oxygen mass fraction contour obtained from numerical simulation of two phase flow (N = 26, T = 353.15 K).

Fig. 7. Flow distribution patterns under two phase and single phase conditions (N = 26, T = 353.15 K).

The lower value of F2 shows a better flow distribution among the individual cells. 6.1. Effect of phase change The partial pressure of the water vapor is a function of temperature. When the temperature of the mixture decreases, the partial pressure of water vapor exceeds the local saturation pressure and condensation takes place. Due to the temperature difference between the operating temperature and the ambient, the temperature of the mixture decreases as the gas moves toward the end of the manifold. The contour of the oxygen mass fraction obtained from the simulation of the two phase flow is displayed in Fig. 6. It is found that the beginning and last cells receive more oxygen than the middle cells. The oxygen flow rates at the center of each cell for single phase and two phase flow simulations are presented in Fig. 7. The horizontal coordinate represents the cell number, where no. 1 is in-

dicated as the cell at the feeding end. Results show that the cells receive almost equal oxygen under the single phase conditions. In the two phase simulation, a very small amount of liquid water passed through the cells can immediately change the flow distribution. Therefore, the flow distribution varies by the time when the vapor condensation occurs in the inlet manifold. At t = 14 s, the oxygen flow rate through the few cells (1 to 4 and 20 to 26) is more than other cells (5 to 19) as explained before. The last cell receives the highest oxygen flow. At t = 25.2 s, a very small amount of liquid water passed through the last cell. As a result, the last cell receives the lowest oxygen mass flow rate and flow in the other cells increases to maintain the constant total flow rate. The significant difference between the flow distribution patterns, under single phase and two phase flow conditions, is in agreement with the experimental finding of Kandlikar et al. [46]. ¯˙ ) obFig. 8 compares the normalized mass flow rate (m˙ i /m tained from the single phase simulation, which is the ratio of the mass flow rate at the center of each cell to the average

8

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

Fig. 8. Normalized mass flow rate under two phase and single phase conditions (N = 26, T = 353.15 K).

Fig. 9. water droplets formation on the manifold top wall (N = 26, T = 353.15 K).

Table 4 Effect of phase change on the flow uniformity.

F1 F2

6.2. Effect of water removal process

Single phase

Two-phase, t = 14 s

Two-phase, t = 25.2s

0.0011 5.02e-5

0.0142 8.15e-4

0.0278 9.74e-4

mass flow rate entering each cell, with the normalized oxy¯˙ oxy ) and normalized total mass flow rate gen flow rate (m˙ oxy /m (m˙ oxy+vapor /m¯˙ oxy+vapor ) obtained from the two phase simulation. Comparison between Fig. 8 and Fig. 7 shows that the total mass distribution among the individual cells in the simulations of two phase and single phase flow is approximately the same. Due to the vapor condensation, some cells receive more oxygen and less water vapor under the two-phase condition. The two factors of F1 and F2 for Fig. 7 are calculated and listed in Table 4. The mass flow rates of oxygen through the unit cell are used to calculate the flow uniformity parameters. The results show that F1 and F2 in the two phase simulation are 25 and 19 times more than the single phase simulation.

The contours of liquid water volume fraction obtained from numerical simulation of two phase flow at different times are displayed in Fig. 9 and Fig. 10. The inlet manifold is enlarged in Fig. 9 and shows water droplets produced by condensation of water vapor on the manifold walls. The drag force, created by the flow passing the water droplet, pushes the droplets toward the end of the manifold and droplets join together in the corner of the inlet manifold. Due to gravity, most of the accumulated water descends to the last cell and discharges from the stack by the outlet manifold. The water removal process from the stack is illustrated in Fig. 10. At t = 24 s, liquid water divide into two parts: most of the water descends into the last cell due to the effect of gravity, while the rest stick on the top right corner of the inlet manifold due to wall adhesion. This phenomenon can lead to blockage of the bipolar plate channels by slug formation. At t = 24.6 s, drained water from the cell 26 sticks to the top wall of the outlet manifold. Due to a large amount of liquid water, some liquid water detaches from the top wall of the outlet manifold and reaches the bottom wall (t = 25.2 s). At this time, the small droplets join together at

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

9

Fig. 10. water removal process (N = 26, T = 353.15 K).

the top right corner of the inlet manifold and start to enter into the cell 26 at t = 25.6 s. The oxygen flow rates in individual cells are plotted in Fig. 11 according to the results shown in Fig. 10. To better show the flow distribution, the left vertical coordinate represents the oxygen flow rates through cells 1 to 25 and the right vertical coordinate indicates the oxygen flow rate through the cell 26. As seen in Fig. 10, the oxygen flow rate in cell 26 is not reduced to zero and indicates a semi-slug flow [46,47]. The flow distribution parameters for selected simulation times are listed in Table 5. The efficiency of a stack depends on the performance of all individuals cells. In the water draining process, the oxygen flow rate in cell 26 drops sharply from 0.0359 gr/s to 0.00236 gr/s. As a result, the value of F1 increases to 0.9355. Thus, the performance of a stack decreases significantly. The value of F2 increases from 5.02e-5 in the single phase simulation to 3.66e-2 in the water removal process in the two phase flow simulation.

Table 5 Effect of water draining process on the flow uniformity.

F1 F2

t = 23.8s

t = 24s

t = 24.6s

t = 25.2s

t = 25.6s

0.0458 1.99e-3

0.8025 3.11e-2

0.9355 3.66e-2

0.0278 9.74e-4

0.4624 1.75e-2

6.3. Effect of the operating temperature The water vapor concentration in the mixture decreases when the operating temperature decreases. Also, the heat transfer rate from the manifold walls reduces as the temperature difference between the ambient and the inlet gas decreases. As a result, the rate of condensation decreases. The maximum local volumetric phase change rate (m˙ ) on the manifold top wall decreases from 10.77 kg/m3 .s to 6.95 kg/m3 .s as the operating temperature changes

10

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

Fig. 11. Oxygen flow rate through the unit cell at different time (N = 26, T = 353.15 K).

Fig. 12. liquid water production on the top wall of the inlet manifold (N = 26). Table 6 Effect of operating temperature on the flow uniformity. time

t=8s

Top (K) F1 F2

343.15 0.0116 4.9e-4

t=9s 353.15 0.0189 8.4e-4

343.15 0.0253 1.3e-3

t = 11.6 s 353.15 0.0141 8.8e-4

from 353.15 K to 343.15 K. The surface integral of volumetric phase change rate over the manifold top wall is shown in Fig. 12. At the beginning of the simulations, there is no liquid water on the manifold. For this reason, the condensation rate has the highest value. As the operating temperature rises, the droplets grow faster and liquid water removes from the stack earlier. The liquid water de-

343.15 0.0201 1.0e-3

t = 14.2 s 353.15 0.0216 9.8e-4

343.15 0.0153 9.3e-4

353.15 0.0151 1.0e-3

scends into the last cell at t = 23.9 s and t = 30.4 s for the operating temperature of 353.15 K and 343.15 K, respectively. The solution is unsteady and the flow distribution changes with different amounts of water flowing into different cells. As a result, the two parameters of F1 and F2 vary with the simulation time. The values of F1 and F2 are tabulated in Table 6 for various simulation

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

Fig. 13. Contours of liquid water volume fraction obtained from simulation of a 39-cell stack (N = 39, T = 353.15 K).

11

12

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

Fig. 14. Oxygen flow rate through cells for different times (N = 39, T = 353.15 K).

times. The flow uniformity parameters are influenced by two factors: inlet mass flow rate and amount of water in the cells. The water vapor mass fraction increases with the operating temperature, causing a higher mass flow rate and more flow non-uniformity. Table 6 shows that the values of non-uniformity parameters for the operating temperature of 343.15 K are sometimes higher due to the presence of liquid water in the cells. The values of F1 and F2 at t = 9 s show that the amount of liquid water flowing in the cells has a dominant effect on the flow uniformity parameters. 6.4. Effect of number of cells When the number of cells increases, the maximum local condensation rate remains constant, but the heat transfer area increases. The surface integral of vapor condensation rate over the manifold top wall shows that the rate of vapor condensation increases from 0.0158 kg/m.s to 0.021 kg/m.s as the number of cells changes from 26 to 39. The contour of liquid water volume fraction obtained from numerical simulation of a 39-cell stack is displayed in Fig. 13. As seen in Fig. 13, the initial droplets coalesce into larger ones (t = 27 s). The larger droplet in the middle of the manifold top wall moves backward and sweeps up small droplets. Thus, the droplet radius continuously increases with time (t = 29.4 s). When the gravity force overcomes the adhesion force (t = 29.7 s), the droplet detaches from the manifold wall and sink into seven cells. Flow distribution for various simulation times is plotted against the number of cells in Fig. 14. In this case, the mass flow rates in seven cells (10 to 16) drop dramatically, and the mass flow rates in other cells increase. Fig. 14 shows the blockage of the cells by slug formation in the cells and resultant severe flow maldistribution. It has been also found from Fig. 14 that the oxygen flow rate in 17.95% of cells (10 to16) reduced from 17.92% to 4.56% in the water removal process. The values of F1 and F2 for different times are calculated and listed in Table 7. By comparing Table 5 (t = 24.6 s) with Table 7 (t = 30 s), it can be seen that the value of F1 for the 26-cells and 39-cells stacks is approximately equal, while the value of F2 increases by more than 50%. Because of the F1 depends on the maximum and minimum values of mass flow rates among cells. Therefore, two cells with the maximum and minimum mass flow

Table 7 Flow distribution parameters for the 39-cell stack.

F1 F2

t = 29.7s

t = 29.8s

t = 29.9s

t = 30s

0.0377 1.37e-3

0.799 3.97e-2

0.859 5.62e-2

0.942 5.32e-2

rates are used to calculate the F1 . Mass flow rates passing through all cells in a stack are used in the calculation of the F2 . Thus the F2 is a better criterion to determine the severity of maldistribution in a fuel cell stack. 7. Conclusion In this paper, two phase flow in the manifold of a stack is simulated to investigate the effect of the water vapor condensation on the flow distribution. A new solver is developed in an open source software, OpenFOAM. The accuracy of the solution is verified by comparing the results of the present study with the experimental and numerical data found in the literature. First, the flow in the bipolar plate of a single fuel cell is simulated to determine the pressure drop through the gas flow channels. Then, flow channels are replaced by straight parallel channels filled with porous media. The numerical results show that the cells receive almost equal oxygen in single phase simulation. In the two phase simulation, the oxygen concentration in the mixture increases near the walls, due to the condensation of the water vapor. Therefore, the beginning and last cells receive more oxygen than the middle cells. The nonuniformity index and the root mean square value of the deviation from the mean flow rate are used to measure the flow maldistribution. Both parameters increase significantly in the two phase simulation. The condensed water forms small droplets on the inside surface of the manifold, and the small droplets coalesce into larger ones. When the gravity overcomes adhesion, the large droplet descends into cells and discharges from the stack. A large amount of liquid water in the cells forms slug and causes a significant increase in the flow maldistribution. The non-uniformity index increases from 0.0011 in the single flow simulation to 0.94 in the water removal process. Almost all liquid water drains through the

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

last cell in the simulation of the 26-cell stack while mass flow rate through seven cells affected by water removal process in the simulation of the 39-cell stack. A comparison between two flow distribution parameters under different conditions shows that the root mean square value of the deviation is a better criterion for measurement of the flow maldistribution. 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. CRediT authorship contribution statement Ahmad Rezaei Sangtabi: Methodology, Software, Formal analysis, Investigation, Data curation, Writing - original draft. Ali Kianifar: Writing - review & editing, Visualization, Supervision, Project administration. Ebrahim Alizadeh: Conceptualization, Resources, Funding acquisition. References [1] J. Wang, J. Yan, J. Yuan, B. Sundén, On flow maldistribution in PEMFC stacks, Int. J. Green Energy. 8 (2011) 585–606, doi:10.1080/15435075.2011.576288. [2] M. Abdollahzadeh, A.A. Ranjbar, Q. Esmaili, (1D + 1D) approach mathematical modeling of two phase multicomponent transport flow in PEMFC, Russ. J. Electrochem. 48 (2012) 1187–1196, doi:10.1134/S1023193512050023. [3] B.H. Lim, E.H. Majlan, W.R.W. Daud, M.I. Rosli, T. Husaini, Numerical analysis of flow distribution behavior in a proton exchange membrane fuel cell, Heliyon 4 (2018) e00845, doi:10.1016/j.heliyon.2018.e00845. [4] J.J. Baschuk, X. Li, Modelling of polymer electrolyte membrane fuel cell stacks based on a hydraulic network approach, Int. J. Energy Res. 28 (2004) 697–724, doi:10.1002/er.993. [5] J.-.H. Koh, H.-.K. Seo, C.G. Lee, Y.-.S. Yoo, H.C. Lim, Pressure and flow distribution in internal gas manifolds of a fuel-cell stack„ J. Power Sources 115 (2003) 54–65, doi:10.1016/S0378-7753(02)00615-8. [6] G. Karimi, J.J. Baschuk, X. Li, Performance analysis and optimization of PEM fuel cell stacks using flow network approach, J. Power Sources 147 (2005) 162– 177, doi:10.1016/j.jpowsour.2005.01.023. [7] S. Maharudrayya, S. Jayanti, A.P. Deshpande, Flow distribution and pressure drop in parallel-channel configurations of planar fuel cells, J. Power Sources 144 (2005) 94–106, doi:10.1016/j.jpowsour.2004.12.018. [8] Pressure drop and flow distribution in parallel-channel configurations of fuel cells: U-type arrangement, Int. J. Hydrogen Energy. 33 (2008) 6339–6350, doi:10.1016/j.ijhydene.2008.08.020. [9] Y. Qin, Q. Du, M. Fan, Y. Chang, Y. Yin, Study on the operating pressure effect on the performance of a proton exchange membrane fuel cell power system, Energy Convers. Manag. 142 (2017) 357–365, doi:10.1016/j.enconman.2017.03. 035. [10] M. Sajid Hossain, B. Shabani, C.P. Cheung, Enhanced gas flow uniformity across parallel channel cathode flow field of proton exchange membrane fuel cells, Int. J. Hydrogen Energy. 42 (2017) 5272–5283, doi:10.1016/j.ijhydene.2016.11. 042. [11] R. Mustata, L. Valiño, F. Barreras, M.I. Gil, A. Lozano, Study of the distribution of air flow in a proton exchange membrane fuel cell stack, J. Power Sources 192 (2009) 185–189, doi:10.1016/j.jpowsour.2008.12.083. [12] S.Y. Kim, W.N. Kim, Effect of cathode inlet manifold configuration on performance of 10-cell proton-exchange membrane fuel cell, J. Power Sources 166 (2007) 430–434, doi:10.1016/j.jpowsour.2006.12.104. [13] C.-.H. Chen, S.-.P. Jung, S.-.C. Yen, Flow distribution in the manifold of PEM fuel cell stack, J. Power Sources 173 (2007) 249–263, doi:10.1016/j.jpowsour. 20 07.05.0 07. [14] J. Lebæk, M.B. Andreasen, H.A. Andresen, M. Bang, S.K. Kær, Particle image velocimetry and computational fluid dynamics analysis of fuel cell manifold, J. Fuel Cell Sci. Technol. 7 (2010) 031001, doi:10.1115/1.3206697. [15] J. Lebæk, M. Bang, S.K. Kær, Flow and pressure distribution in fuel cell manifolds, J. Fuel Cell Sci. Technol. 7 (2010) 061001, doi:10.1115/1.4001319. [16] H.-.H. Liu, C.-.H. Cheng, K.-.L. Hsueh, C.-.W. Hong, Modeling and design of airside manifolds and measurement on an industrial 5-kW hydrogen fuel cell stack, Int. J. Hydrogen Energy. 42 (2017) 19216–19226, doi:10.1016/j.ijhydene. 2017.06.057. [17] J. Dong, X. Xu, B. Xu, CFD analysis of a novel modular manifold with multistage channels for uniform air distribution in a fuel cell stack, Appl. Therm. Eng. 124 (2017) 286–293, doi:10.1016/j.applthermaleng.2017.06.030. [18] J.M. Jackson, M.L. Hupert, S.A. Soper, Discrete geometry optimization for reducing flow non-uniformity, asymmetry, and parasitic minor loss pressure drops in Z-type configurations of fuel cells, j, Power Sources 269 (2014) 274–283, doi:10.1016/j.jpowsour.2014.06.136.

13

[19] E. Alizadeh, M. Khorshidian, S.H.M. Saadat, S.M. Rahgoshay, M. Rahimi-Esbo, The experimental analysis of a dead-end H 2 /O 2 PEM fuel cell stack with cascade type design, Int. J. Hydrogen Energy. 42 (2017) 11662–11672, doi:10. 1016/j.ijhydene.2017.03.094. [20] A. Ramiar, A.H. Mahmoudi, Q. Esmaili, M. Abdollahzadeh, Influence of cathode flow pulsation on performance of proton exchange membrane fuel cell with interdigitated gas distributors, Energy 94 (2016) 206–217, doi:10.1016/j.energy. 2015.10.110. [21] S.M. Rahgoshay, A.A. Ranjbar, A. Ramiar, E. Alizadeh, Thermal investigation of a PEM fuel cell with cooling flow field, Energy 134 (2017) 61–73, doi:10.1016/ j.energy.2017.05.151. [22] M. Bilgili, M. Bosomoiu, G. Tsotridis, Gas flow field with obstacles for PEM fuel cells at different operating conditions, Int. J. Hydrogen Energy. 40 (2015) 2303– 2311, doi:10.1016/j.ijhydene.2014.11.139. [23] S. Ebrahimi, B. Ghorbani, K. Vijayaraghavan, Optimization of catalyst distribution along PEMFC channel through a numerical two-phase model and genetic algorithm, Renew. Energy. 113 (2017) 846–854, doi:10.1016/j.renene.2017. 06.067. [24] J.H. Jo, W.T. Kim, Numerical simulation of water droplet dynamics in a right angle gas channel of a polymer electrolyte membrane fuel cell, Int. J. Hydrogen Energy. 40 (2015) 8368–8383, doi:10.1016/j.ijhydene.2015.04.122. [25] Y. Qin, Y. Yin, K. Jiao, Q. Du, Effects of needle orientation and gas velocity on water transport and removal in a modified PEMFC gas flow channel having a hydrophilic needle, Int. J. Energy Res. (2018), doi:10.1002/er.4116. [26] A. Soopee, A.P. Sasmito, T. Shamim, Water droplet dynamics in a dead-end anode proton exchange membrane fuel cell, Appl. Energy. 233–234 (2019) 300– 311, doi:10.1016/j.apenergy.2018.10.001. [27] Y. Hou, G. Zhang, Y. Qin, Q. Du, K. Jiao, Numerical simulation of gas liquid twophase flow in anode channel of low-temperature fuel cells, Int. J. Hydrogen Energy. 42 (2017) 3250–3258, doi:10.1016/j.ijhydene.2016.09.219. [28] Y. Ding, X.T. Bi, D.P. Wilkinson, Numerical investigation of the impact of twophase flow maldistribution on PEM fuel cell performance, Int. J. Hydrogen Energy. 39 (2014) 469–480, doi:10.1016/j.ijhydene.2013.10.047. [29] J. Kim, G. Luo, C.-.Y. Wang, Modeling two-phase flow in three-dimensional complex flow-fields of proton exchange membrane fuel cells, J. Power Sources 365 (2017) 419–429, doi:10.1016/j.jpowsour.2017.09.003. [30] E. Alizadeh, M. Rahimi-Esbo, S.M. Rahgoshay, S.H.M. Saadat, M. Khorshidian, Numerical and experimental investigation of cascade type serpentine flow field of reactant gases for improving performance of PEM fuel cell, Int. J. Hydrogen Energy. 42 (2017) 14708–14724, doi:10.1016/j.ijhydene.2017.04.212. [31] P. Ribeirinha, M. Abdollahzadeh, J.M. Sousa, M. Boaventura, A. Mendes, Modelling of a high-temperature polymer electrolyte membrane fuel cell integrated with a methanol steam reformer cell, Appl. Energy. 202 (2017) 6–19, doi:10.1016/j.apenergy.2017.05.120. [32] A. Ghanbarian, M.J. Kermani, J. Scholta, M. Abdollahzadeh, Polymer electrolyte membrane fuel cell flow field design criteria – Application to parallel serpentine flow patterns, Energy Convers. Manag. 166 (2018) 281–296, doi:10.1016/j. enconman.2018.04.018. [33] T. Wilberforce, O. Ijaodola, F.N. Khatib, E.O. Ogungbemi, Z. El Hassan, J. Thompson, A.G. Olabi, Effect of humidification of reactive gases on the performance of a proton exchange membrane fuel cell, Sci. Total Environ. 688 (2019) 1016– 1035, doi:10.1016/j.scitotenv.2019.06.397. [34] A. Kolahan, E. Roohi, M.-.R. Pendar, Wavelet analysis and frequency spectrum of cloud cavitation around a sphere, Ocean Eng 182 (2019) 235–247, doi:10. 1016/j.oceaneng.2019.04.070. [35] M.-.R. Pendar, E. Roohi, Investigation of cavitation around 3D hemispherical head-form body and conical cavitators using different turbulence and cavitation models, Ocean Eng 112 (2016) 287–306, doi:10.1016/j.oceaneng.2015.12. 010. [36] P.K. Jithesh, A.S. Bansode, T. Sundararajan, S.K. Das, The effect of flow distributors on the liquid water distribution and performance of a PEM fuel cell, Int. J. Hydrogen Energy. 37 (2012) 17158–17171, doi:10.1016/j.ijhydene.2012.08.058. [37] J... Brackbill, D... Kothe, C. Zemach, A continuum method for modeling surface tension„ J. Comput. Phys. 100 (1992) 335–354, doi:10.1016/0021-9991(92) 90240-Y. [38] N. Khajeh-Hosseini-Dalasm, K. Fushinobu, K. Okazaki, Phase change in the cathode side of a proton exchange membrane fuel cell„ J. Power Sources 195 (2010) 7003–7010, doi:10.1016/j.jpowsour.2010.04.089. [39] F. Barbir, in: Fuel Cell Modeling, in: PEM Fuel Cells, Elsevier, 2013, pp. 217–263, doi:10.1016/B978- 0- 12- 387710- 9.0 0 0 07-2. [40] M.-.R. Pendar, E. Roohi, Cavitation characteristics around a sphere: an LES investigation, Int. J. Multiph. Flow. 98 (2018) 1–23, doi:10.1016/j. ijmultiphaseflow.2017.08.013. [41] Y. Ding, H.T. Bi, D.P. Wilkinson, Three dimensional numerical simulation of gas–liquid two-phase flow patterns in a polymer–electrolyte membrane fuel cells gas flow channel„ J. Power Sources. 196 (2011) 6284–6292, doi:10.1016/j. jpowsour.2011.03.100. [42] J.H. Kim, W.T. Kim, Numerical investigation of gas-liquid two-phase flow inside PEMFC gas channels with rectangular and trapezoidal cross sections, Energies 11 (2018) 1403, doi:10.3390/en11061403. [43] A. Rezaei Sangtabi, A. Ramiar, A.A. Ranjbar, M. Abdollahzadeh, A. Kianifar, Influence of repetitive laser pulse energy depositions on supersonic flow over a sphere, cone and oblate spheroid, Aerosp. Sci. Technol. 76 (2018) 72–81, doi:10.1016/j.ast.2018.02.005. [44] B. van Leer, Towards the ultimate conservative difference scheme. V. a secondorder sequel to Godunov’s method, J. Comput. Phys. 32 (1979) 101–136, doi:10. 1016/0021-9991(79)90145-1.

14

A.R. Sangtabi, A. Kianifar and E. Alizadeh / International Journal of Heat and Mass Transfer 151 (2020) 119471

[45] X. Zhu, P.C. Sui, N. Djilali, Dynamic behaviour of liquid water emerging from a GDL pore into a PEMFC gas flow channel, J. Power Sources 172 (2007) 287– 295, doi:10.1016/j.jpowsour.2007.07.024. [46] Z. Lu, S.G. Kandlikar, C. Rath, M. Grimm, W. Domigan, A.D. White, M. Hardbarger, J.P. Owejan, T.A. Trabold, Water management studies in PEM fuel cells,

Part II: ex situ investigation of flow maldistribution, pressure drop and twophase flow pattern in gas channels, Int. J. Hydrogen Energy. 34 (2009) 3445– 3456, doi:10.1016/j.ijhydene.2008.12.025. [47] E. Grolman, J.M.H. Fortuin, and flow patterns in inclined gas-liquid pipe flow, Exp. Therm. Fluid Sci. 15 (1997) 174–182, doi:10.1016/S0894-1777(97)0 0 021-6.