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

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.


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.

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.
Energies 2018, 11, 2529 where ε LJ is the depth of the potential well, σ LJ is the finite distance at which the inter-particle potential is zero, r LJ is the distance between the particles.

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.
where ε is the depth of the potential well, is the finite distance at which the inter-particle potential is zero, 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   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 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 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.
where ε is the depth of the potential well, is the finite distance at which the inter-particle potential is zero, 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   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 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 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.
where ε is the depth of the potential well, is the finite distance at which the inter-particle potential is zero, is the distance between the particles.Figure 2 shows the modeling of the CCL used in this research.In Figure 2a 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 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 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: 1.
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.

2.
If Pt spheres are coated on the surface of a carbon sphere, then this carbon sphere surface is marked as green.

3.
Since only the process of the oxygen reduction reaction is investigated here, dry-air is applied to the PEM fuel cell model.

4.
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.5.
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 Å.
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.

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 (H 3 O + ).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.

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: where, %Pt is the weight percentage of Pt in Pt-C, N Pt is the amount of Pt, m Pt 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.

Reaction Distance
The electrochemical reaction at the cathode catalyst layer is expressed as Equation ( 6): 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−Pt , and the reaction distance between oxygen and protons, R O−H , are determined by Equations ( 7) and ( 8): 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 Å.
Energies 2018, 11, x FOR PEER REVIEW 6 of 17 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, , and the reaction distance between oxygen and protons, , are determined by Equation ( 7) and Equation ( 8): 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.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.

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, 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 = 4Nq (10) 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: where A is the area of the cathode catalyst layer.

Potentials
The real operating voltage of a fuel cell could be written by subtracting the various overpotential losses from the thermodynamically predicted voltage: where η act , η ohmic , and η conc are activation overpotential loss, ohmic overpotential loss, and concentration overpotential loss, respectively.η act 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): where j 0 is in mA•cm −2 and the absolute temperature T is in K.
η conc 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: where j L is the limiting current density.
Energies 2018, 11, 2529 8 of 18 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 ASR m as:

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 scales with the thickness of the electrolyte should be as thin as possible.To verify the above model, the value of used is 0.47 Ω • cm , suggested by Marr et al. [30].The cell voltage curve of the present model drops sharply at first.The reason for this is that when current density is = 0 mA • cm 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 .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 Several simulations were performed with the above model.The area-specific resistance ASR m scales with the thickness of the electrolyte should be as thin as possible.To verify the above model, the value of ASR 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 Pt = 0.4 mg•cm −2 .The effective mass loading of Pt shown in Figure 5 is m Pt = 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 thermo .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.
Energies 2018, 11, x FOR PEER REVIEW 9 of 17 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.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 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  Each current density value corresponds to an output cell voltage value when using the KMC method.The non-zero minimum current density 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 ); carbon spheres without Pt (  Each current density value corresponds to an output cell voltage value when using the KMC method.The non-zero minimum current density 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  Each current density value corresponds to an output cell voltage value when using the KMC method.The non-zero minimum current density 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  Each current density value corresponds to an output cell voltage value when using the KMC method.The non-zero minimum current density 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 ).The units of length and radius are in Å.Each current density value corresponds to an output cell voltage value when using the KMC method.The non-zero minimum current density j min 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.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.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 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 = 0.120 mg•cm −2 , = 0.250 mg•cm −2 , = 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 = 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.During this movement duration, no oxygen reduction reaction occurs.Therefore, the current density will be zero and the output cell voltage will be E thermo 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.
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 = 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.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/cm 2 has low performance compared to 0.25 mg/cm 2 , because of the lower reaction surface area, A.
Figure 9 shows the reaction surface area per unit volume of the CCL ( ) and the fitting curve of catalyst surface area per unit mass of the catalyst ( ), both as a function of Pt weight percentage.When the weight percentage of Pt was in the range of 10-90%, 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, , 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, , is slightly decreased, and then both of the surfaces areas and 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/cm 2 has low performance compared to 0.25 mg/cm 2 , because of the lower reaction surface area, A.

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: Energies 2018, 11, 2529 16 of 18 1.
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.2.
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.

3.
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.

Figure 2
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.
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.(a) (b)

Figure 2 .
Figure 2. The modeling of catalyst layers: carbon spheres coated by Pt ( ); carbon spheres only ( ).All the units are in Å.(a) 100% carbon particles attached by Pt; (b) 90% carbon particles attached by Pt.

Figure 2 .
Figure 2. The modeling of catalyst layers: carbon spheres coated by Pt ( 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.(a) (b)

Figure 2 .
Figure 2. The modeling of catalyst layers: carbon spheres coated by Pt ( ); carbon spheres only ( ).All the units are in Å.(a) 100% carbon particles attached by Pt; (b) 90% carbon particles attached by Pt.

Figure 2 .
Figure 2. The modeling of catalyst layers: carbon spheres coated by Pt ( ); carbon spheres only ( ).All the units are in Å.(a) 100% carbon particles attached by Pt; (b) 90% carbon particles attached by Pt.

Figure 3 .
Figure 3. Electrochemical reaction processes inside the cathode catalyst layer.

2. 5 .
Current Density and Potentials 2.5.1.Current Density Current expresses the rate of charge transfer when electrons are consumed by electrochemical reactions at the cathode catalyst layer.According to the Faraday's law,

Figure 3 .
Figure 3. Electrochemical reaction processes inside the cathode catalyst layer.

Figure 4 .
Figure 4. Flow chart of the present developed kinetic Monte Carlo (KMC) model.
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 = 0.4 mg • cm .The effective mass loading of Pt shown in Figure 5 is = 0.4 mg • cm .Thus, the required number of Pt-C spheres is 1800.

Figure 4 .
Figure 4. Flow chart of the present developed kinetic Monte Carlo (KMC) model.

Figure 6 .
Figure 6.3-D snapshots of the microstructure of CCL: Pt-C spheres (

Figure 7 .
Figure 7.The comparison of present KMC model and analytical results.

Figure 7
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 Figure8a-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 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 Figure8a, 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 = 0.120 mg•cm −2 , = 0.250 mg•cm −2 , = 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 = 0.250 mg•cm −2 .The average current density is calculated to better observe the differences among the various cases as shown in Figure8a.

Figure 7 .
Figure 7.The comparison of present KMC model and analytical results.

Figure 7
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 O 2 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 thermo 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.

Figure 9 .
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 ) as a function of Pt weight percentage, (---); experimental data by E-TEK, (-*-).

Figure 9 .
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, (-*-).

Table 1 .
The diffusivities and random vector magnitudes of proton, oxygen and water molecules.

Table 2 .
Computation time for each case