Next Article in Journal
Application and Visualization of Fluorescent-Tagged Antiscalants in Electrodialysis Processing of Aqueous Solutions Prone to Gypsum Scale Deposition
Previous Article in Journal
Performance Comparison of Polymeric and Silica-Based Multi-Bed Pervaporation Membrane Reactors during Ethyl Levulinate Production
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Microstructure Reconstruction and Multiphysics Dynamic Distribution Simulation of the Catalyst Layer in PEMFC

1
State Key Laboratory of Advanced Technology for Materials Synthesis and Processing, Wuhan University of Technology, Wuhan 430070, China
2
Foshan Xianhu Laboratory of the Advanced Energy Science and Technology Guangdong Laboratory, Xianhu Hydrogen Valley, Foshan 528200, China
3
School of Automotive Engineering, Wuhan University of Technology, Wuhan 430070, China
4
State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China
*
Author to whom correspondence should be addressed.
Membranes 2022, 12(10), 1001; https://doi.org/10.3390/membranes12101001
Submission received: 5 September 2022 / Revised: 25 September 2022 / Accepted: 11 October 2022 / Published: 14 October 2022

Abstract

:
Due to the complexity of both material composition and the structure of the catalyst layer (CL) used in the proton-exchange membrane fuel cell (PEMFC), conjugated heat and mass transfer as well as electrochemical processes simultaneously occur through the CL. In this study, a microstructure model of CL was first reconstructed using images acquired by Nano-computed tomography (Nano-CT) of a real sample of CL. Then, the multiphysics dynamic distribution (MPDD) simulation, which is inherently a multiscale approach made of a combination of pore-scale and homogeneous models, was conducted on the reconstructed microstructure model to compute the corresponded heat and mass transport, electrochemical reactions, and water phase-change processes. Considering a computational domain with the size of 4 um and cube shape, this model consisting of mass and heat transport as well as electrochemical reactions reached a stable solution within 3 s as the convergence time. In the presence of sufficient oxygen, proton conduction was identified as the dominant factor determining the strength of the electrochemical reaction. Additionally, it was concluded that current density, temperature, and the distribution of water all exhibit similar distribution trends, which decrease from the interface between CL and the proton-exchange membrane to the interface between CL and the gas-diffusion layer. The present study not only provides an in-depth understanding of the mass and heat transport and electrochemical reaction in the CL microstructure, but it also guides the optimal design and fabrication of CL components and structures, such as improving the local structure to reduce the number of dead pores and large agglomerates, etc.

1. Introduction

Renewable energy sources such as solar energy and wind energy are unstable and intermittent during generation, and thus these valuable electric energies are difficult to apply continuously and stably. The employment of PEMFC combined with hydrogen production from water electrolysis may greatly improve the utilization rate and stability of renewable energy [1]. The PEMFC consists of a bipolar plate (BPP), proton-exchange membrane (PEM), catalyst layer (CL), microporous layer (MPL), and gas-diffusion layer (GDL). The role of PEM is to conduct protons from the anode to the cathode, and the most commonly used PEM is Nafion series membranes which are perfluorinated sulfonic acid membranes containing both hydrophobic PTFE backbone and hydrophilic sulfonic acid group [2,3]. Due to its important functions and complicated structure, CL has already been one of the most important components in PEMFC [4]. However, it is difficult to characterize the heat and mass transport as well as electrochemical reaction processes in the CL via visualization and experimental in-situ methods. Therefore, establishing numerical microstructure models of CL is of great significance to understand and capture the complex processes in it.
Nano-Computed Tomography (Nano-CT) technique can provide lossless observation images of three-dimensional microstructures. Shi et al. [5] used Nano-CT to characterize the three-dimensional micro-crack structure of low-volatile coal and simulated the permeability characteristics of single-phase water flow in the micro-crack network using COMSOL. Zhao et al. [6] studied the pore structure characteristics of coal via Nano-CT, discussed the computational fluid dynamics simulation method based on the micro-nano pore structure, and then carried out the permeability simulation of the reconstructed microstructure. The results showed that the permeability of coals is low and anisotropic in all three different directions. Cetinbas et al. [7] reconstructed the PEMFC cathode catalyst layer (CCL) containing high surface-area porous carbon (HSC) using Nano-CT and simulated the effects of HSC morphology, agglomerate structure and relative humidity on the local oxygen transport resistance, concluding that smaller size or lower aspect ratio agglomerates facilitate better O2 transfer to the catalyst surface at high relative humidity. Braaten et al. [8] studied the degradation process of the PEMFC CCL by Nano-CT characterization and discovered that the combination of HSC support and PtCo catalyst was able to limit the loss of ECSA and Pt band formation in the CL better than the combination of Vulcan and pure Pt. Cui et al. [9] used Nano-CT to reconstruct the PEMFC CCL to study how various dimensions carbon materials affect the electrochemical performance of catalyst coated membrane (CCM). The outcomes showed that increasing the electrochemical performance of CCM by using carbon nanotubes in the CCL with low Pt loading is a successful strategy.
In addition to the experimental characterization of Nano-CT for microscopic reconstruction of CLs, numerical simulations are effective methods to study the structure and composition on the performance of PEMFCs. The conventional macroscopic homogeneous model can study the effect of flow channel structures and humidification conditions on cell performance. However, the relationship between Pt loading, Nafion content and cell performance is not considered in the macroscopic model [10,11,12]. The agglomerate model establishes the relationship between the content of CL and the cell performance. Sun et al. [13] built a two-dimensional CL polymer electrolyte coated catalyst particle agglomerate model which was used for fuel cell cathode half-cell simulations. The results showed that the catalyst utilization in the cathode agglomerates was extremely low at high current densities due to the limited catalysts. By modeling CL agglomerate and one-dimensional PEMFC model, Sohn et al. [14] carried out numerical simulations of CCL and concluded that cell with smaller pore in CCL exhibits more concentration loss at high current densities. Hosseini et al. [15] performed simulations of two-dimensional open-anode PEMFC utilizing an agglomerate model. According to the simulation results, a rise in relative humidity promotes water’s diffusion from the cathode to the anode. Cosse et al. [16] proposed and validated a transient PEMFC agglomerate model for predicting the peak current of double-layer capacitor, and the simulation results showed that the gas supply stoichiometry has little effect on the absolute peak of the short-circuit current. Yang et al. [17] investigated the effect of CL on the PEMFC cold start using an agglomerate model and concluded that a larger I/C decreases CL porosity and agglomerate pore size, which significantly reduces the critical ice fraction of cold start failure.
Although the agglomerate model considers the relationship between CL component content and cell performance, the effects of CL pore structure and solid phase structure on cell performance and lifetime have not been deeply studied. Therefore, it is significant to establish a mesoscopic scale model based on experimental characterization data to study the distribution of Multiphysics and enhance the performance of PEMFCs. In recent decades, the study of mesoscopic scale reconstruction [18,19] and simulation [20,21,22,23,24,25,26] in PEMFCs mainly focuses on the GDL. This is because there are nanoscale pore structure and complex electrochemical reaction processes in CL, the resolution of CL reconstruction is required to be higher, and the simulation process is more complicated. Babu et al. [27] used Nano-CT to reconstruct Pt-free catalysts and simulated the CL proton conduction under different Nafion content. Compared with the experimental results, it was shown that higher Nafion content could easily result in CL flooding and then affect proton conduction, which ultimately leads to the decline of fuel cell performance under high current density. Cetinbas et al. [28] used Nano-CT to reconstruct the agglomerate structure with secondary pores in CL, meanwhile, the microstructure with primary pores, Carbon and Pt particle distribution was numerically reconstructed based on Transmission Electron Microscopy (TEM) and Scanning electron microscope (SEM) characterization images. The two reconstructed structures were combined to form a computational domain, and the electron, proton and oxygen transport properties were investigated by using STAR CCM+. Satjaritanun et al. [29] investigated the transport properties in different scale components of the PEMFC through a hybrid model. Firstly, the microstructures of the gas-diffusion media were obtained by Nano-CT, then the transport phenomena in the membrane electrode assembly (MEA) was studied via Lattice Boltzmann Agglomeration Method. Meanwhile, the data calculated by Lattice Boltzmann Agglomeration was exchanged with the flow channel, and conventional CFD method was used to simulate the transport phenomena in flow channel. Finally, the factors of affecting the cell performance were studied.
In this study, the three-dimensional microstructure of the CL was reconstructed by Nano-CT experimental data, and the pore structural characteristics and gas-diffusion media tortuosity were analyzed. The Knudsen effect was considered to correct the transport properties, and a pore-scale model for heat and mass transport as well as electrochemical processes was built. For current research status, i.e., that Nano-CT cannot easily distinguish the Pt, carbon particles, and Nafion, the solid phase region was treated with a homogeneous method, and the multiphysics dynamic distribution characteristics of the CL were studied. The comprehensive simulation methods and results presented in this study not only demonstrate the importance of transport properties in CL, but also help in advancing the in-depth understanding of electrochemical reaction processes in PEMFC.

2. Microstructure Reconstruction and Physical Properties Correction of CL

2.1. Sample Preparation

In this paper, the MEA, consisted of a Nafion 211 membrane and 60% mass fraction Pt/C industrial catalyst, which was produced by Wuhan University of Technology Hydrogen Power Technology Co. Ltd. The Pt loading and active area of the MEA were 0.4 mg/cm2 and 25 cm2. The treatments of MEA were as follows: firstly, the GDL was peeled off from both sides of the MEA and then the CCM was cut into a 2 mm × 2 mm square. Secondly, the treated CCM sample was immersed in saturated cesium sulfate solution (CsSO4) with deionized water for 78 h to exchange H+ with Cs+. Finally, to avoid sulfate deposition, the samples were shocked with deionized water for 5 min and dried to remove excess water.

2.2. Nano-CT Imaging

Nano-CT imaging was performed at the 4W1A Topography Imaging Experiment Station of Beijing Synchrotron Radiation Facility. The sample was analyzed through Nano-CT with a resolution of 50 nm in the Zernike phase contrast mode of 8 keV and the exposure time was 10 s. A series of images were obtained when the sample was rotated at 0.3° between −70° and +70°, then the 2D images were corrected and 3D reconstruction was processed by Xradia XMReconstructor software. Finally, a computational domain of 4 μm × 4 μm × 3 μm was selected from the reconstructed microstructure, and histogram equalization and normalization were performed on it to enhance the contrast. The volume was binarized into solid and pore spaces using artificial thresholding, as shown in Figure 1; there were some isolated pore structures called dead pores.

2.3. Structural Characterizations of Reconstructed Model

The porosity of the reconstructed CL shown in Figure 1 was 32.4%, as obtained by Avizo software. It is consistent with the results obtained by the mercury intrusion method (MIP) experimental method, which indicates the feasibility and validity of the reconstruction method. Figure 2a,b show the pore size distribution and tortuosity distribution of the CL, respectively. The average pore size and average tortuosity of the reconstructed CL were 136 nm and 2.1, which is similar to ref. [28].
The three parameters of transport properties can be used to correct the CL physical properties in the subsequent section.

2.4. Correction of Transport Parameters in CL

The relative gas-diffusion coefficient D i R in the CL is the ratio of effective diffusion coefficient D i e f f to the molecular diffusion coefficient D i b u l k as listed in Table 1. The Knudsen effect is neglected, it can be expressed in terms of porosity ε and tortuosity τ as follows [30]:
D i R = D i e f f D i b u l k = ε τ
where the subscript i can represent the gas species of air in the cathode: oxygen, nitrogen, and water vapor at the cathode.
The gas-diffusion coefficient D i in the CL is related to the Knudsen diffusion coefficient D i k and the molecular diffusion coefficient D i b u l k , it can be solved by the Bosanquet equation [31]:
1 D i = 1 D i k + 1 D i b u l k
When the Knudsen effect is considered, the relative diffusion coefficient Equation (1) is modified as follows:
D i R = D i e f f D i b u l k = ε τ 1 ( 1 + D i b u l k D i k )
D i k = 2 3 ( 8 R T π M i ) 1 / 2 r
where R is the universal gas constant, T is the temperature, r is the average pore radius, M i is the relative molecular mass.
The effective electrical conductivity κ e l e e f f and effective proton conductivity κ i o n e f f of the reconstructed CL are solved as follows:
κ e l e e f f = σ e l e V e l e / τ e l e
κ i o n e f f = σ i o n V i o n / τ i o n
where σ e l e is the electron conductivity, σ i o n is the proton conductivity, V e l e and V i o n are the carbon volume fraction and ionomer volume fraction of the CL, respectively. τ e l e and τ i o n are the tortuosity of carbon and ionomer, respectively.
Based on the above equations, the gas diffusivity and conductivity in the reconstructed CL were corrected. Figure 3a shows the relationship between relative diffusion coefficient and porosity in the CL with and without Knudsen diffusion. The relative diffusion coefficient in the CL corrected by porosity and tortuosity is smaller than that corrected by Bruggmen. Additionally, the relative diffusion coefficient of the CL considering Knudsen diffusion was the smallest; Knudsen diffusion had a great effect on the relative diffusion coefficient. The above results of the relative diffusion coefficient correction are similar to ref. [35].
Figure 3b–d show the relationship between the effective diffusion coefficients and porosity in the CL for oxygen, water vapor, and nitrogen with and without Knudsen diffusion, respectively. The gas’s effective diffusion coefficient in the CL increased with the porosity, which is due to gas transport resistance decreasing with the porosity. Additionally, the gas’s effective diffusion coefficient in the CL when Knudsen diffusion is considered was smaller than those without Knudsen diffusion. In addition, the effective diffusion coefficient of the gas species obtained from the Bruggmen equation was relatively high. Figure 3e,f show the relationship between effective electron/proton conductivity and carbon/ionomer volume fraction, respectively. Both the effective electron and proton conductivity in the CL increased with the carbon/ionomer volume fraction. This is because more carbon particles and ionomer provide more pathways for electrons and protons to transport, respectively.

3. Mathematical Analysis of Reconstructed CL

3.1. Geometric Models

The above three-dimensional reconstructed microstructure of CL was imported into the Multiphysics simulation software COMSOL as shown in Figure 4. The pore phase and solid phase of the CL are represented in red and blue, respectively.
In this study, Nafion, C and Pt were not distinguished in the solid phase; they were divided according to the proportion obtained by the MEA experimental characterization.

3.2. Mathematical Model

There are complicated processes such as heat and mass transport, phase change, and electrochemical reaction in the CL of PEMFC. To simplify the calculation, the following assumptions were made in this study [36]:
(1) All the gas species are ideal gases.
(2) The water generated by the electrochemical reaction in the CL is ionomer water, and then diffuses toward the area of low concentration. The ionomer water can convert to water vapor when it reaches the pore/solid phase interface.
(3) The liquid water at the boundary can be completely removed in time.
(4) The heat, mass transport and electrochemical reactions are calculated using a homogenous model due to the difficulty in distinguishing various components of the solid phase and the amount of calculation in this study.
Based on the above assumptions, the mathematical model in this study includes the following governing equations [31,37]:
The mass conservation equation:
t ( ( 1 s l q ) ρ g ) + ( ρ g u ) = S m a s s
where s l q , ρ g , u and S m a s s represent the liquid water saturation, mixed gas density, mixed gas velocity vector and mass source term, respectively.
Inertial and viscous forces were ignored because they have little effect on performance in the model; the momentum conservation equation was simplified to:
  ( 1 s l q ) u = K g μ g p g
where K g , μ g , and p g represent mixed gas permeability, viscosity, and gas pressure, respectively
The energy conservation equation is:
( c p ρ e f f T ) t + ( u c p ρ e f f T ) = ( K e f f T ) + S Q
where c p , K e f f , ρ e f f , and S Q represent the constant pressure heat capacity, temperature, effective thermal conductivity, effective density, and energy source terms, respectively.
The component conservation equation is:
t ( ( 1 s l q ) ρ g Y i ) + ( ρ g u g Y i ) = ( ρ g D i e f f Y i ) + S i
where Y i , D i e f f , and S i represent the component concentration, the effective diffusion coefficient of the species i, and the component source term, respectively. The values of D i e f f are corrected after considering the Knudsen effect according to the method in the previous section.
The liquid water conservation equation is:
( s l q ρ l q ) t + ( K l q μ g K g μ l q ρ l q u g ) = ( ρ l q D l q s l q ) + S l q
where ρ l q and μ l q represent the density and viscosity of liquid water, K l q represents the liquid phase permeability, and K g represents the gas phase permeability. The liquid water diffusion coefficient D l q after considering the effect of capillary pressure difference can be expressed as:
D l q = K l q μ l q d p c d s l q
p c a = σ cos θ ( ε K 0 ) 0.5 × [ 1.42 ( 1 s l q ) 2.12 ( 1 s l q ) 2 + 1.26 ( 1 s l q ) 3 ] θ < 90
p c a = σ cos θ ( ε K 0 ) 0.5 × [ 1.42 s l q 2.12 s l q 2 + 1.26 s l q 3 ] θ > 90
where p c a , θ , and σ are capillary pressure, contact angle, and surface tension, respectively.
Charge conservation equation:
( κ e l e e f f ϕ e l e ) + S e l e = 0
( κ i o n e f f ϕ i o n ) + S i o n = 0
where ϕ e l e is the electrical potential, ϕ i o n is the proton potential, κ e l e e f f is the effective electrical conductivity, κ i o n e f f is the effective proton conductivity, and S e l e and S i o n are the electron current source term and the proton current source term, respectively. Only the electrochemical reaction of cathode was considered in this study:
S e l e = J c A c
S i o n = J c A c
J c = j c , r e f v ( C O 2 C O 2 , r e f ) γ × ( e α F R T η c )
where J c is the cathode electrochemical reaction rate, A c is the active specific surface area, j c , r e f v is the cathode reference exchange current density, C O 2 , r e f is the oxygen reference concentration, γ is the concentration index, η c is the cathode overpotential, α is the transport coefficient, and F is the Faraday constant. The subscript c indicates the cathode. The phase change equation can be referred to the literature [38].

3.3. Boundary Conditions and Initial Values

The model parameters and boundary conditions used in the model are listed in Table 2 and Table 3, respectively.
Since the computational domain of this study is part of the CL, the setting of the boundary conditions was determined by the PEMFC operating conditions and the location. In this study, the boundary conditions are given according to the literature [35] and were set as initial values as shown in Table 3.

3.4. Numerical Calculation Method

The number of grids in this model is 840,000, it has little effect on the calculation results when the grid number changes by 10%. The physical properties of the CL and the gas-diffusion coefficient were corrected after considering the tortuosity, porosity, and Knudsen effect. Finally, COMSOL Multiphysics 5.3a software was used to solve the governing equations of dynamic distribution of Multiphysics.

4. Results and Discussion

4.1. Oxygen Concentration Distributions

Figure 5a,b show the molar concentration distribution of oxygen in the pore phase at 1 s and 3 s, respectively. In the CL, diffusion is regarded as the only mode of oxygen transport. The oxygen concentration gradually decreased from the CL/GDL interface to the CL/PEM interface due to the consumption of electrochemical reaction. The difference of oxygen concentration between the two interfaces was 0.4 mol/m3 at 1 s and 0.6 mol/m3 at 3 s, respectively. In particularly, there are some dead pores in the computational domain where oxygen cannot enter as shown in Figure 5a,b. Therefore, the O2 concentration in dead pores is constant to zero, and no electrochemical reaction occurs in there. Figure 5c shows the average O2 molar concentration along the CL thickness direction at different moments. The difference between the maximum and minimum values of O2 molar concentration in the CL gradually became larger over time. The differences were 0.35 mol/m3, 0.4 mol/m3, 0.47 mol/m3, 0.6 mol/m3, and 0.61 mol/m3 at 0.5 s, 1 s, 2 s, 3 s, and 4 s, respectively. At 3 s and 4 s, the difference was almost unchanged, indicating that the electrochemical reaction in CL tended to be stable.

4.2. Current Density Distributions

Figure 6 is the contour of the current density distribution calculated by Equation (19) at six different moments. With the progress of the reaction, the electrochemical reaction gradually proceeded from the CL/PEM interface to the CL/GDL interface, and the current density always decreased from the CL/PEM interface to the CL/GDL interface. Because the oxygen concentration was 10.1 mol/m3 and water vapor concentration was 0 mol/m3 at the initial moment in this study, it is difficult for the ionomer in the CL to conduct protons; the oxygen reduction reaction starts at the CL/PEM interface where protons are sufficient. The water generated in CL gradually increased and diffused from the CL/PEM interface to the CL/GDL interface, which led to the increase of proton conductivity in the solid phase. However, there was a gradient of current density from the CL/PEM interface to the CL/GDL interface. This phenomenon indicates that the proton conduction is the dominant factor that affects the electrochemical reaction when the oxygen is sufficient, which are consistent with ref. [36]. The current density distribution in the calculation domain during 3–4 s did not change obviously, and it was in a dynamic equilibrium state.
In addition, the current density was larger in the solid phase region near the pore phase. This is because the oxygen diffusion coefficient in the CL pore phase is much larger than that in the solid phase, and there is sufficient oxygen in the solid phase region near the pore, which leads to the more intense electrochemical reaction and larger current density at this pore/solid interface. To sum up, the oxygen diffusion and the intensity of the electrochemical reaction are all restricted in the large solid phase due to the lack of sufficient contact area with the pores.

4.3. Temperature Distributions

Figure 7a shows the average temperature distribution of the CL/GDL interface, the CL/PEM interface, and the CL middle section (along the Z-axis direction) within 4 s. The temperature increased over time for all three different sections. Additionally, the temperature of the CL/PEM interface was the highest, followed by the CL middle section, and the lowest temperature was at the CL/GDL interface. After 3 s, the temperature of all the three interfaces did not increase, indicating that the electrochemical reaction became stable. Figure 7b shows the contour of temperature distribution on different sections along the CL thickness direction when the reaction was stable, which is similar to the current density distribution in Figure 6. Besides, the temperature was higher in the solid phase region near the pore. The overall temperature in the CL increased by 0.6 K from the beginning of the reaction to stability.

4.4. Liquid Water and Membrane Water Distribution

Figure 8a shows the profiles of the average liquid water saturation calculated by Equation (11) along the Z-axis at different moments. The liquid water saturation gradually decreased from the CL/PEM interface to the CL/GDL interface. This is because the electrochemical reaction firstly generated water at the CL/PEM interface and the water gradually diffused to the GDL/CL interface. The average liquid water saturation in the CL gradually increased with time, and the liquid water saturation profiles at 3 s and 4 s basically overlapped, indicating that the reaction was basically stable at 3 s. Additionally, the distribution of liquid water saturation is consistent with the distribution of current density. Figure 8b shows the distribution of the average liquid water saturation in the y = 0 section when the reaction was stable, which is consistent with the curve at 3 s in Figure 8a.
Figure 8c shows the profiles of the average ionomer water content in CL over time, assuming that the initial ionomer water content λ is 6. The water generated by electrochemical reaction diffused toward the region with a lower water concentration. Because the electrolyte in the CL absorbs water vapor and reaches saturation, the ionomer water did not increase any more after 3 s, and the stable average ionomer water content in the CL was around 17.5.
Figure 8d is the iso-surface of ionomer water along the Z-axis direction of CL when the reaction was stable, which is the same as the distribution of liquid water at the same time. It can be found that the distribution of ionomer water is hilly and undulating on any iso-surface due to the complex pore structure, the strength of electrochemical reactions, and the wettability of the ionomer inside the CL.
The distribution of oxygen, current density, temperature, and water in the computational domain tended to be stable after 3 s. Considering the dimensions of the computational domain, this transition time from the beginning of reaction to the steady state is basically consistent with the conclusion of ref. [40], which indicates the feasibility and validity of the study results.

5. Conclusions

In this study, a multiscale approach was developed to incorporate the effective transport properties computed from 3D reconstructed CL with a macroscopic model created in COMSOL Multiphysics 5.3a. In the first step, a three-dimensional microstructure geometry of CL was reconstructed using Nano-CT images of a commercial sample as well as its experimental data. In the following step, the distribution of pore sizes and tortuosity was analyzed. In the pore phase, the transport properties were first corrected by the Knudsen effect, and then the COMSOL Multiphysics 5.3a was used to simulate the dynamic processes of heat, mass transport, and electrochemical reaction in the CL. The results are listed as follows:
(1) It can be seen that the relative diffusion coefficient and effective diffusion coefficient of the gas in CL were smaller after considering Knudsen diffusion, indicating that Knudsen diffusion had a significant effect on the gas’s effective diffusion coefficient.
(2) Due to the difficulties in distinguishing multiple components of CL in the solid phase and the desire to reduce the computational cost, a Multiphysics simulation method was proposed in which the mesoscopic-scale structure was combined with the pore-scale model and the homogeneous model.
(3) It is evident that current density, temperature, and water distribution all exhibited similar distribution trends, decreasing from the CL/PEM interface to the CL/GDL interface.
(4) Under the introduced boundary and operating conditions, the mass transport, heat transport, and electrochemical reactions reached a stable solution after 3 s as the convergence time. Upon completion of this 3 s, oxygen concentration, current density, temperature, and water content no longer changed.
The simulation results help to understand the mass and heat transport and electrochemical reaction in the CL microstructure. It has guiding significance for the optimal design and fabrication of porous structures, such as improving the local structure to reduce the number of dead pores and large agglomerates, etc. Future research should focus on developing more precise reconstruction methods to realize the differentiation of the CL’s components, and then simulate and study the influence of components and structure on the electrochemical performance of the CL.

Author Contributions

Conceptualization, Z.Z. and H.Z.; methodology, Z.Z., H.S., X.Y., P.J., R.C. and H.Z.; writing—original draft, Z.Z. and H.S.; Writing—review & editing, H.B.H. and H.Z.; Software, H.S., X.Y., P.J. and R.C. Formal analysis, Z.Z. and H.S.; funding acquisition, H.Z. and M.P.; supervision, H.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China grant number (22179103, 21676207) and the Foshan Xianhu Laboratory Open Fund Key Project grant number (XHD2020-002).

Data Availability Statement

Data is contained within the article.

Acknowledgments

We thank the National Natural Science Foundation of China (22179103, 21676207) and the Foshan Xianhu Laboratory Open Fund Key Project (XHD2020-002) for financial support.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

D i R relative gas-diffusion coefficient (m2/s) σ i o n proton conductivity (S/m)
D i e f f effective diffusion coefficient (m2/s) τ e l e tortuosity of carbon
D i b u l k molecular diffusion coefficient (m2/s) τ i o n tortuosity of ionomer
D i gas-diffusion coefficient (m2/s) ρ density (kg/m3)
D i k Knudsen diffusion coefficient (m2/s) μ dynamic viscosity (kg/(m s))
R universal gas constant (J/(mol K)) ρ e f f effective density
T Temperature (K) σ surface tension
r average pore radius (m) θ contact angle
M i relative molecular mass (kg/mol) ϕ e l e electronic potential (V)
V e l e carbon volume fraction ϕ i o n ionic potential (V)
V i o n ionomer volume fraction γ concentration index
s l q liquid water saturation η c cathode overpotential (V)
u mixed gas velocity vector (m/s) α transport coefficient
S source term (kg/(m3 s), W/m3…)Subscripts
p g gas pressure (Pa) i gas species
c p constant pressure heat capacity (J/(mol K)) e l e electron
K e f f effective thermal conductivity (W/(m K)) i o n proton
Y i gas species mass fraction l q liquid water
K intrinsic permeability (m2) g gas phase
p c a capillary pressure (Pa) c cathode
J c cathode electrochemical reaction rate (A/m2) H 2 O water
A c active specific surface area (m−1) O 2 oxygen
C O 2 , r e f oxygen reference concentration (mol/m3) N 2 nitrogen
Greek letters
ε porosity
τ tortuosity
κ e l e e f f effective electrical conductivity (S/m)
κ i o n e f f effective ionic conductivity (S/m)
σ e l e electron conductivity (S/m)

References

  1. Sun, C.; Zhang, H. Review of the Development of First-Generation Redox Flow Batteries: Iron-Chromium System. ChemSusChem 2022, 15, e202101798. [Google Scholar] [CrossRef]
  2. Sun, C.; Negro, E.; Nale, A.; Pagot, G.; Vezzù, K.; Zawodzinski, T.A.; Meda, L.; Gambaro, C.; Di Noto, V. An Efficient Barrier toward Vanadium Crossover in Redox Flow Batteries: The Bilayer [Nafion/(WO3)x] Hybrid Inorganic-Organic Membrane. Electrochim. Acta 2021, 378, 138133. [Google Scholar] [CrossRef]
  3. Dong, P.; Xie, G.; Ni, M. Improved Energy Performance of a PEM Fuel Cell by Introducing Discontinuous S-Shaped and Crescent Ribs into Flowing Channels. Energy 2021, 222, 119920. [Google Scholar] [CrossRef]
  4. Majlan, E.H.; Rohendi, D.; Daud, W.R.W.; Husaini, T.; Haque, M.A. Electrode for Proton Exchange Membrane Fuel Cells: A Review. Renew. Sustain. Energy Rev. 2018, 89, 117–134. [Google Scholar] [CrossRef]
  5. Shi, X.; Pan, J.; Pang, L.; Wang, R.; Li, G.; Tian, J.; Wang, H. 3D Microfracture Network and Seepage Characteristics of Low-Volatility Bituminous Coal Based on Nano-CT. J. Nat. Gas Sci. Eng. 2020, 83, 103556. [Google Scholar] [CrossRef]
  6. Zhao, Y.; Sun, Y.; Liu, S.; Chen, Z.; Yuan, L. Pore Structure Characterization of Coal by Synchrotron Radiation Nano-CT. Fuel 2018, 215, 102–110. [Google Scholar] [CrossRef]
  7. Cetinbas, F.C.; Ahluwalia, R.K.; Kariuki, N.N.; De Andrade, V.; Myers, D.J. Effects of Porous Carbon Morphology, Agglomerate Structure and Relative Humidity on Local Oxygen Transport Resistance. J. Electrochem. Soc. 2020, 167, 013508. [Google Scholar] [CrossRef]
  8. Braaten, J.P.; Ogawa, S.; Yarlagadda, V.; Kongkanand, A.; Litster, S. Studying Pt-Based Fuel Cell Electrode Degradation with Nanoscale X-Ray Computed Tomography. J. Power Sources 2020, 478, 229049. [Google Scholar] [CrossRef]
  9. Cui, L.; Zhang, J.; Wang, H.; Lu, S.; Xiang, Y. The Effects of Different Dimensional Carbon Additives on Performance of PEMFC with Low-Pt Loading Cathode Catalytic Layers. Int. J. Hydrogen Energy 2021, 46, 15887–15895. [Google Scholar] [CrossRef]
  10. Fan, L.; Zhang, G.; Jiao, K. Characteristics of PEMFC Operating at High Current Density with Low External Humidification. Energy Convers. Manag. 2017, 150, 763–774. [Google Scholar] [CrossRef]
  11. Chen, H.; Liu, B.; Zhang, T.; Pei, P. Influencing Sensitivities of Critical Operating Parameters on PEMFC Output Performance and Gas Distribution Quality under Different Electrical Load Conditions. Appl. Energy 2019, 255, 113849. [Google Scholar] [CrossRef]
  12. Wang, Y.; Sun, Z.Y.; Yang, L. Enhancement Effects of the Obstacle Arrangement and Gradient Height Distribution in Serpentine Flow-Field on the Performances of a PEMFC. Energy Convers. Manag. 2022, 252, 115077. [Google Scholar] [CrossRef]
  13. Sun, W.; Peppley, B.A.; Karan, K. An Improved Two-Dimensional Agglomerate Cathode Model to Study the Influence of Catalyst Layer Structural Parameters. Electrochim. Acta 2005, 50, 3359–3374. [Google Scholar] [CrossRef]
  14. Sohn, Y.-J.; Yim, S.-D.; Park, G.-G.; Kim, M.; Cha, S.-W.; Kim, K. PEMFC Modeling Based on Characterization of Effective Diffusivity in Simulated Cathode Catalyst Layer. Int. J. Hydrogen Energy 2017, 42, 13226–13233. [Google Scholar] [CrossRef]
  15. Hosseini, M.; Afrouzi, H.H.; Arasteh, H.; Toghraie, D. Energy Analysis of a Proton Exchange Membrane Fuel Cell (PEMFC) with an Open-Ended Anode Using Agglomerate Model: A CFD Study. Energy 2019, 188, 116090. [Google Scholar] [CrossRef]
  16. Cosse, C.; Schumann, M.; Grumm, F.; Becker, D.; Schulz, D. Numerical Investigation of PEMFC Short-Circuit Behaviour Using an Agglomerate Model Approach. Energies 2020, 13, 4108. [Google Scholar] [CrossRef]
  17. Yang, L.; Cao, C.; Gan, Q.; Pei, H.; Zhang, Q.; Li, P. Revealing Failure Modes and Effect of Catalyst Layer Properties for PEM Fuel Cell Cold Start Using an Agglomerate Model. Appl. Energy 2022, 312, 118792. [Google Scholar] [CrossRef]
  18. Kotaka, T.; Tabuchi, Y.; Mukherjee, P.P. Microstructural Analysis of Mass Transport Phenomena in Gas Diffusion Media for High Current Density Operation in PEM Fuel Cells. J. Power Sources 2015, 280, 231–239. [Google Scholar] [CrossRef]
  19. Zenyuk, I.V. Gas-Diffusion-Layer Structural Properties under Compression via X-Ray Tomography. J. Power Sources 2016, 328, 364–376. [Google Scholar] [CrossRef] [Green Version]
  20. Chen, W.; Jiang, F. Impact of PTFE Content and Distribution on Liquid–Gas Flow in PEMFC Carbon Paper Gas Distribution Layer: 3D Lattice Boltzmann Simulations. Int. J. Hydrogen Energy 2016, 41, 8550–8562. [Google Scholar] [CrossRef]
  21. Satjaritanun, P.; Weidner, J.W.; Hirano, S.; Lu, Z.; Khunatorn, Y.; Ogawa, S.; Litster, S.E.; Shum, A.D.; Zenyuk, I.V.; Shimpalee, S. Micro-Scale Analysis of Liquid Water Breakthrough inside Gas Diffusion Layer for PEMFC Using X-Ray Computed Tomography and Lattice Boltzmann Method. J. Electrochem. Soc. 2017, 164, E3359–E3371. [Google Scholar] [CrossRef]
  22. Jinuntuya, F.; Whiteley, M.; Chen, R.; Fly, A. The Effects of Gas Diffusion Layers Structure on Water Transportation Using X-Ray Computed Tomography Based Lattice Boltzmann Method. J. Power Sources 2018, 378, 53–65. [Google Scholar] [CrossRef] [Green Version]
  23. Yu, J.; Froning, D.; Reimer, U.; Lehnert, W. Apparent Contact Angles of Liquid Water Droplet Breaking through a Gas Diffusion Layer of Polymer Electrolyte Membrane Fuel Cell. Int. J. Hydrogen Energy 2018, 43, 6318–6330. [Google Scholar] [CrossRef]
  24. Cetinbas, F.C.; Ahluwalia, R.K.; Shum, A.D.; Zenyuk, I.V. Direct Simulations of Pore-Scale Water Transport through Diffusion Media. J. Electrochem. Soc. 2019, 166, F3001–F3008. [Google Scholar] [CrossRef]
  25. Zhu, L.; Zhang, H.; Xiao, L.; Bazylak, A.; Gao, X.; Sui, P.-C. Pore-Scale Modeling of Gas Diffusion Layers: Effects of Compression on Transport Properties. J. Power Sources 2021, 496, 229822. [Google Scholar] [CrossRef]
  26. Shangguan, X.; Li, Y.; Qin, Y.; Cao, S.; Zhang, J.; Yin, Y. Effect of the Porosity Distribution on the Liquid Water Transport in the Gas Diffusion Layer of PEMFC. Electrochim. Acta 2021, 371, 137814. [Google Scholar] [CrossRef]
  27. Babu, S.K.; Chung, H.T.; Zelenay, P.; Litster, S. Resolving Electrode Morphology’s Impact on Platinum Group Metal- Free Cathode Performance Using Nano-CT of 3D Hierarchical Pore and Ionomer Distribution. ACS Appl. Mater. Interfaces 2016, 8, 32764–32777. [Google Scholar] [CrossRef]
  28. Cetinbas, F.C.; Ahluwalia, R.K.; Kariuki, N.; De Andrade, V.; Fongalland, D.; Smith, L.; Sharman, J.; Ferreira, P.; Rasouli, S.; Myers, D.J. Hybrid Approach Combining Multiple Characterization Techniques and Simulations for Microstructural Analysis of Proton Exchange Membrane Fuel Cell Electrodes. J. Power Sources 2017, 344, 62–73. [Google Scholar] [CrossRef] [Green Version]
  29. Satjaritanun, P.; Hirano, S.; Zenyuk, I.V.; Weidner, J.W.; Tippayawong, N.; Shimpalee, S. Numerical Study of Electrochemical Kinetics and Mass Transport inside Nano-Structural Catalyst Layer of PEMFC Using Lattice Boltzmann Agglomeration Method. J. Electrochem. Soc. 2020, 167, 013516. [Google Scholar] [CrossRef]
  30. Inoue, G.; Yokoyama, K.; Ooyama, J.; Terao, T.; Tokunaga, T.; Kubo, N.; Kawase, M. Theoretical Examination of Effective Oxygen Diffusion Coefficient and Electrical Conductivity of Polymer Electrolyte Fuel Cell Porous Components. J. Power Sources 2016, 327, 610–621. [Google Scholar] [CrossRef]
  31. Zhu, L.; Wang, S.; Sui, P.-C.; Gao, X. Multiscale Modeling of an Angled Gas Diffusion Layer for Polymer Electrolyte Membrane Fuel Cells: Performance Enhancing for Aviation Applications. Int. J. Hydrogen Energy 2021, 46, 20702–20714. [Google Scholar] [CrossRef]
  32. Sone, Y.; Ekdunge, P.; Simonsson, D. Proton Conductivity of Nafion 117 as Measured by a Four-Electrode AC Impedance Method. J. Electrochem. Soc. 1996, 143, 1254–1259. [Google Scholar] [CrossRef]
  33. Probst, N.; Grivei, E. Structure and Electrical Properties of Carbon Black. Carbon 2002, 40, 201–205. [Google Scholar] [CrossRef]
  34. Vetter, R.; Schumacher, J.O. Free Open Reference Implementation of a Two-Phase PEM Fuel Cell Model. Comput. Phys. Commun. 2019, 234, 223–234. [Google Scholar] [CrossRef]
  35. Lange, K.J.; Sui, P.-C.; Djilali, N. Pore Scale Simulation of Transport and Electrochemical Reactions in Reconstructed PEMFC Catalyst Layers. J. Electrochem. Soc. 2010, 157, B1434. [Google Scholar] [CrossRef]
  36. Zhang, J.; Yang, W.; Xu, L.; Wang, Y. Simulation of the Catalyst Layer in PEMFC Based on a Novel Two-Phase Lattice Model. Electrochim. Acta 2011, 56, 6912–6918. [Google Scholar] [CrossRef]
  37. Zhang, H.; Xiao, L.; Chuang, P.-Y.A.; Djilali, N.; Sui, P.-C. Coupled Stress–Strain and Transport in Proton Exchange Membrane Fuel Cell with Metallic Bipolar Plates. Appl. Energy 2019, 251, 113316. [Google Scholar] [CrossRef]
  38. Cai, Y.; Fang, Z.; Chen, B.; Yang, T.; Tu, Z. Numerical Study on a Novel 3D Cathode Flow Field and Evaluation Criteria for the PEM Fuel Cell Design. Energy 2018, 161, 28–37. [Google Scholar] [CrossRef]
  39. Parthasarathy, A.; Dav, B.; Srinivasan, S.; Appleby, A.J.; Martin, C.R. The Platinum Microelectrode/Nafion Interface: An Electrochemical ImpedanceSpectroscopicAnalysisof Oxygen Reduction Kinetics and Nation Characteristics. J. Electrochem. Soc. 1992, 139, 8. [Google Scholar]
  40. Wang, Y.; Wang, C.-Y. Dynamics of Polymer Electrolyte Fuel Cells Undergoing Load Changes. Electrochim. Acta 2006, 51, 3924–3933. [Google Scholar] [CrossRef]
Figure 1. CL microstructure reconstruction. (a) Pore phase; (b) solid phase.
Figure 1. CL microstructure reconstruction. (a) Pore phase; (b) solid phase.
Membranes 12 01001 g001
Figure 2. (a) Pore size distribution and (b) tortuosity distribution of CL.
Figure 2. (a) Pore size distribution and (b) tortuosity distribution of CL.
Membranes 12 01001 g002
Figure 3. The relationship between various variables and porosity: (a) relative diffusion coefficient, (b) effective oxygen diffusivity, (c) effective water vapor diffusivity, (d) effective nitrogen diffusivity, (e) the relationship between the effective electron conductivity and carbon volume fraction, and (f) the relationship between the effective proton conductivity and ionomer volume fraction.
Figure 3. The relationship between various variables and porosity: (a) relative diffusion coefficient, (b) effective oxygen diffusivity, (c) effective water vapor diffusivity, (d) effective nitrogen diffusivity, (e) the relationship between the effective electron conductivity and carbon volume fraction, and (f) the relationship between the effective proton conductivity and ionomer volume fraction.
Membranes 12 01001 g003aMembranes 12 01001 g003b
Figure 4. Geometric model of microstructure of CL.
Figure 4. Geometric model of microstructure of CL.
Membranes 12 01001 g004
Figure 5. Distribution of oxygen molar concentration at different moments: (a) 1 s, (b) 3 s, and (c) oxygen molar concentration along the thickness of CL at different moments.
Figure 5. Distribution of oxygen molar concentration at different moments: (a) 1 s, (b) 3 s, and (c) oxygen molar concentration along the thickness of CL at different moments.
Membranes 12 01001 g005aMembranes 12 01001 g005b
Figure 6. Distribution of current density in CL at different moments.
Figure 6. Distribution of current density in CL at different moments.
Membranes 12 01001 g006aMembranes 12 01001 g006b
Figure 7. (a) Profiles of temperature at CL/GDL interface, CL/PEM interface, and CL middle section at different moments. (b) Temperature distribution of different sections along the CL thickness direction when the reaction is stable.
Figure 7. (a) Profiles of temperature at CL/GDL interface, CL/PEM interface, and CL middle section at different moments. (b) Temperature distribution of different sections along the CL thickness direction when the reaction is stable.
Membranes 12 01001 g007
Figure 8. (a) The average liquid water saturation along the CL thickness direction at different moments. (b) Distribution of the average liquid water saturation on the y = 0 μm section when the reaction is stable. (c) Profiles of average ionomer water in the CL over time. (d) The iso-surface map of ionomer water along the thickness direction of the CL.
Figure 8. (a) The average liquid water saturation along the CL thickness direction at different moments. (b) Distribution of the average liquid water saturation on the y = 0 μm section when the reaction is stable. (c) Profiles of average ionomer water in the CL over time. (d) The iso-surface map of ionomer water along the thickness direction of the CL.
Membranes 12 01001 g008
Table 1. Transport parameters used in the model [32,33,34].
Table 1. Transport parameters used in the model [32,33,34].
ParametersValues
D O 2 b u l k (cm2/s)0.129
D H 2 O b u l k (cm2/s)0.166
D N 2 b u l k (cm2/s)0.161
σ i o n (S/cm)0.02
σ e l e (S/cm)10.0
Table 2. Model parameters [35,39].
Table 2. Model parameters [35,39].
ParametersValue
j c , r e f v (A/cm2)1.659 × 10−6
α 0.61
C O 2 , r e f (mol/cm3)40.96 × 10−6
γ 1.0
F (C/mol)96,487.0
R (J/(mol K))8.314
Table 3. Boundary conditions and initial values.
Table 3. Boundary conditions and initial values.
VariablesPEM SideGDL Side
Oxygen concentration (mol/cm3)/10.1 × 10−6
Liquid water saturation/0
Water vapor concentration (mol/cm3)/0
Proton potential (V)0.9080.9
Electron potential (V)1.7381.73
Temperature (K)353.0353.15
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhan, Z.; Song, H.; Yang, X.; Jiang, P.; Chen, R.; Harandi, H.B.; Zhang, H.; Pan, M. Microstructure Reconstruction and Multiphysics Dynamic Distribution Simulation of the Catalyst Layer in PEMFC. Membranes 2022, 12, 1001. https://doi.org/10.3390/membranes12101001

AMA Style

Zhan Z, Song H, Yang X, Jiang P, Chen R, Harandi HB, Zhang H, Pan M. Microstructure Reconstruction and Multiphysics Dynamic Distribution Simulation of the Catalyst Layer in PEMFC. Membranes. 2022; 12(10):1001. https://doi.org/10.3390/membranes12101001

Chicago/Turabian Style

Zhan, Zhigang, Hao Song, Xiaoxiang Yang, Panxing Jiang, Rui Chen, Hesam Bazargan Harandi, Heng Zhang, and Mu Pan. 2022. "Microstructure Reconstruction and Multiphysics Dynamic Distribution Simulation of the Catalyst Layer in PEMFC" Membranes 12, no. 10: 1001. https://doi.org/10.3390/membranes12101001

APA Style

Zhan, Z., Song, H., Yang, X., Jiang, P., Chen, R., Harandi, H. B., Zhang, H., & Pan, M. (2022). Microstructure Reconstruction and Multiphysics Dynamic Distribution Simulation of the Catalyst Layer in PEMFC. Membranes, 12(10), 1001. https://doi.org/10.3390/membranes12101001

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop