Using MRT lattice Boltzmann method to simulate gas flow in simplified catalyst layer for different inlet–outlet pressure ratio

International Journal of Heat and Mass Transfer 88 (2015) 122–132

Contents lists available at ScienceDirect

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

Using MRT lattice Boltzmann method to simulate gas flow in simplified catalyst layer for different inlet–outlet pressure ratio Yuan Gao Clean Energy Automotive Engineering Centre & School of Automotive Studies, Tongji University, China

a r t i c l e

i n f o

Article history: Received 9 September 2014 Received in revised form 2 April 2015 Accepted 11 April 2015

Keywords: PEM fuel cells Pore-scale modelling Catalyst layers Lattice Boltzmann model Flow rate Permeability

a b s t r a c t Gas flow in catalyst layer is important as it is the place where the electrochemical reactions take place. But it is challenging because the size of the pores is less than one micron and gas flow in the catalyst layer is no longer continuum, the Knudsen number cannot be neglected. As a result, the wall–gas collision must be re-considered, meaning that in the macroscopic models, the absolute permeability of the catalyst layer for different gases varies. To investigate this, we will simplify the pore geometry in catalyst layer into a bundle of tubes whose diameters can be derived from the pore-size distribution of the 3D (focused ion beam) FIB images. A model for gas flow in each tube is then simulated, the flow rate in the small pores have high resistance than big one, large difference start after pore size bigger than 100 nm. Use the pore volume percentage with different pore size diameter and porosity, we have calculate the average cross section area of the catalyst layer and calculate the average flow rate based on them which is used for the permeability calculation later. This permeability value can reflect the permeability in the true catalyst. The results show that permeability of catalyst layer is not a constant but varies with Knudsen number, meaning that the permeability calculation in the micro/nano scale is depend on the Knudsen number. Assuming a constant permeability for all the gases, as used in the available fuel cell models in literature, could give rise to significant errors. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Fuel cells based on proton exchange membrane (PEM) and fuelled by hydrogen and air have many attractive features, including high power density, rapid start-up and high efficiency. It has received much attention as one of the leading candidates for the power sources of stationary, mobile and portable device [1]. The catalyst layer, probably, the most important part of the fuel cell is adjacent to the gas diffusion layer. It is the layer where electrochemical reactions take place on the surface of the catalyst [2,3] and hence plays a central role in cell performance [4]. It is believed that the associated species transport in the catalyst layer is the key limitation in cell performance and hence the permeability prediction of catalyst layer, especially in the cathode part is key parameter in the recent research [5]. Increasing catalyst layer efficiency and durability play a important role in PEM fuel cell design and manufacturing [6,7]. In catalyst layer, the averaged pore size ranges from 60 nm to 188 nm which depend on the technologies used for visualisation [8,9], and fluid–wall collision could have a considerable impact on fluid

E-mail address: [email protected] http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.04.031 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.

properties. What characterises the relative dominance of fluid– fluid collision and fluid–wall collision is the Knudsen number [10], which is defined as the ratio between an averaged distance that a gas molecule travels between two consecutive collisions with other gas molecules and the characteristic size of the domain in which gas flows. As the Knudsen number increases, the impact of fluid–wall collision becomes increasingly important, and the dynamical property of the fluid is no longer able to be described by a single viscosity. Such flow is often called micro-flow, and in fuel cells it occurs in the catalyst layers. There has been an increased interest in micro-flow over the past decade, the lattice Boltzmann method (LBM) has also been applied to such flows [11–15]. Micro-flow cannot be described by the continuum approach such as the Navier–Stokes equation and some extra processes need to be considered [16]. The most important one is rarefaction, which is attributed to the original characteristics of fluid itself. The fluid–wall collision frequency of micro-flow is comparable with fluid–fluid collision frequency. As such, the fluid velocity at the wall surface can no longer be assumed to be non-slip, and as a result, it cannot be given a priori value. Instead, the slip velocity is part of the problem to be solved in micro-flow. Traditional modelling of micro fluid flow is based on molecular dynamics (MD). In MD modelling, the fluid is treated as a collection

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132

of particles, whose number should be sufficient to yield average-meaningful results to characterise the flow at a scale at which flow phenomena are measurable. The MD models are computational costly. Over the last few years, efforts have been made to develop the lattice Boltzmann model as an alternative to simulate fluid flow at nano scales. For example, Nie et al. [17] developed a LB model to simulate compressible flow in 2D micro-channels, finding that the LB model can capture fluid behaviours such as velocity slip, nonlinear pressure distribution along the channel and the dependence of mass flow rate on Knudsen number. In their model, the impact of Knudsen number was accounted for in the viscosity, which in turn was solved by allowing the relaxation-time to change with the Knudsen number. Review on micro-fluid modelling has been given by several researchers [16,18–22] and the lattice Boltzmann models for simulating micro-flows were also improved [17,23–25]. When gas flows in micro-pores, because of the gas slippage on the solid boundaries, the frictional dragging force is reduced. Such reduced dragging force depends on Knudsen number; the higher the Knudsen number, the less the dragging force is. In terms of the impact of the Knudsen number on gas flow in porous media, this means that the medium permeability is no longer the property of the medium, but also depends on the Knudsen number. This has been experimentally proven for H2, N2 and CO2 flow through a porous medium [25,26]. Shen [11] extended the work of Nie et al. [17], and compared their simulated velocity and pressure distributions in micro-channel with that obtained from the MD model. In their lattice Boltzmann (LB) model simulations, Shen used bounce-back method to treat the solid walls and the extrapolation scheme to treat the inlet and outlet boundaries of the channel. The predicted flow rate was found to be in good agreement with the results obtained from other methods with various Knudsen numbers (0.0194, 0.194 and 0.388). The results of Shen showed that the LB model of Nie is feasible to simulate gas flow only when the Knudsen number is small or moderate, and that it will give rise to considerable errors as the Knudsen number becomes large. Lee and Lin [27] proposed a second order definition of Knudsen number and a wall equilibrium boundary condition for simulating gas flow in micro-channels using the LB model; their tested examples show that the slip velocity simulated by the LB model was in close agreement with Arkilic’s prediction [28], indicating that their definition of the Knudsen number and the wall equilibrium boundary conditions are more meaningful than those used in previous studies [17,23,25]. The catalyst layer is an important component in PEM fuel cells as it is the place where all electrochemical reactions take place. However, since the catalyst layer is only a few hundreds of nano-metre thick, it is difficult to measure and quantify its ability to conduct gases. Therefore, in most macroscopic modelling, the catalyst layer was treated as an infinitely thin layer which is solved as prescribed flux boundary to the gas diffusion layer. The flux is described by the Butler–Volmer equation which relates the flux rate to over-potential and local concentration of the reactants [29]. This is an over-simplification as the local concentration in the catalyst layer depends not only on reactant concentration at the interface between the catalyst layer (CL) and the Gas diffusion layer (GDL), but also on how the reactants move into the catalyst layer. Over the past few years, there has been an increased interest to explicitly include the catalyst layer in macroscopic modelling. This, however, needs to know the effective permeability and effective diffusion coefficients of the catalyst layer to conduct gases, driven by both pressure gradient and molecular diffusion. Several stochastically models have been developed to investigate gas flow in catalyst layer, and Lange et al. [30] compared their impact on the effective transport properties. However, in these models, the gas

123

flow was assumed to be diffusion-dominated and the flow in cathode involves oxygen, hydrogen and water vapour. The movement of each species is described by the Stefan–Maxwell equations with a modified diffusion coefficient to account for the Knudsen diffusion due to gas–wall collision. In a recent work, Wu [31] used pore-network models to simulate oxygen flow in catalyst layer at the cathode by assuming that oxygen flow was driven by molecular and Knudsen diffusion only. The pore geometry in catalyst layer is heterogeneous with the size of its pores ranging in nano-metres. As a result, the Knudsen number is also heterogeneous with smallest pore having the highest Knudsen number. How such a heterogeneous Knudsen number affect the effective transport and reactive ability of the catalyst has attracted attention. For example, Siddique and Liu [32] investigated the impact of heterogeneous Knudsen number on the effective oxygen diffusion based on a numerically re-constructed catalyst layer, finding that the commonly used Bruggemann relation considerably underestimated the increase of effective diffusion coefficient with porosity. A typical catalyst layer consists of carbon to conduct electron and offer support to platinum nano-particles; ionometer binding to provide paths for proton conduction, and pores to provide path ways for gases to flow. To investigate how the carbon affects gas flow and proton conduction, Lange et al. [33] numerically generated a catalyst layer assuming that the carbon are uniform spheres. Their result showed that effective oxygen diffusion coefficient depends not only on porosity, but also on how the carbon spheres are packed. Chiavazzo and Asinari [34] illustrate a complex geometry reconstruction based on statistical analysis of images (Scanning electron microscope SEM) to be used in LB simulations and carried out on the thermal conductivity. In their group recent research, Salomov and Chiavazzo [35] construct a stochastic model for the catalyst layer (CL) based on clusterization of carbon particles, and investigate the effect of catalyst particle distribution inside the CL on the degree of clusterization and on the microscopic fluid flow. In general, the work on catalyst layers assumed that the bulk velocity of the mixture of the gases is zero, and that gas flow was dominated by molecular and Knudsen diffusion. The underlined assumption is that there is no friction between gas molecules and solid surface. This is a simplification. The purpose of this paper is using lattice Boltzmann method to investigate the impact of nano-pores on the ability of catalyst layer to conduct gases when gas flow is also driven by pressure gradient, in which the dragging force of the solid wall to the gas flow varies with the Knudsen number. And we proposed a new simplified catalyst structure to let the permeability calculation of catalyst is linked with the pore size distribution of the true catalyst layer.

2. Methodology 2.1. An approach to simplify the catalyst layer The pores in catalyst layer (CL) are multi-modal with the pore-size distribution centred around 50–100 nm. The size of the carbon grains in CL is in the order of 10 nm, which can agglomerate to 100 nm clusters. The catalysts deposit on the surface of carbon grains, and their size is in the order of 1 nm [36]. A typical PEFC catalyst layer, is consisting of carbon supported platinum and ionomer, was used in this study with a platinum loading of 0.45 mg=cm2 (0.21 lm3 =lm2 ) and a porosity of 40%. Fig. 1(a) shows a SEM image of a typical catalyst layer and Fig. 1(b) shows a catalyst layer cross section. Fig. 2 shows a 3D reconstruction of a mesoporous CL which acquired by a focused ion beam tomography at resolution of 15

124

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132

Fig. 1. (a) A SEM image of a cathode catalyst layer (b) SEM image of a catalyst layer cross section, the ellipse drawn shown in the image is Pt particle.

Fig. 2. FIB/SEM image of the catalyst layer at resolution of 15 nm.

nanometres, The dimensions of CL is 1  1  1 lm. This image has a 5  5  10 nm pixel size in x, y and z directions. In this image, the dark blue is for pores and light blue is for solid and red is for platinum particles larger than 5 nm. Such a high resolution imaging requires careful sample preparation. The sample preparation involves: deposition of protective Pt layer; rough milling or pocket milling; milling or recognising a fiducial mark; slice and view;

image alignment through fiducial mark and image binarisation with 3D reconstruction. Fig. 3 shows the pore-size distribution for the image shown in Fig. 2. This distribution is acquired from the analysis of the digital model generated using FIB nano-tomography process. It is evident that 69% of the pores has diameter in the range of 20–60 nm, and 25% in 60–140 nm. The reconstructed CL also reveals that the characteristic mean pore radius is 55 nm and the characteristic mean pore–pore lengths are 55 nm to 78 nm, respectively. Directly simulating gas flow through the nano-scale pores shown in Fig. 2 is difficult since the Knudsen number varies not only from pores to pores, but also with flow directions. To simplified the structure of the true catalyst layer in Fig. 2, we idealised the complicated pore geometry into a bundle of tubes. The 3D real catalyst layer in Fig. 4(a) is acquired by focused ions beam/scanning electron microscopy (FIB/SEM) technology, it shows the agglomerates which is platinum supported by carbon grains and these carbon grains are bound by an ionomer. In manufacturing, the carbon grains tend to aggregate, making the catalyst layer a bi-mode porous medium with micropores inside the agglomerates and macropores between [37]. Fig. 4(b) is the 3D binary tomography image which acquired form true catalyst layer, where the black is the solid (agglomerate, ionomer and carbon) and the white are pores inside. We choose x–z plane of Fig. 4(b) for example, you will found that have many pores in the plane, assume the pores are consist of circle and ellipse. It is simplify the micropores and macropores inside of the CL, macropores provide pathways for oxygen to diffuse from the gas diffusion layer to the catalyst layer, the micropores inside the agglomerates transport the oxygen to the catalyst where the electrochemical reaction takes place.

Fig. 3. Pore size distribution of a CL.

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132

125

Fig. 4. From 3D FIB/SEM image (a) to binary FIB image (b) 3D solid and pore image (c) 2D image in X–Z plane (d) tubes in the 3D model.

Based on these macropores and micropores, simplify the void space in catalyst layer into bundle shapes of tubes which shown in Fig. 4(d) in that gas flow. In order to make the simplified tubes reconstruction structure more close to the true catalyst layer, the diameter of these pores are all from the FIB/SEM nanotomography reconstruction which shown in Fig. 3, it is possible to elucidate certain transport mechanism and transport properties in true CL. Knudsen diffusion is a mode of transport where the likelihood of molecular–wall interactions exceed that of inter-molecular interactions due to the restrictive diameter of the pore wall and can occur in system where the mean pore diameter is of the order of 10 nm. This characteristic mean pore diameter of 55 nm confirms that the Knudsen regime is dominant in CL structures. For a given gas, each tube has a specified Knudsen number that are affected by pore diameter in Fig. 3. Simulating gas flow in tube with the specified Knudsen numbers could reveal how the heterogeneous Knudsen number affects the ability of the catalyst layer conduct gases. Subsequently, we can know how the pore diameter influence the route of fluid flow and used for the optimal structure design.

2.2. Lattice Boltzmann modelling The lattice Boltzmann method is able to handle extremely complex geometries and is a promising method for flow calculations on porous media, we have already successfully use it into the permeability calculation of gas diffusion layer [38]. But for flow in the catalyst layer, the wall–gas collision frequency is comparable to gas–gas collision frequency. As such, there is a slip on the solid wall, and the dynamic property of the gas depends not only on its viscosity, but also on the Knudsen number. In LB models, this is solved by accounting for the impact of Knudsen number into

the relaxation time parameter and modifying the boundary treatment. 2.2.1. The relationship between relaxation time and Knudsen number In kinematic theory, the mean free path of a gas molecule is proportional to its dynamic viscosity and inversely proportional to its pressure; that is given by Li [39]:

rffiffiffiffiffiffiffiffiffi

k¼

g pRT p

2

ð1Þ

where k is the mean free path, g ¼ ql is the dynamic viscosity, p ¼ qRT is gas pressure, R is the gas constant and T is temperature. For flow in a 3D tubes of catalyst layer is driven by pressure gradient. Because of the pressure change, the Knudsen number will be increase with the pressure decreases. If the pressure at the outlet of the tube is known, the value of the Knudsen number at the outlet can be determined by:

K n0 ¼

pffiffiffiffi k g p ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Hz Hz 2q2 RT

ð2Þ

where q is density, Hz is the diameter of the tube in the through-plane direction. When fluid flows over a solid surface, a Knudsen layer exists near the solid surface with its thickness comparable to the mean free path. If the Knudsen number is not negligible, the effective Knudsen number is defined as follows [12,14]:

K ne ¼ K n wðK n Þ

ð3Þ

where K ne is the effective Knudsen number. Assuming the Knudsen number at the outlet is K n0 , K n is the local Knudsen number. The local Knudsen number in the Z direction can be expressed as follows:

126

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132

Fig. 5. Lattice velocity direction of the LBM D3Q19.

K n ðzÞ ¼ K n0  Pout =PðzÞ

ð4Þ

where P out is the prescribed pressure at the outlet. K n0 is the export Knudsen number that calculate by Eq. (2), PðzÞ is the changing local pressure in through-plane (z) direction. The value of w in Eq. (3) can be described as follows [40]:

wðK n Þ ¼ 1 þ 0:5½ða  1Þea þ ðb  1Þeb  a2 Ei ðaÞ  b2 Ei ðbÞ

ð5Þ

R1

where a ¼ x=k, b ¼ ðHz  xÞ=k, and Ei ðxÞ ¼ 1 t 1 ext dt x ¼ 1; 2; . . . Hz is the exponential integral function, which considers the variation of the effective Knudsen number in both stream-wise and span-wise directions. When gas transport through the macro-pores (where the Knudsen number is low K n < 0:1), the relaxation parameter s can   calculated from kinematic viscosity t ¼ c2s s  12 dt . But if the gas transport through the nano-pores (big Knudsen number K n > 0:1), the relaxation time from the viscosity is not adapted in the LB model. Recent studies showed that most LB model for micro-flows are inadequate when K n > 0:1 [12] due to their inaccuracy to capture the Knudsen layer near the solid surface. To improve the LB model for simulating flow with high Knudsen number, this work use effective relaxation time as detailed below. For gaseous micro-flow, the effective relaxation parameter s in the lattice Boltzmann equation approximation is modified to consider the gas compressibility, and the effective Knudsen number is calculated from [41]

pffiffiffi pffiffiffiffi s ¼ 0:5 þ 6Ny  K ne = p

ð6Þ

where Ny ¼ L=dy is the lattice number in the direction of the characteristic length, L is the length of the tube, dy is the side length of the voxels, and K ne is the local effective Knudsen number. The relaxation time in Eq. (6) will be used in LB model when the gas transport through the nano-pores. 2.2.2. Evolution equation In the Multiple relaxation time (MRT) LBM model [42], the evolution of the distribution functions f a ðxi ; tÞ is described below which has been proved to be able to model gas flows with Knudsen numbers in Eq. (3):

  eq f a ðxi þ ca Dt; t þ DtÞ  f a ðxi ; tÞ ¼ S f a ðxi ; tÞ  f a ðxi ; tÞ

ð7Þ

where Dt is a time step during which the particle moves from one lattice into another, S ¼ M 1 ^SM is a matrix in which M is the transform matrix that transforms the distribution functions into moments, and ^S is a diagonal matrix to perform the relaxation collision in moment space. For the D3Q19 LB model used in this paper, Fig. 5 shows the lattice structure of the D3Q19 LB model, in which the physical space is divided into regular lattice and the velocity space at which the particles move is discretised into a finite set of lattice velocities. ^ are as follows: The terms of S

^S ¼ diagðs0 ;s1 ;s2 ;s3 ;s4 ;s5 ;s6 ;s7 ;s8 ;s9 ;s10 ;s11 ;s12 ;s13 ;s14 ;s15 ;s16 ;s17 ;s18 Þ ð8Þ in which

s0 ¼ s3 ¼ s5 ¼ s7 ¼ 0 s1 ¼ s2 ¼ s9 ¼ s10 ¼ s11 ¼ s12 ¼ s13 ¼ s14 ¼ s15 ¼ 1=s s4 ¼ s6 ¼ s8 ¼ s16 ¼ s17 ¼ s18 ¼ 8ðð2  1=sÞ=ð8  1=sÞÞ where s is a relaxation time which acquired from Eq. (6). Transform the distribution functions on the right-hand side of Eq. (7) into moments format gives:

f a ðxi þ ca Dt; t þ DtÞ ¼ f a ðxi ; tÞ  M1 ^Sðma ðxi ; tÞ  meq a ðxi ; tÞÞ

ð9Þ

where Dt is a time step during which the particle moves from one eq lattice into another, ma ¼ M  f a is the moment, and meq a ¼ M  fa is the value of ma when the system is in equilibrium. M is a matrix to transform the particle distribution functions into moments, and its details is given in [43,44]. For the D3Q19 model used in this paper, the lattice velocity ca is given as follows:

8 a¼0 > < ð0;0; 0Þc=Dt ca ¼ ð1; 0; 0Þc=Dt; ð0;1; 0Þc=Dt; ð0;0;1Þc=Dt a ¼ 1;. .. ;6 > : ð1; 1;0Þc=Dt; ð1; 0; 1Þc=Dt; ð0; 1; 1Þc=Dt a ¼ 7;. .. ;18 ð10Þ where c is the side length of each cubic lattice. The equilibrium distribution functions for the D3Q19 mode are given by:

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132



 3 u2 ¼ xa q 1  a¼0 2 c2s " # ca  u 9 ðca  uÞ2 3 u2 a ¼ 1; . . . ; 18 ¼ xa q 1 þ 3 2 þ  cs 2 2 c2s c4s

eq

fa

eq

fa

The combination of the bounce-back and specular-reflection boundary treatment is.

ð11Þ

where x0 ¼ 1=3; x1;...;6 ¼ 1=18; x7;...;18 ¼ 1=36, c2s ¼ c2 =3Dt 2 and the macroscopic density q and bulk fluid velocity u in the above equations are related to the distribution functions through 19 X fa

q¼

a¼1

qu ¼

ð12Þ

19 X f a ca

Implementation of the above model involves a collision step and a streaming step to advance one time step. Since the size of the pores in catalyst layers is in nano scale, it is very difficult to measure the ability of the catalyst layer to conduct gases. The model presented in this work provides a way to qualitatively investigate the average flow phenomena in every nano-pores. An important application of the digital three dimensional models which is able to closely represent the structural features of a sample is for determining structure properties of the CL. The digital model reconstructed has been used to find that the porosity of the CL is 40%. The analysis of the digital model generate from the FIB nano-tomography shows that the pore size distribution peaks at 20 nm, 40 nm and 60 nm, the reconstructed digital model of the CL also reveals that the characteristic mean pore radii is 55 nm and the characteristic mean pore–pore length is 78 nm. The average volume of pore in true catalyst layer is calculate by



i¼1

Di 2

p

 wi Va

2

 78

¼ 40%

ð13Þ

where Di is the diameter of the ith tube which is from Fig. 3, wi is the percentage of the pore volume. V a is the average cross-section area of pores which will be used into the LB modelling. We calculate the fluid flow in z direction which is through-plane direction and the average velocity in this direction is qn which calculated from

qn ¼

c3

P

i ui;zn

ð14Þ

Va

in which Lx ; Ly ; Lz are the lengths of the simulated domain in the x, y and z directions respectively, ui;zn are the bulk fluid rate components through the z direction in ith voxel of nth simplified tube. All the simulations were run using a Dell Precision T5400 workstation with two Intel CPU E5420. In the simulations, a steady flow was adjudged using the parameter X where

P

X¼

f a ðx; tÞ ¼ rfa ðx; tÞ þ ð1  rÞf a0 ðx; tÞ

ð16Þ

where 0 6 r 6 1 is a coefficient that weights the portion of the contribution of the no-slip bounce-back and specular reflection, respectively. a is the unknown point, a  is the point on the diagonal and a0 is the point on the horizontal. As an illustration example, we take the X–Y plane as an example. The particle distribution functions are computed as follows after the bounce-back treatment:

f 8 ðx; tÞ ¼ rf10 ðx; tÞ þ ð1  rÞf 7 ðx; tÞ f 7 ðx; tÞ ¼ rf9 ðx; tÞ þ ð1  rÞf 8 ðx; tÞ

a¼1

Pn

127

i juj ðxi tÞ

 uj ðxi t  1Þj < 106 i juj ðxi tÞj

P

ðj ¼ x; y; zÞ

ð15Þ

f 1 ðx; tÞ ¼ f 2 ðx; tÞ f 10 ðx; tÞ ¼ rf8 ðx; tÞ þ ð1  rÞf 9 ðx; tÞ

ð17Þ

f 9 ðx; tÞ ¼ rf7 ðx; tÞ þ ð1  rÞf 10 ðx; tÞ f 4 ðx; tÞ ¼ f 3 ðx; tÞ The value of r in the above equation can be chosen as follows [45]. For the MRT model

2c2 r ¼ ppffiffiffi þ c1 6

ð18Þ

In the MRT-LBE, it is only depended on the gas–solid interaction parameters c1 ¼ 1:11 and c2 ¼ 0:61, these parameter have already validate in [45]. We assume that the gas flow is driven by a pressure drop between the inlet and outlet. The pressure at the inlet and outlet are Pin and P out respectively; P r ¼ Pin =P out is pressure ratio. When the Knudsen number can be neglected, the pressure is linearly distributed along the tube:

z Pl ðzÞ ¼ Pin þ ðPout  P in Þ L

ð19Þ

With the impact of Knudsen number, the pressure distribution is P(z) which means that the pressure distribution in Z directions in Fig. 5, and the difference between P(z) and Pl(z) is defined as pressure deviation, which, after normalisation, is,

P ¼ ½PðzÞ  Pl ðzÞ=P out

ð20Þ

The analytical solution for the normalised pressure deviation along the channel is given as follows [46].

P ðzÞ ¼ 8K n0 þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð8K n0 Þ2 þ ð1 þ 16K n0 Þz þ ðP2r þ 16K n0 Pr Þð1  zÞ

P ¼ ½P ðzÞ  Pl ðzÞ=Pout ð21Þ where z ¼ z=L. 3. Numerical results 3.1. Model validation

Flow is assumed to have reached steady state when the tolerance X < 106 is satisfied. 2.2.3. Fluid–solid boundary after modifying For the gas transport in the nano-pores, the gas velocity at the solid wall is not zero but depends on the Knudsen number. It thus needs a special treatment. For the wall boundary, the bounce-back method is usually used with a specular reflection to reflect that the solid wall is no longer a non-slip boundary [23]. For the D3Q19 MRT LB model used here, the lattice velocity direction of the LBM D3Q19 is shown in Fig. 5. The whole domain is divided into 200  200  200 in the x, y and z directions.

Gas flow in tubes with different diameter is shown in Fig. 3. The simulations are for gas flow in both the cathode and anode. As such, we consider hydrogen, oxygen and water vapour. The associated Knudsen number for H2, O2, and H2O (vapour) at the outlet is 0.1, 0.054472 and 0.062413, respectively, and they are calculated by Eq. (2). Fig. 6(a)–(c) compare the analytical normalised pressure deviation distributions across the tube in Fig. 4 with the pressure deviation distribution for the three gases under different pressure drops P r . The prediction results are fit well with Arkilic’s analytical solutions [28]. From Fig. 6, we found some points: (1) for hydrogen,

128

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132

Fig. 6. Nonlinear pressure profiles for: (a) hydrogen with K n0 = 0.10 and pressure ratio from 1.2 to 2.0; (b) oxygen with K n0 = 0.055 and pressure ratio from 1.2 to 2.0; (c) water vapour K n0 = 0.062413, pressure ratio from 1.2 to 2.0 and (d) pressure ratio at 2.0 with different Knudsen numbers.

oxygen and water vapour, the pressure gradient increase, the pressure deviation increase also. Like, for hydrogen, when pressure gradient increase from 1.2 to 2.0, the deviation is from below 0.01 to over 0.07. For other two gases is similar. (2) For each gas, the small pressure gradient will let the prediction results more close to the theoretical one. Observe the hydrogen one, when the pressure gradient equal to 1.2, the prediction results is nearly cover the analytical one, but have a large gap at the 2.0 gradient. The condition is same for oxygen and vapour. (3) with the increase of Knudsen number (from oxygen to hydrogen), the prediction results is more close to the theoretical value. Compare Fig. 4(a) with Fig. 4(b), at the same pressure gradient equal to 2.0, the deviation of hydrogen is more close to analytical resolution compare to oxygen. Fig. 4(d) is take out the prediction results of three gas at the 2.0 pressure gradient to shown, the curve clear show that the large pressure deviation is happened for the small Knudsen number. The pressure of oxygen is close to 0.09 and hydrogen is just over 0.07. Fig. 7 shows the change of local Knudsen number K n along the tube under different pressure ratios for the three gases. Some points conclude: (1) Through the flow direction, from inlet to outlet, the Knudsen number is creased. Like, when gas flow in the z directions (through-plane), the Knudsen number is increase from 0.026 to 0.055 for oxygen from 0.031 to 0.0624 for vapour from 0.05 to 0.1 for hydrogen. (2) Compare three gases, the hydrogen have the highest number of local Knudsen number, and have the biggest increase through the flow direction. (3) For the high pressure gradient, it have a slower local Knudsen number increase. For hydrogen, the local Knudsen number at 1.2 pressure drop have more steeper curve compare to 2.0 pressure drop.

3.2. Relationship between flow rate and pressure gradient P Under a given pressure gradient, i ui;zn is the flow rate for nth tubes (these tubes have different diameter which is shown in Fig. 3) when fluid transport in z direction, which reflects the resistance of the material to gas flow. In the pore size distribution direction in Fig. 8, we found that the low flow rate happened at the pore size between 10 nm and 70 nm and because here have high resistance than other size. For the small pore part where have the high Knudsen number, the complex structure will play an important part role, and also the effects of rarefaction are play an important role with the increase of Knudsen number. The flow rate increase quickly when pore size from 110 to 130. The pore size is nearly two times than previous, but the flow rate is nearly four times than 10 to 70 parts. These help us to know more about the internal transfer characteristic for each different tube here. Through analysis these structure and transfer character, it is help to improve the pores size distribution in the catalyst layer and enhance the flow rate, let more reaction gas enter into the pores of agglomerate to improve the catalyst utilisation. The variation trend is nearly same for the two different gases (oxygen and hydrogen), but we found that the oxygen have the large flow rate for the same pore size compare with hydrogen. Due to our simplified model regards the CL is consisted of different tubes with various diameters, and these diameters were measured from a true CL sample. For different cylindrical pore structures with various diameters, the fluid velocity shows a parabolic profile which like the Poiseuille flow. The laminar poiseuille flow, it has a exact analytic solution for the cylindrical pipes as follow [47]:

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132

129

Fig. 7. Change of local Knudsen number K n of the three gases along the tube under different P r .

Fig. 8. The calculated gas flow rate under different pressure drops in the nano-tubes of catalyst layer (a) hydrogen, and (b) oxygen.

Q¼

pR4 @p 8g @x

ð22Þ

where g is the dynamic viscosity, R is the radius of the tube and @p is @x the pressure gradient, the permeability K is computed by applying Darcy’s law qx ¼  Kg

dp ; dx

the relationship between permeability and

pore diameter is given by:

K¼

pR4 8A

ð23Þ

where A is the cross-sectional area, Fig. 9 shows calculated permeability from our simplified LB model and predicted permeability from Eq. (23). We calculated the permeability with different average diameter from Fig. 3, Good agreement with the theoretical

130

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132

Fig. 9. Simulated permeability by the lattice Boltzmann method (red dots) and the theoretical predictions (blue solid dot line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

prediction with our simulated values, the small discrepancy at small diameter comes from the discretization effect. To validate the permeability calculation with other prediction results, we use a cubic computational box made of 200  200  200 nodes and porosity is close to 36%. The PEFC-CL is in reality a heterogeneous structure with a range of nano-scale features, the carbon grains which are typically of 10 nm orders can agglomerate to 100 nm order, and catalytic deposits on the carbon grains that are typically of 1 nm order. When pressure gradient used is equal to 2:0  106 l:u:, the permeability is calculated and compare with Uktam’s results. The results are in the same order

Table 1 Compare the sensitivity of the permeability to the cluster size in the CL. LBM results

U.R. Salomov results Ref. [35]

Resolution Average size of cluster Lattice pressure gradient

200  200  200 10–100 nm

200  200  200 40 nm

2:0  106 l:u:

2:51256  106 l:u:

Lattice length

1:0  108 m 0.07441 mD

1:0  108 m 0.107451 mD

Permeability

stage, and have small difference that might because of the structure difference or the pressure gradient small difference (see Table 1). Once the flow rates in all the tube are calculated, the permeability is calculated from the pressure and flow data [48]

K¼

2lqn pout p2in  p2out

ð24Þ

where l is the kinematic viscosity of the gas, pin pout are inlet pressure and outlet pressure, and qn is the averaged flow rate which is acquired from Eq. (14). The gas permeability which depended on the simplified tube model of catalyst layer is calculated by Eq. (24), theunit of the permeability is milidarcy (MD). In Fig. 10, As the Knudsen number increases, the rarefaction effect strengthen will increase, and hence the gas permeability increases. And also, at relative large pressure gradient point, compare the hydrogen with the oxygen, the increase of flow rate with pressure drop is nonlinear, and also non-linearity increases with the Knudsen number. The gas permeability profiles show an obvious nonlinear behaviour due to the slip velocity and nonlinear drag-force influence.

Fig. 10. Gas permeability comparisons with different pressure gradient.

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132

The unit of Permeability is Darcy, one Darcy is about 9:869  1013 m2 , and the calculated permeability can be converted from lattice unit (l.u) to Darcys using the following expression

dx2 K ½Darcy ¼ K ½l:u 0:9869

131

Johnson Matthey Fuel Cells Ltd., Saati Group Inc., Technical Fiber Products Ltd and academic partners Liverpool University, Loughborough University and Brimingham University. References

ð25Þ

where dx is the lattice spacing in microns. 4. Conclusions and discussion In this paper, we have presented a simplified model to simulate gas flow in nano-pores of catalyst layer, in an attempt to understand the flow characteristic when gas–wall collision frequency is comparable to gas–gas collision frequency for gas flow in the FEM catalyst layers. The complicated pore geometry in the catalyst layer is simplified into a bundle of tubes and the number and the diameters of the tubes are derived from the porosity and pore-size distribution of true catalyst layer. The lattice Boltzmann method is then presented to simulate gas flow in each tube. Summing up the flow rates in all the tubes could reveal how the nano-pores affect the gas flow behaviour in the catalyst layer. The impact of the Knudsen number on the dynamic properties of the gases is adsorbed in the relaxation time parameters, and the non-slip boundary was solved by a modified bounce–back method to yield a slip velocity at the gas–wall interface. The simulated results were validated against analytical solutions, and the comparison shows good agreement. The model was then applied to simulate gas flow in the simplified catalyst layer. The results reveal that, when the Knudsen number becomes significant, the flow rate and permeability through the catalyst layer is no longer linearly increasing with pressure gradient. With an increase in Knudsen number, the increase of flow rate with pressure gradient becomes increasing nonlinear. The flow rate change with the pore size distribution can reflect the resistance of the material which helpful to the structure improving. This new simplified tube model explore the internal characteristics of catalyst layer and shed some insight into gas flow in the catalyst layers, which is difficult to implementation by experiment or recent modelling simulation due to the complex structure and nano sized pore size. Although the catalyst layer is idealised, the results still has some implications. For example, it revealed that the permeability of the catalyst layer for hydrogen differs from its permeability for oxygen and water vapour due to the effect of gas-wall collision. Therefore, in macroscopic modelling the permeability of the catalyst layer for different gases and water should take different values. For the different pressure drop in the flow direction, the growth gradient of the flow rate has much different with the pore size distribution. For hydrogen and oxygen, the oxygen is more affected by the pore size; the flow rate will change more with pore size for the same pressure drop, especially for the low pressure drop. Conflict of interest None declared. Acknowledgements This research is partly supported by the UK Technology Strategy Board (TSB Project no.: TP/6/S/K3032 H) and partly supported by the Fundamental Research Funds for the Central Universities (Grant No. 1700219139). We acknowledge the kind support from industrial partners AVL List GmbH, Intelligent Energy Ltd.,

[1] M.F. Mathias et al., in: W. Vielstich, A. Lamm, H.A. Gasteiger (Eds.), Handbook of Fuel Cells-Fundamentals, Technology and Application, vol. 3, John Wiley and Sons, Chichester, England, 2003. [2] C.J. Zhong et al., Nanostructured catalysts in fuel cells, Nanotechnology 21 (2010) 1–20. [3] F. Barbir, PEM Fuel Cells: Theory and Practice, Elsevier Academic Press, 2005. [4] S.S. Zhang et al., A review of platinum-based catalyst layer degradation in proton exchange membrane fuel cells, J. Power Sources 194 (2009) 588–600. [5] G.Q. Wang, P.P. Mukherjee, C.Y. Wang, Optimization of polymer electrolyte fuel cell cathode catalyst layers via direct numerical simulation modeling, Electrochim. Acta 52 (2007) 6367–6377. [6] X. Yu, J.L. Yuan, B. Sunden, Review on the properties of nano-/microstructure in the catalyst layer of PEMFC, in: J. Fuel Cell Sci. Technol. 8 (2011). [7] R. Balgis et al., Nanostructured design of electrocatalyst support materials for high-performance PEM fuel cell application, J. Power Sources 203 (2012) 26– 33. [8] S. Thiele, R. Zengerle, C. Ziegler, Nano-morphology of a polymer electrolyte fuel cell catalyst layer-imaging, reconstruction and analysis, Nano Res. 4 (2011) 849–860. [9] W.K. Epting, J. Gelb, S. Lister, Resolving the three-dimensional microstructure of polymer electrolyte fuel cell electrodes using nanometer-scale X-ray computed tomography, Adv. Funct. Mater. 22 (2012) 555–560. [10] G. Arya, H.C. Chang, E.J. Maginn, Molecular simulations of Knudsen wall-slip: effect of wall morphology, Mol. Simul. 29 (10–11) (2003) 697–709. [11] C. Shen et al., Examination of the LBM in simulation of microchannel flow in transitional regime, Microscale Thermophys. Eng. 8 (2004) 423–432. [12] Z. Guo, T.S. Zhao, Y. Shi, Physical symmetry, spatial accuracy and relaxation time of the lattice Boltzmann equation for micro gas flows, J. Appl. Phys. 99 (7) (2006) 074903. [13] F. Verhaeghe, L.-S. Luo, B. Blanpain, Lattice Boltzmann modeling of microchannel flow in slip flow regime, J. Comput. Phys. 228 (2009) 147–157. [14] Z. Guo et al., Discrete effects on boundary conditions for the lattice Boltzmann equation in simulating microscale gas flows, Phys. Rev. E 76 (2007) 074903. [15] X.L. Liu, Z. Guo, A lattice Boltzmann study of gas flows in a long micro-channel, Comput. Math. Appl. 65 (2013) 186–193. [16] G. Karniadakis, A. Beskok, N. Aluru, Microflows and Nanoflows: Fundamentals and Simulations, Springer, New York, 2005. [17] X. Nie, G.D. Doolen, S. Chen, Lattice Boltzmann simulations of fluid flow in MEMs, J. Stat. Phys. 107 (2002) 279–289. [18] H. Bruus, Theoretical Microfluidics, Oxford University Press, Oxford, 2008. [19] N. Kockmann, Transport Phenomena in Micro Process Engineering, Springer, Berlin, 2008. [20] T.R. Dietrich, Microchemical Engineering in Practice, Wiley, Oxford, 2009. [21] G.M. Whitesides, The origins and the future of microfluidics, Nature 442 (7101) (2006) 368–373. [22] J. Zhang, Lattice Boltzmann method for microfluidics: models and applications, Nanofluid 10 (2011) 1–28. [23] C. Lim, C. Shu, X.D. Niu, Application of lattice Boltzmann method to simulate microchannel flows, Phys. Fluids 14 (7) (2002) 2299–2308. [24] G. Tang, W. Tao, Y. He, Gas flow study in MEMs using lattice Boltzmann method, in: 1st International Conference on Microchannels and Minichannels, ASME, New York, 2004. [25] G. Tang, W. Tao, Y. He, Three dimensional lattice Boltzmann model for gaseous flow in rectangular microducts and microscale porous media, J. Appl. Phys. 97 (10) (2005). [26] J.L. Ge, Porous Flow Mechanics in Gas and Oil Reservoir, Petroleum Industry Press, Beijing, 1982. [27] T. Lee, C.L. Lin, Rarefaction and compressibility effects of the lattice Boltzmann equation method in a gas microchannel, Phys. Rev. E. 71 (4) (2005) 046706. [28] E.B. Arkilic, M.A. Schmidt, Gaseous slip flow in long microchannels, J. Microelectromech. Syst. 6 (2) (1997) 167–178. [29] T. Berning, D.M. Lu, N. Djilali, Three dimensional computational analysis of transport phenomena in a PEM fuel cell, J. Power Sources 106 (2002) 284–294. [30] K.J. Lange, P.C. Sui, N. Djilali, Determination of effective transport properties in a PEMFC catalyst layer using different reconstruction algorithms, J. Power Sources 208 (2012) 354–365. [31] R. Wu et al., Pore network modeling of cathode catalyst layer of proton exchange membrane fuel cell, Int. J. Hydrogen Energy 37 (2012) 11255–11267. [32] N.A. Siddique, F.Q. Liu, Process base reconstruction an simulation of a three dimensional fuel cell catalyst layer, Electrochim. Acta 55 (19) (2010) 5357– 5366. [33] K.J. Lange, P.C. Sui, N. Djilali, Pore scale modeling of a proton exchange membrane fuel cell catalyst layer: effect of water vapor and temperature, J. Power Sources 196 (6) (2011) 3195–3203. [34] E. Chiavazzo, P. Asinari, Reconstruction and modeling of 3D percolation networks of carbon fillers in a polymer matrix, Int. J. Therm. Sci. 49 (2010) 2272–2281.

132

Y. Gao / International Journal of Heat and Mass Transfer 88 (2015) 122–132

[35] U.R. Salomov, E. Chiavazzo, P. Asinari, Pore scale modeling of fluid flow through gas diffusion and catalyst layers for high temperature proton exchange membrane (HT-PEM) Fuel cells, Comput. Math. Appl. 67 (2014) 393–411. [36] M. Uchida et al., Investigation of the microstructure in the catalyst layer and effects of both perfluorosulfonate ionomer and PTFE-loaded carbon on the catalyst layer of polymer electrolyte fuel cells, J. Electrochem. Soc. 142 (1995) 4143–4149. [37] N.A. Siddique, F.Q. Liu, Process based reconstruction and simulation of a threedimensional fuel cell catalyst layer, Electrochim. Acta 55 (2010) 5357–5366. [38] Y. Gao et al., Modeling fluid flowing the gas diffusion layers in pemfc using the multiple relaxation-time lattice Boltzmann method, Fuel Cells 12 (3) (2012) 365–381. [39] Q. Li et al., Lattice Boltzmann modeling of microchannel flows in the transition flow regime, Microfluid. Nanofluid. 10 (2011) 607–618. [40] Z. Guo, C.G. Zheng, B. Shi, Lattice Boltzmann equation with multiple effective relaxation times for gaseous microscale flow, Phys. Rev. E 77 (2008) 036707. [41] X.L. Liu, Z. Guo, A lattice Boltzmann study of gas flows in a long micro-channel, Comput. Math. Appl. 65 (2) (2013) 186–193.

[42] D. d’Humieres et al., Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. London Ser. A – Math. Phys. Eng. Sci. 360 (1792) (2002) 437–451. [43] P.K. Sinha, P. Halleck, C.Y. Wang, Quantification of liquid water saturation in a PEM fuel cell diffusion medium using X-ray microtomography, Electrochem. Solid State 9 (7) (2006) A344–A348. [44] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, 2001. [45] Z.H. Chai et al., Gas glow through square arrays of circular cylinders with Klinkenberg effect: a lattice Boltzmann study, Commun. Comput. Phys. 8 (2010) 1052–1073. [46] C.I. Weng, W.L. li, C.C. Hwang, Gaseous flow in microtubes at arbitrary Knudsen numbers, Nanotechnology 10 (1999) 373–379. [47] Y. Keehm, Computational rock physics: transport properties in porous media and applications, in: Department of Geophysics, Stanford University, 2003. [48] A.E. Scheidrgger, The Physics of Flow Through Porous Media, University of Toronto Press, Toronto, 1972.