Investigation of current density spatial distribution in PEM fuel cells using a comprehensively validated multi-phase non-isothermal model

International Journal of Heat and Mass Transfer 150 (2020) 119294

Contents lists available at ScienceDirect

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

Investigation of current density spatial distribution in PEM fuel cells using a comprehensively validated multi-phase non-isothermal model Guobin Zhang a,b, Jingtian Wu b, Yun Wang b,∗, Yan Yin a, Kui Jiao a,∗ a b

State Key Laboratory of Engines, Tianjin University, 135 Yaguan Road, Tianjin, 300350, China Renewable Energy Resources Laboratory, Department of Mechanical and Aerospace Engineering, University of California, Irvine, CA 92697-3975, USA

a r t i c l e

i n f o

Article history: Received 4 August 2019 Revised 27 December 2019 Accepted 29 December 2019

Keywords: Proton exchange membrane fuel cell Multi-phase model Temperature Validation Current density distribution

a b s t r a c t In this study, a three-dimensional (3D) multi-phase non-isothermal model of proton exchange membrane (PEM) fuel cell is present to investigate the current density spatial distributions under different output voltages, temperatures, current densities and relative humidities (RH). Two sets of experimental data are selected for validation, including those from the Los Alamos National Laboratory in the United States and the University of Waterloo in Canada. Reasonable agreements are achieved between the model prediction and experimental measurements, indicative of the validity of this 3D model. In addition, it is found that under a higher RH of 50%, most electric current is produced near the cathode inlet, which is primarily due to the local availability of abundant oxygen and hence small transport polarization. Low current occurs near the outlet of the cathode air flow, as a result of oxygen consumption by the oxygen reduction reaction (ORR). Under a lower RH of 25%, the high current density region shifts to the middle of the fuel cell, which is primarily attributed to hydration of the dry membrane by water production. In our case study, the operating temperature has little impact on the current density distribution. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction Nowadays, the carbon-based energy sources, such as: coals, crude oils and natural gasses, are still the major portion of energy consumption in the world, which cause severe pollution and global warming due to the release of carbon dioxide (CO2 ) [1]. Proton exchange membrane (PEM) fuel cell, whose by-product is only water, is widely regarded as a promising candidate in future automobile application thanks to its zero emissions, fast dynamic response and high power density [2–4]. Modeling plays an important role in fuel cell design and development for commercialization [5,6]. A multidimensional multi-physics non-isothermal model predicts the distributions of oxygen, vapor, liquid water, temperature, and current density in a fuel cell, which is of vital importance to the performance and durability improvement as well as cost reduction. In PEM fuel cell, the electrochemical and transport phenomena are complex, which involve gas and liquid two-phase flow, gas species transport, water condensation and evaporation, membrane water absorption and release, proton and electron transport

∗

Corresponding authors. E-mail addresses: [email protected] (Y. Wang), [email protected] (K. Jiao).

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

and heat transfer, etc. In addition, these transport phenomena mentioned above are always highly coupled [7,8]. For example, the temperature distribution will greatly influence local water content, meanwhile, the local high current density leads to high heat production rate and increase local temperature. Therefore, developing reliable analytical or numerical models has been an important way to predict PEM fuel cell performance and understand these highly-coupled transport phenomena. In order to take these transport phenomena into account, a series of conservation equations, including mass, momentum, energy, and species, need to be included in the model. Three-dimensional (3D) multi-phase description is necessary for high-fidelity prediction in comparison with other reduced-dimension models [9], e.g. 0D [10], 1D [11] and 2D [12,13] models, owing to inclusion of the realistic geometry structure [14–16] and liquid water impacts. One important aspect of modeling is validation, which compares the computational prediction with experimental data. The polarization curve is widely used to validate fuel cell models [17,18]. For more reliable validation, two or more polarization curves at different conditions can be used for comparison [19,20]. Penga et al. [21] compared two polarization curves for constant and spatially variable temperature heat removal rates at the outer surface of BP with corresponding experimental data. In addition, the polarization curve, along with three elementary voltage losses,

2

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

Nomenclature A a aPt C Cp D Eagg EW F h H I iref 0 J Jion J0ref K k kc M m mC mPt N P Pc R RO2 RatioPtC ragg RH S Sp

S s T t u Yi

Area (m2 ) Water activity The specific platinum surface area per unit catalyst volume (m2 m−3 ) Molar concentration (mol m−3 ) Heat capacity (J kg−1 K−1 ) Diffusion coefficient (m2 s−1 ) Effectiveness factor of agglomerate Equivalent weight of membrane The Faraday’s constant (96,487 C mol−1 ) Latent heat of water (J mol−1 ) Henry’s coefficient (Pa m3 mol−1 ) Current density (A m−2 ) Exchange current density (A m−2 ) Volumetric reaction rate (A m−3 ) Ionic current density (A m−2 ) Reference volumetric current density (A m−3 ) Intrinsic permeability (m2 ) Relative permeability or thermal conductivity (W m−1 K−1 ) Electrochemical reaction constant (s−1 ) Molar mass (kg mol−1 ) Mass flow rates (kg s−1 ) Carbon loading (kg m−2 ) Platinum loading (kg m−2 ) Number of agglomerates per unit volume of CL (m−3 ) Pressure (Pa) Capillary pressure (Pa) Universal gas constant (J mol−1 K−1 ) Oxygen reduction reaction rate (mol m−3 s−1 ) Platinum ratio Agglomerate radius (m) Relative humidity Source term (kg m−3 s−1 or mol m−3 s−1 ) or reaction surface area per unit platinum mass (m2 kg−1 ) Source term of membrane water caused by pressure difference (mol m−3 s−1 ) Entropy change (J mol−1 K−1 ) Liquid saturation Temperature (K) Time (s) Velocity (m s−1 ) Gas species mass fraction

Greek letters α Transfer coefficient γ Phase change coefficient (s−1 ) δ Thickness (m) ɛ Porosity ɛPt Effective platinum surface ratio ɛPt/C Volume fraction of carbon supporting platinum particle η Overpotential (V) θ Contact angle (°) κe Electric conductivity (S m−1 ) κ ion Ionic conductivity (S m−1 ) λ Membrane water content μ Dynamic viscosity (kg m−1 s−1 ) ξ Stoichiometric ratio ρ Density (kg m−3 ) σ Surface tension coefficient (N m−1 )

ϕ ϕe ϕ ion χ ω ωagg

Thiele’s modulus for spherical agglomerate Electric potential (V) Proton potential (V) An auxiliary variable Electrolyte volume fraction Electrolyte volume fraction in agglomerate

Subscripts and superscripts 0 Standard state a Anode act Activation state agg Agglomerate c Cathode CL Catalyst layer d Membrane water d-v Membrane water to vapor e Electrical eq Equivalent state eff Effective GDL Gas diffusion layer g Gas phase H Henry’s law H2 Hydrogen H2 O Water i Gas species ion Ionic l Liquid phase m Mass mem Membrane mw Membrane water N Diffusion in electrolyte O2 Oxygen out Outlet Pt Platinum ref Reference state rev Reversible sat Saturation state T Temperature u Momentum v-l Vapor to liquid phase

i.e. activation, Ohmic and concentration losses, has been used in validation. Hao et al. [22] validated the polarization curve and Ohmic loss for their model at the platinum loading ranging from 0.025 to 0.2 mg cm−2 . Zhang et al. [15] compared two polarization curves at different electrolyte volume fractions and polarization curve together with Ohmic loss with experimental data selected from two independent groups. However, only comparing the polarization curve or elementary voltage loss is inadequate for 3D model validation, and the local quantities may be very different even under similar polarization curves [23]. For 3D models, more comprehensive validation is required, which compares the spatial distributions of key variables between model prediction and experimental data. The current density distribution is one of the most important parameters to characterize fuel cell performance, which provides local reaction activity and electrochemical consequence resulting from local reactant concentrations, temperature, liquid water, and materials. In the past two decades, the measurement of current density distribution has been achieved by segmented cell techniques [24], such as printed circuit board (PCB) [25], resistors network [26] and Hall effect sensors [27], which require to modify cell components including current collector, bipolar plate (BP) or gas diffusion layer (GDL). Besides, some non-invasive approaches, e.g. magnetotomography [28], have

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

3

Table 1 Governing equations. Description Mass (gas phase)

Momentum (gas phase)

Gas species Liquid saturation Liquid pressure Membrane water content Electric potential Ionic potential Energy

Equations

Components

∂ (ε (1 − s )ρg ) + ∇ · (ρg ug ) = Sm ∂t        ug ug T ρg ug ug ∂  ρg ug  +∇ · = −∇ Pg + μg ∇ · ∇ +∇ 2 2 ∂ t ε (1 − s ) ε ( 1 − s ) ε ( 1 − s ) ε (1 − s )   u  2 g − μg ∇ ∇ · + Su 3 ε (1 − s )   ∂ (ε (1 − s )ρgYi ) + ∇ · (ρg ugYi ) = ∇ · ρg Deff i ∇ Yi + Si ∂t ∂ (ρ s ) + ∇ · (ρl ul s ) = Sl , ul = ug ∂t l   ∂ Kk (ρ ε s ) = ∇ · ρl l ∇ Pl + Sl ∂t l μl  J  ρ   ρmem ∂ (ωλ ) + ∇ · nd ion = mem ∇ · Deff d ∇λ + Smw EW ∂ t F EW   0 = ∇ · κeeff ∇ ϕe + Se  eff  ∇ ϕion + Sion 0 = ∇ · κion

Flow fields, GDLs, MPLs and CLs

Flow fields, GDLs, MPLs and CLs

Flow fields, GDLs, MPLs and CLs Flow fields GDLs, MPLs and CLs Membrane and CLs BPs, GDLs, MPLs and CLs Membrane and CLs

∂ (ε sρlCp,l T + ε (1 − s )ρgCp,g T ) + ∇ · (ε sρlCp,l ul T + ε (1 − s )ρgCp,g ug T ) = ∇ · (keff ∇ T ) + ST ∂t

been proposed without segmentation of the original cell structure. These experimental studies not only improve understanding of internal operating condition in PEM fuel cell but also provide important data for model validation. As for the validation of 3D models, Carnes et al. [29] validated their 3D multi-phase model in term of the current density distribution at the outer surface of the cathode bipolar plate (BP). The experimental data were collected from a single 50 cm2 cell with 10 × 10 segmented current collectors. Similar validation of current density spatial distribution was performed by Ju et al. [30]. Li et al. [31] conducted the current density validation in the in-plane direction between channel and land and achieved good agreement at a submillimeter resolution for both air and oxygen cathode flows. In addition, liquid water distribution was also used in validation. Hao et al. [32] compared the 3D multi-phase model predictions with experimental data for a large-scale cell with 25 × 4 segmented current collectors, including the current density, temperature and liquid water distributions in the in-plane direction. Iranzo et al. [33] validated their 3D multi-phase model using the in-plane liquid water distribution obtained by neutron radiography. In addition, comparison of the predicted through-plane and in-plane liquid distributions and neutron/X-Ray radiography data were conducted by Wang and Chen [34,35]. Though detailed validation has been attempted by several groups, few were reported to compare with experimental data of the current density distribution from various groups and for different fuel cell designs. In this study, a 3D multi-phase model of PEM fuel cell is comprehensively validated against the experimental current density distributions from two well-known research groups with one fuel cell of 50 cm2 and the other of 40 cm2 . The current density distributions at different output voltages, temperatures, current densities and humidification conditions are favorably compared with the experimental data. 2. Numerical model

All the components

and liquid water equations are separated with the condensation/evaporation given as a source term for each of the equations [9]. Specifically, mass and momentum equations describe the gas flows. Gaseous reactant transport is described by the species equation. For the liquid phase in porous components such as GDL, MPL, and CL, its flow is primarily driven by the capillary pressure, which is described by the liquid pressure equation. Besides, an individual governing equation is added for the water content in membrane and CL considering that its transport characteristic is different from vapor and liquid water [4]. The electric and ionic potential equations describe the electron and proton transports, respectively. In the energy equation, the fuel cell heat generation is added as the source term. The governing equations and source terms are shown in Tables 2 and 3, respectively, with the model parameters listed in Table 4. It should be noted that the present model is a transient model, but only steady simulations are conducted in this study by neglecting the transient terms in all the governing equations. 2.2. Agglomerate model for CLs In order to take the complex structures of CL [36] into account, the spherical agglomerate description is adopted, which expresses the volume fraction of carbon supporting platinum (ɛPt/C , conducting electrons), electrolyte volume fraction (ω, conducting proton) and porosity (ɛCL , for mass transport) [37] by:

εPt/C = ω= N=

mPt



δCL

1

ρPt

+

1 − RatioPt/C 1 RatioPt/C ρC



(1)

3 4 3 π N ragg ωagg + (ragg + δm )3 − ragg , 3 4 3

π

εPt/C (1 − ωagg )

3 ragg

(2)

2.1. Governing equations

εCL = 1 − εPt/C − ω

The 3D non-isothermal multi-phase model is based on first principle conservation laws, which are similar to many other 3D modeling studies [5,6,21–23,32]. It consists of a series of governing equations, as listed in Table 1. The two-fluid model is adopted to describe the gas and liquid two-phase flow, in which the gas

where mPt (kg m−2 ) is the platinum loading, δ CL and δ m (m) the thickness of CL and electrolyte film in CLs, RatioPt/C = mPt /(mPt + mC ) the platinum ratio, N(m−3 ) the number of agglomerates per unit volume of CL, ragg (m) the agglomerate radius, and ωagg the electrolyte volume fraction in agglomerate.

(3)

4

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

Table 2 Physical and transport properties. Parameters

Expression

Effective electric conductivity in MPL (S m Effective electric conductivity in CL (S m

−1

−1

)

)

Effective electric conductivity in GDL (S m−1 )

Diffusion coefficient of membrane water (m2 s−1 )

κeeff = κe (1 − ε )1.5 κeeff = κe εPt/C 1.5 ⎧   3ε −0.016 ⎪ exp(0.367(1 − ε )) κeeff,GDL ⎨1 − 0.962(1 − ε ) 3 − (1 − ε ) =   3ε κe,GDL ⎪ ⎩1 − 0.962(1 − ε )−0.007 exp(0.889(1 − ε )) 3 − (1 − ε )  −2346  ⎧ 0<λ<3 3.1 × 10−7 λ(exp(0.28λ ) − 1 ) exp ⎪ ⎪ T ⎪ ⎪ ⎪  −2346  ⎨ −8 3 ≤ λ < 17 Dd = 4.17 × 10 λ(161 exp(−λ ) + 1 ) exp T ⎪ ⎪  0.15    ⎪ ⎪ λ − 2.5 ⎪ ⎩4.1 × 10−10 λ λ ≥ 17 1 + tanh

in − plane through − plane

1.4

25

Gas diffusivities (m2 s−1 )

κion = (0.5139λ − 0.326)exp[1268(1/303.15 − 1/T )]      (ωagg − 1 ) δm 1/3 eff κion , χ = min 0, = (1 − εCL ) 1 + κion + (1 − ωagg ) − 1 3 ragg (1 + δm /ragg + χ ) = ε 1.5 Di Deff i

Relative permeability

kg = ( 1 − s )3.0 , kl = s3.0

Equilibrium membrane water content

λeq =

Water activity in CL

a = 2s + CH2 O RT /Psat

Water saturation pressure (Pa)

log10

Proton conductivity (S m

−1

)

Effective proton conductivity in CL (S m−1 )



0.043 + 17.81a − 39.85a2 + 36.0a3 14.0 + 1.4(a − 1.0)

 P  sat 101325

= −2.1794 + 0.02953(T − 273.15) − 9.1837 × 10−5 (T − 273.15)2 + 1.4454 × 10−7 (T − 273.15)3



Oxygen diffusion coefficient in electrolyte (m2 s−1 ) Henry’s law

Pc = σ cos θ

Leverett-J function



1 RT ln PHin2 + ln POin2 2F 2  T − 273.15  − 1.6461 × 10−10 λ0.708 + 5.2 × 10−10 DO2 ,N = 1.3926 × 10−10 λ0.708 exp 106.65 PH2 PO2 , COH2 = CHH2 = HH2 HO2 Erev = 1.229 − 0.9 × 10−3 (T − 298.15) +

Reversible voltage (V)

0≤a≤1 1
 ε 0 . 5 K



J (s ), J (s ) =

1.42(1 − s ) − 2.12(1 − s ) + 1.26(1 − s ) 1.42s − 2.12s2 + 1.26s3 2

3

θ < 90◦ θ > 90◦

Table 3 Source terms of the conservation equations. Source terms

Expressions

Hydrogen (kg m−3 s−1 )

SH2 = −MH2 Ja /2F, ACL

Oxygen (kg m−3 s−1 )

SO2 = −MO2 Jc /4F, CCL

Water vapor (kg m

−3

s

−1

SH2 O =

)

⎧ ⎨−Sv−l

GDL, MPL ACL CCL

−Sv−l + Sd−v v−l + Sd−v + MH2 O Jc /2F

⎩−S

Liquid water (kg m−3 s−1 )

Sl = Sv−l , flow fields, GDLs, MPLs and CLs

Gas mixture (kg m−3 s−1 )

Sg =

Water evaporation and condensation (kg m−3 s−1 )

Sv−l =

Electric and ionic potential (A m−3 )

Se =

Membrane water content (mol m−3 s−1 )

Smw =

⎧ ⎨−Sv−l + Sd−v − Ja MH2 /2F − Jc MO2 /4F −S + S ⎩−Sv−l d−v v−l



γv−l ε (1 − s)(CH2 O − Csat )MH2 O γl−v ε s(CH2 O − Csat )MH2 O   −Ja Jc



Membrane water absorption and release (kg m−3 s−1 ) Membrane water caused by pressure (mol m

Energy (W m−3 )

−3

s

−1

ACL CCL Flow fields, GDLs and MPLs

)

ACL S = CCL ion

−Sd−v /MH2 O − Sp −Sd−v /MH2 O + Sp

Ja −Jc

CH2 O >Csat CH2 O
ACL CCL

ACL CCL

Sd−v = γd−v ρmem /EW(λ − λeq )MH2 O

ρl Kmem (Pla − Plc ) μl MH2 O δmem δCL ⎧ ∇ϕe 2 κeeff ⎪ ⎪ ⎪ ⎪∇ϕ 2 κ eff + S h ⎪ ⎪ e v−l e ⎪ ⎪ ⎪ ⎪ ⎨Ja |ηa | + ∇ϕe 2 κ eff + ∇ϕ 2 κ eff + Ja Sa T + (S − S )h ion v−l d−v e act ion 2F ST = Sc T ⎪ c 2 eff 2 eff ⎪ + J + | η | + ∇ϕ  κ + ∇ϕ  κ ( S − S J ⎪ c e c ion v−l d−v )h e act ion ⎪ 4F ⎪ ⎪ ⎪∇ϕion 2 κ eff ⎪ ion ⎪ ⎪ ⎩

Sp =

Sv−l h

BPs GDLs, MPLs ACL CCL Membrane Flow fields

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

5

Table 4 Fuel cell and model parameters. Parameter

Value

Permeability (m2 ) (GDL, MPL, CL, Membrane) Effective platinum surface ratio Agglomerate radius in CL (μm) Electrolyte fraction in agglomerate Reference hydrogen concentration (mol m−3 ) Reference oxygen concentration (mol m−3 ) Transfer coefficient Dry membrane density (kg m−3 ) Membrane equivalent weight (kg mol−1 ) Water condensation and evaporation rate (s−1 ) Latent heat of water evaporation (J mol−1 ) Liquid density (kg m−3 ) Henry’s constant (Pa m3 mol−1 ) Platinum density (kg m−3 ) Carbon density (kg m−3 ) Entropy change (J mol−1 K−1 ) Faraday’s constant (C mol−1 ) Universal gas constant (J mol−1 K−1 ) Surface tension coefficient (N m−1 ) Membrane water absorption and release rate (s−1 ) Thermal conductivity (W m−1 K−1 ) (BP, MPL, CL, Membrane) GDL thermal conductivity (W m−1 K−1 )

2.0e-12, 1.0e-12, 1.0e-13, 2.0e-20 0.5 1.0 0.226 56.4 3.39 0.5 1980 1.1 100 44,900 970 Hydrogen: 4.56e3, Oxygen: 2.8e4 21,450 1800 Anode: 130.68, Cathode: 32.55 96,487 8.314 0.0625 1.3 20, 1, 1, 0.95 In-plane: 21, Through-plane: 1.7

The anode and cathode reaction rates in CLs, Ja and Jc (A m−3 ), are given by:



eff iref 0,a iT,a aPt,a

Ja =



HH2 CHref2

 2F α a

× exp −

Jc = 4F

PO2 HO2





iT ,a

0.5

PH2

  2F 1 − α  ( a) a a ηact − exp − ηact

RT

RT

1 (ragg + δm )δm + Eagg kc (1 − εCL ) aagg ragg DO2 ,N

1

1 = exp −1400 − T 353.15

−1 (5)

4F ( 1 − ε )

ref CL CO2



Eagg =

1



RT

CO2 =

  1 − α 4F  ( c) c c ηact − exp − ηact RT

1

−

1



iref (A 0

1 T

−

(7)

 (8)

kc DO2,N ωagg 1.5

iT,c = exp −7900

(9) 1 353.15





mPt

δCL

S

(11)

S = 227.79Ratio3Pt/C − 158.57Ratio2Pt/C − 201.53RatioPt/C + 159.5 ×103



(12) S(m2

kg−1 )

(13)

2F CH2 ρgc Iξc Aact 4F CO2 Pga,out + Pga − RHa P sat RT  0.21 Pgc,out + Pgc − RHc P sat



RHa,c P RT MiCi Yi =  MiCi

RT sat

= RHa,cCHsat2 O

(14) (15) (16) (17) (18)

In addition, the constant temperatures are specified at all the surrounding surfaces of the computational domain. For the solution of electric potential equation, the reference voltage (0 V) is implemented at the end surface of cathode BP, while at the anode, the constant current density [20,39] or overpotential [40] is defined for the potentiostatic or galvanostatic mode, respectively. 2.4. Numerical implementation

where denotes the exchange current density, H(Pa m3 mol−1 ) Henry’s coefficient, and ηact (V) activation overpotential. The subscripts a and c denote anode and cathode, respectively. The effective reaction surface area per unit platinum mass (aeff , m2 m−3 ) Pt is calculated as:

εPt

ρga Iξa Aact

(10)

m−2 )

aeff Pt = εPt aPt =

mc =

(6)

ϕ tanh (3ϕ ) 3ϕ  0.5

ragg 3

ma =

CHa,2cO =

 4α F c

× exp −

For gas flow, the mass flow rate (m, kg s−1 ) together with species mass fraction and flow pressure are specified at the inlet and outlet [38], respectively,

CH2 =



eff iref 0,c aPt,c iT,c

kc =

ϕ=

(4)

2.3. Boundary conditions

where ɛPt and are the effective platinum surface ratio and the reaction surface area, respectively.

The 3D multi-phase model was implemented in the commercial CFD software, ANSYS FLUENT, through its user-defined function (UDF) interface. The grid independence test was conducted to ensure the current computational domain is adequate for fuel cell study in our previous work [15]. For example, two different grid cells in the in-plane direction, 168 × 100 and 210 × 125 were tested and there were negligible differences in the main simulation results, e.g. performance, oxygen concentration and liquid saturation. Totally 3.06 and 11.66 million cells were employed in the computational domains of PEM fuel cell (I) and (II), respectively. Simulations were conducted at a small workstation with 24 processors (Intel (R) Xeon (R) CPU E5-2620 v3 @ 2.40 GHz) and 32GB DDR RAM. In general, about several to tens thousands iteration steps are essential to obtain convergence results.

6

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

3. Results and discussion 3.1. Polarization curve comparison In our previous work, the polarization curve and Ohmic loss calculated by this model have been compared with experimental data (Fig. 2 in Ref. [15]). Incorporation of a CL agglomerate model, species transport equations, and flow equations make it possible to predict the concentration and kinetic losses accurately [41,42] and agrees better with experimental results. In fact, the more accurate prediction of concentration loss is a significant advantage over early-stage 3D models [6,43]. 3.2. Current density distribution validation In this section, the predicted current density distribution was first compared with the experimental data reported by the Los Alamos National Laboratory in the United States using their fuel cell configuration [29]. The schematic of the fuel cell is shown in Fig. 1 (PEM fuel cell I). Experimental data are available for the current density distributions at different current densities, temperatures and relative humidities. And the operation conditions and geometry parameters are shown in Table 5. Fig. 2 compares two

Fig. 2. Comparison of the polarization curves (353.15 K) for PEM fuel cell (I) [29].

Fig. 1. Schematic of PEM fuel cell (I) [29] and (II) [45].

Table 5 Fuel cell geometry and validation parameters. Parameters 2

MEA area (cm ) Flow field dimension (mm × mm) Channel height (mm) Channel width (mm) Manifold width (mm) BP height (mm) Land width (mm) Thickness (Membrane, ACL, CCL, MPL, GDL) (μm) Stoichiometry ratio Operation pressure (atm) Relative humidity Operating temperature (K) Porosity (GDL, MPL) Electrical conductivity (S m−1 ) (MPL, CL) GDL electrical conductivity (S m−1 ) Contact angle (GDL, MPL, CL) (°) Platinum loading in CLs (mg cm−2 ) Pt/C ratio in CLs Electrolyte film thickness in CLs (nm)

PEM fuel cell (I)

PEM fuel cell (II)

50.41 71 × 71 1.0 1.0 3.0 3.0 1.03 18.0, 7.0, 12.0, 50.0, 200.0 Anode: 1.2, Cathode: 2.0 2.72 0.5 or 0.25 353.15 or 333.15 0.6, 0.6 5000, 5000 In-plane: 5000, Through-plane: 500 110, 120, 98 Anode: 0.2, Cathode: 0.4 0.3 150

43.56 66 × 66 1.0 1.0 – 2.0 1.0 53.0, 10.0, 10.0, -, 230.0 Anode: 1.2, Cathode: 2.0 1.49 1.0 338.15 0.75, -, 3000 In-plane: 943.42, Through-plane: 59.44 120, -, 100 Anode: 0.5, Cathode: 0.5 0.35 120

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

7

Fig. 3. Comparison of the predicted current density distributions with experimental data (353.15 K, 50% RH) under three current densities for PEM fuel cell (I) [29].

polarization curves at different relative humidities (353.15 K, 25% RH, 50% RH) between simulation and experiment. The error bars in Fig. 2 represent ± 5% error and the mean absolute error values for 25% and 50% RHs are 20.6 and 9.8 mV, respectively, indicating that the agreements are good. In our study, we also had a few computational grids across the membrane to better capture the local

Ohmic resistance. Again, all the three physical dimensions and important physical-chemical processes are taken into account in the simulation, which contributes to the successful validations. Comparison of the current density distributions at 0.1, 1.0 and 1.2 A cm−2 are shown in Fig. 3, which correspond to the three data points in Fig. 2 (353.15 K, 50% RH). The current density

8

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

Fig. 4. (a): Comparison of the local current density from inlet to outlet (or from top to bottom) and (b): absolute error under three current densities for PEM fuel cell (I) [29] (353.15 K, 50% RH).

Fig. 5. (a): Oxygen molar concentration; (b): liquid water saturation; (c): membrane water content; and (d): temperature distribution in the middle plane of cathode CL (353.15 K, 50% RH, 1.2 A cm−2 ) for PEM fuel cell (I).

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

distribution in experiment was obtained by 10 × 10 segmented cathode BPs. For comparison, the simulation results are coarsened to 10 × 10 data points by averaging the current density over each segment. It shows that experimental and simulation results indicate similar distributions. Specifically, the current density at the top area is much higher than that at the bottom, and the lowest current density is present near the air outlet, as shown in Fig. 4(a), in which the row number in Fig. 4 is consistent with that in Fig. 3. It can be seen in Fig. 4(b) that all the local absolute errors are within 0.15 A cm−2 , and the mean absolute error (MAE) values are 0.017, 0.052 and 0.078 A cm−2 for 0.1, 1.0 and 1.2 A cm−2 , respectively. The mean relative absolute error values are 16.6%, 5.4% and 6.3% for 0.1, 1.0 and 1.2 A cm−2 , respectively. Meanwhile, both experimental data and simulation results indicate that the local current density difference increases with the increment of current density. This is mainly caused by the oxygen consumption due to the ORR and hence decreasing oxygen concentration down the air flow path, as shown in Fig. 5(a). In this study, the local concentration loss is clearly predicted as shown in Fig. 3 despite that there is no obvious concentration loss in the overall performance (Fig. 2). However, there is little difference between simulation and experimental results. The simulation results show that the highest current density occurs near the cathode inlet, which is contrary to the experimental result. This difference is probably due to the

9

difference in the inlet arrangement between the simulation and experiment. In experiment, see Fig. 1 in Ref. [29], the inlet flow is perpendicular to the BP when flowing to the fuel cell, while the inlet flow is modified to be parallel to the BP in the simulation. In addition, there are assumptions in the present model, such as no two-phase flow in the channel. Again, the main trends of the current density distribution change are similar for the simulation and experiment. In addition, the increased liquid saturation along the air flow path (Fig. 5(b)), as a result of water production, also contributes to this observed trend of the current density distribution because liquid water hampers oxygen diffusion toward the CL. Fig. 5(c) displays the membrane water content distribution, which is similar as that of liquid saturation. As for the temperature distribution in Fig. 5(d), it shows that the temperature under the manifold area is high. This is because only a small portion of heat produced in fuel cell is removed by the gas flow [44]. The majority of the heat produced under the manifold has to be transported to the neighbor ribs (or lands) via the GDL for removal, which has a long transport path and consequent high resistance. To show the effects of relative humidity and operation temperature, Figs. 6 and 7, respectively, compare the spatial distributions of current density between the model prediction and experiment. Meanwhile, Fig. 8 shows the comparison of the local current

Fig. 6. Comparison of the predicted current density distributions with experimental data [29] (353.15 K, 1.0 A cm−2 ) under (a): 50% RH and (b): 25% RH for PEM fuel cell (I).

10

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

Fig. 7. Comparison of the predicted current density distributions at two temperatures with experimental data [29] (50% RH) for PEM fuel cell (I).

Fig. 8. (a): Comparison of the local current density from inlet to outlet (top to bottom) with experimental data [29] and (b): absolute error under 50% and 25% RHs for PEM fuel cell (I) (353.15 K, 1.0 A cm−2 ).

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

density from inlet to outlet (or from top to bottom) and the corresponding absolute errors. Specifically, the MAE values are 0.052 and 0.042 A cm−2 for 50% and 25% RHs, respectively, and the maximum local absolute error is lower than 0.15 A cm−2 . The mean relative absolute error values are 5.4% and 4.1% for 50% and 25% RHs, respectively. At high humidification conditions, the oxygen is a main factor affecting the current density distribution. For a lower relative humidity (25% RH), as shown in Figs. 6 and 8, the high current density region shifts downstream to the middle region of fuel cell. This is likely controlled by the Ohmic polarization, which is large near the inlet due to the local dry membrane as a result of the 25% RH inlet flow. Fig. 9(a) displays the water content distributions at the middle plane of the membrane for 25% and 50% RH, respectively. It is shown that for the former a much larger area of inlet region or cathode upstream region is subjected to dry member even though water is produced in the cathode CL. The majority

11

of the top region shows a membrane water content as low as 6, thus a large membrane resistance and Ohmic loss are present locally, leading the observed low current density for 25% RH as shown in Fig. 6(b). In this case study, the two operating temperatures showed similar trend of the current density distribution variation, as shown in Fig. 7. At the same RH, the water is more likely to exist in the vapor state under higher temperature (Table 6), which is favorable for membrane hydration due to high water activity (Fig. 9(b)). At the same time, higher water vapor concentration at higher temperature dilutes oxygen [39], which will increase local oxygen transport polarization. However, decrease of the liquid water level under high temperature will benefit oxygen transport. Thus, the temperature’s impacts are complex and need a 3D study to precisely capture. In the next, the predicted current density distribution was compared with the experimental data reported by the University of

Fig. 9. Water content distributions in the middle plane of membrane at different relative humidities (a) and temperatures (b) for PEM fuel cell (I).

12

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294 Table 6 Comparison of some variables at different temperatures for PEM fuel cell (I). Parameters

333.15 K

353.15 K

Ohmic loss caused by membrane and CLs Concentration loss in cathode CL Oxygen concentration in cathode CL Water vapor concentration in cathode CL Liquid saturation in cathode CL

0.0592 V 0.4688 V 3.347 mol m−3 19.359 mol m−3 26.84%

0.0446 V 0.4647 V 2.870 mol m−3 26.289 mol m−3 20.14%

Fig. 10. Comparison of the predicted current density distribution with experimental data [45] at two output voltages (338.15 K, 100% RH) for PEM fuel cell (II).

Waterloo in Canada using their fuel cell configuration [45], which contains a three-path serpentine flow field and a MEA area of 40 cm2 . The schematic of the fuel cell is shown in Fig. 1 (PEM fuel cell II). In addition, the anode and cathode flow fields are the same and arranged in co-flow, i.e. the anode and cathode inlets or outlets are put on the same sides of the fuel cell. Fig. 10 shows comparison of the predicted current density distribution with experiment PEM fuel cell II. The operation conditions and geometry parameters are shown in Table 5. It can

be seen that the predicted current density distributions agree reasonably with the experimental data for the two operating voltages. In both cases, the peak electric current is present near the cathode inlet primarily due to the local abundant oxygen availability and hence small transport polarization. Low current density occurs near the cathode outlet, caused by oxygen consumption by the oxygen reduction reaction (ORR). This shows that the oxygen concentration is a major factor affecting the current density spatial distribution.

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

4. Conclusions In this study, the current density distributions predicted by a three-dimensional (3D) multi-phase model of PEM fuel cell were validated against the experimental data from the Los Alamos National Laboratory in the United States (PEM fuel cell I) and the University of Waterloo in Canada (PEM fuel cell II). The complex structure of catalyst layer (CL) was described using the spherical agglomerate sub-model, which was incorporated in the 3D multi-phase fuel cell model. For PEM fuel cell I, we first validated the polarization curves at different relative humidities, and then compared the current density distributions under different current densities, temperature and relative humidities (RH). The maximum local absolute error between simulation and experimental results is lower than 0.15 A cm−2 , and the mean absolute error (MAE) values are 0.017, 0.052, 0.078 and 0.042 A cm−2 for 0.1, 1.0, 1.2 A cm−2 (353.15 K, 50% RH) and 1.0 A cm−2 (353.15 K, 25% RH), respectively. In addition, the mean relative absolute error values are 16.6%, 5.4%, 6.3% and 4.1 for 0.1, 1.0 and 1.2 A cm−2 (353.15 K, 50% RH) and 1.0 A cm−2 (353.15 K, 25% RH), respectively. It was found that the high current region moves downstream from the inlet area to the middle of fuel cell when reducing RH from 50% to 25%. Simulation indicated that the inlet flow under 25% RH operation causes a much larger area of the inlet region to dryness with a membrane water content as low as 6, which is attributed to the shift in the current density. In our case study, the two operating temperatures showed similar trend of the current density distribution variation. For PEM fuel cell II, the comparison shows reasonably good agreement between model prediction and experiment. High current density was present near the cathode inlet and low current density appeared near the outlet, as a result of oxygen consumption by the ORR. This validation study indicated that the 3D multi-phase model is reliable to predict the current density distributions for PEM fuel cells. Declaration of Competing Interests The authors declare that they have no conflict of interest. CRediT authorship contribution statement Guobin Zhang: Methodology, Software, Validation, Investigation, Writing - original draft, Visualization. Jingtian Wu: Software, Investigation, Data curation. Yun Wang: Conceptualization, Resources, Writing - review & editing, Funding acquisition. Yan Yin: Writing - review & editing, Funding acquisition. Kui Jiao: Conceptualization, Resources, Writing - review & editing, Funding acquisition. Acknowledgements This work is supported by the National Key Research and Development Program of China (Grant No. 2018YFB0105601), and the China-UK International Cooperation and Exchange Project (Newton Advanced Fellowship) jointly supported by the National Natural Science Foundation of China (Grant No. 51861130359) and the UK Royal Society (Grant No. NAF\R1\180146). Y. Wang also thanks Shanghai Everpower Technologies Ltd for partial financial support. References [1] Y. Wang, K.S. Chen, S.C. Cho, PEM Fuel cells: Thermal and Water Management Fundamentals, Momentum Press, 2013. [2] Y. Wang, D.F.R. Diaz, K.S. Chen, Z. Wang, X.C Adroher, Materials, technological status, and fundamentals of PEM fuel cells-a review, Mater. Today (2019) [in press], doi:10.1016/j.mattod.2019.06.005. [3] J. Larminie, A. Dicks, S McDonald M, Fuel Cell Systems explained. Chichester, UK: J. Wiley, 2003.

13

[4] K. Jiao, X. Li, Water transport in polymer electrolyte membrane fuel cells, Prog. Energy Combust. Sci. 37 (3) (2011) 221–291. [5] Y. Cai, Z. Fang, B. Chen, et al., Numerical study on a novel 3D cathode flow field and evaluation criteria for the PEM fuel cell design, Energy 161 (2018) 28–37. [6] F. Cao T, T. Mu Y, J. Ding, et al., Modeling the temperature distribution and performance of a PEM fuel cell with thermal contact resistance, Int. J. Heat Mass Transf. 87 (2015) 544–556. [7] J. Shen, Z. Tu, H. Chan S, Enhancement of mass transfer in a proton exchange membrane fuel cell with blockage in the flow channel, Appl. Therm. Eng. 149 (2019) 1408–1418. [8] Q. Meyer, K. Ronaszegi, B. Robinson J, et al., Combined current and temperature mapping in an air-cooled, open-cathode polymer electrolyte fuel cell under steady-state and dynamic conditions, J. Power Sources 297 (2015) 315–322. [9] G. Zhang, K. Jiao, Multi-phase models for water and thermal management of proton exchange membrane fuel cell: a review, J. Power Sources 391 (2018) 120–133. [10] R. Wang, G. Zhang, Z. Hou, et al., Comfort index evaluating the water and thermal characteristics of proton exchange membrane fuel cell, Energy Convers. Manage. 185 (2019) 496–507. [11] Y. Jiang, Z. Yang, K. Jiao, et al., Sensitivity analysis of uncertain parameters based on an improved proton exchange membrane fuel cell analytical model, Energy Convers. Manage. 164 (2018) 639–654. [12] B. Wang, H. Deng, K. Jiao, Purge strategy optimization of proton exchange membrane fuel cell with anode recirculation, Appl. Energy 225 (2018) 1–13. [13] B. Wang, K. Wu, Z. Yang, et al., A quasi-2D transient model of proton exchange membrane fuel cell with anode recirculation, Energy Convers. Manage. 171 (2018) 1463–1475. [14] G. Zhang, L. Fan, J. Sun, et al., A 3D model of pemfc considering detailed multiphase flow and anisotropic transport properties, Int. J. Heat Mass Transf. 115 (2017) 714–724. [15] G. Zhang, B. Xie, Z. Bao, et al., Multi-phase simulation of proton exchange membrane fuel cell with 3D fine mesh flow field, Int. J. Energy Res. 42 (15) (2018) 4697–4709. [16] G. Zhang, X. Xie, B. Xie, et al., Large-scale multi-phase simulation of proton exchange membrane fuel cell, Int. J. Heat Mass Transf. 130 (2019) 555–563. [17] J. Zhou, D. Stanier, A. Putz, et al., A mixed wettability pore size distribution based mathematical model for analyzing two-phase flow in porous electrodes ii. Model validation and analysis of micro-structural parameters, J. Electrochem. Soc. 164 (6) (2017) F540–F556. [18] Y. Wang, Modeling of two-phase transport in the diffusion media of polymer electrolyte fuel cells, J. Power Sources 185 (1) (2008) 261–271. [19] L. Fan, G. Zhang, K Jiao, Characteristics of PEMFC operating at high current density with low external humidification, Energy Convers. Manage. 150 (2017) 763–774. [20] G. Zhang, K. Jiao, Three-dimensional multi-phase simulation of PEMFC at high current density utilizing Eulerian-Eulerian model and two-fluid model, Energy Convers. Manage. 176 (2018) 409–421. [21] P. Ž, T. I, F. Barbir, Computational fluid dynamics study of PEM fuel cell performance for isothermal and non-uniform temperature boundary conditions, Int. J. Hydrogen Energy 41 (39) (2016) 17585–17594. [22] L. Hao, K. Moriyama, W. Gu, et al., Modeling and experimental validation of PT loading and electrode composition effects in PEM fuel cells, J. Electrochem. Soc. 162 (8) (2015) F854–F867. [23] Q. Tao W, H. Min C, L. Liu X, et al., Parameter sensitivity examination and discussion of PEM fuel cell simulation model validation: part I. current status of modeling research and model development, J. Power Sources 160 (1) (2006) 359–373. [24] C. Pérez L, L. Brandão, M. Sousa J, et al., Segmented polymer electrolyte membrane fuel cells—A review, Renew. Sustain. Energy Rev. 15 (1) (2011) 169–185. [25] J. Shan, R. Lin, X. Chen, et al., EIS and local resolved current density distribution analysis on effects of MPL on PEMFC performance at varied humidification, Int. J. Heat Mass Transf. 127 (2018) 1076–1083. [26] S. Yoshioka, A. Yoshimura, H. Fukumoto, et al., Development of a PEFC under low humidified conditions, J. Power Sources 144 (1) (2005) 146–151. [27] D. Liang, Q. Shen, M. Hou, et al., Study of the cell reversal process of large area proton exchange membrane fuel cells under fuel starvation, J. Power Sources 194 (2) (2009) 847–853. [28] H. Hauer K, R. Potthast, T. Wüster, et al., Magnetotomography—a new method for analysing fuel cell performance and quality, J. Power Sources 143 (1–2) (2005) 67–74. [29] B. Carnes, D. Spernjak, G. Luo, et al., Validation of a two-phase multidimensional polymer electrolyte membrane fuel cell computational model using current distribution measurements, J. Power Sources 236 (2013) 126–137. [30] H. Ju, Y. Wang C, S. Cleghorn, et al., Nonisothermal modeling of polymer electrolyte fuel cells I. experimental validation, J. Electrochem. Soc. 152 (8) (2005) A1645–A1653. [31] J. Li, Y. Wang C, A. Su, Prediction and experimental validation of in-plane current distribution between channel and land in a PEFC, J. Electrochem. Soc. 155 (1) (2008) B64–B69. [32] L. Hao, K. Moriyama, W. Gu, et al., Three dimensional computations and experimental comparisons for a large-scale proton exchange membrane fuel cell, J. Electrochem. Soc. 163 (7) (2016) F744–F751. [33] A. Iranzo, P. Boillat, F. Rosa, Validation of a three dimensional PEM fuel cell CFD model using local liquid water distributions measured with neutron imaging, Int. J. Hydrogen Energy 39 (13) (2014) 7089–7099.

14

G. Zhang, J. Wu and Y. Wang et al. / International Journal of Heat and Mass Transfer 150 (2020) 119294

[34] Y. Wang, S. Chen K, Through-plane water distribution in a polymer electrolyte fuel cell: comparison of numerical prediction with neutron radiography data, J. Electrochem. Soc. 157 (12) (2010) B1878–B1886. [35] Y. Wang, S. Chen K, Effect of spatially-varying gdl properties and land compression on water distribution in PEM fuel cells, J. Electrochem. Soc. 158 (11) (2011) B1292–B1299. [36] J. Zhao, Catalyst Layers in Polymer Electrolyte Membrane Fuel Cells: Formation, Characterization and Performance, University of Waterloo, 2019. [37] L. Xing, X. Liu, T. Alaje, et al., A two-phase flow and non-isothermal agglomerate model for a proton exchange membrane (PEM) fuel cell, Energy 73 (2014) 618–634. [38] G. Zhang, K. Jiao, R. Wang, Three-Dimensional simulation of water management for high-performance proton exchange membrane fuel cell, SAE Int. J. Alternative Powertrains 7 (2018–1–1309) (2018). [39] G. Zhang, H. Yuan, Y. Wang, et al., Three-dimensional simulation of a new cooling strategy for proton exchange membrane fuel cell stack using a non-isothermal multiphase model, Appl. Energy 255 (2019) 113865.

[40] L. Fan, Z. Niu, G. Zhang, et al., Optimization design of the cathode flow channel for proton exchange membrane fuel cells, Energy Convers. Manage. 171 (2018) 1813–1821. [41] K. Broka, P. Ekdunge, Modelling the PEM fuel cell cathode, J. Appl. Electrochem. 27 (3) (1997) 281–289. [42] M. Moein-Jahromi, J. Kermani M, Performance prediction of PEM fuel cell cathode catalyst layer using agglomerate model, Int. J. Hydrogen Energy 37 (23) (2012) 17954–17966. [43] M. Moore, P. Wardlaw, P. Dobson, et al., Understanding the effect of kinetic and mass transport processes in cathode agglomerates, J. Electrochem. Soc. 161 (8) (2014) E3125–E3137. [44] G. Zhang, G. Kandlikar S, A critical review of cooling techniques in proton exchange membrane fuel cell stacks, Int. J. Hydrogen Energy 37 (3) (2012) 2412–2429. [45] I. Alaefour, G. Karimi, K. Jiao, et al., Experimental study on the effect of reactant flow arrangements on the current distribution in proton exchange membrane fuel cells, Electrochim. Acta 56 (5) (2011) 2591–2598.