Lattice Boltzmann simulation of oxygen diffusion and electrochemical reaction inside catalyst layer of PEM fuel cells

Lattice Boltzmann simulation of oxygen diffusion and electrochemical reaction inside catalyst layer of PEM fuel cells

International Journal of Heat and Mass Transfer 143 (2019) 118538 Contents lists available at ScienceDirect International Journal of Heat and Mass T...

7MB Sizes 0 Downloads 33 Views

International Journal of Heat and Mass Transfer 143 (2019) 118538

Contents lists available at ScienceDirect

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

Lattice Boltzmann simulation of oxygen diffusion and electrochemical reaction inside catalyst layer of PEM fuel cells Hao Deng a, Yuze Hou a, Wenmiao Chen b, Fengwen Pan b, Kui Jiao a,⇑ a b

State Key Laboratory of Engines, Tianjin University, 135 Yaguan Rd, Tianjin 300350, China Weichai Power Co. Ltd., 197A Fushou St. E., Weifang 261016, China

a r t i c l e

i n f o

Article history: Received 5 June 2019 Received in revised form 5 August 2019 Accepted 6 August 2019 Available online 17 August 2019 Keywords: Catalyst layer Oxygen transport Electrochemical reaction Lattice Boltzmann Stochastic algorithm

a b s t r a c t An in-depth understanding of pore-scale transport process and electrode reaction kinetics is critical to ease the sluggish oxygen reduction reaction of proton exchange membrane fuel cells. In this study, the effects of catalyst layer microstructure on the oxygen diffusion and electrochemical reaction consumption are analyzed, based on the developed 3D single-phase lattice Boltzmann model and optimized stochastic reconstruction algorithm. The results suggest that increasing the ionomer content and platinum loading are beneficial to increasing the active sites, while the oxygen transport along the thickness direction becomes more difficult due to the narrowed pores, so the platinum particles in the back region have less participation in the reaction, which causes the significant catalyst waste and limits the further cost reduction and improvement of power density. Improving the carbon particle diameter ensures the enlarged pores and more effective Pt particles covered by ionomer, thus enhances the oxygen supply and electrochemical reaction rate. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Sluggish electrode reaction kinetics in the cathode is a critical factor to limit the further performance improvement of proton exchange membrane (PEM) fuel cells, as the reaction rate of oxygen reduction is almost 1/6–1/4 compared to that of hydrogen oxidation. To ease the oxygen reduction reaction (ORR) in PEM fuel cells, high loading use of platinum (Pt) based catalyst is essential owing to its high chemical stability and catalytic activity, while cost and durability issues also rise because Pt is expensive and easily poisoned, which accounts for almost half of all PEM fuel cell manufacturing costs [1]. Considerable efforts were devoted to the cost reduction of Ptbased catalyst and its durability improvement, by developing new electrode materials and structures [2–4]. On the one hand, non-noble metal catalysts have been studied to partially or fully replace Pt, such as the transition metal catalysts, carbon materials with unique nanostructures and nitrogen-doped carbon structures, even though their current activity and durability are still not enough especially in the cathode. On the other hand, new nanostructured catalysts and supports are also exploited to extend the active surface area and enhance the chemical stability, such as

⇑ Corresponding author. E-mail address: [email protected] (K. Jiao). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118538 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

the Pt alloying with non-precious metals, carbon nanotubes and wires. Because catalyst layer (CL) should simultaneously provide effective path ways for gas reactants, electrons and protons, improving the understanding of transport processes and electrode reaction kinetics are expected to guide the designs of electrode structure. On account of the opaque nature of traditional CL materials, gas reactants and liquid water transports in an operating fuel cell are the challenging phenomena to study using in-situ visualization techniques, thus appropriate numerical methods may be the promising tool. The traditional macroscopic numerical models treat CL as the isotropic porous media, by characterizing the structure and transport properties with the statistical or empirical parameters, such as porosity, ionomer content, effective diffusion coefficient, permeability and thermal conductivity, and the relationship between reaction rate and activation overpotential is established to use the assumed reference current density and reaction transfer coefficient, then the macroscopic physical quantity distributions are obtained by solving the conservation equations [5–7]. However, the hypothesis of isotropic material ignores the various microstructure characteristics of porous media, so it is unable to explore the coupling effects of microstructure, transport process and macroscopic performance. On this basis, a series of improved macroscopic CL models have been developed, such as the thin-film, discrete-volume and agglomerate models. Harvey

2

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

et al. [8] compared three models using the finite volume method, and the results show that only the agglomeration model, which is based on the more realistic CL structure, is capable to capture the ‘‘oxygen starvation” phenomenon at high current density. Mesoscopic lattice Boltzmann (LB) method is another promising approach to explore the transport processes and electrode reaction kinetics at pore-scale, due to its natural background of microscopic particle, ease dealing of complicated porous boundaries and extensible parallel computation [9–13]. At present, the LB studies inside the porous layers of fuel cells can be divided into three categories: one is to predict the structure and physical properties of porous media, such as the permeability, effective mass diffusivity, thermal conductivity and mass transfer coefficient [14–19]; the second is to study the phenomena related to gas species and liquid water transports [20–25]; and the third is the research on electrode reaction kinetics [26–35]. Debates about whether the LB method is suitable for CL simulation are existed before 2010, which are focused on how to deal with the electrochemical reaction. In recent years, LB method is recognized as a powerful tool evidenced by the increased research studies in this aspect. The coupling of electrochemical reaction is realized by adding the additional source term at the right side of LB evolution equation, and the source term is proportional to the reaction rate. The relationship between reaction rate and local oxygen concentration is determined by ButlerVolmer equation, which is usually reflected as the change of the concentration distribution function at the active lattice points [32,36]. Although LB method has obvious advantages in dealing with complex structural boundaries, many previous LB studies in CL are based on the simplified electrode structure. In 2010, Ostadi et al. [26] developed a LB model to study the multi-component gas transport in the cathode, and CL structure is brought into the LB study for the first time. However, their study simplified CL into a thin film to define the flow boundary, and assumed that the current density distribution in the whole CL was uniform. Chen et al. [27] also simplified CL into the interface at the bottom side of gas diffusion layer (GDL) in their earliest model, and in the subsequent work [28], CL is replaced by a regular structure with spacing distribution for pore and solid structure. In the recent work [29], they reconstructed 3D CL nanostructures using the Quartet Structure Generation Set (QSGS) method, which can freely control the content and distribution of carbon support, ionomer and catalyst, thus the porosity, pore distribution and specific surface area are adjustable. However, the CL reconstruction in their work assumed carbon support is formed by the accumulation of carbon lattices, rather than the spherical structure, so the resulted pore size and distribution still lags behind the realistic microstructure. The LB simulation results [29] show that the effective diffusion coefficient of gas species in CL significantly depend on the microstructure, and the simulated tortuosity is larger than the value obtained by Bruggeman approximation. The effective proton conductivity obtained from LB study [29] is close to that calculated by Lange et al. [37], but is much lower than the experimental result [38], which implies that except ionomer, there may be other available paths for proton transport in CL, such as in the liquid water [39]. Molaeimanesh et al. [30–32] also devoted considerable efforts to the pore-scale modeling in CL, and their most recent work [30] reconstructed CL with randomly distributed ellipsoids, then a 3D agglomeration model is established based on the LB and finite difference methods, to study the concentration distribution of gas species in the pore, ion potential distribution in the ionomer and current density distribution at the interface between membrane and CL under different CL micro-morphologies. However, the existence of Pt particle and detailed agglomeration structure are neglected in their study. Chen et al. [33] recently reconstructed the single carbon particle covered by ionomer and the dispersed

Pt is at the smaller lattice resolution of 0.5 nm, then the effects of microstructures and reactive transport conditions on oxygen transport are investigated at such a tiny computational domain. On the whole, CL reconstruction remains to be improved, to more accurately incorporate the effect of dissimilar electrode microstructures on the catalyst activity and effective transport of gas reactant and proton. In this study, the optimized stochastic model is developed to reconstruct carbon support microstructure with the spherical particles and the lattice spatial resolution is improved to 2 nm, then more realistic 3D CL micromorphology is obtained by considering carbon support, ionomer and Pt particle simultaneously. A 3D single-phase LB model is developed to study the pore-scale oxygen transport process and electrochemical reaction, and the effects of ionomer amount, Pt loading and carbon particle size are investigated. 2. Reconstruction of CL microstructure Compared with GDL and micro porous layer (MPL), the reconstruction of CL with multi-scale structure is more complicated, as it consists of three different phases: porous phase (or gas phase), mixed electronic phase (carbon and Pt particles) and ionomer phase. CL thickness and diameters of single Pt and carbon particle are about 5–10 lm, 2–5 nm and 10–100 nm, respectively, and the diameters of shaped pores usually below than 100 nm. CL is often partially simplified in the literature to lower the computation of reconstruction and flow simulation, specifically: a. Assuming CL is composed of porous and solid phases (carbon particle, Pt powder and ionomer are considered as a whole), which simplifies the code and lowers the computational cost with the sacrifice of details in microstructure characteristics. b. More elaborate treatment is to take account into the ionomer phase, while carbon and Pt phases are regarded as one without distinction. Such a treatment is able to consider the gas species and proton transports in the ionomer and effect of ionomer distribution on the pore structure, with neglecting the surface morphologies of catalyst. In this study, based on the previous studies reported by Siddique and Liu [40] and Chen et al. [29], CL microstructures are reconstructed with considering Pt particle, carbon support, ionomer and shaped pores simultaneously, and the reconstruction process is described below. 2.1. Calculation of control parameters The first step of reconstruction is to calculate volume fractions of carbon, Pt, ionomer (eC , ePt , eion ) and Pt loading cPt (mg cm2) using the predefined statistics, which are used as the control parameters of following steps. The predefined statistics include porosity e, thickness d (cm), Pt/C mass ratio hPt=C and ionomer content hion (mass ratio of ionomer to electrode), and their relationships are given as follows:

eC ¼ ð1  eÞ

1=qC   hPt=C =qPt þ 1=qC þ 1 þ hPt=C hion =ðqion ð1  hion ÞÞ

ePt ¼ ð1  eÞ

ð1Þ

h =q  Pt=C Pt  hPt=C =qPt þ 1=qC þ 1 þ hPt=C hion =ðqion ð1  hion ÞÞ ð2Þ 

eion ¼ ð1  eÞ

 1 þ hPt=C hion =ðqion ð1  hion ÞÞ   hPt=C =qPt þ 1=qC þ 1 þ hPt=C hion =ðqion ð1  hion ÞÞ ð3Þ

cPt ¼ ePt qPt d

ð4Þ

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

where qPt , qC and qion (mg cm3) are densities of Pt, carbon and ionomer, respectively.

3

achieved with QSGS method, and assuming the lattice surrounded by multiple growth nucleuses automatically converts into growth phase is helpful to achieve the effective connection of ionomer.

2.2. Generation of carbon phase 3. Model formulation The growth process of carbon phase is specially treated to improve the connectivity, which is also reported by Chen et al. [29]. In the first treatment, the growth of carbon lattice is limited in six adjacent directions (up, down, left, right, forward and back) without considering the diagonal direction, and in the second, the lattice around by two or more carbon points will automatically convert into carbon phase. Chen et al. [29] assumed that the carbon phase grows from the accumulation of carbon points with diameter of 5 nm, which neglects the spherical morphology of carbon particles. In this study, the carbon support structure is reconstructed starting from the spherical particles with the improved spatial resolution of 2 nm, and the reconstruction process is depicted as follows: (1) Generating a new carbon sphere: generate a carbon sphere at the random location with the diameter randomly varied in a specified range. (2) Growth nucleus judgement: generate a random number in [0, 1] and assume the initial distribution probability of growth nucleus is P, if the random number is less than P, the carbon sphere generated in (1) is regarded as a new growth nucleus. Then the position of carbon sphere needs to be calibrated, if there is overlap between this new carbon sphere and any previous one, the carbon sphere will be regenerated back to step (1) until the carbon sphere does not intersect with all other carbon spheres. This treatment is expected to obtain more even distribution of carbon support structure rather than clustered together. (3) Overlap judgment: if the carbon sphere generated in (1) is not a new growth nucleus, the overlap judgement needs to be conducted. If the new one overlaps with other carbon spheres, go on to the next step, if not, regenerate the carbon sphere randomly and make the judgment in step (2). (4) Overlap ratio judgment: define the maximum allowable overlap ratio between carbon spheres as r, and calculate the volume fraction of new carbon sphere occupied by other carbon spheres, if the volume fraction is not greater than r, proceed to the next step, if yes, regenerate the carbon sphere randomly and go back to steps (2) and (3) for judgment. (5) Generation confirmation: when step (2) or steps (3) and (4) simultaneously meet the requirements, the generated carbon sphere is regarded as the effective one, and all lattices belong to the region of new carbon sphere are converted into carbon phase lattices. (6) Volume fraction check: check if the volume fraction of carbon phase reaches the predefined target, if not, repeat steps (1)–(5). 2.3. Generation of Pt and ionomer phases The growth process of Pt particles is almost consistent with the steps described by QSGS method [29], except that the growth position is limited in the lattice points adjacent to carbon phase, which improves the effective contact between Pt particle and carbon support. In the simplest treatment for ionomer phase, ionomer thickness is assumed to be equal at the surface of carbon support, which has been widely adopted in most of agglomerate models [41] and a few of pore-scale models [37,42]. However, the uniform assumption of ionomer thickness may overestimate the diffusion coefficient of oxygen and underestimate the proton conductivity [37]. The non-uniform distribution of ionomer at Pt/C surface can be

3.1. Physical problems For the current LB model, it is difficult to consider all the transport processes and electrochemical reaction simultaneously, due to the limitations of numerical method and computational ability. The results presented by Chen et al. [29] show that the variations of electric and ionic potentials are insignificant in such a small computational domain, so their distributions are approximately uniform, and the overpotential can be fixed in simulations. Based on the reconstructed microstructures, a single-phase LB model is developed in coupling the electrochemical reaction. 3.2. Macroscopic physical-electrochemical model Without considering the transport of water, electron and proton, and based on the assumption of isothermal operating condition, only the transport and electrochemical reaction related to oxygen are studied. The gas convection-diffusion equation is generally described as:

@C þ r  ðCuÞ ¼ r  ðDrC Þ þ S @t

ð5Þ

where the four terms in above equation are transient, convection, diffusion and sour terms, respectively. The convection term can be neglected in nano-scale CL, due to the insignificant flow velocity and dominate role of diffusion transport for gas species, thus the above governing equation can be rewritten as:

  @C O2 ¼ r  DO2 rC O2 þ SO2 @t

ð6Þ

The transport paths of oxygen in the cathode CL include: diffusion from the pore to ionomer surface, dissolution and transport in the ionomer, which can be described as follows:

  @C O2 ;p ¼ r  DO2 ;p rC O2 ;p þ SO2 ;p @t

ð7Þ

RT C O2 ;e ¼ H C O2 ;p

ð8Þ

  @C O2 ;e ¼ r  DO2 ;e rC O2 ;e þ SO2 ;e @t

ð9Þ

Similar to the jump phenomenon of liquid saturation at the interface of different porous layers [7,43], the oxygen concentration also jumps at the two sides of ionomer surface, due to the significant difference of oxygen diffusivity in the pore and ionomer, which can be described by Henry’s Law. For CL, Henry’s law stipulates that the amount of oxygen dissolved in the ionomer is proportional to the partial pressure of oxygen in the adjacent pores, that is, the relationship of oxygen concentrations at the two sides of interface is linear (Eq. (8)), where R (J mol1 K1) is the ideal gas constant, T (K) is the fixed temperature, and Henry constant H is given as follows:

  666 H ¼ 1:01325  101 exp þ 14:1 T

ð10Þ

Owing to the Henry’s law, the conservation equations (7) and (9) can be rewritten as the united formula as in Eq. (6), which means only a set of lattices in LB model is required to describe the evolution of oxygen transport in different phase regions, thus

4

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

reduces the computation. In this study, the realistic oxygen concentration in the pores is stipulated as the direct solving variable, and the oxygen concentration in the ionomers can be obtained by inverse calculation according to Henry’s law. Specifically, the diffusion coefficients and source terms in the pore and ionomer should be defined respectively in different regions. The macro numerical method usually treats gas as the continuous medium, while this assumption is not available to CL due to its nanoscale microstructure. In the flow problems, it is generally judged whether the fluid is suitable for the continuous medium hypothesis based on the Knudsen number (the ratio of the mean molecular free path k of gas to the characteristic length L). When Knudsen number is less than 0.01, the gas can be regarded as the continuous medium, and if not, considering the Knudsen diffusion is more accurate to describe the gas transport process. For oxygen diffusion in the CL pores, the characteristic length L is adopted as the pore diameter, which is typically less than 100 nm. The mean molecular free path k of gas satisfies the following formula:

kT k ¼ pffiffiffi 2 2p d p

ð11Þ

Table 1 Oxygen diffusivity and source term in different regions [44–47]. Parameter Oxygen diffusivity in the pore

Equation DO2 ;p ¼



Unit 1 DO2 ;b

þ DO1;Kn

1

m2 s1

[44]

2

1:5

ðT=293:2Þ [45] DO2 ;b ¼ 0:22  104 p=1:01310 5 ð Þ q ffiffiffiffiffiffiffiffiffi ffi d DO2 ;Kn ¼ 3p p8RT [46] MO 2

Oxygen diffusivity in the ionomer Source term

m2 s1

DO2 ;e ¼ 1010  ð0:1543ðT  273Þ  1:65Þ [47] SO2 ¼

H n i  RT DL 4F ; 0;

mol m3 s1

in an active site cell in other cells

where k ¼ 1:38  1023 J K1 is the Boltzmann constant, T (K) and p (Pa) are the ambient temperature and pressure, respectively, and d (m) is the effective molecule diameter. Take T = 333.15 K, p = 151950 Pa, d ¼ 3:46  1010 m, then the mean molecular free path of oxygen is calculated as k ¼ 5:691  108 m, which is the same order to CL pore diameter and Knudsen number is close to 1. Therefore, the bulk and Knudsen diffusions are simultaneously considered in the calculation of oxygen diffusivity (Table 1, [44– 47]), and the pore diameter dp (m), which is used to calculate the Knudsen diffusion coefficient DO2 ;Kn (m2 s1), varies at different positions. In addition to the diffusion transport, another major process related to oxygen is the electrochemical reaction consumption, that is, oxygen dissolved in the ionomer transports to the surface of Pt particle, and reacts with the proton to produce water. As given in Table 1, the source term of oxygen in the non-activated region is 0, and the consumption rate in the activated region is related to the current density and number of activated lattices. The activated region is defined as the three-phase interface, and should satisfy that the lattice is in the ionomer region, and at least one adjacent lattice is occupied by Pt. The n in the source term represents the number of adjacent lattices occupied by Pt, as n increases, the catalysis is improved and oxygen consumption is faster. DL (m) is the physical length represented by a lattice cell, and the current density i (A m2) is given by the simplified Tafel equation:

i ¼ iref

  C O2 ;e aF exp Du RT C O2 ;ref

ð12Þ

where iref (A m2) is the reference current density, C O2 ;ref (mol m3) is the reference oxygen concentration, a is the reaction transfer coefficient, F (C mol1) is the Faraday constant, R (J mol1 K1) is the ideal gas constant, T (K) is the operating temperature, and Du Table 3 Base and calculated conditions.

Table 2 Model and operating parameters. Parameter

Value

Unit

Temperature, T Pressure, P

333.15

K Pa

Oxygen inlet concentration, C O2 ;in Activation overpotential, Du Reference current density, iref Reference oxygen concentration, C O2 ;ref Reaction transfer coefficient, a

1:5  1:0135  10 10.28 0.6 0.01 40.96 0.592

5

mol m3 V A m2 mol m3 /

Parameter

Base condition

Calculated conditions

Mass ratio of ionomer to electrode Pt loading

0.3

0.2, 0.3, 0.4

0.10 mg cm2

Carbon particle diameter

20–32 nm (100%)

Mass ratio of Pt to carbon

0.8

0.06, 0.08, 0.10, 0.12 mg cm2 20–32 nm (100%) 32–40 nm (100%) 12–20 nm (50%) + 32– 40 nm (50%) 12–20 nm (100%) 0.8

Fig. 1. Computational domain and boundary conditions.

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

where x is the space coordinate, t is the time, s is the dimensionless discrete time, a is the discrete velocity direction, ea is the discrete velocity in the a direction, xa is the weight coefficient in the a direction, and SO2 represents the oxygen source term. The lattice structure of D3Q19 is adopted in this study, and the weight coefficients of 19 directions are x0 ¼ 1=3, x16 ¼ 1=18 and x718 ¼ 1=36, and the discrete velocity ea and equilibrium distribution function g eq a ðx; tÞ satisfy:

(V) is the activation overpotential, and more details are given in Table 2. 3.3. Lattice Boltzmann model Corresponding to the oxygen governing equation (Eq. (6)) in Section 3.2, the LB evolution equation describes the spatialtemporal variation of distribution function for oxygen concentration, which is given as follows:

1 g a ðx þ ea Dt;t þ DtÞ  g a ðx;tÞ ¼  ðg a ðx;tÞ  g eq a ðx; tÞÞ þ xa SO2

s

2

0 1 0

0 1

0

0

1 1 0 1 1

0

1 0 1 1

6 e ¼ 40 0 1 0

0

1

0 0

0

0

0 1

1 0 1 1

0

0

0

ð13Þ

1

1 1

1 1

5

0

1

0

0

1

1 1

1 1 1

0

0

1

0

3

7 1 5 1

Fig. 2. Reconstructed CL microstructures with different mass ratios between ionomer and electrode: (a) hion = 0.2; (b) hion = 0.3; (c) hion = 0.4.

ð14Þ

6

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

g eq a ðx; tÞ ¼ xa C

ð15Þ P

where the oxygen concentration C ¼ g a , and the relationship between dimensionless relaxation time s and lattice diffusivity D of oxygen is given as:

D ¼ c2s ðs  0:5Þ

ð16Þ

As suggested in Section 3.2, when calculating the Knudsen diffusion coefficient of oxygen in the pores, the effective pore diameter dp is required. The calculation of dp in the lattice computational domain is organized as this process: first, determine whether the current lattice point is a pore lattice point, and if yes, start the diameter calculation, for lattice structure of D3Q19, there are 18 adjacent lattice points around a pore lattice point. Second, expand the current lattice point along the 18 adjacent directions until the non-pore lattice point is encountered, calculate the expanded distance in each direction to represent the pore radius in this direction, then the effective pore diameter at this lattice coordinate is calculated as:

2 dp ¼

18 P

a¼1

the model’s efficiency, accuracy and stability. For the research on oxygen transport and electrochemical reaction, there are several important unit conversion coefficients: length scale C l (m), concentration scale C c (mol m3), diffusivity scale C D (m2 s1), time scale C t (s), and source term scale C S (mol m3 s1), specifically, (1) Length scale C l : As shown in Fig. 1, the physical dimension of computational domain is selected as 200 nm  200 nm  2000 nm (the first two numbers represent cross-sectional dimensions, and the third number is thick-

Ra

18

ð17Þ

3.4. Unit conversion and boundary condition In the LB simulations, choosing appropriate unit conversion between lattice units and physical properties is critical to improve

Fig. 3. (a) Distribution probability of pore diameter and (b) number of effective Pt particles along the thickness direction with different mass ratios between ionomer and electrode (hion =0.2, 0.3, 0.4).

Fig. 4. Variations of average oxygen concentration in the (a) pore and (b) ionomer as flow time evolves, and (c) oxygen consumption rate along the thickness direction at equilibrium state with different mass ratios between ionomer and electrode (hion =0.2, 0.3, 0.4).

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

ness). To take account into the micromorphology of carbon support, ionomer and Pt powder, the length scale C l ¼ 2  109 m is adopted to spatial discretization, corresponding to the lattice size of 100  100  1000 . Note that when the length scale C l has been determined, the unit conversion of length in lattice and real systems is established, and thus the consistency of Knudsen number is enforced   P KnL ¼ dkp;LL ¼ dkp;LL CCl ¼ dkp;P ¼ KnP . l

(2) Concentration scale C c : C c can be defined as the arbitrary value (take C c ¼ 1 mol m3 in this study), because it only affects the source term scale, and is not related to the other three conversion coefficients. (3) Diffusivity scale C D : C D is related to the dimensionless relaxation time s (Eq. (16)) and thus affects the model stability. s should be carefully chosen, and first, select a s that is small enough and meet the stability requirement, corresponding to the minimum value of oxygen diffusivity in the whole

7

computational domain, and then C D is calculated. Second, calculate the oxygen diffusivity of each lattice point in the pore or ionomer, then inversely calculate the s by using C D and Eq. (16), to establish the correspondence between s and local oxygen diffusivity. According to Table 1, the oxygen diffusivity in the ionomer is minimum and the value is calculated as 7:608  1010 m2 s1 (60 °C), corresponding to s = 0.556 and then diffusivity scale can be calculated as C D ¼ 4:076  108 m2 s1 . (4) Time scale C t : C t is related to length and diffusivity scales, and satisfies C t ¼ C l 2 =C D . Therefore, time scale is calculated as C t ¼ 9:814  1011 s based on the determined C l and C D . (5) Source term scale C S : The unit of oxygen source term in Eq. (6) is mol m3 s1, which implies that C S is related to concentration and time scales, and satisfies C S ¼ C c =C t . Source term scale is calculated as C S ¼ 1:019  1010 mol m3 s1 based on the determined C l and C D .

Fig. 5. Oxygen concentration distribution at equilibrium state with different mass ratios between ionomer and electrode: (a) hion = 0.2; (b) hion = 0.3; (c) hion = 0.4.

8

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

For the boundary definition, periodic boundary is applied to the four sides of computational domain, and the internal solid structure (Pt and carbon lattices) is defined as the half-bounce back boundary. The flow boundaries are applied at both sides, the left side is the oxygen concentration inlet (acts as the interface between GDL or MPL and CL), and the right side is set as open boundary (represents a section inside the CL). At the left concentration inlet (x = 0), the oxygen concentration is fixed to C O2 ;in ¼ 0:21P=RT (details of operating conditions are given in Table 2), thus the unknown distribution functions at the inlet are obtained by:

20 nm, corresponding to the fact that an increase in the ionomer content will cause the loss of pore size. In addition to the reduction of pore size, the variation of ionomer content also affects the number of effective Pt particles, and ‘‘effective” means the Pt particle is partially or fully covered by ionomer. Fig. 3(b) shows the number variation of effective Pt particles along the thickness direction, and the results show that when hion ¼ 0:2; 0:3; 0:4, the number of effective Pt particles fluctuates around 4000, 4500 and 5000, respectively, which indicates that more Pt particles are covered as the ionomer content increases.

2 g 1 ð0; y; zÞ ¼ 18 C O2 ;in  g 4 ð0; y; zÞ 2 g 7 ð0; y; zÞ ¼ 36 C O2 ;in  g 10 ð0; y; zÞ 2 g 8 ð0; y; zÞ ¼ 36 C O2 ;in  g 11 ð0; y; zÞ

g 13 ð0; y; zÞ ¼ g 14 ð0; y; zÞ ¼

2 C 36 O2 ;in 2 C 36 O2 ;in

ð18Þ

 g 16 ð0; y; zÞ  g 17 ð0; y; zÞ

At the last layer of lattice (x = N  1), oxygen concentration in the pores at this interface satisfies the Neumann open boundary, thus the unknown distribution functions at the right boundary are given as follows:

g 4 ðN  1; y; zÞ ¼ g 4 ðN  2; y; zÞ g 10 ðN  1; y; zÞ ¼ g 10 ðN  2; y; zÞ g 11 ðN  1; y; zÞ ¼ g 11 ðN  2; y; zÞ

ð19Þ

g 16 ðN  1; y; zÞ ¼ g 16 ðN  2; y; zÞ g 17 ðN  1; y; zÞ ¼ g 17 ðN  2; y; zÞ In this study, the self-developed LB code, which is based on the single relaxation time model introduced in Section 3.3, is parallelized with the compute unified device architecture (CUDA) platform and NVIDIA Titan X GPUs. 4. Results and discussion Based on the reconstructed CL microstructures and developed single-phase LB model, the effects of three critical structural parameters on micromorphology, oxygen transport and electrochemical reaction are studied, which includes ionomer and electrode mass ratio (hion ), Pt loading (cPt ) and carbon particle diameter (DC ). The basic and operating conditions are given in Table 3. 4.1. Ionomer and electrode mass ratio In this section, CL microstructures with different mass ratios of ionomer to electrode are reconstructed (Fig. 2), and other structural parameters used in the reconstruction are consistent with the basic condition in Table 3. When the Pt loading cPt and Pt/C mass ratio hPt=C are kept constant, it can be known from Eqs. (1)– (4) that the variation of the mass ratio between ionomer and electrode hion does not change the volume fraction of carbon and Pt phases (eC , ePt ), and only affects the volume fraction of ionomer phase eion , resulting in the varied total porosity e. when hion ¼ 0:2; 0:3; 0:4, it can be calculated that eC is equal to 0.34710, ePt is equal to 0.0233, eion are 0.14056, 0.24083 and 0.37472, and e are 0.48904, 0.38877 and 0.25488, respectively. As the ionomer content increases, the remaining pore space will shrink and more Pt particles are covered (Fig. 2). Fig. 3(a) shows the distribution probabilities of pore diameter in 10 numerical spaces: 1–10, 11–20, 21–30, 31–40, 41–50, 51–60, 61–70, 71–80, 81–90 and 91–100 nm. For all three cases, the probability of pore diameter larger than 60 nm is close to 0. The higher the ionomer content, the greater the probability that pore diameter is less than

Fig. 6. Oxygen concentration along the thickness direction at different time moments with different mass ratios between ionomer and electrode: (a) hion = 0.2; (b) hion = 0.3; (c) hion = 0.4.

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

Based on the computational domains obtained from CL reconstruction, the effect of ionomer content on oxygen transport and electrochemical reaction are further investigated. Fig. 4 shows the variations of average oxygen concentration in the CL pore

9

and ionomer as time evolves, as well as the oxygen consumption rate along the thickness direction at the equilibrium state. It can be seen from Fig. 4(a) and (b) that the average oxygen concentrations in the pore and ionomer show synchronous change for all

Fig. 7. Reconstructed CL microstructures with different Pt loadings: (a) cPt = 0.06; (b) cPt = 0.08; (c) cPt = 0.10; (d) cPt = 0.12 mg cm2.

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

4.2. Pt loading Carbon-supported Pt is the typical structure of the catalyst powder, by dispersing Pt particles at the surface of carbon support through physical adsorption. The mass ratio of Pt to carbon hPt=C in the same batch of catalyst powder is generally fixed, and does not change significantly as the amount of catalyst powder varies, which only affects the CL thickness. Based on the stochastic reconstruction model, the CL microstructures with four different Pt loadings (cPt =0.06, 0.08, 0.10, 0.12 mg cm2) are reconstructed (Fig. 7),

and the thickness of 2000 nm is fixed in four computational domains, as well as same mass ratios of ionomer to electrode and Pt to carbon (hion and hPt=C ). According to Eqs. (1)–(4), the change in Pt loading simultaneously affects the volume fractions of ionomer, carbon and Pt phases (eion , eC , ePt ) and total porosity e. WhencPt = 0.06, 0.08, 0.10, 0.12 mg cm2, it can be calculated that eC are 0.20821, 0.27762, 0.34710 and 0.41643, ePt are 0.01398, 0.01864, 0.02330 and 0.02796, eion are 0.14455, 0.19272, 0.24083 and 0.28908, and e are 0.63326, 0.51102, 0.38877 and 0.26653, respectively. As cPt increases and to maintain the unchanged hion , more pore space will be occupied by carbon-supported Pt and ionomer, resulting in a decrease of CL porosity. Fig. 8(a) shows that most of pores are smaller than 20 nm when cPt ¼ 0:12 mg cm2 , while as Pt loading decreases, these small pores are enlarged and even a fair amount of them are close to 100 nm when cPt ¼ 0:06 mg cm2 . Although the probability of single Pt particle being covered by ionomer is close in four cases due to the fixed hion , the total amount of effective Pt particles still increases with increment of Pt loading (Fig. 8(b)). The difference in pore diameter and amount of effective Pt particles, directly affects the oxygen diffusion and electrochemical reaction rate (Fig. 9). The reduced Pt loading empties more pore space for oxygen transport, which enhances the Knudsen diffusion and improves the effective diffusion coefficient of oxygen, corresponding to a considerable improvement of oxygen concentration in equilibrium state (Fig. 9(a) and (b)). Moreover, the increment of oxygen concentration at low Pt loading is more pronounced,

60

γPt = 0.06 mg cm-2 50

γPt = 0.08 mg cm-2 γPt = 0.10 mg cm-2

40

Frequency, %

three cases. At the beginning stage, the internal oxygen concentration is small, thus oxygen rapidly diffuses into the CL interior under the drive of large concentration gradient. As the overall oxygen concentration increases, the oxygen diffusion gradually slows down. The consistent variation trend of oxygen concentration in the pore and ionomer is attributed to the Henry’s law (Eq. (8)), which defines the linear relationship of gas concentration at the two sides of interface. The fluctuations of oxygen concentration curves in the middle and second half are caused by the irregular pore structure, which leads to the uneven transmission of oxygen in the plane and through the plane, and this effect becomes significant as the oxygen diffusion gradually slows down. Moreover, higher ionomer content lowers the average oxygen concentration, because the increased ionomer occupies more pore space, resulting in the weakened Knudsen diffusion and thus the reduced effective diffusion coefficient of oxygen (Table 1). Fig. 5 clearly shows that oxygen transport becomes more difficult at higher electrolyte content, corresponding to the smaller penetration distance in the thickness direction. Fig. 6 gives the comparison of oxygen penetration distance and concentration along the thickness direction at different time moments, the sparser the curve is, the easier for oxygen to transport through this position, resulting in the faster growing oxygen concentration. The stagger phenomenon is observed in the oxygen reaction rate curves as shown in Fig. 4(c). At the positions that thickness is smaller than 700 nm, higher ionomer content accelerates the oxygen reaction consumption, while for positions that thickness is larger than 700 nm, this trend is reversed. This phenomenon may be explained by the fact that although lower ionomer content is beneficial to improving the average oxygen concentration in the ionomer and thus is favorable to ORR, oxygen consumption rate is not determined solely by the electrochemical reaction rate. According to the definition of oxygen source term in Table 1, oxygen reaction consumption rate is also affected by the number of adjacent Pt particles, that is, the trigger condition for oxygen consumption is the presence of Pt around ionomer lattice point, and the more Pt lattice points, the stronger the catalytic effect and the faster the oxygen consumption. For the case with higher ionomer content, more Pt particles can be covered by the ionomer, representing the high potential catalyst utilization. At the position close to the oxygen concentration inlet, the concentration difference for three cases is not obvious due to the facile diffusion (Fig. 6). Therefore, the number of effective Pt particles is dominated for oxygen consumption, and higher ionomer content improves the oxygen consumption rate. As the thickness position advances, the difference in oxygen concentration gradually manifests (Fig. 6) and the effect of oxygen concentration becomes important, resulting in a reversal of the oxygen consumption rate. It can be seen that at the position close to outlet, the oxygen reaction rate drops rapidly to zero because oxygen concentration approaching zero for hion ¼ 0:4 case, indicating that in this region, although the case with high ionomer content can cover more Pt particles, the catalysts are not effectively utilized due to the insufficient oxygen.

γPt = 0.12 mg cm-2

30 20 10 0 0

20

(a)

40

60

80

100

Pore diameter, nm

8000

γPt = 0.12 mg cm-2

7000

Covered Pt number

10

6000

γPt = 0.10 mg cm-2 5000 4000

γPt = 0.08 mg cm-2

3000 2000

γPt = 0.06 mg cm-2

1000 0

(b)

400

800

1200

1600

2000

Thickness, nm

Fig. 8. (a) Distribution probability of pore diameter and (b) number of effective Pt particles along the thickness direction with different Pt loadings (cPt =0.06, 0.08, 0.10, 0.12 mg cm2).

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

11

10

O2 concentration in pore, mol m-3

γPt = 0.06 mg cm-2 8

6

γPt = 0.08 mg cm-2

4

γPt = 0.10 mg cm-2

2

γPt = 0.12 mg cm-2

0 0.0

0.2

0.6

0.8

1.0

0.8

1.0

Time, ms

(a)

O2 concentration in ionomer, mol m-3

0.4

γPt = 0.06 mg cm-2

0.20

0.16

γPt = 0.08 mg cm-2 0.12

γPt = 0.10 mg cm-2

0.08

γPt = 0.12 mg cm-2

0.04

0.00 0.0

(b)

0.2

0.4

0.6

Time, ms

Fig. 9. Variations of average oxygen concentration in the (a) pore and (b) ionomer as flow time evolves, and (c) oxygen consumption rate along the thickness direction at equilibrium state with different Pt loadings (cPt =0.06, 0.08, 0.10, 0.12 mg cm2).

which is attributed to the increased pores with larger size and thus the significant speed up of oxygen transport along the thickness direction. This phenomenon is also observed in Fig. 10, according to the intrusion distance of oxygen along the thickness direction at different time moments. At low Pt loading, the spacing of curves is sparser, indicating that oxygen is easier to reach and pass through the current position. As the Pt loading increases, the reduced effective diffusion coefficient is insufficient to support oxygen reaching the back region. The significant difference in diffusion capability for four Pt loading cases, resulting in the distinction of their concentration gradients along the thickness direction, thus as Pt loading increases, the uniformity of concentration distribution along the thickness direction decreases.

Fig. 10. Oxygen concentration along the thickness direction at different time moments with different Pt loadings: (a) cPt = 0.06; (b) cPt = 0.08; (c) cPt = 0.10; (d) cPt = 0.12 mg cm2.

The probability of Pt particles covered by ionomer along the thickness direction, is almost the same for four Pt loading cases,

12

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

the uniformities of oxygen consumption rate are determined by the difference in oxygen concentration. It can be seen from Fig. 9 (c) that decreasing Pt loading is beneficial to improving the uniformity of oxygen consumption rate along the thickness direction,

and when cPt ¼ 0:06 mg cm2 , the curve fluctuates along the horizontal line. Similar to the results of ionomer content cases, the curves of oxygen consumption rate with different Pt loadings also intersect. In the area close to oxygen concentration inlet, the

Fig. 11. Reconstructed CL microstructures with different diameter ranges of carbon particle: (a) 12–20 nm (100%); (b) 20–32 nm (100%); (c) 32–40 nm (100%); (d) 12–20 nm (50%) and 32–40 nm (50%).

13

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

4.3. Carbon particle diameter In this section, four CL microstructures are reconstructed with different carbon particle diameters (Fig. 11): 12–20 nm (100%), 20–32 nm (100%), 32–40 nm (100%), 12–20 nm (50%) and 32– 40 nm (50%), and the number in parentheses represents the percentage of carbon particles in this diameter range. The thickness of 2000 nm, mass ratios of ionomer to electrode hion and Pt to carbon hPt=C , and Pt loading cPt are fixed in four cases, resulting in the consistent volume fractions of ionomer, carbon and Pt phases (eion , eC , ePt ), as well as the total porosity e based on Eqs. (1)–(4). Fig. 12(a) gives the distribution probability of pore diameter, and the results show that almost all pore diameters are less than 30 nm in the case of 12–20 nm (100%). As the carbon particle

diameter increases, the shaped pore size gradually migrates to the increased diameter range. For the case of 32–40 nm (100%), a small amount of pore diameter is even close to 80 nm, and for the case of 12–20 nm (50%) and 32–40 nm (50%), the pore distribution is close to the case of 32–40 m (100%), while the blending of carbon particles with smaller size will lower the number of large pores. It can be seen from Fig. 13 that increasing carbon particle diameter is beneficial to improving the oxygen concentration in both the ionomer and pore, because the increased large pores enhance the Knudsen diffusivity of oxygen and transport becomes more easier, which can also be observed from the curves of concentration

4

O2 concentration in pore, mol m-3

oxygen concentration difference in four cases is insignificant, and the amount of effective Pt particles plays a dominate role, explains why the oxygen consumption with high Pt loading is faster. As the thickness position moves backward, the difference in oxygen concentration and its effect on the electrochemical reaction rate become significant. The high Pt loading condition is dragged down by the reduced oxygen concentration, thus the oxygen consumption rate gradually decreases and is even exceeded by low Pt loading case, which greatly causes the waste of Pt catalyst. The highest value of oxygen consumption rate is observed at the Pt loading of cPt ¼ 0:08 mg cm2 , indicating that both the enhancement of oxygen diffusion and effective reaction area are critical to the improvement of cell performance.

3

DC=12~20 nm (100%)

2

DC=20~32 nm (100%) DC=32~40 nm (100%)

1

DC=12~20 nm (50%) 32~40 nm (50%)

0 0.0

0.2

0.4

0.6

0.8

1.0

Time, ms

(a) O2 concentration in ionomer, mol m-3

0.10

0.08

0.06 DC=12~20 nm (100%)

0.04

DC=20~32 nm (100%) DC=32~40 nm (100%)

0.02

DC=12~20 nm (50%) 32~40 nm (50%)

0.00 0.0

0.2

0.4

0.6

0.8

1.0

Time, ms

(b) O2 consumption rate, 108 mol m-3 s-1

6 5

DC=12~20 nm (100%) Mean Value 1.166 DC=20~32 nm (100%) 1.581 D =32~40 nm (100%)

4

DC=12~20 nm (50%)

1.934

32~40 nm (50%)

1.720

t=0.981 ms

C

3 2 1 0 0

(c) Fig. 12. (a) Distribution probability of pore diameter and (b) number of effective Pt particles along the thickness direction with different diameter ranges of carbon particle (12–20 nm (100%); 20–32 nm (100%); 32–40 nm (100%); 12–20 nm (50%) and 32–40 nm (50%)).

400

800

1200

1600

2000

Thickness, nm

Fig. 13. Variations of average oxygen concentration in the (a) pore and (b) ionomer as flow time evolves, and (c) oxygen consumption rate along the thickness direction at equilibrium state with different diameter ranges of carbon particle (12–20 nm (100%); 20–32 nm (100%); 32–40 nm (100%); 12–20 nm (50%) and 32–40 nm (50%)).

14

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

0.25

DC=12~20 nm (100%) 0.20

0.15

s

1m

0.10

,

88

0.1 43 0.24 m 2m s s

0.05

0.00 0

400

8 0.9

0.5

0.3

800

(a)

1200

1600

2000

Thickness, nm 0.25

DC=20~32 nm (100%) 0.20

0.15

s

1m

0.10

,

8 0.9

88

0.05

0.5

, 90

0.1 43 0.24 ms 2 m s

s m 44 04 0. 7 ms 24 0.0 ms 8e-3 4.93

O2 concentration in ionomer, mol m-3

, 90

s m 44 04 s 0. m 7 24 0 .0 s -3 m 38e 4.9

O2 concentration in ionomer, mol m-3

distribution along the thickness direction at different time moments (Fig. 14). The average oxygen concentration in the case

of 12–20 nm (50%) and 32–40 nm (50%) is close to the case of 32– 40 nm (100%), indicating that further increasing carbon particle diameter has no obvious benefit for oxygen transport when a certain number of large pores is formed. And note that the adjustment of carbon particle diameter only affects the pore size and will not create additional space, because the Pt loading, and mass ratios of ionomer to electrode and Pt to carbon are fixed in four cases, explaining why the oxygen concentration difference in this section is not as obvious as the previous group cases. The average number of effective Pt particles along the thickness direction in four cases are 4459.9, 4655.2, 4795.5 and 4512.0, respectively, and satisfies: 32–40 nm (100%) > 20–32 nm (100%) > 12–20 nm (50%) and 32–40 nm (50%) > 12–20 nm (100%) (Fig. 12(b)). In summary, smaller carbon particle can extend the surface area to carry Pt particle and ionomer, which means catalyst is difficult to be covered under the same ionomer content and Pt loading. This result indicates that in order to obtain sufficient three-phase interfaces and continuous proton transport channels, the ionomer content should increase as the carbon particle diameter decreases, which is also reported in the experiment literature conducted by Soboleva et al. [48]. Moreover, the interlaced phenomenon is not observed for the curves of oxygen consumption rate in four cases (Fig. 13(c)), which is attributed to the enlarged pores and more effective coverage of catalyst as the carbon particle diameter increases, the former enhances the oxygen transport and the latter brings more reactive sites, and both of them are beneficial to improving the oxygen consumption rate.

0.3

5. Conclusion

0.00 0

400

800

1200

1600

2000

Thickness, nm 0.25

DC=32~40 nm (100%) 0.20

0.15

81

ms

.9 ,0

0.10

88

0.05

0.5

, 90

0.1

s m 44 04 s 0. 7m 24 0.0 ms -3 38e 4.9

O2 concentration in ionomer, mol m-3

(b)

43

ms

0.24

0.3

2m

s

0.00 0

800

1200

1600

2000

Thickness, nm 0.25

DC=12~20 nm (50%) 32~40 nm (50%)

0.20

0.15

81

0.10

ms

.9 ,0

8

.58

,0

0.05

0.00

(d)

400

s m 44 04 s 0. m 7 24 0.0 ms -3 38e 4.9

O2 concentration in ionomer, mol m-3

(c)

0

400

0.1 43 0.242 ms ms

800

1200

0 .39

0

1600

2000

Thickness, nm

Fig. 14. Oxygen concentration along the thickness direction at different time moments with different diameter ranges of carbon particle: (a) 12–20 nm (100%); (b) 20–32 nm (100%); (c) 32–40 nm (100%); (d) 12–20 nm (50%) and 32–40 nm (50%).

In this study, different catalyst layer microstructures, with the varied mass ratios of ionomer to electrode, platinum (Pt) loadings and carbon particle diameters, are reconstructed based on the improved stochastic algorithm, which considers Pt particle, carbon support and ionomer simultaneously, and the pore size distribution and active sites number are analyzed. A single-phase lattice Boltzmann (LB) model is then applied to study the effect of micromorphology on the oxygen transport process and electrochemical reaction. The results suggested that the average oxygen concentrations in the pore and ionomer show synchronous change, that is, increase rapidly at the beginning stage and then gradually slow down, which are attributed to the evolution of concentration gradient and Henry’s law. Decreasing the ionomer content and Pt loading is beneficial to the oxygen transport along the thickness direction, due to the enlarged pore space and enhanced Knudsen diffusion, resulting in the increased oxygen concentrations in the pore and ionomer, as well as the longer oxygen penetration distance. However, the resulted reduction of effective Pt particles complicates the oxygen consumption rate, in the regions close to oxygen inlet, the oxygen consumption rate is dominated by the number of effective Pt particles, due to the facile oxygen diffusion, thus oxygen consumes slower for the cases of lower ionomer content and Pt loading. As the thickness position advances, the oxygen concentration rapidly decreases, resulting in a reversal of the oxygen consumption rate. The interlaced phenomenon is not observed in carbon particle diameter cases, because increasing the carbon particle diameter improves both the pore space and Pt particle coverage. The presented results in this study are helpful to understanding the transport process and reaction kinetics inside CL and are expected to guide the electrode microstructure designs. Declaration of Competing Interest The authors declared that there is no conflict of interest.

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

15

Acknowledgements

References

This research is supported by the National Key Research and Development Program of China (Grant No. 2017YFB0102703), 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), and the National Natural Science Foundation of Tianjin (China) for Distinguished Young Scholars (Grant No. 18JCJQJC46700).

[1] Fuel Cell Technical Team Roadmap. U.S. Drive, 2013. [2] E.H. Majlan, D. Rohendi, W.R.W. Daud, T. Husaini, M.A. Haque, Electrode for proton exchange membrane fuel cells: a review, Renew. Sustain. Energy Rev. 89 (2018) 117–134. [3] D. Banham, S. Ye, K. Pei, J. Ozaki, T. Kishimoto, Y. Imashiro, A review of the stability and durability of non-precious metal catalysts for the oxygen reduction reaction in proton exchange membrane fuel cells, J. Power Sources 285 (2015) 334–348. [4] M. Kiani, J. Zhang, Y. Luo, C. Jiang, J. Fan, G. Wang, J. Chen, R. Wang, Recent developments in electrocatalysts and future prospects for oxygen reduction reaction in polymer electrolyte membrane fuel cells, J. Energy Chem. 27 (2018) 1124–1139. [5] K. Jiao, X. Li, Water transport in polymer electrolyte membrane fuel cells, Prog. Energy Combust. Sci. 37 (2011) 221–291. [6] H. Wu, A review of recent development: transport and performance modeling of PEM fuel cells, Appl. Energy 165 (2016) 81–106. [7] 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. [8] D. Harvey, J.G. Pharoah, K. Karan, A comparison of different approaches to modelling the PEMFC catalyst layer, J. Power Sources 179 (2008) 209–219. [9] G.R. Molaeimanesh, H.S. Googarchin, A.Q. Moqaddam, Lattice Boltzmann simulation of proton exchange membrane fuel cells – a review on opportunities and challenges, Int. J. Hydrogen Energy 41 (2016) 22221–22245. [10] Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, Q. Liu, Lattice Boltzmann methods for multiphase flow and phase-change heat transfer, Prog. Energy Combust. Sci. 52 (2016) 62–105. [11] L. Chen, Q. Kang, Y. Mu, Y.L. He, W.Q. Tao, A critical review of the pseudopotential multiphase lattice Boltzmann model: methods and applications, Int. J. Heat Mass Transf. 76 (2014) 210–236. [12] A. Xu, L. Shi, T.S. Zhao, Accelerated lattice Boltzmann simulation using GPU and OpenACC with data management, Int. J. Heat Mass Transf. 109 (2017) 577– 588. [13] A. Xu, W. Shyy, T.S. Zhao, Lattice Boltzmann modeling of transport phenomena in fuel cells and flow batteries, Acta Mech. Sin. 33 (2017) 555–574. [14] L. Hao, P. Cheng, Capillary pressures in carbon paper gas diffusion layers having hydrophilic and hydrophobic pores, Int. J. Heat Mass Transf. 55 (2012) 133–139. [15] W.Z. Fang, Y.Q. Tang, L. Chen, Q.J. Kang, W.Q. Tao, Influences of the perforation on effective transport properties of gas diffusion layers, Int. J. Heat Mass Transf. 126 (2018) 243–255. [16] M. Espinoza-Andaluz, M. Andersson, B. Sunden, Comparing through-plane diffusibility correlations in PEFC gas diffusion layers using the lattice Boltzmann method, Int. J. Hydrogen Energy 42 (2017) 11689–11698. [17] J. Liu, S. Shin, S. Um, Comprehensive statistical analysis of heterogeneous transport characteristics in multifunctional porous gas diffusion layers using lattice Boltzmann method for fuel cell applications, Renewable Energy 139 (2019) 279–291. [18] S. Shin, A.R. Kim, S. Um, Computational prediction of nanoscale transport characteristics and catalyst utilization in fuel cell catalyst layers by the lattice Boltzmann method, Electrochim. Acta 275 (2018) 87–99. [19] A. Xu, T.S. Zhao, L. Shi, J.B. Xu, Lattice Boltzmann simulation of mass transfer coefficients for chemically reactive flows in porous media, J. Heat Transfer 140 (2018), 052601-1–8. [20] F. Jinuntuya, M. Whiteley, R. Chen, A. Fly, The effects of gas diffusion layers structure on water transportation using X-ray computed tomography based Lattice Boltzmann method, J. Power Sources 378 (2018) 53–65. [21] D.H. Jeon, H. Kim, Effect of compression on water transport in gas diffusion layer of polymer electrolyte membrane fuel cell using lattice Boltzmann method, J. Power Sources 294 (2015) 393–405. [22] D.H. Jeon, Effect of channel-rib width on water transport behavior in gas diffusion layer of polymer electrolyte membrane fuel cells, J. Power Sources 423 (2019) 280–289. [23] G.R. Molaeimanesh, M.H. Akbari, Role of wettability and water droplet size during water removal from a PEMFC GDL by lattice Boltzmann method, Int. J. Hydrogen Energy 41 (2016) 14872–14884. [24] A.H. Kakaee, G.R. Molaeimanesh, M.H. Elyasi Garmaroudi, Impact of PTFE distribution across the GDL on the water droplet removal from a PEM fuel cell electrode containing binder, Int. J. Hydrogen Energy 43 (2018) 15481–15491. [25] K.N. Kim, J.H. Kang, S.G. Lee, J.H. Nam, C.J. Kim, Lattice Boltzmann simulation of liquid water transport in microporous and gas diffusion layers of polymer electrolyte membrane fuel cells, J. Power Sources 278 (2015) 703–717. [26] H. Ostadi, K. Jiang, P.D. Prewett, Micro/nano X-ray tomography reconstruction fine-tuning using scanning electron microscope images, Micro Nano Lett. 3 (2008) 106–109. [27] L. Chen, H.B. Luan, Y.L. He, W.Q. Tao, Pore-scale flow and mass transport in gas diffusion layer of proton exchange membrane fuel cell with interdigitated flow fields, Int. J. Therm. Sci. 51 (2012) 132–144. [28] L. Chen, H. Luan, Y. Feng, C. Song, Y.L. He, W.Q. Tao, Coupling between finite volume method and lattice Boltzmann method and its application to fluid flow and mass transport in proton exchange membrane fuel cell, Int. J. Heat Mass Transf. 55 (2012) 3834–3848.

Appendix A. Nomenclature C unit conversion coefficient; concentration (mol m3) cs lattice sound speed d diameter (m) D diffusivity (m2 s1) e discrete velocity F Faraday constant (C mol1) g distribution function H Henry constant i current density (A m2) k Boltzmann constant (J K1) L length (m) p pressure (Pa) R ideal gas constant (J K1 mol1); radius S source term (mol m3 s1) t time (s) T temperature (K) u Velocity (m s1) w weight coefficient x, y, z spatial position Greek letters a reaction transfer coefficient; c Pt loading (mg cm2) d thickness (cm) e volume fraction h mass ratio k mean molecular free path (m) q density (mg cm3) s relaxation time u overpotential (V) x weight coefficient Subscripts and superscripts c concentration C carbon D diffusivity e electrolyte eq equilibrium in inlet ion ionomer Kn Knudsen l length L lattice unit O2 oxygen p pore P physical unit Pt platinum Pt/C platinum to carbon ref reference S source term t time

16

H. Deng et al. / International Journal of Heat and Mass Transfer 143 (2019) 118538

[29] L. Chen, G. Wu, E.F. Holby, P. Zelenay, W.Q. Tao, Q. Kang, Lattice Boltzmann pore-scale investigation of coupled physical-electrochemical processes in C/Pt and non-precious metal cathode catalyst layers in proton exchange membrane fuel cells, Electrochim. Acta 158 (2015) 175–186. [30] G.R. Molaeimanesh, M.A. Bamdezh, M. Nazemian, Impact of catalyst layer morphology on the performance of PEM fuel cell cathode via lattice Boltzmann simulation, Int. J. Hydrogen Energy 43 (2018) 20959–20975. [31] G.R. Molaeimanesh, M.H. Akbari, Agglomerate modeling of cathode catalyst layer of a PEM fuel cell by the lattice Boltzmann method, Int. J. Hydrogen Energy 40 (2015) 5169–5185. [32] G.R. Molaeimanesh, M.H. Akbari, A pore-scale model for the cathode electrode of a proton exchange membrane fuel cell by lattice Boltzmann method, Korean J. Chem. Eng. 32 (2015) 397–405. [33] L. Chen, R. Zhang, P. He, Q. Kang, Y.L. He, W.Q. Tao, Nanoscale simulation of local gas transport in catalyst layers of proton exchange membrane fuel cells, J. Power Sources 400 (2018) 114–125. [34] D. Zhang, Q. Cai, S. Gu, Three-dimensional lattice-Boltzmann model for liquid water transport and oxygen diffusion in cathode of polymer electrolyte membrane fuel cell with electrochemical reaction, Electrochim. Acta 262 (2018) 282–296. [35] X. Li, J. Cai, F. Xin, X. Huai, J. Guo, Lattice Boltzmann simulation of endothermal catalytic reaction in catalyst porous media, Appl. Therm. Eng. 50 (2013) 1194– 1200. [36] M.R. Kamali, S. Sundaresan, H.E.A. Van den, J.J.J. Gillissen, A multi-component two-phase lattice Boltzmann method applied to a 1-D Fischer-Tropsch reactor, Chem. Eng. J. 207 (2012) 587–595. [37] K.J. Lange, P.C. Sui, N. Djilali, Pore scale modeling of a proton exchange membrane fuel cell catalyst layer: Effects of water vapor and temperature, J. Power Sources 196 (2011) 3195–3203.

[38] Y. Liu, M.W. Murphy, D.R. Baker, W. Gu, C. Ji, J. Jorne, H.A. Gasteiger, Proton conduction and oxygen reduction kinetics in PEM fuel cell cathodes: effects of ionomer-to-carbon ratio and relative humidity, J. Electrochem. Soc. 156 (2009) B970–980. [39] P. Choi, N.H. Jalani, R. Datta, Thermodynamics and proton transport in Nafion II. Proton diffusion mechanisms and conductivity, J. Electrochem. Soc. 152 (2005) E123–130. [40] N.A. Siddique, F. Liu, Process based reconstruction and simulation of a threedimensional fuel cell catalyst layer, Electrochim. Acta 55 (2010) 5357–5366. [41] W. Sun, B.A. Pepple, K. Karan, An improved two-dimensional agglomerate cathode model to study the influence of catalyst layer structural parameters, Electrochim. Acta 50 (2005) 3359–3374. [42] S.H. Kim, H. Pitsch, Reconstruction and effective transport properties of the catalyst layer in PEM fuel cells, J. Electrochem. Soc. 156 (2009) B673–681. [43] H. Deng, D. Wang, X. Xie, Y. Zhou, Y. Yin, Q. Du, K. Jiao, Modeling of hydrogen alkaline membrane fuel cell with interfacial effect and water management optimization, Renewable Energy 91 (2016) 166–177. [44] W.G. Pollard, R.D. Present, On gaseous self-diffusion in long capillary tubes, Phys. Rev. 73 (1948) 762–774. [45] E.L. Cussler, Diffusion: Mass Transfer in Fluid Systems, third ed., Cambridge University, New York, 1997. [46] R.E. Cunningham, R.J.J. Williams, Diffusion in Gases and Porous Media, Plenum, New York, 1980. [47] K.J. Lange, P.C. Sui, Pore scale simulation of transport and electrochemical reactions in reconstructed PEMFC catalyst layers, J. Electrochem. Soc. 157 (2010) B1434–1442. [48] T. Soboleva, K. Malek, Z. Xie, T. Navessin, S. Holdcroft, PEMFC catalyst layers: the role of micropores and mesopores on water sorption and fuel cell activity, ACS Appl. Mater. Interfaces 3 (2011) 1827–1837.