Next Article in Journal
A Multidisciplinary Approach for the Development of Smart Distribution Networks
Next Article in Special Issue
Quantitative Comparison of Vernier Permanent-Magnet Motors with Interior Permanent-Magnet Motor for Hybrid Electric Vehicles
Previous Article in Journal
Estimation of the Relative Arrival Time of Microseismic Events Based on Phase-Only Correlation
Previous Article in Special Issue
Steering Stability Control for a Four Hub-Motor Independent-Drive Electric Vehicle with Varying Adhesion Coefficient
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Simulation of the Oxygen Reduction Reaction (ORR) Inside the Cathode Catalyst Layer (CCL) of Proton Exchange Membrane Fuel Cells Using the Kinetic Monte Carlo Method

Department of Mechanical Engineering, University of Nevada, Las Vegas, NV 89154, USA
*
Author to whom correspondence should be addressed.
Energies 2018, 11(10), 2529; https://doi.org/10.3390/en11102529
Submission received: 19 August 2018 / Revised: 15 September 2018 / Accepted: 20 September 2018 / Published: 21 September 2018
(This article belongs to the Collection Electric and Hybrid Vehicles Collection)

Abstract

:
In this paper, a numerical model of the kinetic Monte Carlo (KMC) method has been developed to study the oxygen reduction reaction (ORR) that occurs inside the cathode catalyst layer (CCL). Firstly, a 3-D model of the CCL that consists of Pt and carbon spheres is built using the sphere packing method; secondly, an efficient procedure of the proton-oxygen reaction process is developed and simulated. In the proton-oxygen reaction process, all of the continuous movements of protons and oxygen are considered. The maximum reaction distance is determined to be 8 Å. The input pressures of protons and oxygen are represented by the number of spheres of the species. The value of the current density is calculated based on the amount of reaction during the interval time. Indications are that the results of the present model match reasonably well with the published results. A new way to apply the KMC method in the proton exchange membrane fuel cell (PEMFC) research field is developed in this paper.

1. Introduction

Proton exchange membrane fuel cells (PEMFCs) are considered as one of the best fuel cells to provide power in the future, due to its having zero emissions, a rapid start-up time, and a high energy density at low operating temperatures [1]. PEMFCs have also attracted much attention as a potential replacement to traditional engines to power both mobile and stationary devices. The cost of the catalyst layer (CL) of PEM fuel cells is one of the crucial obstacles in fuel cell commercialization. Improving the platinum-use efficiency through optimizing the CL structure, and replacing the platinum with a cheap and highly reactive metal catalyst are the two technological ways that have been used to reduce the costs of the CL over the past decade [2,3,4]. Either technological method requires catalyst and void space in PEM fuel cells to provide continuous pathways for oxygen and protons to have oxygen reduction reaction (ORR) occurred at the cathode catalyst layer (CCL).
Studies of PEMFCs have been developed from 1-D to 3-D over the last two decades. In a review, Weber et al. [5] concluded the main methods for fuel cell research. Up to today, most of the fuel cell research has focused on macroscopic models. As the oxygen reduction reaction (ORR) occurs inside the catalyst layers, which are made of porous media, two important factors must be understood before further study. First, the porous media structure of the cathode catalyst layer needs to be understood; second, the behaviors of the ORR happening inside the layer need to be known [6].
Since the physical structure of the CCL is on a microscopic scale, numerical modeling has been found advantageous in practice to evaluate the effect of the CCL structure on PEMFC performance [7,8]. However, volumetric average values have been used in the numerical modeling to predict the performance of PEM fuel cells, which cannot explicitly explain the CCL structure [3,9,10]. Sphere packing and the kinetic Monte Carlo (KMC) method are ways to create a porous media structure of a CCL and to best simulate the reaction occurring inside the CCL. The KMC method can simulate geometries on a microscopic scale, and has been widely applied to the research of PEM fuel cells. Water distribution within a gas diffusion layer (GDL) was obtained from a Monte Carlo (MC) simulation by Seidenberger et al. [11,12]. Zhang et al. [13] studied the nanostructured catalyst layer and catalyst utilization in PEMFCs using the MC method. In this research, the authors constructed a lattice model for the catalyst layer to understand the composition-structure-performance relations in PEMFC. Each “cubic” element generated in the authors’ models represented an individual species, such as Pt, carbon, ionomer, and pore. Quiroga et al. [14] applied the MC method to investigate the electrochemical reactivity of materials in electrodes for a Pt-based PEM fuel cell. The authors reported on the oxygen transport at the microscopic interface between the active material and electrolyte. Gao et al. [3] investigated six CL structures using the MC method. The CL structures and their ability to diffuse oxygen determined the distribution of the dissolved oxygen. The differences of cell voltages among the six CL models were very small. Zhang et al. [15] constructed a 3-D microstructure of a solid oxide fuel cell electrode. Based on the researchers’ construction of the electrode, the authors estimated several properties of the solid oxide fuel cell, such as area specific resistance and thermal aging stability [15,16,17,18]. Most modern catalyst layer models are based on the fact that the oxygen needs to dissolve in the ionomer film first, before moving to reactive sites [19]. In the published articles, the researchers all discussed overall fuel cell performance, including the changes of output cell voltage and current density, as well as the location and distribution of water in a fuel cell. However, none of them clearly explained the electrochemical processes in the CCL, specifically, the motions of oxygen and protons.
In this paper, a systematic study on the PEMFC cathode catalyst layer was investigated using the kinetic Monte Carlo (KMC) method. With the consideration of the movements of oxygen and protons, the aims of this paper are to calculate the cell voltage and current density, and to investigate the PEM fuel cell performance. The model presented in this paper was validated by a comparison between the results obtained from the KMC model and available published simulation/experimental data.

2. Modeling and Theoretical Analysis

2.1. Model of the CCL

The PEMFC structure is schematically shown in Figure 1. Protons will be carried as hydronium through the membrane to the cathode catalyst layer (CCL), reacting with oxygen molecules and producing water as product.
In many papers, researchers have done experiments and used commercial computational fluid dynamics (CFD) software to study the performance of PEMFCs, which can show the results of the overall performance of fuel cells on the macro-scale. However, as the reaction happens within the catalyst layer of the cathode in a very short time, it is important to study the reaction on molecule-scale, in order to better understand the reaction behavior inside the catalyst layer. In previous work, Soontrapa and Chen [6] built a catalyst layer through the sphere random packing (SRP) method. SRP method is an arrangement of non-overlapping spheres within a containing space. The positions of spheres are generated randomly, while the distance between any two spheres is controlled by the Lennard-Jones (L-J) potential. The equation of Lennard-Jones potential is given in Equation (1). The spheres generated were mono-sized and represented carbon spheres. In this paper, the SRP method was used to generate a new CCL model with Pt coated on carbon spheres.
V L J = 4 ε L J [ ( σ L J r L J ) 12 ( σ L J r L J ) 6 ]
where ε L J is the depth of the potential well, σ L J is the finite distance at which the inter-particle potential is zero, r L J is the distance between the particles.
Figure 2 shows the modeling of the CCL used in this research. In Figure 2a, 100% of the carbon spheres are coated by Pt (the green spheres), and in Figure 2b only 90% of the carbon spheres are coated by Pt. All the self-developed codes used in this research were written and developed in Mathematica 10.1.
The van der Waals radii of carbon and Pt are 1.7 Å and 2.1 Å, respectively. If a CCL with a thickness of 10 µm is created and used, then over 2 × 10 23 carbon molecules need to be packed; this would require tremendous time to compute, which is not practical in numerical modeling and simulation. Therefore, the movements and reaction behaviors of proton and oxygen molecules were studied with the radius of 1.2 Å and 1.5 Å, respectively. If the size of the computational model is too large compared to the sizes of protons and oxygen, the amount of required proton and oxygen molecules would be huge, and this could cause a computing time issue.
To simplify the computation, five assumptions are made:
  • A rectangular container, shown in Figure 2, is created and used in this paper. The thickness, height, and width of the container are 50 Å, 100 Å, and 100 Å, respectively. Therefore, the total number of spheres can be reduced to 2000.
  • If Pt spheres are coated on the surface of a carbon sphere, then this carbon sphere surface is marked as green.
  • Since only the process of the oxygen reduction reaction is investigated here, dry-air is applied to the PEM fuel cell model.
  • The contact points where Pt is in contact with electrolytes, carbon, and gas are called the triple-phase boundary (TPB), and thus can be considered on the surface of Pt-carbon spheres, where the oxygen reduction reaction can occur.
  • The movements of electrons are not considered in the present study because of the following two reasons: (i) The electrons are moving very fast. Electrons move very fast only through the external wire or external solid conductive phase. (ii) The number of electrons in the external solid conductive phase is very high. It is very convenient to obtain electrons for the ORR.
There is no overlap between any two spheres. Because the sizes of the carbon and Pt spheres in a CCL are variable, the radii of the catalyst spheres in this paper follow the Gaussian distribution given in Equation (2). The minimum radius should be equal to the Van der Waals radius of carbon, which is 1.7 Å.
f ( r | r ¯ , σ 2 ) = 1 2 σ 2 π exp ( ( r r ¯ ) 2 2 σ 2 )
PEMFCs use hydrogen as fuel and oxygen as an oxidant. In this study, protons and oxygen are continuously supplied to the cathode catalyst layer from the interface (left side, as shown in Figure 2) of membrane and CCL, and the interface (right side, as shown in Figure 2) between CCL and GDL, respectively.

2.2. Theoretical Analysis

Under the STP condition, the volume of one mole of gas is about 22.4 L and consists of about 6.022 × 10 23 molecules. Therefore, about 27 molecules are in a cube with the volume of 1000   nm 3 . Considering the cost of PEMFC application, air is used instead of pure oxygen. Thus, about 6 oxygen molecules are in a 1000   nm 3 container of air at 1 atm, and 3, 9, 12, and 15 oxygen molecules are in the air at 0.5 atm, 1.5 atm, 2 atm, and 2.5 atm, respectively. Therefore, different inlet pressures of oxygen at the cathode gas flow channel can be considered in terms of the input number of oxygen molecules. Both oxygen and protons are supplied to the CCL continuously. Accordingly, the oxygen molecules and protons from the last iteration keep moving, while a certain number of new oxygen and proton molecules will enter the CCL per new iteration. For the spheres moving out of the model, the periodic boundary condition (PBC) was applied in this simulation. The number of total calculation iterations would be more than 100 million. To save calculation time, in each step, 24 protons were supplied to the CCL, and the simulation will run 200 times steps.
For the KMC method, the interaction energy between two particles is considered. This interaction energy can be calculated using simple pair potentials, such as the More or Lennard-Jones potentials [20]. In this study, the Lennard-Jones potential between each two molecules will be calculated and the reaction at the CCL can be determined.
The diffusivities of protons, oxygen and water for a Nafion 117 membrane are ( 0.6 ~ 5.8 ) × 10 6   cm 2   s 1 , ( 5.27 ~ 5.35 ) × 10 6   cm 2   s 1 , and ~ 2.65 × 10 6   cm 2   s 1 [21,22], respectively. In this simulation, the movements of proton, oxygen, and water molecules are based on the diffusivities and random vectors generated by the Monte Carlo technique. The values of diffusivities of proton/oxygen are 5.8 / 5.27 × 10 6   cm 2   s 1 in this study. The values are close because hydrogen loses electrons at the anode catalyst layer; then, they cross the membrane in the form of hydronium (H3O+). The size of hydronium is larger than protons and closer to oxygen. For simplification, in this paper protons were simulated instead of hydronium, which is reasonable for the electrochemical reaction occurring at the cathode side.
The magnitudes of the random vectors for proton, oxygen and water molecules are listed in Table 1. The magnitude of the random vector for oxygen is 50% of the van der Waals radius of oxygen. The reason for this is that the size of an oxygen molecule is very close to that of a carbon molecule. This moving algorithm has to make sure these oxygen molecules consecutively move through the CCL and do not jump through the void space. If this magnitude is too large then there is a chance that the oxygen spheres will not move forward following a proper path, but jump to the final location from the initial [23]. This is the same reason for choosing the vectors of proton and water spheres. The magnitude of a random vector of proton is 100,000 times the proton’s radius, as the diffusivities of protons and oxygen in the Nafion 117 membrane are very close, and the atomic radius of a proton is 0.8775 femtometer [24]. The value of the diffusion coefficient of water is about half that of oxygen. Therefore, the magnitude of the random vector of water is 25% of its radius.

2.3. Mass Loading of Pt

The mass loading of Pt is one of the key factors for PEM fuel cell performance. The amount of Pt can affect the amount of proton-oxygen reaction. Equations (3) through (5) are given to calculate the mass loading percentage or the number of Pt in the Pt-C particles:
% P t = N P t r ¯ P t 3 ρ P t N P t r ¯ P t 3 ρ P t + ( N N P t ) r ¯ C 3 ρ C
N P t = % P t   N   r ¯ C 3   ρ C r ¯ P t 3   ρ P t   ( 1 % P t ) + % P t   r ¯ C 3   ρ C
m P t = 4 π r ¯ P t 3 ρ P t N P t / 3 [ 0.5 ( 1 ε ) 3 / 2 a ] 2
where, % P t is the weight percentage of Pt in Pt-C, N P t is the amount of Pt, m P t is the effective mass loading of Pt, N is the total amount of Pt and C, r ¯ is the average radius, ρ is the density, a is the length of the model along y and z directions, ε is the porosity. The factor of 0.5 is applied to balance the fact that oxygen molecules would diffuse in or out of the catalyst layer randomly.

2.4. Reaction Distance

The electrochemical reaction at the cathode catalyst layer is expressed as Equation (6):
4 H + + 4 e + O 2 2 H 2 O
Figure 3 indicates the electrochemical reaction processes inside the CCL. As a matter of fact, the electrochemical reaction occurs on the triple phase boundary (TPB). The distance between protons and oxygen should be close enough to react. Water [25] suggested that the reaction distance should be 1 Å, for the bond length of hydrogen in a water molecule that is less than 1 Å. In addition, when the Lennard-Jones (L-J) potential of protons and oxygen reaches the minimum, this distance between protons and oxygen is about 3.2 Å. Equation (6) indicates that for each reaction, 4 protons and 1 oxygen molecule are needed. However, it is difficult to have 4 protons at a distance of 3.2 Å from the center of the oxygen molecule at the same time. To save computational time and satisfy the minimum image conventions, Smit [26] and Frenkel et al. [27] suggested that the L-J potential can be truncated at a cut-off distance of 2.5 σ , where σ is the distance at which the L-J potential between two particles is zero. Another reason for the reaction is the collision between oxygen molecules and Pt spheres. Therefore, as shown in Figure 3, the collision distance between oxygen and Pt, R O P t , and the reaction distance between oxygen and protons, R O H , are determined by Equations (7) and (8):
R O P t 1.05 ( r O 2 + r P t )
R O H 2.5   σ O H
The maximum reaction distance between oxygen and protons are determined by Equation (8) in the article, which is 8 Å. The electrochemical reaction will not occur when the distance between oxygen and protons is greater than 8 Å, because the Lennard-Jones potential will not be calculated at any distance greater than 8 Å.
Since the abundant electrons move very fast only through the external solid conductive phase, it is much easier to capture electrons around the reaction sites. Therefore, the movements of electrons are not considered in the present study.

2.5. Current Density and Potentials

2.5.1. Current Density

Current i expresses the rate of charge transfer when electrons are consumed by electrochemical reactions at the cathode catalyst layer. According to the Faraday’s law,
i = dQ / τ
where dQ is the charge in Coulomb (C) during the time interval τ in seconds (s). In the KMC method, the time interval, τ , can be calculated by the jump frequency, ω . Since the jump distance for oxygen is S O = 0.75 Å, and the oxygen mass transfer coefficient is K O = 2.80 × 10 4   cm · s 1 [22], it indicates that the real time interval has the value of τ = ω 1 = S O / K O 2.7 × 10 5   s . The following equation gives dQ as:
dQ = 4 N q
where q is the charge per electron, and N is the amount of electrochemical reactions occurring in time interval τ . The value of the reaction number N is monitored and obtained by simulation.
Thus, current density j can be calculated by the following equation:
j = i / A
where A is the area of the cathode catalyst layer.

2.5.2. Potentials

The real operating voltage of a fuel cell could be written by subtracting the various overpotential losses from the thermodynamically predicted voltage:
V = E t h e r m o η a c t η o h m i c η c o n c
where η a c t , η o h m i c , and η c o n c are activation overpotential loss, ohmic overpotential loss, and concentration overpotential loss, respectively. η a c t is sacrificed to overcome the activation barrier associated with the electrochemical reaction, and yields Equation (13). In Equation (13), j 0 is the exchange current density correlated with the experimental data of Parthasarathy et al. [28], with a formulation in Equation (14):
η a c t = R T α n F ln j 0 + R T α n F ln j
log 10 ( j 0 ) = 3.507 4001 T
where j 0 is in mA · cm 2 and the absolute temperature T is in K.
η c o n c is expressed in Equation (15) as shown below [29]. Equation (15) shows that concentration affects fuel cells’ Nernst voltage and reaction rates. The real reversible thermodynamic voltage and the reaction kinetics of a fuel cell are both determined by the reactant and product concentrations at the reaction sites:
η c o n c = ( R T n F ) ( 1 + 1 α ) ln j L j L j
where j L is the limiting current density.
The charges of protons do not transport through the membrane with a frictionless process. The cell voltage will have loss due to membrane resistance. The ohmic loss is defined from the current density and the area-specific resistance A S R m as:
η o h m i c = j   A S R m

3. Results and Discussion

The main flowchart of the present developed KMC model is shown in Figure 4.
Several simulations were performed with the above model. The area-specific resistance A S R m scales with the thickness of the electrolyte should be as thin as possible. To verify the above model, the value of A S R m used is 0.47 Ω · cm 2 , suggested by Marr et al. [30]. Figure 5 depicts the cell voltage curves among the results calculated by the present developed KMC model and the simulation model obtained by Marr et al. [30], with m P t = 0.4   mg · cm 2 . The effective mass loading of Pt shown in Figure 5 is m P t = 0.4   mg · cm 2 . Thus, the required number of Pt-C spheres is 1800. The cell voltage curve of the present model drops sharply at first. The reason for this is that when current density is j = 0   mA · cm 2 at the beginning, the oxygen reduction reaction (ORR) has not occurred at the cathode side. Thus, the fuel cell works as a capacitor for a very short time. Therefore, the cell voltage is equal to the theoretical reversible voltage or E t h e r m o . Once the ORR occurs at the cathode, the minimum amount of reaction is 1. This means that one oxygen and four protons appear within a reaction distance range during a specific time of τ , and one oxygen reduction reaction happens. Thus, the cell voltage drops down to the turning point immediately, due to various overpotential losses. Figure 5 shows that the results obtained from the present developed KMC model match reasonably well with the previously published results.
One of the interesting things about this research is computation time. This self-developed code was programmed in Mathematica, and the computation time is obtained and shown in Table 2 for each case.
Figure 6 shows snapshots of the processes of the ORR inside the microstructure of the CCL at various time steps for the case of 6 oxygen molecules supplying into the CCL per step. The total amount of Pt-C spheres and carbon spheres are both 1000. With a continuous supply of protons (red spheres) and oxygen (blue spheres), proton spheres keep diffusing into the CCL from the left side, while oxygen molecules move into CCL from the right side of the layer. Because of their smaller size and larger diffusivity, protons move faster than oxygen molecules in the CCL, and then at about 70 τ they start to encounter each other, at the right side of the CCL, about 20 Å from the oxygen-supplied side. Once four protons and one oxygen are within the reaction distance, the oxygen reduction reaction (ORR) occurs and two water molecules are generated. A few of the water molecules (yellow spheres) can be seen in Figure 6b. When the encountering of protons and oxygen increases, more water molecules are produced inside the CCL, as can be seen in Figure 6c,d.
Each current density value corresponds to an output cell voltage value when using the KMC method. The non-zero minimum current density j m i n can be calculated by Equation (9) through Equation (11) when N = 1. The j-V curve of the present developed KMC model is shown in Figure 7. Figure 7 shows that the output cell voltage decreases when the current density increases. The KMC results obtained from the present model are displayed as the “square”. The results obtained by the KMC method can be considered as discrete values. The reason for this is that the amounts of the electrochemical reaction occurred at each time step are not the same. The more reactions are occurred, the higher the current density is obtained. However, no matter what the amounts of reaction are for a certain case, the j-V curve could always be the same. However, the j-V curve is not a continuous curve, but shown as discrete points obtained from the KMC simulation, such as the “square” dot displayed in Figure 7. To validate the trace of the j-V curve obtained from the present developed KMC model, the analytical cell voltage is calculated using Equations (9) through (16) and the comparison of results is shown in Figure 7. The agreement of j-V curves is reasonably well matched.
Figure 7 only gives the overall relationship between the current density and the cell voltage, while Figure 8 shows that the cell voltage and current density change in the present developed KMC model per interval of time. The current densities and the cell voltages of the PEM fuel cells at various mass loadings of Pt are shown in Figure 8a–d, for the cases of 6 and 12 O2 diffusing in the CCL per interval time, respectively. Oxygen and protons must diffuse in the CCL before reacting. It takes about 65 τ –75 τ for protons and oxygen to move to the reaction site before reaching the reaction distance. During this movement duration, no oxygen reduction reaction occurs. Therefore, the current density will be zero and the output cell voltage will be E t h e r m o due to the capacitor mechanism. The amount of proton-oxygen reaction will be captured during each moving interval, and the values of the current density and output cell voltage jumped to corresponding points.
In Figure 8a, the red-squared line, blue-squared line, and the black-squared line represent the current density in the case of 6 O2 diffusing in the CCL per interval time with m P t = 0.120 mg·cm−2, m P t = 0.250 mg·cm−2, m P t = 0.350 mg·cm−2, respectively. The red-squared line jumped back to zero because no proton-oxygen reaction occurred. This is because the effective mass loading of Pt is low and does not have enough reaction sites for oxygen and protons to encounter. The number of proton-oxygen reactions increases when increasing the effective mass loading of Pt. However, the existing Pt will occupy the space and narrow down the movement pathway of protons and oxygen. Therefore, higher current density will appear in the case of m P t = 0.250 mg·cm−2. The average current density is calculated to better observe the differences among the various cases as shown in Figure 8a. It is clear that the higher current density will cause more overpotential losses. Even though the case of m P t = 0.250 mg·cm−2 has the highest current density, the average output cell voltage is the lowest, as shown in Figure 8b. The results shown in Figure 8c,d are similar to those shown in Figure 8a,b. The average current density and output cell voltage in the case of m P t = 0.250 mg·cm−2 are higher than those of the other two cases, but the output cell voltage difference between the blue line and black line is very small, only 0.003 V. Another important factor shown in Figure 8c,d is the amount of oxygen diffusing to the CCL. Increasing the inlet channel pressure would increase the amount of oxygen and increase the number of reactions.
The cell voltages displayed in Figure 8 show that the catalyst loading of 0.35 mg/cm2 has low performance compared to 0.25 mg/cm2. Another crucial factor for oxygen reduction reaction (ORR) is reaction surface area, A, which is determined by the Pt loading, the thickness of cathode catalyst layer (CCL), and the catalyst surface area per unit mass of catalyst, As. The catalyst surface area per unit mass of catalyst, As, is the experimental data provided by E-TEK [31], shown in the following Table 3. Equations (17) and (18) are used to calculate the reaction surface area, A:
A = A s · m P t / δ
m P t δ = ( 1 ε ) / ( 1 ρ P t + 1 % P t % P t ρ C )
where m P t is the platinum mass loading per unit area of the CCL, δ is the thickness of the CCL, % P t is the weight percentage of Pt in the total weight of solid spheres including Pt spheres and carbon spheres. ρ P t and ρ C are the densities of Pt and carbon, respectively. According to the experimental data shown in Table 3, A s can be described as the polynomial fitting function of %Pt:
A s = 73164.7 x 8 + 355613 x 7 719382 x 6 + 784759 x 5 499180 x 4 + 186732 x 3 39233.1 x 2 + 3883.54 x 1
where x represents %Pt.
Figure 9 shows the reaction surface area per unit volume of the CCL ( A ) and the fitting curve of catalyst surface area per unit mass of the catalyst ( A s ), both as a function of Pt weight percentage. When the weight percentage of Pt was in the range of 10–90%, A s decreased when the weight percentage of Pt in Pt-C is increased. It was found that the Pt at 40% of weight of Pt-C spheres is larger than those at the 20% and 10% of weight of Pt-C [32]. The size of Pt particles at the 40% Pt-C is almost twice as large as that of the 20% Pt-C, and large particle size corresponds to a smaller surface area [33]. However, the amount of Pt particles increased when the weight percentage of Pt is increased and then the reaction surface area, A , is increased to promoting the reaction occurred in CCL. With the accumulation of Pt spheres in the CCL when keeps increasing the weight percentage of Pt spheres, the reaction surface area, A , is slightly decreased, and then both of the surfaces areas A and A s are increased because of the increased amount of Pt spheres. This has caused that the amount of total ORR is slightly decreased to around 90% weight percentage of Pt in Pt-C. Therefore, catalyst loading of 0.35 mg/cm2 has low performance compared to 0.25 mg/cm2, because of the lower reaction surface area, A.

4. Conclusions

Enormous efforts have been made over the past few years to figure out the electrochemical reaction processes in PEM fuel cells. Since the dimensions of the cathode catalyst layer of a PEM fuel cell are at a microscopic scale, the electrochemical reaction processes are very difficult to measure. Traditional experimental methods and numerical methods, such as using CFD software or partial differential equations (PDEs) or ordinary differential equations (ODEs), can easily draw out the performance of PEM fuel cells in macroscopic form. However, the transport processes of oxygen and protons, the oxygen reduction reaction processes, and other electrochemical reaction processes involved still need to be focused on, especially at the microscopic scale. The cathode catalyst layer of PEM fuel cells, as the thinnest layer located between the gas diffusion layer and membrane, is always neglected when applying experimental methods or traditional numerical methods. The Monte Carlo method has been proven to be one of the most powerful and reliable methods on the microscopic scale. The MC method builds up a bridge connecting different scale sizes and brings possibility to the microscopic world.
This study proposes a numerical model to estimate microstructure features and PEMFC performance. The comparison between the results obtained from the present developed kinetic Monte Carlo (KMC) model and available published works and data lends support to the reliability of this developed KMC model in simulating and studying oxygen reduction reactions inside of a cathode catalyst layer. The main conclusions are summarized as follows:
  • Protons and oxygen are continuously supplied to the CCL. At each interval time, τ , proton and oxygen spheres diffuse into the CCL through interfaces between membranes and the CCL, and between the CCL and GDL, respectively.
  • The maximum reaction distance is 8 Å, determined by the Lennard-Jones potential. This reaction distance ensures that one oxygen sphere and four proton spheres will meet and the ORR will occur.
  • The current density is calculated by the amount of ORR during each interval time, τ . The amount of ORR per interval time is monitored; thus, the total charge transferred is obtained to calculate the current density.
The procedures in this paper explicitly give a new evaluation method when using the KMC method to obtain more detailed information about the performance of the CCL in PEMFCs. Our goal in fuel cell research is to fully understand how the performance of PEM fuel cell affected by water generated in the CCL. This study shows the procedures to simulate that water molecules are generated. However, the water cannot be fully removed from the PEM fuel cell. At low temperature, ice will form from the remaining water. The ice formation using the KMC method is currently under developed and studied to better understand how to control water management in PEM fuel cell research. The present developed KMC method and results will be useful for further fuel cell research.

Author Contributions

Conceptualization, Y.-T.C. and B.B.; Methodology, Y.-T.C. and B.B.; Software, B.B.; Validation, B.B.; Analytical Analysis, B.B.; Investigation, Y.-T.C. and B.B.; Resources, Y.-T.C.; Data Curation, B.B.; Writing—Original Draft Preparation, B.B.; Writing—Review & Editing, Y.-T.C.; Visualization, B.B.; Supervision, Y.-T.C.; Project Administration, Y.-T.C.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclatures

Areaction surface area per unit volume of the CCL, cm2/cm3
Ascatalyst surface area per unit mass of the catalyst, cm2/g
a The length along y and z direction, cm
D i j Diffusivity, cm 2 · s 1
m P t The effective mass loading of Pt, mg · cm 2
NThe number of particles
r The radius of the particle, cm
δ The thickness of the CCL, cm
ε The porosity
ρ The density, g · cm 3
Subscripts
c carbon
o oxygen
p proton
Pt platinum
w water

References

  1. Ji, M.; Wei, Z. A review of water management in polymer electrolyte membrane fuel cells. Energies 2009, 2, 1057–1106. [Google Scholar] [CrossRef]
  2. Park, J.C.; Park, S.H.; Chung, M.W.; Choi, C.H.; Kho, B.K.; Woo, S.I. Optimization of catalyst layer composition for PEMFC using grapheme-based oxygen reduction reaction catalyst. J. Power Sources 2015, 286, 166–174. [Google Scholar] [CrossRef]
  3. Gao, Y.; Zhang, X. Geometrical structures of catalyst layer and their impact on oxygen reduction in proton exchange membrane fuel cell. Electrochim. Acta 2016, 218, 101–109. [Google Scholar] [CrossRef]
  4. Zhu, S.; Wang, S.; Jiang, L.; Xia, Z.; Sun, H.; Sun, G. High Pt utilization catalyst prepared by ion exchange method for direct methanol fuel cells. Int. J. Hydrogen Energy 2012, 37, 14543–14548. [Google Scholar] [CrossRef] [Green Version]
  5. Weber, A.Z.; Newman, J. Modeling Transport in Polymer-Electrolyte Fuel Cells. Chem. Rev. 2004, 104, 4679–4726. [Google Scholar] [CrossRef] [PubMed]
  6. Soontrapa, K.; Chen, Y. Mono-sized sphere packing algorithm development using optimized Monte Carlo technique. Adv. Powder Technol. 2013, 24, 955–961. [Google Scholar] [CrossRef]
  7. Owejan, J.P.; Owejan, J.E.; Gu, W.B. Impact of Platinum loading and Catalyst Layer Structure on PEMFC Performance. J. Electrochem. Soc. 2012, 160, F824–F833. [Google Scholar] [CrossRef]
  8. Yoon, W.; Weber, A.Z. Modeling Low-Platinum-Loading Effects in Fuel-Cell Catalyst Layers. J. Electrochem. Soc. 2011, 158, B1007–B1018. [Google Scholar] [CrossRef]
  9. Moein-Jahromi, M.; Kermani, M.J. Performance prediction of PEM fuel cell cathode catalyst layer using agglomerate model. Int. J. Hydrogen Energy 2012, 37, 17954–17966. [Google Scholar] [CrossRef]
  10. Weber, A.Z.; Borup, R.L.; Darling, R.M.; Das, P.K.; Dursch, T.J.; Gu, W.; Harvey, D.; Kusoglu, A.; Litster, S.; Mench, M.M.; et al. A Critical Review of Modeling Transport Phenomena in Polymer-Electrolyte Fuel Cells. J. Electrochem. Soc. 2014, 161, F1254–F1299. [Google Scholar] [CrossRef] [Green Version]
  11. Seidenberger, K.; Wilhelm, F.; Schmitt, T.; Lehnert, W.; Scholta, J. Estimation of water distribution and degradation mechanisms in polymer electrolyte membrane fuel cell gas diffusion layers using a 3D Monte Carlo model. J. Power Sources 2011, 196, 5317–5324. [Google Scholar] [CrossRef]
  12. Seidenberger, K.; Wilhelm, F.; Haußmann, J.; Markötter, H.; Manke, I.; Scholta, J. Grand canonical Monte Carlo study on water agglomerations within a polymer electrolyte membrane fuel cell gas diffusion layer. J. Power Sources 2013, 239, 628–641. [Google Scholar] [CrossRef]
  13. Zhang, J.; Cao, P.; Xu, L.; Wang, Y. Modeling nanostructured catalyst layer in PEMFC and catalyst utilization. Front. Chem. Eng. China 2011, 5, 297–302. [Google Scholar] [CrossRef]
  14. Quiroga, M.A.; Franco, A.A. A Multi-Paradigm Computational Model of Materials Electrochemical Reactivity for Energy Conversion and Storage. J. Electrochem. Soc. 2015, 162, E73–E83. [Google Scholar] [CrossRef] [Green Version]
  15. Zhang, Y.; Ni, M.; Yan, M.; Chen, F. Thermal aging stability of infiltrated solid oxide fuel cell electrode microstructures: A three-dimensional kinetic Monte Carlo simulation. J. Power Sources 2015, 299, 578–586. [Google Scholar] [CrossRef]
  16. Zhang, Y.; Wang, Y.; Wang, Y.; Chen, F.; Xia, C. Random-packing model for solid oxide fuel cell electrodes with particle size distributions. J. Power Sources 2011, 196, 1983–1991. [Google Scholar] [CrossRef]
  17. Zhang, Y.; Sun, Q.; Xia, C.; Ni, M. Geometric Properties of Nanostructured Solid Oxide Fuel Cell Electrodes. J. Electrochem. Soc. 2013, 160, F278–F289. [Google Scholar] [CrossRef] [Green Version]
  18. Zhang, Y.; Ni, M.; Xia, C. Microstructural Insights into Dual-Phase Infiltrated Solid Oxide Fuel Cell Electrodes. J. Electrochem. Soc. 2013, 160, F834–F839. [Google Scholar] [CrossRef]
  19. 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]
  20. Claassens, C.H.; Terblans, J.J.; Hoffman, M.H.; Swart, H.C. Kinetic Monte Carlo simulation of monolayer gold film growth on a graphite substrate. Surf. Interface Anal. 2005, 37, 1021–1026. [Google Scholar] [CrossRef]
  21. Zawodzinski, T.A.; Neeman, M.; Sillerud, L.O.; Gottesfeld, S. Determination of water diffusion coefficients in perfluorosulfonate ionomeric membranes. J. Phys. Chem. 1991, 95, 6040–6044. [Google Scholar] [CrossRef]
  22. Chae, K.; Choi, M.; Ajayi, F.; Park, W. Mass transport through a proton exchange membrane (nafion) in microbial fuel cells. Energy Fuels 2008, 22, 169–176. [Google Scholar] [CrossRef]
  23. Soontrapa, K. Study of Water Transport Phenomena on Cathode of PEMFCs Using Monte Carlo Simulation. Ph.D. Thesis, University of Nevada, Las Vegas, NV, USA, 2014. [Google Scholar]
  24. Mohr, P.; Taylor, B.; Newell, D. CODATA recommended values of the fundamental physical constants: 2010. Rev. Mod. Phys. 2012, 84, 1527–1605. [Google Scholar] [CrossRef] [Green Version]
  25. Franks, F. Water: A Matrix of Life; The Royal Society of Chemistry: London, UK, 2000. [Google Scholar]
  26. Smit, B. Phase diagrams of Lennard-Jones fluids. J. Chem. Phys. 1992, 96, 8639–8640. [Google Scholar] [CrossRef] [Green Version]
  27. Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: London, UK, 2002; Volume 638. [Google Scholar]
  28. Parthasarathy, A.; Srinivasan, S.; Appleby, J. Temperature dependence of the electrode kinetics of oxygen reduction at the platinum/Nafion interface—A microelectrode investigation. J. Electrochem. Soc. 1992, 139, 2530–2537. [Google Scholar] [CrossRef]
  29. O’Hayre, R.; Cha, S.W.; Colella, W.; Prinz, F.B. Chapter 5: Fuel Cell Mass Transport. In Fuel Cell Fundamentals; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2016; pp. 167–202. [Google Scholar]
  30. Marr, C.; Li, X. Composition and performance modelling of catalyst layer in a proton exchange membrane fuel cell. J. Power Sources 1999, 77, 17–27. [Google Scholar] [CrossRef]
  31. K-TEK. Gas Diffusion Electrodes and Catalyst Materials; Catalogue: Natick, MA, USA, 1995. [Google Scholar]
  32. Ticianelli, E.A.; Derouin, C.R.; Redondo, A.; Srinivasan, S. Methods of advance technology of proton exchange membrane fuel cells. J. Electrochem. Soc. 1988, 135, 2209–2214. [Google Scholar] [CrossRef]
  33. Qi, Z.; Kaufman, A. Low Pt loading high performance cathodes for PEM fuel cells. J. Power Sources 2003, 113, 37–43. [Google Scholar] [CrossRef]
Figure 1. Schematic structure of proton exchange membrane fuel cell (PEMFC).
Figure 1. Schematic structure of proton exchange membrane fuel cell (PEMFC).
Energies 11 02529 g001
Figure 2. The modeling of catalyst layers: carbon spheres coated by Pt ( Energies 11 02529 i001); carbon spheres only ( Energies 11 02529 i002). All the units are in Å. (a) 100% carbon particles attached by Pt; (b) 90% carbon particles attached by Pt.
Figure 2. The modeling of catalyst layers: carbon spheres coated by Pt ( Energies 11 02529 i001); carbon spheres only ( Energies 11 02529 i002). All the units are in Å. (a) 100% carbon particles attached by Pt; (b) 90% carbon particles attached by Pt.
Energies 11 02529 g002
Figure 3. Electrochemical reaction processes inside the cathode catalyst layer.
Figure 3. Electrochemical reaction processes inside the cathode catalyst layer.
Energies 11 02529 g003
Figure 4. Flow chart of the present developed kinetic Monte Carlo (KMC) model.
Figure 4. Flow chart of the present developed kinetic Monte Carlo (KMC) model.
Energies 11 02529 g004
Figure 5. Comparisons of cell voltage curves between the present developed KMC model simulation and the published data. m P t = 0.4   mg · cm 2 , A S R m = 0.47   Ω · cm 2 [30].
Figure 5. Comparisons of cell voltage curves between the present developed KMC model simulation and the published data. m P t = 0.4   mg · cm 2 , A S R m = 0.47   Ω · cm 2 [30].
Energies 11 02529 g005
Figure 6. 3-D snapshots of the microstructure of CCL: Pt-C spheres ( Energies 11 02529 i003); carbon spheres without Pt ( Energies 11 02529 i004); oxygen spheres ( Energies 11 02529 i005); proton spheres ( Energies 11 02529 i006); water spheres ( Energies 11 02529 i007). The units of length and radius are in Å.
Figure 6. 3-D snapshots of the microstructure of CCL: Pt-C spheres ( Energies 11 02529 i003); carbon spheres without Pt ( Energies 11 02529 i004); oxygen spheres ( Energies 11 02529 i005); proton spheres ( Energies 11 02529 i006); water spheres ( Energies 11 02529 i007). The units of length and radius are in Å.
Energies 11 02529 g006
Figure 7. The comparison of present KMC model and analytical results.
Figure 7. The comparison of present KMC model and analytical results.
Energies 11 02529 g007
Figure 8. Current densities and cell voltages: (a,b) for the case of 6 O2 diffusing into cathode catalyst layer (CCL) per time step; (c,d) for the case of 12 O2 diffusing into CCL per time step.
Figure 8. Current densities and cell voltages: (a,b) for the case of 6 O2 diffusing into cathode catalyst layer (CCL) per time step; (c,d) for the case of 12 O2 diffusing into CCL per time step.
Energies 11 02529 g008aEnergies 11 02529 g008b
Figure 9. Reaction surface area per unit volume of the CCL ( A ) as a function of Pt weight percentage, (---); fitting curve of catalyst surface area per unit mass of the catalyst ( A s ) as a function of Pt weight percentage, (---); experimental data by E-TEK, (-*-).
Figure 9. Reaction surface area per unit volume of the CCL ( A ) as a function of Pt weight percentage, (---); fitting curve of catalyst surface area per unit mass of the catalyst ( A s ) as a function of Pt weight percentage, (---); experimental data by E-TEK, (-*-).
Energies 11 02529 g009
Table 1. The diffusivities and random vector magnitudes of proton, oxygen and water molecules.
Table 1. The diffusivities and random vector magnitudes of proton, oxygen and water molecules.
Species D i j ( c m 2 / s ) Vector FormulaVector Magnitude (Å)
H + 5.8 × 10 6 100,000 r p 0.88
O 2 5.27 × 10 6 0.5 r o 0.75
H 2 O 2.65 × 10 6 0.25 r w 0.35
Table 2. Computation time for each case
Table 2. Computation time for each case
m P t   ( m g · c m 2 ) Number of Oxygen Supplying into CCL Per Step
612
0.1231 h45 h
0.2537 h50 h
0.3542 h57 h
0.4049 h66 h
0.5061 h77 h
Table 3. Representative catalyst surface areas [31].
Table 3. Representative catalyst surface areas [31].
Percentage of Pt on Carbon Black (%Pt) A s   ( m 2 / g )
10%140
20%112
30%88
40%72
60%32
80%11
100%28

Share and Cite

MDPI and ACS Style

Bai, B.; Chen, Y.-T. Simulation of the Oxygen Reduction Reaction (ORR) Inside the Cathode Catalyst Layer (CCL) of Proton Exchange Membrane Fuel Cells Using the Kinetic Monte Carlo Method. Energies 2018, 11, 2529. https://doi.org/10.3390/en11102529

AMA Style

Bai B, Chen Y-T. Simulation of the Oxygen Reduction Reaction (ORR) Inside the Cathode Catalyst Layer (CCL) of Proton Exchange Membrane Fuel Cells Using the Kinetic Monte Carlo Method. Energies. 2018; 11(10):2529. https://doi.org/10.3390/en11102529

Chicago/Turabian Style

Bai, Baosheng, and Yi-Tung Chen. 2018. "Simulation of the Oxygen Reduction Reaction (ORR) Inside the Cathode Catalyst Layer (CCL) of Proton Exchange Membrane Fuel Cells Using the Kinetic Monte Carlo Method" Energies 11, no. 10: 2529. https://doi.org/10.3390/en11102529

APA Style

Bai, B., & Chen, Y. -T. (2018). Simulation of the Oxygen Reduction Reaction (ORR) Inside the Cathode Catalyst Layer (CCL) of Proton Exchange Membrane Fuel Cells Using the Kinetic Monte Carlo Method. Energies, 11(10), 2529. https://doi.org/10.3390/en11102529

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