Selection and peer-review under the responsibility of the scientific committee of the CEN2022.
Applied Energy Symposium 2022: Clean Energy towards Carbon Neutrality (CEN2022) April 23-25, 2022, Ningbo, China
Paper ID: 141
Pore-scale Investigation of a Gas Diffusion Layer Model Coupled Mass Transfer and Two-phase Flow: Effects of the Compression Ratio
Lingyi Guo1, Li Chen 1*, Wen-Quan Tao 1
1 Key Laboratory of Thermo-Fluid Science and Engineering of MOE, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi, 710049, China (*Corresponding Author)
ABSTRACT
Reducing water flooding and enhancing oxygen transport in the gas diffusion layer (GDL) of proton exchange membrane fuel cells are of great importance for optimizing cell performance. In this study, a pore- scale model based on the lattice Boltzmann method is proposed, which considers two-phase flow, oxygen diffusion and electrochemical reaction. The pore-scale model is then adopted to explore multiphase reactive transport processes in the GDL. Effects of compression on the liquid water saturation and current density are explored. It is found that, by increasing the compression ratio from 0.2 to 0.4, the current density increases by 4.7%, with obvious relieve of the water flooding under the rib and increment of the number of reaction sites under the gas channel. Besides, the results demonstrate that while reducing the total saturation in the GDL is important, decreasing the local saturation near the interface of the microporous layer/GDL is also crucial for enhancing cell performance. The coupled method can directly reveal the effects of GDL compression on water flooding and oxygen reactive transport, which provides a new insight for optimizing cell performance.
Keywords: proton exchange membrane fuel cell, liquid water, gas diffusion layer, oxygen reactive transport, compression ratio.
INTRODUCTION
The proton exchange membrane fuel cell (PEMFC) is a typical device that directly converts chemical energy to electrical energy. As one of the key components in PEMFC, the gas diffusion layer (GDL) has been studied extensively to improve the cell performance [1]. One of the crucial tasks for improving GDL performance is to enhance the oxygen transport as well as the liquid water migration [2, 3].
Due to the assembly pressure of the fuel cell stack, GDL is clamped together with other components, resulting in structure deformation and change of permeability and pore size. Excessive compression will
increase mass transfer resistance and change the holdup of water in the GDL, resulting in poor performance.
During the past decade, many studies have been conducted to investigate the effects of GDL compression on water transport and mass transfer [4-6].
In all the pore-scale numerical studies about effects of GDL compression, the liquid water into the GDL was assumed with constant pressure or constant velocity boundary condition, neglecting the fact that such boundary condition with constant value cannot be maintained as the liquid water will hinder the oxygen reactive transport and in turn the formation rate of liquid water will be reduced. Consequently, it’s necessary to develope a method coupled the oxygen reaction and diffusion with the air-water two-phase flow. To the best of our knowledge, studies related to this point are rare.
Only Zhang et al. [7] performed such pore-scale numerical study.
In this paper, a pore-scale model coupling the air- water two-phase flow with the oxygen transport is developed to investigate effects of the compression ratio (CR) of GDL. The GDL domain as well as the coupled processes are introduced in Section 2. In Section 3, effects of compression ratio (CR) are explored, focusing on the liquid water migration and average saturation, oxygen concentration, and current density. Finally, in Section 4 important conclusions are drawn.
METHOD AND MODEL
2.1 The computational domain and physicochemical processes
The computational domain is shown in Fig. 1 with size of W×H=2048×486 µm and space resolution of 1 µm, which contains a gas channel (GC) with width of 1024 µm, two ribs with width of 512 µm, a microporous layer (MPL) and a GDL. Along the in-plane (IP) direction, the GDL domain is divided into three parts: two compressed regions under the rib and one uncompressed region under the GC.
According to Tötzke et al. [8], excessive compression deforms the uniform structure of the GDL, bending the surface of the uncompressed region into an arc shape and reducing the thickness of the compressed region.
Such arc shape means that the middle uncompressed region maintains the original thickness of L, while the parts under the rib are compressed to thickness of l0, as shown in Fig. 1. Compression ratio (CR) is employed to measure the GDL bending, which is defined as the ratio of the reduced thickness (l0) versus the original GDL (L) [9]
=L l− 0
CR L (1)
Based on the above description, the thickness of the GDL (l(x)) at the location x can be described by a quadratic function
0 2 0
0
( 512)( 1536) 512 x 1536 (x) 512
l 0 x<512 or 1536 x 2048
l L x x l
l
− − − + ≤ ≤
=
< < ≤
（ ）
（ ）
(2)
Accordingly, the local GDL porosity (ε(x)) is calculated based on the original (or maximum) porosity (ε0) of 0.8 [9]
ε( )x =ε0l x( )
L (3)
Based on the above porosity distribution, the GDL is reconstructed by randomly placing solid particles of 7 µm radius until the local porosity ε(x) is satisfied.
Under the GDL bottom is the MPL with a thickness of 50 µm and the porosity of 0.65. Note that the MPL intrudes into the GDL with thickness of 10 µm, to reflect the integration of the actual MPL/GDL interfaces, as shown in the zoom-in image B of Fig. 1. Besides, considering the micropores and particles in the MPL are extremely small compared with that in the GDL [4], , solid particles of 2 µm radius are randomly placed in the MPL until the prescribed porosity of 0.65 are satisfied.
Meanwhile, as shown in the zoom-in image C in Fig.
1, totally 12 cracks are considered and evenly distributed
across the MPL, with the same width as 10 µm. The size of these cracks is determined based on the experimental results in Ref. [12]. The red particles and cavity beneath the MPL are the reaction boundary with thickness of 11 µm.
In the above 2D domain, coupled processes, including two-phase flow, oxygen mass transport and the chemical reaction, are described as follows. The processes are described as follows. Oxygen with prescribed concentration at the GC top wall diffuses through the GC, GDL, MPL and reaches the reaction boundary where the oxygen reduction reaction occurs.
The water generated by the reaction, which is assumed to be liquid phase, accumulates inside the bottom cavities until the cavity is fully filled, then invades into the GDL from random sites. The invading rate of liquid water follows neither constant pressure condition nor constant velocity condition. Instead, the rate is calculated each time step according to the real oxygen consumption rate, and thus the generating rate is not constant and actually is dynamically changed. The details of determining the oxygen consumption rate and the liquid water invading rate are discussed as follows.
In practice, after entering and further transporting through the MPL, the oxygen will diffuse into the CL.
Then the oxygen transfers in the macrospores among carbon agglomerates inside the CL, and is further consumed at the local structure around a carbon particle.
The transport resistance in the CL is considered, leading to the following reactive boundary condition [10]
= ∂ =
∂ + +
oxygen
oxygen GDL oxygen
CL CL,local elec
1 1
( )
J D C C
n R R
k
(4)
where J is the flux, C is the local concentration, and n is the normal direction of reactive surface. R is the transport resistance, and RCL is calculated as δCL/DCL
(1×10-5 m/3×10-5 m2∙s), respectively, with δ (1×10-5 m) as the thickness and D (3×10-5 m2∙s) as the effective diffusivity. RCL,local (13 s∙m-1) is the local transport resistance as discussed above. Based on the Butler- Volmer equation, the reaction rate constant kelec is expressed as follows [10]
(
α)
α η η
=
oxygen elec 0
oxygen,ref
exp(- )-exp(- 1- ) 4
C F F
k i
FC RT RT (5)
where i0 (0.01 A∙m-2) stands for the reference exchange current density, and Coxygen,ref (40.96 mol∙m-3) is the reference oxygen concentration (40.96 mol∙m-3).
α (0.61), F, R, T (353 K), and η (1 V) are the transfer coefficient, Faraday constant, universal gas constant, temperature, and overpotential, respectively.
Fig. 1. The two-dimensional computational domain and boundary conditions
GC210 μm
H=486 μm
0
Periodic Periodic
W=2048 μm
Rib
GDL MPL Coxygen=7.68 mol∙m-3& No-slip
Reaction boundary & No-slip Oxygen Water
5 μm y
x MPL 50 µm
Boundary
Crack 10 μm
L l0
Compressed surface Compressed surface
10 µm A
B
C
Based on Eq. (4), the water generating rate can be determined
water =0.5oxygen
J J (6)
The boundary conditions of the above processes are as follows. The periodic boundary condition is adopted for both mass transport and two-phase flow at the left and right boundaries. For the mass transport, constant oxygen concentration is prescribed at the GC top wall, reactive boundary condition (Eq. (4)) is adopted at the surface of the bottom red particles and cavity, and non- flux boundary condition is employed for the interior fluid-solid interface. For the air-water two-phase flow, non-slip boundary condition is applied for the fluid-solid interfaces. The contact angle for the carbon particles is θ=120°.
Detailed information about the physicochemical model and the lattice Boltzmann method can be found in our previous work [11].
RESULTS AND DISCUSSION
In this section, the pore-scale model is used to investigate effects of the CR on multiphase oxygen reactive transport.
3.1 Liquid water dynamic behaviors
Figs. 2 (a)-(c) show the dynamic behaviors of the liquid water in the domain. Three cases with CR as 0.2, 0.3 and 0.4 are studied, and herein are called as compression-0.2, compression-0.3 and compression-0.4.
In Fig. 2(a), for the compression-0.2 case at t=150000, the liquid water mainly enters the GDL through the cracks in the MPL from the domain bottom. The cracks as main paths for water migration indicate a typical behavior of capillary fingering. The liquid water clusters, marked from 1 to 12 in Fig. 2(a), accumulate at the MPL/GDL interface after breakthrough the MPL. At t=250000, it can be observed that the cluster 11 firstly
touched the bottom of the rib. Subsequently, it transports along the rib bottom, which is consistent with the description of Sakaida et al. [4]. Simultaneously, the coalescence between the water clusters frequently occurs at the interface of GDL/MPL, such as that between clusters 1-2, that between clusters 6-7-8, and between clusters 10-11-12. These coalescent clusters are called sub-clusters, marked from (a)-(c) in Fig. 2(a). Due to large local transport resistance, it is noteworthy that the aggregated water clusters at the MPL/GDL interface flow back into the micropores inside the MPL and cause the local water flooding, as shown in Fig. 2(a) with green arrow. Afterward, cluster 4 and sub-clusters (b) under GC continuously breakthrough at t=4100000 and 820000, and form droplets A and B. Subsequently, the breakthrough droplet A is snapped off when it grows up and touches the hydrophilic rib at t=1140000. The snap- off phenomenon, which means the disconnection of the water clusters, happens because as the local pore size suddenly increases, the capillary force and local resistance decreases, leading to the acceleration of the liquid water movement. However, the supply rate of the liquid from the underlying clusters cannot support such a suddenly increasing flow rate. As a result, droplet A is separated from cluster 4 and adheres to the hydrophilic surface of rib.
The dynamic behaviors of the liquid water in compression-0.3 and compression-0.4 cases are similar.
Thus, they are not repeated here for the purpose of brevity. Here effects of compression are further discussed. Firstly, it is obvious that in all the cases, the liquid water clusters breakthrough at the middle of the GC/GDL interface, where the local porosity is the biggest.
The biggest local porosity is caused by the compression of the rib. Thus, setting a hydrophilic block at the middle of the top GC to guide the breakthrough droplet flow out of the GDL may alleviate the water flooding and enhance
Fig. 2. Evolution of the liquid water distribution in the domain with different CR. (a) Compression-0.2; (b) compression-0.3; (c) compression-0.4 (dt=1×10-8 s, dx=1×10-6 m).
the cell performance. Secondly, as shown in Figs. 2(a)-(c), the raise of the CR makes more liquid water backflow from the GDL to the MPL, especially under the compressed surface (marked by the red arrows). The higher CR decreases the porosity differences between MPL and part of GDL under the compressed surface, and impedes amount of liquid water flow into the GDL.
3.2 Oxygen concentration distribution
Corresponding to the results of the two-phase flow, the time variations of the oxygen concentration for compression-0.2 case, are displayed in Fig. 3. Due to the consumption of oxygen at the reaction boundary, it is obvious that the oxygen concentration decreases from the top of the GC to the bottom of the MPL. Besides, due to the obstruction of the rib, the oxygen concentration along the IP direction is also different. In Fig. 3, the oxygen concertation of GDL under the rib is much lower than that under the GC.
The above characters of the oxygen concertation reflect the effects of structure on the mass transfer.
However, due to the coupled process of the two-phase flow and mass transfer, the dynamic behaviors of the water also impact a lot. The first one is the deformation and movement of the liquid water clusters around the rib.
As shown in the zoom-in image in Fig. 3 at t=900000, droplet A grows up at the site near the lateral wall of the
rib and adheres to it. The adherent droplet A and water cluster 4 block the oxygen transfer from the GC to the GDL under the rib. Subsequently, droplet A is snapped off at t=1140000 and adheres to the surface of the hydrophilic rib, reopening the direct pathway for mass transfer under the rib. Fig. 4 counts the water-free points on the reactive boundary (called reaction sites) and their corresponding local current density (Ilocal) at IP direction.
It can be observed from Fig. 4(a) that Ilocal at x=300-600 increases (from t=9000000 to t=1100000), reflecting effects of the hydrophilic rib reported in Ref [12].
Secondly, the coalescence of water clusters at the MPL/GDL interface increases the oxygen concentration in the GDL and simultaneously isolates more reaction sites on the pathway of mass transfer. For instance, on the one hand, the coalescence at the MPL/GDL interface will block the pathway of oxygen diffusion and increase the oxygen concentration around the reaction sites. It is shown in Fig. 3(a), the oxygen concentration in the zoom- in images increases from 6.34 mol m-3 at t=250000 to 6.80 mol m-3 at t=410000. On the other hand, the
coalescence at the MPL/GDL interface simultaneously Fig. 3. Evolutions of the oxygen concentration
distributions in the compression-0.2 case.
Fig. 4. The local current density along the IP direction for (a)compression-0.2 case at t=900000 and 1100000, (b) the cases with different CRs at t=1100000
reduces the reaction sites. Based on Eq. (4), the increment of the oxygen concentration and decrement of reaction sites are reverse variables to control the average current density. As a result, the local saturation at the MPL/GDL interface is a significant evaluation criterion to judge the relative effects of the coalescence.
3.3 Liquid water saturation and current density
The local saturation (Slocal, y) along the TP direction (y direction) is plotted in Fig. 5. Slocal, y is defined as the ratio between the local volume of liquid water to the entire volume of the void space at a certain y position. In Fig.
5(a), the curve of Slocal,y of compression-0.2 case increases with time proceeds and eventually converges to a stable curve. The can be divided into two parts: the part of the MPL and GDL. Along the TP direction, the Slocal, y of the MPL part is almost identical except for some small fluctuations. At the GDL part, there is a maximum Slocal,y
at the MPL/GDL interface, reflecting the water accumulation caused by the suddenly increased porosity from MPL to GDL. Then Slocal,y decreases to the interface of GC/GDL.
The CR has affected Slocal, y a lot. As shown in Fig. 5(b), the average Slocal,y in the MPL increases with CR increment. The increment of the Slocal,y inside the MPL part can be ascribed to the accumulation of the liquid water at the MPL/GDL interface, which is described in
section 3.2. However, in the GDL part, with the CR increasing, there is no such clear tendency of Slocal,y
among the three cases.
The time variations of the total saturation (Stotal) and the current density (Iave) are presented in Fig. 6. Stotal is the ratio between the volume of liquid water and that of the total pore volume. As the time proceeds, Stotal
gradually increases, and the increasing rate is reduced.
The maximum value of Stotal is obtained at the breakthrough moment. The maximum Stotal of compression-0.3 case is the lowest, followed by compression-0.4 case and compression-0.2 case. It is expected that the low Stotal inside the GDL will provide more oxygen transfer pathways, increasing the effective diffusion coefficient. However, the compression-0.3 case with the lowest Stotal does not have the highest Iave in Fig.
6. The maximum Slocal, y at the MPL/GDL interface explains that counterintuitive association. The zoom-in image in Fig. 5 displays the Slocal, y for the compression-0.4 case owns the minimum saturation peak, followed by compression-0.3 and compression-0.2 cases. The maximum local saturation at the MPL/GDL interface indicates higher mass transport resistance in the corresponding region, and the local current density in Fig.
4(b) also proves it.
In Fig. 4(b), it is obvious that the Ilocal under the rib is almost at the same level as CR varies, indicating that the mass transfer under the rib is less affected by CR. On the contrary, the Ilocal under the channel side is different. Due to the less accumulation of liquid water at the MPL/GDL interface under the channel, the compression-0.4 case has more reaction sites in the corresponding location compared with the other two cases. As a result, the Iave
of compression-0.4 case is the highest. In conclusion, it can be found that effects of CR on the GDL under the rib are inferior to that under the channel. Therefore, it is
Fig. 5. Slocal,y along TP direction for cases. Curves with time proceeds for (a) compression-0.2. (b) Comparison
of Slocal,y between different CRs at t=1100000.
Fig. 6. The Stotal and Iave for the cases with different CRs
necessary to build an effective pathway for liquid water removal and keep mass transfer under the channel clear.
CONCLUSION
Relieving the water flooding and enhancing oxygen transport in the GDL are of great importance for optimizing cell performance. To investigate effects of the GDL compression, in this study, a pore-scale model based on the lattice Boltzmann method is proposed, in which the coupled processes of mass transfer, two-phase flow and electrochemical reaction are considered. It is found that the lowest total liquid water saturation in the GDL does not necessarily generate the highest current density. The local distribution of the liquid water is also of great importance. Particularly, the case with lower local saturation near the MPL/GDL interface may own better cell performance, even if in this case the total saturation in the GDL is not the lowest. Moreover, compared with the local current density under the rib, the higher local current density under the channel directly displays the importance of building an effective pathway for draining liquid water and enhancing the mass transfer under the channel. The relative increment of CR is a possible method to realize both. With an increased compression ratio from 0.2 to 0.4, the current density increases by 4.7%, and the water flooding under the rib is obviously relieved and simultaneously increases the number of the reaction sites under the channel. From the results of the coupled model, focusing on the local parameters variation of the GDL and changing the local structures at the corresponding region give a new insight to improve the cell performance.
ACKNOWLEDGEMENT
We acknowledge the support of National key research and development project (2021YFB4001701) and the Fundamental Research Funds for the Central Universities.
REFERENCE
[1] Xing L, Shi WD, Su HN, Xu Q, Das PK, Mao BD, et al.
Membrane electrode assemblies for PEM fuel cells: A review of functional graded design and optimization.
Energy. 2019;177:445-64.
[2] Fadzillah DM, Rosli MI, Talib MZM, Kamarudin SK, Daud WRW. Review on microstructure modelling of a gas diffusion layer for proton exchange membrane fuel cells.
Renewable & Sustainable Energy Reviews.
2017;77:1001-9.
[3] Chen Q, Niu Z, Li H, Jiao K, Wang Y. Recent progress of gas diffusion layer in proton exchange membrane fuel
cell: two-phase flow and material properties. Int J Hydrogen Energy. 2021;46:8640-71.
[4] Sakaida S, Takahashi Y, Tanaka K, Konno M.
Improvement of polymer electrolyte fuel cell performance by gas diffusion layer with relatively hydrophilic surface.
[5] Simon C, Hasche F, Gasteiger HA. Influence of the Gas Diffusion Layer Compression on the Oxygen Transport in PEM Fuel Cells at High Water Saturation Levels. Journal of the Electrochemical Society. 2017;164:F591-F9.
[6] Toetzke C, Gaiselmann G, Osenberg M, Arlt T, Markoetter H, Hilger A, et al. Influence of hydrophobic treatment on the structure of compressed gas diffusion layers. Journal of Power Sources. 2016;324:625-36.
[7] Zhang D, Cai Q, Gu S. 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.
2018;262:282-96.
[8] Toetzke C, Gaiselmann G, Osenberg M, Bohner J, Arlt T, Markoetter H, et al. Three-dimensional study of compressed gas diffusion layers using synchrotron X-ray imaging. J Power Sources. 2014;253:123-31.
[9] Chi PH, Chan SH, Weng FB, Su A, Sui PC, Djilali N. On the effects of non-uniform property distribution due to compression in the gas diffusion layer of a PEMFC. Int J Hydrogen Energy. 2010;35:2936-48.
[10] Ferreira RB, Falcao DS, Oliveira VB, Pinto AMFR.
Experimental study on the membrane electrode assembly of a proton exchange membrane fuel cell:
effects of microporous layer, membrane thickness and gas diffusion layer hydrophobic treatment. Electrochim Acta. 2017;224:337-45.
[11] Zhou C, Guo L, Chen L, Tian X, He T, Yang Q. Pore- scale modeling of air-water two phase flow and oxygen transport in gas diffusion layer of proton exchange membrane fuel cell. Energies. 2021;14.
[12] Chen L, Luan H-B, He Y-L, Tao W-Q. Pore-scale flow and mass transport in gas diffusion layer of proton exchange membrane fuel cell with interdigitated flow fields. International Journal of Thermal Sciences.
2012;51:132-44.
