A Molecular Model of PEMFC Catalyst Layer: Simulation on Reactant Transport and Thermal Conduction

Minimizing platinum (Pt) loading while reserving high reaction efficiency in the catalyst layer (CL) has been confirmed as one of the key issues in improving the performance and application of proton exchange membrane fuel cells (PEMFCs). To enhance the reaction efficiency of Pt catalyst in CL, the interfacial interactions in the three-phase interface, i.e., carbon, Pt, and ionomer should be first clarified. In this study, a molecular model containing carbon, Pt, and ionomer compositions is built and the radial distribution functions (RDFs), diffusion coefficient, water cluster morphology, and thermal conductivity are investigated after the equilibrium molecular dynamics (MD) and nonequilibrium MD simulations. The results indicate that increasing water content improves water aggregation and cluster interconnection, both of which benefit the transport of oxygen and proton in the CL. The growing amount of ionomer promotes proton transport but generates additional resistance to oxygen. Both the increase of water and ionomer improve the thermal conductivity of the C. The above-mentioned findings are expected to help design catalyst layers with optimized Pt content and enhanced reaction efficiency, and further improve the performance of PEMFCs.


Introduction
Proton exchange membrane fuel cells (PEMFCs) have attracted attention from energy devices such as portable, mobile, and stationary devices because it helps effective reductions of energy shortage and environment pollution [1][2][3][4]. However, the high price of PEMFCs components, especially the platinum (Pt) group catalyst becomes one of the bottlenecks that limit their commercial development. Numerous strategies and approaches have been tried to reduce the catalyst platinum loading, which is still hard to meet the Department of Energy (DOE) target of a cathode Pt loading of 0.1 mg/cm 2 [5]. In pursuit of reducing platinum loading and raising power efficiency, many scholars have found that the resistance of oxygen transmission is becoming particularly obvious under low platinum loading and high current density, which greatly limits the oxygen reduction reaction (ORR) and hence reduces the current density. At present, the studies about reactant transport and reaction efficiency in CL are in the primary stage. Gasteiger et al. [6] quantified the activation loss and voltage loss and proposed two theoretical methods to test the performance of catalysts. Some of the tested non-Pt catalysts have a 2.5-4 times enhancement of mass activity. Wonseok et al. [7] studied the influence of a layer of a thin film on the surface of catalysts on oxygen transport and battery performance using an agglomerate model. The review of Weber et al. [8] summarized the sharp resistance increase when Pt capacity was less than 0.1 mg/cm 2 from the experimental and numerical model results. In general, the local resistance of oxygen is the sum of many parts when thronging the multi-components of the catalyst layer. However, the mathematical model and experimental measurement method for the quantitative and qualitative analyses have been not unified due to the different morphologies, structures, and complex micro-interaction [9][10][11][12].
As for numerical simulation research, analytical models [13,14], mesoscopic methods such as lattice Boltzmann method (LBM) [15,16], microscopic methods such as molecular dynamics (MD) [17][18][19][20][21][22][23][24][25][26][27] simulation, and density functional theory (DFT) [28][29][30][31] are usually used to construct the CL structure and study the mechanisms of reactant transport, reaction, and heat conduction. For example, Liang et al. [14] proposed the fractal theory of porous media to quantify the effective electrolyte diffusivity in porous media with consideration of the electrical double layer (EDL) effects. Chen et al. [15,16] proposed a "watermelon model" containing primary and secondary pores to simulate the transport under different conditions. Rao et al. [18,19] established the molecular dynamics model of polymer membrane to study the proton and thermal conductivity. Feng et al. investigated the thermal-mechanical properties in [20] and gave the thermal expand coefficients with different Pt-C combination styles. She also studied the multicomponent interaction and mechanical capacity in [21] and obtained the stress-strain characteristics with various water contents. Fan et al. [22] expounded that over 90% of O 2 transportation occurred at the upper corner and edge regions rather than the faces of Pt because of the different arrangement tightness of perfluorosulfonic. He also studied the effects of side chain length on the thermal conductivity of the perfluorosulfonic acid (PFSA) membrane in another paper [26]. Zhao et al. [27] used the oxidized graphite to design a three-phase microenvironment in CL by MD simulation and experiments. Taeyoon et al. [29] calculated the adsorption characteristics and electron transfer of the catalysts with optimized geometry using the density functional theory (DFT) with carbon nanotube added in the CL. However, there are some deficiencies in the model establishment and parameter calculation based on the above research analysis. Firstly, some micro-structure are not built corresponding to the real components such as the carbon pores. Secondly, the parameter study is not so comprehensive from the molecular-level. Moreover, there is a paucity of research about the thermal conductivity in the CL using MD simulations. The MD simulation method make the parameter study enable profiting from the modifiable number of molecules and state of operation. For the catalyst layer, the parameter study contains multiple structural parameters and operating parameters. In this study, two typical parameters-water content and ionomer-to-carbon ratios-are used to investigate their effect on transport and thermal conduction in CL. In this study, a molecular model containing the multicomponent cathode catalyst layer is built and the corresponding equilibrium MD and nonequilibrium MD simulations with different water contents and ionomer-to-carbon ratios are conducted. The radial distribution functions (RDFs), mean square displacements (MSDs), and water contour distribution are calculated to illuminate the relationships between each component. The thermal conduction characteristics are obtained via the given heat flux and temperature gradients. This work is expected to help understand the multicomponent nanostructures, transport properties, and thermal features of the catalyst layer in PEMFCs.

Model and Dynamics Simulation Details
The cathode CL is composed of carbon-supported Pt particles, water, hydronium (proton), oxygen, and ionomer. Figure 1 shows the schematic construction of the PEMFC catalyst layer between the proton exchange membrane (PEM) and gas diffusion layer (GDL). In the CL the interactive multicomponent mixture gathers around the catalysts after oxygen transport through the GDL, and hydronium ions through the membrane. After a complex transport process at the three-phase interface of the reactants, the reaction takes place near the catalyst reactive place and produces water. In order to have a comprehensive reflection of the complex microstructures and the multi-components reactions in the catalyst layer, a molecular model of CL is built as shown in Figure 2. The diameter of carbon particles is 5 nm, and each of them contains 8000 carbon atoms. The pores inside the carbon particles, which are the so-called primary pores and the secondary pores, have been set to make the results closer to reality. A total of 4-6 Pt particles with 1.2 nm diameter are dispersed Membranes 2021, 11, 148 3 of 15 randomly on the carbon particle surface and inner pore. The mass fraction of Pt is 25 wt% in the Pt-C particle system. in the catalyst layer, a molecular model of CL is built as shown in Figure 2. The diameter of carbon particles is 5 nm, and each of them contains 8000 carbon atoms. The pores inside the carbon particles, which are the so-called primary pores and the secondary pores, have been set to make the results closer to reality. A total of 4-6 Pt particles with 1.2 nm diameter are dispersed randomly on the carbon particle surface and inner pore. The mass fraction of Pt is 25 wt% in the Pt-C particle system.  In addition to the Pt-C particles, the ionomer is an important part of CL owing to its proton conduction effects. Although some new short-chain polymers have been used in the catalyst layer such as adding aromatic groups. In addition, it is noted that some properties would change like ionic conductivity, oxygen resistance, and thermal conductivity with different polymers. This study focuses on the reactant transport and thermal conduction at the specific polymer. Hence, Nafion polymer, which is the most commonly used polymer in an experimental and commercial application, is adopted in this work. The chemical structure of the ionomer monomer from the Nafion polymer is demonstrated in Figure 3. In this study, m = 1, n = 2, x = 6, and y = 1, and each ionomer particle in the catalyst layer, a molecular model of CL is built as shown in Figure 2. The diameter of carbon particles is 5 nm, and each of them contains 8000 carbon atoms. The pores inside the carbon particles, which are the so-called primary pores and the secondary pores, have been set to make the results closer to reality. A total of 4-6 Pt particles with 1.2 nm diameter are dispersed randomly on the carbon particle surface and inner pore. The mass fraction of Pt is 25 wt% in the Pt-C particle system.  In addition to the Pt-C particles, the ionomer is an important part of CL owing to its proton conduction effects. Although some new short-chain polymers have been used in the catalyst layer such as adding aromatic groups. In addition, it is noted that some properties would change like ionic conductivity, oxygen resistance, and thermal conductivity with different polymers. This study focuses on the reactant transport and thermal conduction at the specific polymer. Hence, Nafion polymer, which is the most commonly used polymer in an experimental and commercial application, is adopted in this work. The chemical structure of the ionomer monomer from the Nafion polymer is demonstrated in Figure 3. In this study, m = 1, n = 2, x = 6, and y = 1, and each ionomer particle In addition to the Pt-C particles, the ionomer is an important part of CL owing to its proton conduction effects. Although some new short-chain polymers have been used in the catalyst layer such as adding aromatic groups. In addition, it is noted that some properties would change like ionic conductivity, oxygen resistance, and thermal conductivity with different polymers. This study focuses on the reactant transport and thermal conduction at the specific polymer. Hence, Nafion polymer, which is the most commonly used polymer in an experimental and commercial application, is adopted in this work. The chemical structure of the ionomer monomer from the Nafion polymer is demonstrated in Figure 3. In this study, m = 1, n = 2, x = 6, and y = 1, and each ionomer particle contains 10 monomers. To investigate the influence of water content and ionomer quantities on the reactants transport and thermal conductivity, the various water contents are set as 1, 4, 7, and 10, and ionomer-to-carbon ratios are 0.4, 0.8, 1.2, and 1.6. As for the determination of I/C ratios, the commercial catalyst ionomer content is about 10-30%. In addition, this content can reach 50% to the maximum in some reported experimental studies. After unit conversion, the I/C ratios of 0.4, 0.8, 1.2, and 1.6 used in this paper are 24%, 39%, 49%, and 56% ionomer content. The I/C = 0.4-1.2 are within the scope of the commercial and experimental application while preparing the catalyst layer. Moreover, in order to obtain the influence rules with a slightly excessive ionomer, the I/C= 1.6 are also set as the simulation cases. The detailed molecule numbers of each component in the simulations are listed in Tables 1 and 2. contains 10 monomers. To investigate the influence of water content and ionomer quantities on the reactants′ transport and thermal conductivity, the various water contents are set as 1, 4, 7, and 10, and ionomer-to-carbon ratios are 0.4, 0.8, 1.2, and 1.6. As for the determination of I/C ratios, the commercial catalyst ionomer content is about 10-30%. In addition, this content can reach 50% to the maximum in some reported experimental studies. After unit conversion, the I/C ratios of 0.4, 0.8, 1.2, and 1.6 used in this paper are 24%, 39%, 49%, and 56% ionomer content. The I/C = 0.4-1.2 are within the scope of the commercial and experimental application while preparing the catalyst layer. Moreover, in order to obtain the influence rules with a slightly excessive ionomer, the I/C= 1.6 are also set as the simulation cases. The detailed molecule numbers of each component in the simulations are listed in Tables 1 and 2.    NIonomer  10  20  30  40  NPt-C  3  3  3  3  NO2  200  200  200  200  NH2O  1200  1200  1200  1200  NH3O+  200  200  200  200 All simulations in this work are conducted using Materials Studio 2017 molecular modeling software (Accelrys, San Diego, CA, USA), and the COMPASS force field is used, which can calculate uniformly the organic molecular system and the inorganic molecular system covering all the adopted molecules in the simulation.
In the COMPASS force field, total energy between the atoms contain the valance energy and the non-bond energy [32], which can be expressed as follows: The valance energy contained the bond item, angle item, torsion item, and out-ofplane vibration item in Equation (2). The non-bond energy contained the Van der Waals item and the Columbia item in Equation (3).  All simulations in this work are conducted using Materials Studio 2017 molecular modeling software (Accelrys, San Diego, CA, USA), and the COMPASS force field is used, which can calculate uniformly the organic molecular system and the inorganic molecular system covering all the adopted molecules in the simulation.
In the COMPASS force field, total energy between the atoms contain the valance energy and the non-bond energy [32], which can be expressed as follows: The valance energy contained the bond item, angle item, torsion item, and out-ofplane vibration item in Equation (2). The non-bond energy contained the Van der Waals item and the Columbia item in Equation (3).
where b, θ, and Φ represent the bond length, bond angle, and dihedral angle, respectively, and χ is the heterogeneous shape parameter.

Radial Distribution Function
To find the interaction between the multicomponent mixture in the CL, the radial distribution function (RDF) is calculated by Equation (6) as follows: where V represents the volume of the simulation system volume and n B is the number of atom B located in a spherical shell centered on atom A. The shell's volume is 4πr 2 dr. N B is the number of B atoms in the system.

Mean Square Displacement and Diffusion Coefficient
The mean square displacement (MSD) can be expressed as Equation (7).
In addition, the diffusion coefficient D could be obtained by the Einstein formula as follows:

Thermal Conductivity
The thermal conductivity is calculated using the non-equilibrium molecular dynamics (NEMD) [33,34] method in which the energy flux is generated by exchanging the kinetic energy of two particles (in this paper are atoms). After a certain number of times of momentum exchanges, the average heat flux is obtained via the total energy exchange divided by response time and cross-sectional area. As shown in Figure 4, the heat flux is imported Input inward the system from both ends of the model and the temperature gradient is provided via a statistic of the interior temperature after the system reaches stability and then the thermal conductivity is calculated by the Fourier law of heat conduction as follows: where dT/dz is the temperature gradient. Because the direction of the flux is opposite from the gradient, the thermal conductivity is always positive. J Z represents the energy flux in the Z-direction which can be expressed as where ∆t is the interval time of every kinetic energy exchange and ∆E is the generating energy between two layers. A is the area of XY-direction of the model. Factor 2 is due to periodic boundary conditions so the amount of heat flux flows at either side of the model where ∆t is the interval time of every kinetic energy exchange and ∆E is the generating energy between two layers. A is the area of XY-direction of the model. Factor 2 is due to periodic boundary conditions so the amount of heat flux flows at either side of the model is E/2. The three-dimensional period model is 200 Å in the Z direction and 50 × 50 Å in the XY direction (the aspect ratio is 4:1), as recommended in the references. In addition, the model is divided into 40 layers, as shown in Figure 4. The time step is 0.5 fs. First, the geometry optimization of the system is conducted to obtain energy minimization, followed by an equilibration stage with a thermostat for 200 ps. Finally, MD simulations are conducted for 500 ps at 300 K in the NVT ensemble, and the trajectory and dynamical properties are recorded for the thermodynamic analyses. Figure 5 depicts the variation of RDF as a function of water content in the catalyst layer. Specifically, Figure 5a shows flat and wide peaks in all these Pt-C curves, indicating the uniformity of the distribution of Pt in the supported carbon. The minimum concatenate distance between carbon and Pt is set at 2.8 Å in both the inner and outer pores. It is clear that the RDF value first increases and then decreases with the λ ranges from 1 to 10, indicating that the increasing water facilitates the adsorption of Pt in carbon but the excess water reverses this trend. The RDF values shown in Figure 5b indicate that the Pt-H + distribution is not influenced by the water contents. Unlike the H + , Figure 5c shows that the gPt-O(r) increases with increasing λ, which illustrates that water promotes the accumulation of O2 near the catalyst. The RDF of sulfonate-sulfonate represents the distribution homogeneity of the sulfonic acid group, and it can be seen from Figure 5d that with the increase of λ, the short-range effect of sulfonate decreases, indicating that the ionomer uniformity is enhanced. The agglomeration phenomenon is obvious when the water content is very small (λ = 1). The RDF values of Pt-S in Figure 5e show that the water content has little effect on the ionomer distribution around the Pt-C surface. In Figure 5f,g, both of the RDF values of S-H2O and S-H3O + decrease with increasing water content, and the absolute value of gS-H3O+(r) is higher than that of the S-H2O. This indicates that water promotes the aggregation of H + around the acid group, and the degree of aggregation is much higher than that of water, which can be attributed to the ion adsorption and desorption effect of cation and anionic. It is worth noting that the second peaks of S-H3O + RDFs are the hydrogen bond interaction between the sulfur atom and hydronium ion. Therefore, when λ = 1, the peak is missing because of the lack of hydronium ion. The gO-O(r) in Figure  5h is small when λ = 1 and becomes higher when the λ continues to increase, indicating that the water content has little effect on the oxygen short-range order except for the state of extreme lack of water. To summarize, the increase of water is beneficial to the Pt-C combination and the accumulation of O2 near the catalyst. It also reduces the aggregation effect of water and hydronium ion around the sulfonic acid ion. The promotion effects are in accordance with Feng et al simulation results [21] and further lift the current density as experimental data showed in Carcadea et al. [35]. The three-dimensional period model is 200 Å in the Z direction and 50 × 50 Å in the XY direction (the aspect ratio is 4:1), as recommended in the references. In addition, the model is divided into 40 layers, as shown in Figure 4. The time step is 0.5 fs. First, the geometry optimization of the system is conducted to obtain energy minimization, followed by an equilibration stage with a thermostat for 200 ps. Finally, MD simulations are conducted for 500 ps at 300 K in the NVT ensemble, and the trajectory and dynamical properties are recorded for the thermodynamic analyses. Figure 5 depicts the variation of RDF as a function of water content in the catalyst layer. Specifically, Figure 5a shows flat and wide peaks in all these Pt-C curves, indicating the uniformity of the distribution of Pt in the supported carbon. The minimum concatenate distance between carbon and Pt is set at 2.8 Å in both the inner and outer pores. It is clear that the RDF value first increases and then decreases with the λ ranges from 1 to 10, indicating that the increasing water facilitates the adsorption of Pt in carbon but the excess water reverses this trend. The RDF values shown in Figure 5b indicate that the Pt-H + distribution is not influenced by the water contents. Unlike the H + , Figure 5c shows that the g Pt-O (r) increases with increasing λ, which illustrates that water promotes the accumulation of O 2 near the catalyst. The RDF of sulfonate-sulfonate represents the distribution homogeneity of the sulfonic acid group, and it can be seen from Figure 5d that with the increase of λ, the short-range effect of sulfonate decreases, indicating that the ionomer uniformity is enhanced. The agglomeration phenomenon is obvious when the water content is very small (λ = 1). The RDF values of Pt-S in Figure 5e show that the water content has little effect on the ionomer distribution around the Pt-C surface. In Figure 5f,g, both of the RDF values of S-H 2 O and S-H 3 O + decrease with increasing water content, and the absolute value of g S-H3O+ (r) is higher than that of the S-H 2 O. This indicates that water promotes the aggregation of H + around the acid group, and the degree of aggregation is much higher than that of water, which can be attributed to the ion adsorption and desorption effect of cation and anionic. It is worth noting that the second peaks of S-H 3 O + RDFs are the hydrogen bond interaction between the sulfur atom and hydronium ion. Therefore, when λ = 1, the peak is missing because of the lack of hydronium ion. The g O-O (r) in Figure 5h is small when λ = 1 and becomes higher when the λ continues to increase, indicating that the water content has little effect on the oxygen short-range order except for the state of extreme lack of water. To summarize, the increase of water is beneficial to the Pt-C combination and the accumulation of O 2 near the catalyst. It also reduces the aggregation effect of water and hydronium ion around the sulfonic acid ion. The promotion effects are in accordance with Feng et al simulation results [21] and further lift the current density as experimental data showed in Carcadea et al. [35].  Figure 6 depicts the variation of RDF as a function of the ionomer-to-carbon ratios in the catalyst layer. Figure 6a indicates that the RDF values of Pt-C slightly increase as the I/C increase and then remain constant when further increases the ionomer, which confirms that an optimized ionomer content can benefit the combination of Pt and carbon. Figure 6b shows the RDF values of Pt-H + are not affected by the I/C, which is similar to the effect of water content. Figure 6c shows that the gPt-O(r) decreases with an increasing ionomer, which illustrates that the ionomer impedes the accumulation of O2 near the catalyst and will further be averse to the reaction. Figure 6d shows that the S-S short-range order becomes higher with the high ionomer content. However, the peak values are not  Figure 6 depicts the variation of RDF as a function of the ionomer-to-carbon ratios in the catalyst layer. Figure 6a indicates that the RDF values of Pt-C slightly increase as the I/C increase and then remain constant when further increases the ionomer, which confirms that an optimized ionomer content can benefit the combination of Pt and carbon. Figure 6b shows the RDF values of Pt-H + are not affected by the I/C, which is similar to the effect of water content. Figure 6c shows that the g Pt-O (r) decreases with an increasing ionomer, which illustrates that the ionomer impedes the accumulation of O 2 near the catalyst and will further be averse to the reaction. Figure 6d shows that the S-S short-range order becomes higher with the high ionomer content. However, the peak values are not multiplied in accordance with the increasing sulfonate group. Figure 6e shows a slight decline of the g Pt-S (r), illustrating that the increase of ionomer will not cluster around the catalysts and therefore it is not advantageous to the proton transport. Figure 6f shows that g S-H2O (r) values are almost the same. The different influences in the sulfonate-water distribution induced by water content and ionomer content can be explained in terms of that the water is adequate. As can be seen in Figure 6g, the values of g S-H+ (r) first increase with the I/C ranging from 0.4 to 1.2 and then decline when the I/C reaches 1.6. This can be attributed to the relative absence of enough hydronium ion. The results of O-O values indicate that little effect on the oxygen short-term order is induced by varying the value of I/C, which is similar to the effect of water content. To summarize, the increase of ionomer can a little bit promote the Pt-C combination. However, it impedes the transportation of oxygen to Pt and proton around the sulfonic group. The rules fit well with Carcadea and Huang et al. experimental results [35,36], in which the current density did not significantly increase when ionomer increase.

Effects of Water Content and Ionomer-to-Carbon Ratio on Diffusion Coefficients
According to Equations (6) and (7), the diffusion coefficients are calculated of various water contents and the results are shown in Figure 7. To have a better understanding of the properties and mechanism of the oxygen and hydronium ion diffusion, the water cluster morphology in the catalyst layer is worth studying. In this paper, the water cluster morphology is given via the water molecular density contours after the dynamics simulation. Moreover, the connecting threshold is set to 2 Å, referring to the 3.5 Å in literature [27] and 1.4 Å in [37]. The water cluster morphology diagrams are demonstrated in Figure 8 and the yellow isosurface represents the water clusters region. In Figure 7, the diffusion coefficient of oxygen first reduces from 2.6 × 10 −6 cm 2 /s to 2.1 × 10 −6 cm 2 /s and then increases to 3.5 × 10 −6 cm 2 /s with increasing water contents. Moreover, the diffusion coefficient of the proton shows an increase of 0.1, 0.25, 0.3, and 0.36 (10 −6 cm 2 /s) with the corresponding λ of 1, 4, 7, and 10, respectively. The D O2 is an order of magnitude larger than the D H+ , which is in good agreement with the D O2 in literature [21] and D H+ in [38]. Due to the measurement difficulties in experiments and modeling error in simulation methods, the same order of magnitude of error level can be acceptable. Figure 8 shows the variations of water clusters form in the CL system. When the λ = 1, the hydronium ions are present in discontinuous spots. According to the RDF results of S-H 3 O + mentioned in Section 3.1, these water clusters tend to spread around the sulfonic acid group. With increasing water content, the clusters volume become grew and interconnected until all the water formed continuous pathways (λ = 10). As the S-H 2 O RDF values display in Figure 5f, the water will spread evenly in the CL system rather than merely around the ionomer when λ = 10. The reason why the D H+ is an order of magnitude smaller than oxygen can be explained as the received stronger interaction from the sulfonate anion. The addition of D H+ is due to that the excessive water declines this kind of interaction as verified in Figure 5f.
indicate that little effect on the oxygen short-term order is induced by varying the value of I/C, which is similar to the effect of water content. To summarize, the increase of ionomer can a little bit promote the Pt-C combination. However, it impedes the transportation of oxygen to Pt and proton around the sulfonic group. The rules fit well with Carcadea and Huang et al. experimental results [35,36], in which the current density did not significantly increase when ionomer increase. content, the clusters volume become grew and interconnected until all the water formed continuous pathways (λ = 10). As the S-H2O RDF values display in Figure 5f, the water will spread evenly in the CL system rather than merely around the ionomer when λ = 10. The reason why the DH+ is an order of magnitude smaller than oxygen can be explained as the received stronger interaction from the sulfonate anion. The addition of DH+ is due to that the excessive water declines this kind of interaction as verified in Figure 5f. these water clusters tend to spread around the sulfonic acid group. With increasing water content, the clusters volume become grew and interconnected until all the water formed continuous pathways (λ = 10). As the S-H2O RDF values display in Figure 5f, the water will spread evenly in the CL system rather than merely around the ionomer when λ = 10. The reason why the DH+ is an order of magnitude smaller than oxygen can be explained as the received stronger interaction from the sulfonate anion. The addition of DH+ is due to that the excessive water declines this kind of interaction as verified in Figure 5f. The variation of the DO2 can be explained in terms of the three diffusion ways of oxygen in the CL system-diffusing in water, diffusing in gas, and diffusing across the gaswater interface. In fact, the "gas" refers to the space that is not occupied by water and ionomer. When the λ is quite small, the diffusion of O2 mainly occurs in the gas region so the diffusion coefficient is relatively large. As the water increase, the area of the gas-water interface increases, which further results in a higher oxygen diffusion resistance, as depicted in Figure 8b,c. Therefore, the DO2 reduces when the λ= 4 and λ = 7. With a continuous increase in the water content (λ = 10), the connected water pathways increases and the area of the gas-water interface becomes smaller, resulting in a remarkable increase in the DO2. In summary, the oxygen diffusion is affected to a great extent by the water cluster states and their connected area. The proton diffusion is only affected by the water content.
The effects of the ionomer-to-carbon ratio on the reactant diffusion coefficients and water cluster morphology are depicted in Figures 9 and 10, respectively. With increasing I/C ratio, the DO2 value declines from 2.2 to 2.1 (10 −6 cm 2 /s), 2.05, and 1.7 (10 −6 cm 2 /s), while the DH+ increases with from 0.05 to 0.08, 0.3, and 0.36 (10 −6 cm 2 /s), respectively. The results imply that the addition of ionomer promotes proton diffusion but impedes the transpor- The variation of the D O2 can be explained in terms of the three diffusion ways of oxygen in the CL system-diffusing in water, diffusing in gas, and diffusing across the gas-water interface. In fact, the "gas" refers to the space that is not occupied by water and ionomer. When the λ is quite small, the diffusion of O 2 mainly occurs in the gas region so the diffusion coefficient is relatively large. As the water increase, the area of the gas-water interface increases, which further results in a higher oxygen diffusion resistance, as depicted in Figure 8b,c. Therefore, the D O2 reduces when the λ= 4 and λ = 7. With a continuous increase in the water content (λ = 10), the connected water pathways increases and the area of the gas-water interface becomes smaller, resulting in a remarkable increase in the D O2 . In summary, the oxygen diffusion is affected to a great extent by the water cluster states and their connected area. The proton diffusion is only affected by the water content.
The effects of the ionomer-to-carbon ratio on the reactant diffusion coefficients and water cluster morphology are depicted in Figures 9 and 10, respectively. With increasing I/C ratio, the D O2 value declines from 2.2 to 2.1 (10 −6 cm 2 /s), 2.05, and 1.7 (10 −6 cm 2 /s), while the D H+ increases with from 0.05 to 0.08, 0.3, and 0.36 (10 −6 cm 2 /s), respectively. The results imply that the addition of ionomer promotes proton diffusion but impedes the transportation of oxygen. The promotion to proton transport is understandable because of the sulfonic acid conduction carrier. The variation of D O2 can also be explained by the water morphology in Figure 10. It can be seen in Figure 10 that with increasing ionomer content, the connected pathways become less, and then the independent clusters occur. According to the S-H 2 O RDF depicted in Figure 6f, these clusters will be center on the sulfonic acid group all the time. It should be noted that the declining extent of D O2 is higher than water, implying a more obvious oxygen impediment effects from ionomer. On the contrary, the promotion effect of ionomers on the proton conductivity is more significant than that of water. In summary, the effects of ionomer increase on oxygen transport will also be influenced by the water formation in CL. picted in Figure 8b,c. Therefore, the DO2 reduces when the λ= 4 and λ = 7. With a continuous increase in the water content (λ = 10), the connected water pathways increases and the area of the gas-water interface becomes smaller, resulting in a remarkable increase in the DO2. In summary, the oxygen diffusion is affected to a great extent by the water cluster states and their connected area. The proton diffusion is only affected by the water content.
The effects of the ionomer-to-carbon ratio on the reactant diffusion coefficients and water cluster morphology are depicted in Figures 9 and 10, respectively. With increasing I/C ratio, the DO2 value declines from 2.2 to 2.1 (10 −6 cm 2 /s), 2.05, and 1.7 (10 −6 cm 2 /s), while the DH+ increases with from 0.05 to 0.08, 0.3, and 0.36 (10 −6 cm 2 /s), respectively. The results imply that the addition of ionomer promotes proton diffusion but impedes the transportation of oxygen. The promotion to proton transport is understandable because of the sulfonic acid conduction carrier. The variation of DO2 can also be explained by the water morphology in Figure 10. It can be seen in Figure 10 that with increasing ionomer content, the connected pathways become less, and then the independent clusters occur. According to the S-H2O RDF depicted in Figure 6f, these clusters will be center on the sulfonic acid group all the time. It should be noted that the declining extent of DO2 is higher than water, implying a more obvious oxygen impediment effects from ionomer. On the contrary, the promotion effect of ionomers on the proton conductivity is more significant than that of water. In summary, the effects of ionomer increase on oxygen transport will also be influenced by the water formation in CL.

Effects of Water Content and Ionomer-to-Carbon Ratio on Thermal Conductivity
The atom momentum exchange NEMD program described in Section 2.2.3 is conducted to calculate the energy flux and temperature gradient, which is further utilized to obtain the thermal conductivity coefficients. Figure 11a shows the thermal conductivity

Effects of Water Content and Ionomer-to-Carbon Ratio on Thermal Conductivity
The atom momentum exchange NEMD program described in Section 2.2.3 is conducted to calculate the energy flux and temperature gradient, which is further utilized to obtain the thermal conductivity coefficients. Figure 11a shows the thermal conductivity of the catalyst layer as a function of water content. With increasing water content, the thermal conductivity coefficients increase gradually from 0.5 to 0.8, 1.06, and 2.5 (W/(m·K)). It should be noted that the result in this study accords with the studies of Kang et al. [39] and Jiang et al. [40]. This increment can be attributed to a huge variation of the thermal conductivity of each component in CL, and the Pt-C system has a higher thermal conductivity value than others. In the periodic volume-fixed computational domain, the volume fraction of water enlarges, while the volume fractions of the other components decrease with increasing λ. The inconsistent amount of these decreases makes the volume fraction of the Pt-C bigger; hence, the thermal conductivity goes up. The result proves that the increasing water is beneficial to the enhancement of the heat conduction of the catalyst layer. The thermal conductivity variations under different I/C ratios are shown in Figure 11b. Similarly, the thermal conductivity coefficient increases gradually from 0.3 to 0.7, 1.06, and 1.6 (W/(m·K)). The effect of I/C ratios on the thermal conductivity is inapparent compared with λ variation, indicating that the promotion effect of ionomer is a slightly weaker than water.
Membranes 2021, 11,148 12 of 14 (a) (b) Figure 11. Thermal conductivity of CL with (a) different water contents and (b) different ionomer-to-carbon ratios. Figure 12 shows the variation of temperature as a function of water content. It is found that the temperature gradient lines are not straight. The bending instances occur at the location where the Pt-C particles are dispersed. Because of the huge difference in thermal conductivity for Pt-C particles, the temperature gradient displays a sudden rise and fall in the whole system with increasing layers. One notable aspect is the effect of particle dispersion on the error bar value. Specifically, the error value becomes higher when the particles move closer to both the hot and cold ends of the model. Moreover, the error bars obtained in Figure 11 are higher than in studies by Zheng et al. [18] and Fan et al. [26], which can be attributed to the existence of carbon-supported Pt in this system model. Figure 11. Thermal conductivity of CL with (a) different water contents and (b) different ionomer-to-carbon ratios. Figure 12 shows the variation of temperature as a function of water content. It is found that the temperature gradient lines are not straight. The bending instances occur at the location where the Pt-C particles are dispersed. Because of the huge difference in thermal conductivity for Pt-C particles, the temperature gradient displays a sudden rise and fall in the whole system with increasing layers. One notable aspect is the effect of particle dispersion on the error bar value. Specifically, the error value becomes higher when the particles move closer to both the hot and cold ends of the model. Moreover, the error bars obtained in Figure 11 are higher than in studies by Zheng et al. [18] and Fan et al. [26], which can be attributed to the existence of carbon-supported Pt in this system model. the location where the Pt-C particles are dispersed. Because of the huge difference in thermal conductivity for Pt-C particles, the temperature gradient displays a sudden rise and fall in the whole system with increasing layers. One notable aspect is the effect of particle dispersion on the error bar value. Specifically, the error value becomes higher when the particles move closer to both the hot and cold ends of the model. Moreover, the error bars obtained in Figure 11 are higher than in studies by Zheng et al. [18] and Fan et al. [26], which can be attributed to the existence of carbon-supported Pt in this system model.

Conclusions
In this study, a molecular model containing Pt, carbon, ionomer, water, oxygen, and hydronium ion is built to investigate the effects of water content and ionomer-to-carbon ratio on the RDFs between components. The reactants′ diffusion properties and the thermal conductivity are analyzed after equilibrium MD and nonequilibrium MD simulations. The main conclusions are drawn as follows.
The increase of water is beneficial to the Pt-C combination and the oxygen accumulation near the catalyst. It also reduces the aggregation effect of water and hydronium ion around the sulfonic acid ion. The oxygen diffusion coefficient firstly grows and then declines, which is dependent closely on the water cluster morphology and connection area.

Conclusions
In this study, a molecular model containing Pt, carbon, ionomer, water, oxygen, and hydronium ion is built to investigate the effects of water content and ionomer-to-carbon ratio on the RDFs between components. The reactants diffusion properties and the thermal conductivity are analyzed after equilibrium MD and nonequilibrium MD simulations. The main conclusions are drawn as follows.
The increase of water is beneficial to the Pt-C combination and the oxygen accumulation near the catalyst. It also reduces the aggregation effect of water and hydronium ion around the sulfonic acid ion. The oxygen diffusion coefficient firstly grows and then declines, which is dependent closely on the water cluster morphology and connection area. In addition, the proton diffusion coefficient increases with increasing water content. The increase of ionomer impedes the transportation of oxygen to Pt and protons around the sulfonic group. The oxygen transport is retarded while the proton diffusion is promoted, which is also closely related to the water cluster morphology and ionic interaction. The thermal conductivity of CL can be enhanced by increasing the water and ionomer contents. On the basis of this conclusion, we suggest that the I/C ratio should be within reasonable bounds, and the best ratio is approximately 0.8 when designing the structure of the catalyst.