Lattice Boltzmann Modeling of the Gas Diffusion Layer of the Polymer Electrolyte Fuel Cell with the Aid of Air Permeability Measurements

Polymer electrolyte fuel cells (PEFCs, PEMFCs) with high efficiency and low environmental impact recently have attracted considerable interest. However, further improvement in performance and reliability is required to realize practical use of PEFCs as future power generation devices. To improve PEFC performance, an appropriate water balance between the water content and product water is a key technology. Loss of water content in the polymer electrolyte membrane decreases proton conductivity, thereby increasing the internal resistance of the cell. A PEFC basically consists of a membrane electrode assembly (MEA), gas diffusion layers (GDLs) and separators having flow fields with flow channels and ribs. The design parameters for the GDL, such as thickness, pore size distribution, and gas permeability play important roles in characterizing the gas flow and water management during PEFC operation[1]. In this chapter, 2D anisotropic modeling of a monolayer of the GDL substrate is carried out by comparing calculated and measured gas permeability with the lattice Boltzmann method (LBM)[2, 3] and through-plane/in-plane gas permeability measurements, respectively.


Introduction
Polymer electrolyte fuel cells (PEFCs, PEMFCs) with high efficiency and low environmental impact recently have attracted considerable interest. However, further improvement in performance and reliability is required to realize practical use of PEFCs as future power generation devices. To improve PEFC performance, an appropriate water balance between the water content and product water is a key technology. Loss of water content in the polymer electrolyte membrane decreases proton conductivity, thereby increasing the internal resistance of the cell. A PEFC basically consists of a membrane electrode assembly (MEA), gas diffusion layers (GDLs) and separators having flow fields with flow channels and ribs. The design parameters for the GDL, such as thickness, pore size distribution, and gas permeability play important roles in characterizing the gas flow and water management during PEFC operation [1]. In this chapter, 2D anisotropic modeling of a monolayer of the GDL substrate is carried out by comparing calculated and measured gas permeability with the lattice Boltzmann method (LBM) [2,3] and through-plane/in-plane gas permeability measurements, respectively.

Lattice Boltzmann method
LBM is a numerical fluid dynamic simulation method that describes macroscopic fluid dynamic phenomena by analyzing the behavior of virtual particles of which fluid is regarded as aggregate. LBM gives simplified kinetic models that incorporate the essential physics of microscopic processes so that the macroscopic averaged properties obey the macroscopic Navier-Stokes equations. Because the conventional Navier-Stokes equation takes long time to calculate and results in poor convergence in porous media, LBM has been developed to take advantage the simplicity of the algorithm and flexibity for complex geometries such as porous media [4]. It is therefore reasonable to apply LBM to fluid flows in the porous structure of the GDLs [5].

Governing equations
LBM analyzes flow by solving the lattice Boltzmann equation (LBE) that describes particle distribution function, which represents flow velocities of virtual particles. Macroscopic parameters such as the flow velocity and pressure are derived from the summation of the moment of the particle velocity.
In general, LBM uses a single relaxation time approximation by the Bhatnager, Gross, Krook (BGK) model [3,6]. Equation 1 presents the Boltzmann equation with the BGK approximation where f represents the distribution function depending on space, x, velocity, v and time, t. f eq is the local equilibrium distribution function, and τ is the relaxation time to local equilibrium. The discrete Boltzmann equation is thus since v-space is discretized by a finite set of particle velocities, v α and associated distribution function f α (x, t).
Discretizing with δt and x + e α δt, Eq. 2 gives The collision-streaming process of the LBM is calculated with the following equations.
Here,f i is the distribution function after the streaming process. The collision process represents the process that the distribution function converges to the equilibrium state, while the streaming process is the process that the virtual particles move to the neighboring sites.
2-dimensional fluid calculation uses the 2D9V model. The virtual particle velocity vector is where C = δx/δt is velocity of a virtual particle.
The local distribution function for the 2D9V model is where ρ is the density per node, u = u x u y T is the fluid velocity, ̟ α is the weighting function expressed as follows.
The distribution function for each direction is The density and macroscopic flow velocity at a node are defined by respectively.
Equations 6 and 19 give the velocity as follows.
In 2D9V model, the sound velocity C s and the pressure p are Kinetic viscosity is expressed as

Boundary conditions
LBM defines the velocity distribution function at boundary from the velocity and pressure to use the boundary condition. In general, the following bounce back boundary condition has been employed.

Half-way wall bounce-back boundary condition
Half-way wall bounce back boundary condition is no-slip boundary condition at a given solid surface as follows [2,3].
where f " α represents the distribution function after the streaming step. Since this boundary condition gives higher precision than the conventional bounce back boundary condition [3], it is employed in the present chapter.

Periodic boundary condition
For large area calculation, periodicity of the solution can be assumed. In this case, the periodic boundary condition is employed along the axis direction. The distribution function is where Nx is the maximum value of x.

Pressure difference boundary condition
The pressure difference boundary condition is applied to a case that there is pressure difference between inlet and outlet while the velocity distribution is the same at the inlet and outlet. The distribution functions at the inlet are assumed as follows [7].
while the distribution functions at the outlet are

Pressure and velocity boundary conditions
Pressure and velocity boundary conditions proposed by Zou and He [8] are also used for the calculation. The case for an inlet at x = 0 is considered for instance here. At the boundary, pressure, that is, ρ in , and u y = 0 are applied. Because f " 2 , f " 3 , f " 4 , f " 6 , f " 7 , f " 9 after streaming step are known, u x , f " 1 , f " 5 , f " 8 are derived as follows. Equation 18 gives while Eqs. 20 and 21 lead to Hence Then the deviation from the equilibrium shall be equal for the distribution function of i = 1,3 as follows to determine the remaining distribution functions.
Thus Eqs. 9 and 11 yield Eqs. 42, 43, and 46 give Thereby u x , f 1 , f 5 , f 8 are determined. On the other hand, the distribution function at the corner should be dealt with in other way. The case for the bottom of the inlet at x = 0, y = 0 is described here for instance. After the streaming step, f " 3 , f " 4 , f " 7 , ρ in is obtained. The no-slip boundary condition gives u x = 0, u y = 0. Thus f 1 , f 2 , f 5 , f 6 , f 8 can be determined as follows.
Eq. 46 provides In a similar manner, Equations 18 and 47 yield From Eq. 44,

LBM binary mixtures with different molecular weights
In this section, LBM binary mixtures with different molecular weights (LBM-BMD) model proposed by Luo and Girimaji [9] and extended by McCracken and Abraham [10] is described. LMB-BMD model consists of LBM and the effect of diffusion. This model deals with two components A and B having different molecular weights. The model can analyze advective flow in addition to diffusion. i and j represent the functions, variables, and constants for the species A and B, respectively.
The equilibrium distribution function is where Ω ii α and Ω ij α are the following self-collision and cross-collision terms for A-A and A-B, respectively. Velocity vector is defined as the similar manner as LBM in the previous section.
where τ i and τ ij D are the relaxation times for the kinetic viscosity, ν i , and diffusion coefficient, In a similar manner as the single component LBM in the previous section, where u i diff,x and u i diff,y are x and y direction of u i diff = u i − u Collision terms are: where The density and flow velocity are derived by The total density and mass averaged velocity are Since the partial pressure is total pressure is The relation between the sound velocities, C i s and C j s is where m i > m j are the molecular weights of A and B, respectively. The kinetic viscosity and diffusion coefficient have the following relations.

Streaming step of species with different velocities
In LBM-BMD model, species A snd B have different velocities. During δt , the species A travel δx, while the species B travel m i /m j δx. Since δx = e i α δt, the species B have different streaming distance. So, the distribution function of B on the nodes should be determined from the interpolation of the distribution functions of surrounding particles. Although McCracken and Abraham proposed a second-order Lagrangian interpolation [10], and Joshi et. al proposed bi-linear interpolation [11], these interpolation methods seem not appropriate for porous structure despite their higher accuracy. Thus linear interpolation of fewer nodes is employed here. The case for α = 1 is depicted in Fig. 4. After the streaming step, the distribution function at (x, y) is determined by the interpolation of the distribution functions at (x − δx, y) and (x, y) before the streaming step.
when (x − δx, y) is an obstacle node,f

GDL models
In the present research, 3D structure of the GDL is projected to 2D structure. So, the following models are created.
Model 1 Cross-section of the carbon fiber is simulated as a circle so that averaged number of the fiber in unit area is the same as that of the actual GDL.  Model 2 Fiber is simulated so that porosity is the same as the actual GDL.
Model 3 Fiber and its cross-section are simulated by so that porosity is the same as the actual GDL. Figure 5 illustrates the GDL models. Through-plane and in-plane anisotropic Darcy coefficients are obtained by LBM employed to these GDL models. These Darcy coefficients are compared with those from the following permeability measurements so that an anisotropic GDL model which agrees the most with the measurements can be found. The GDL model found is then used for the LBM-BMD flow analysis in the GDL with the flow channels and ribs having actual PEFC flow field geomeory under actual operation condition.
The calculation in the present chapter was carried out with a personal computer having Intel Core2 Quad CPU Q6600 2.4GHz and 4GB memory on ASUS P5K motherboard. MatLab (MathWorks, Inc.) was used for the LBM and LBM-BMD calculations. Figure 6 shows a schematic diagram of GDL permeability measurement apparatus. GDL, which was a commercial carbon paper (SIGRACET GDL 24AA, SGL Carbon Inc.) with a thickness of 190 µm, was placed between two cylindrical plates. A soft O-ring was used for gas sealing between the plates. The force required to deform the O-ring was negligible compared with the compression force acting on the GDL. The compression force was controlled using a clamp screw and was measured with a load cell. For air permeability tests, the compression pressure was set at 1 MPa, as measured in a typical PEFC. Fig. 7 presents geometries of the GDL used for the through-plane and in-plane permeability tests [12]. Volumetric air flow rates in through-plane direction, Qth, and in-plane direction, Qin, in the following equations were measured using a mass flow meter (KOFLOC). Pressure drop by the apparatus was compensated beforehand.

Experimental
where k, µ, and r are the Darcy coefficient of the GDL, viscosity of air [13], and radius of the GDL, respectively. P i , P O , δ, and A are inlet and outlet air pressures, thickness of the GDL, and cross-sectional area of air flow, respectively. r i and r O are inner and outer radii of the GDL for the in-plane permeability measurement. The Darcy coefficients in through-plane and in-plane directions are thereby obtained from relations between the flow rates and the pressure difference.   Figure 10 illustrates a parallel-serpentine flow field in the cathode of a PEFC. Flow analyses for the GDL between the flow channels indicated in the blue and red circles. The pressure difference between the channels in the former part is rather smaller than the lattar. This comparison represents that between the parallel and serpentine flow fields.

Results and discussion
The GDL is modeled with the flow channel and rib as presented in Fig. 11.
Cell temperature is 75 • C and water vapor pressure shall be the saturation vapor pressure at 75 • C. Air utilization is 30%. Oxygen and nitrogen partial pressure is assumed to linearly change from the inlet to the outlet, yielding 1.60 and -0.81 kPa, respectively, between the flow channels in the red part at 1.0 A cm −2 . It is assumed that current distribution is uniform and there is no liquid water in the GDL to simplify the calculation. Thickness of the GDL, rib width, and channel width are 190 µm, 0.6 mm, and 0.5 mm, respectively. Viscosities of nitrogen and oxygen of 20.00 × 10 −6 and 23.26 × ×10 −6 Pa·s, respectively [13], and binary diffusion coefficient for nitrogen-oxygen mixture of 2.59 × 10 −5 m 2 s −1 [14] are used for the LBM-BMD calculation. Oxygen flow velocity distribution at 1.0 A cm −2 is presented in Fig. 12. Oxygen amount consumed by the electrochemical reaction is calculated with the Faraday's law.

Channel Channel
Channel Channel Rib Figures 12(a) and (b) depict oxygen flow velocities for the GDLs without and with pressure differences between the flow channels, respectively. Oxygen is transported to the surface of the MEA mainly by diffusion in the case without the pressure differences since there is small forced convection thorough the GDL. On the other hand, oxygen flow velocity is rather larger under the rib in the case with pressure difference that leads to forced convection in the GDL. The forced convection also enhances the discharge of liquid and vapor product water in actual cells. The forced convection in the GDL by the pressure difference between flow channels plays a significant role on the exhaust of the product water and oxygen transport for the interdigitated flow field [15,16].

Conclusion
In this chapter, an anisotropic 2D GDL model is proposed by comparisons of the through-plane and in-plane permeabilities between those obtained by LBM calculation and permeability measurements. The modeled carbon fiber structure agrees well with the actual GDL in terms of the permeability. Moreover, the difference of oxygen flow in GDLs with parallel and serpentine flow channels is visualized with the oxygen-nitrogen two components LBM-BMD calculations using the above anisotropic GDL model. This procedure can be used to optimize the GDL porous structure, flow field patterns, and operation conditions of PEFCs. Liquid-gas two phase modeling [17,18], modeling, modeling of microporous layers [19][20][21], and expansion to 3D modelings are future studies.