J. Phys. Chem. C 2007, 111, 13627-13634
13627
Self-Organization in Catalyst Layers of Polymer Electrolyte Fuel Cells Kourosh Malek,†,‡ Michael Eikerling,*,†,‡ Qianpu Wang,† Titichai Navessin,† and Zhongsheng Liu† Institute for Fuel Cell InnoVation, National Research Council of Canada, 4250 Wesbrook Mall, VancouVer, BC V6T 1W5, Canada, and Department of Chemistry, Simon Fraser UniVersity, Burnaby, BC V5A 1S6, Canada ReceiVed: April 5, 2007; In Final Form: June 10, 2007
Understanding the factors that control microstructure formation in catalyst layers (CLs) of polymer electrolyte fuel cells is of vital importance for improving the operation of these cells. Here, we employ, for the first time, coarse-grained molecular dynamics simulations to perform a structural analysis of the microphase segregation occurring during the fabrication process of CLs. Our mesoscale simulations provide insights into the structural correlations and dynamical behavior of different phases in the catalyst layer composite. This versatile computational study, moreover, rationalizes how the solvent used in catalyst layer fabrication influences the evolution of stable agglomerated conformations. In this realm, we evaluate dispersion media with distinct dielectric properties in view of capabilities for controlling the sizes of carbon/Pt agglomerates and ionomer domains and the resulting pore network topology. These insights are highly valuable for the structural design of catalyst layers with optimized performance and stability.
Introduction Enhancements in the performance and stability of membrane electrode assemblies (MEAs) for polymer electrolyte fuel cells (PEFCs) hinge decisively on the structure-related properties of the catalyst layers on the cathode and anode sides.1-4 Typical catalyst layers (CLs) are fabricated as random heterogeneous composites to meet the multifunctional requirements of transport properties and electrocatalytic activity. They consist of Pt nanoparticles, randomly dispersed on a carbonaceous substrate. Moreover, CLs contain a Nafion ionomer that forms a protonconducting network.5 Carbon particles exhibit a tendency to aggregate, thereby forming “agglomerates” that possess internal porosities due to primary pores. Evidently, the CL structure involves particles and pores that span a wide range of length scales, including Pt particles with sizes of 2-4 nm, carbon particles and primary pores with sizes in the range of 3-10 nm, agglomerates of carbon/Pt with sizes of ∼50 nm, and secondary pores between agglomerates with sizes in the range of 10-50 nm. This implies a hierarchical structural picture for CLs at different levels of complexity.3,6,7 The foremost objective in optimizing the composition and structure of CLs is to provide the most uniform distributions of reactants and reaction rates possible and, thus, attain the highest catalyst utilizations and voltage efficiencies.3-8 This demands well-balanced distributions and percolation properties of the interpenetrating phases of carbon/Pt agglomerates, ionomer domains, and pore spaces that guarantee simultaneous high rates of transport of electrons, protons, and gaseous reactants/product water to or from Pt sites. Recent efforts in improving the fabrication of CLs have explored ways to improve the Pt catalytic utilization by increasing the interfacial area between Pt particles and the ionomer.9-11 This is achieved by mixing the ionomer with dispersed Pt/C catalysts in the ink suspension prior to deposition † ‡
National Research Council of Canada. Simon Fraser University.
to form a CL. The choice of a dispersion medium determines whether the ionomer is to be found in the solubilized, colloidal, or precipitated forms. This influences the microstructure and pore size distribution of the CLs.10,11 The random morphology in turn has a marked impact on the transport and electrochemical properties of the emerging composite medium.3-8,11 In principle, only catalyst particles at Pt/liquid water interfaces are electrochemically active. Electrochemical reactions, thus, occur on the walls of wetted pores inside and between agglomerates as well as at interfaces between Pt and the wetted ionomer. Relative contributions of these distinct types of interfaces, viz., Pt/water in carbon pores and Pt/water in ionomer pores, will depend on the corresponding interfacial areas. Evidently, high catalytic yields hinge on compositions and porous structures with wellattuned pore size distributions and wetting properties of the pore network. To improve the structure-performance relationships of CLs, preparation methods exploit variations in terms of the applicable solvent, particle sizes of the primary carbon powders, wetting properties of the carbon materials, and composition of the catalyst layer ink.9,12-17 These factors determine the complex interactions among carbon/Pt particles, ionomer molecules, and solvent molecules and, thus, control the catalyst layer formation process.18 As demonstrated in numerous experimental and theoretical studies, the microstructure and composition of structural elements at the mesoscopic scale usually referred to as agglomerates, which consist of carbon/Pt particles, are vital for the performance of the layer.3,5,10,11,19 Apart from characterizing the size, internal porosity, and wetting properties of the internal/ external surfaces of such agglomerates, it is critical to understand whether the ionomer is able to penetrate into the primary pores inside the agglomerates and form a proton-conducting, embedded phase therein or whether such a proton-conducting phase can only be formed in the larger void spaces between the agglomerates.20 The latter question is of significant fundamental interest as well as of practical importance for performance
10.1021/jp072692k CCC: $37.00 © 2007 American Chemical Society Published on Web 08/17/2007
13628 J. Phys. Chem. C, Vol. 111, No. 36, 2007 modeling of CLs and PEFCs at the macroscopic device level.2,3,5-8,19,21,22 In general, we expect valuable insights for the advanced design of catalyst layers from understanding the micromorphology of interconnected phases of agglomerates, ionomer, and pore spaces. Various phenomenological models have been employed to rationalize catalyst layer operation.3,5-8,19,22 Less effort has focused on exploring the effects of the microstructure of the CL on transport and reactivity. Concomitantly with the lack of understanding structure formation in CL, there is a lack of systematic techniques to fabricate catalyst layers with controlled microstructures. In experimental characterization, more advanced analytical tools are necessary to directly correlate transport and reaction processes with structural details at the meso- to microscale and down to the atomistic scale. Molecular simulation techniques have emerged as independent tools that could complement experiments in studying formation and characterization of complex materials.23-25 Coarse-grained molecular simulation (coarse-grained molecular dynamics, CG-MD) techniques can describe the system at the micro- to mesolevel, while still being able to capture the morphology at long time and length scales.26-28 The application of CG models, however, requires special care. Due to the requisite reduction in the number of degrees of freedom, CG simulations may not be able to accurately predict physical properties that rely upon time correlation functions (e.g., diffusion).28 Here, we present the first study of microstructure formation in CLs based on mesoscale calculations. In standard experimental procedures, a composite network of ionomer and carbon particles is formed by colloidal dispersion in an organic (alcoholic) solvent.9 It is important to understand the random formation process to devise fabrication techniques that allow control of the structure. Our CG-MD calculations simulate these processes during CL fabrication. The purposes are to understand how the dynamical-structural properties of CLs evolve and how they are correlated after equilibration. We will rationalize the effect of distinct dispersion media in terms of their dielectric properties. Our study provides vital insight into the phase segregation of ionomer and carbon/Pt domains. The structures formed will be evaluated with respect to composition and the sizes of the agglomerates, sizes of the ionomer domains, pore size distribution, and pore network morphology. Computational Methodology A significant number of computational approaches have been employed to understand the phase-segregated morphology and transport properties of water-swollen Nafion membranes.26-35 Because of computational limitations, full atomistic models are not able to probe the random morphology of these systems. However, as demonstrated by these membrane simulations and applications to other random composite media, mesoscale models are computationally feasible to capture their morphology. Several approaches have been used such as lattice Monte Carlo,31 termed cellular automata,32 and coarse-grained mesodynamics based on self-consistent mean field theory.26,27 For Nafion, most of these simulations support the idea that narrow water-filled channels and irregularly shaped, nanometer-size clusters of ionic head groups and water form the protonconducting network that is embedded in the hydrophobic matrix. Structural complexity is even more pronounced in CLs, since they consist of a mixture of Nafion ionomer, Pt clusters supported on carbon particles, solvent, and water. Independent computational strategies are generally needed to simulate the microstructure formation in CLs. For instance, the structure of
Malek et al.
Figure 1. (a) Sketch of the coarse-grained model of Nafion ionomer molecules and carbon particles. The upper image shows a typical folded monomeric sequence of Nafion used in our simulations. The spheres show the coarse-grained beads corresponding to the hydrophobic backbone (red) or hydrophilic side chains (green). Nafion polymer is represented by 20-unit oligomers of length ∼30 nm. The lower image represents the optimized configuration of a coarse-grained model for carbon particles. (b) Interparticle interaction energies, showing the relaxation of the structure to equilibrium, starting from the initial optimized structure. B ) backbone of Nafion, W ) water, and H ) hydronium beads. The equilibration is established after ∼0.08 µsec.
the ionomer phase in CLs cannot be trivially inferred from that in membrane simulations, as there are distinct correlations among Nafion, water, and carbon particles in CLs. Here, we utilize a computational approach based on CG-MD simulations. The coarse-grained model was developed in two major steps. In the first step, we replaced the Nafion chains, water and hydronium molecules, other solvent molecules, and carbon/Pt particles with corresponding spherical beads with a predefined subnanoscopic length scale. In the second step, we specified the parameters of renormalized interaction energies between the distinct beads. Coarse-Graining Model. We considered four main types of spherical beads in our simulations: polar, nonpolar, apolar, and charged beads.36 Clusters including a total of four water molecules or three water molecules plus a hydronium ion were represented by polar beads of radius 0.43 nm. The configuration of a folded Nafion chain is shown in Figure 1a. According to the literature,26 a side chain unit in the Nafion ionomer has a molecular volume of 0.306 nm3, which is comparable to the molecular volume of a four-monomer unit of poly(tetrafluoroethylene) (PTFE) (0.325 nm3). Therefore, each four-monomer unit (-CF2CF2-CF2CF2-CF2CF2-CF2CF2-) in our simulations and each side chain (represented by a charged bead) were coarse-grained as spherical beads of volume 0.315 nm3 (r ) 0.43 nm). The hydrophobic Nafion backbone was replaced by a coarse-grained chain of 20 apolar beads, as illustrated in Figure 1a. Carbonaceous particles can be coarse-grained in various ways, building upon a new technique called multiscale coarse graining (MS-CG).37-39 In this method, the CG potential parameters are systematically obtained from atomistic-level interactions.38 Using this technique, we constructed a model for semispherical carbon particles based on CG sites in the C60 system, Figure 1a. In our simulations, we assumed carbon particles of about 5 nm diameter. Each carbon particle was represented in the coarse-grained model by 60 beads of diameter 0.86 nm. Each bead consisted of eight carbon atoms. The initial spherical structure of the carbon particles was energy minimized to convert to a random structure, which was used for CG-MD
Self-Organization in CLs of PEFCs
J. Phys. Chem. C, Vol. 111, No. 36, 2007 13629
simulations. Pt nanoparticles supported on carbon particles were not considered explicitly, which is an obvious simplification. One limitation of this approach is that it overlooks the effect of Pt particles on the effective wetting properties of carbon particles. We can empirically incorporate this effect by adjusting the interaction parameters between the carbon particles and dispersion medium. CG-MD Simulations. The initial simulation box had a size of 50 × 50 × 50 nm3. It included 52 carbon particles and 72 coarse-grained Nafion oligomers, each consisting of 20 monomers (and therefore 20 side chains). A total of 1440 CG hydronium ions were added for electroneutrality. The number of water beads was fixed at a value corresponding to approximately 13 waters per side chain (i.e., λ ) 13). The Nafion concentration was kept at 32-38 wt %, which is in the range of typical values used in catalyst layer fabrication. To determine the effect of the system size, independent simulations were performed using a box with 2 times larger length by combining eight equilibrated systems from the previous simulation. This led to a total of 71040 CG particles (∼0.7 million effective atoms). Simulations were performed using a modified version of the GROMOS96 force field by adapting the nonbonding and bonding parameters of CG beads.36 In our force field, interactions between beads are divided into nonbonded interactions, acting between any pair of beads that are within a given cutoff radius, and bonded interactions between beads in a Nafion chain connected by chemical bonds. In the case of nonbonded interactions involving charged side chains, a positive charge (+1) and parameters for repulsion and attraction are assigned to each side chain. The bonded interaction consists of bond, angle, and dihedral terms. The force constant of the harmonic bonding potential is 1250 kJ mol-1 nm-2, allowing sufficient deviations from the equilibrium bond length. A harmonic potential of the cosine type is used for the angles, with a force constant of 25 kJ mol-1 rad -2. The dihedrals represent the energy involved in changing the out-of-plane angles. This enhances the flexibility of Nafion and comprises the advantage of our CG-MD methodology over other mesoscale techniques such as dissipative particle dynamics (DPD).28,40 The bond stretching potential and bond angle bending are expressed in terms of simple harmonic potentials. Nonbonded interactions between beads were modeled by pair potentials of the Lennard-Jones (LJ) type:
[( ) ( ) ]
ULJ(r) ) 4ij
σij r
12
-
σij r
6
(1)
where σij is the minimum distance between two beads and ij represents their interactions. In analogy to ref 36, only five levels of interactions were considered in our simulations: attractive (ij ) 5 kJ mol-1), semiattractive (ij ) 4.2 kJ mol-1), intermediate (ij ) 3.4 kJ mol-1), semirepulsive (ij ) 2.6 kJ mol-1), and repulsive (ij ) 1.8 kJ mol-1). The effective bead diameter (σij) is 0.43 nm. The LJ parameters 4ijσij6 and 4ijσij12 can, thus, be calculated depending on the level of interactions. A cutoff of 1.0 nm was used for van der Waals interactions in our simulations. Carbon particles of the VULCAN type are partially hydrophobic (contact angle ∼90° 41) with a semirepulsive (ij ) 2.6 kJ mol-1) interaction with water and Nafion side chains and semiattractive interactions with other carbon particles as well as with Nafion backbones.10,11 In the initial configuration, the centers of all beads were randomly placed on cubic lattice points. First, we optimized the initial structure for 100 ps at T ) 0 °C using 20 fs time
steps in the course of a steep integration procedure. This short energy minimization slightly displaced the overlapped beads of carbon particles, ionomer, water molecules, and hydronium ions in the initial structure. The resulting structure was subsequently embedded in a continuum dielectric medium. In this so-called implicit solvent approach, a continuous medium with relative dielectric permittivity fills all void spaces between the beads. To exploit the effect of the solvent on the evolving CL microstructure, we examined different solvents that represent a wide range of dielectric permittivities, viz., water (W; ≈ 80), isopropyl alcohol (IPA; ≈ 20), and hexane (HX; ≈ 2). Equilibration of the embedded structure was performed by an annealing procedure. The structure was first expanded over a period of 50 ps by gradually increasing the temperature from 295 to 395 K. Then, a short MD simulation was performed for another 50 ps in an NVT ensemble (thermalization), followed by a cooling procedure to 295 K. Following this procedure, we conducted MD simulations in an NPT ensemble. In the latter simulations, the total energy decreased steeply during an initial equilibration period of ∼0.08 µs after which the total energy of the system converged and became stable (∆E < 500 kJ mol-1), as shown in Figure 1b. After the stable density was reached, production runs were continued in the NPT ensemble for up to 0.7 µs for sampling of configurations that were used in the structural analysis. The time step was 0.05 ps. The temperature was controlled by the Berendsen algorithm, which mimics a weak coupling to an external heat bath with given temperature T0. The weak coupling algorithm was applied separately for each component (polymer, carbon particles, water, hydronium ions) with a time constant of 0.1 ps and a temperature of 295 K. During the production run, structures were saved every 500 steps (25 ps) and used for the analysis. We employed a Monte Carlo procedure to determine the size distributions of ionomer and carbon particle domains in our coarse-grained MD simulations.42 Equilibrated morphologies in the presence of three different solvents ( ) 80, 20, and 2) were used as the input configurations. All simulations were carried out using a modified GROMACS package (http://www.gromacs.org).42,43 Visualization was done by using the VMD v1.8.1 commercial package.44 Results and Discussion The main focus of the present study is to understand the effect of interactions in the carbon/Pt-ionomer-solvent system on the segregation of ionomer and carbon/Pt domains in the colloidal dispersion of the catalyst layer ink. We characterize equilibrated mesostructures in terms of phase segregation of carbon, Nafion, and solvent phases, as well as the extent of agglomerate formation. Moreover, for quantitative structural analysis we determined site-site radial distribution functions (RDFs), size distributions of ionomer domains and carbon agglomerates, pore size distributions, pore shapes, and density map profiles. A Monte Carlo procedure was used to determine the size distribution of the ionomer and carbon domains from our CG-MD simulations.43 The radius of a pore was determined by using a procedure based on the CHANNEL algorithm.45 An initial random position was chosen inside the pore network and at any given distance along the pore axis. Pore sizes were then evaluated by calculating the maximum size for a diffusing spherical probe to still fit in the pore without overlap with the van der Waals radii of beads in the pore wall. Figure 2 represents snapshots of the evolution process for a volume element of the CL blend, consisting of a few ionomer clusters and carbon agglomerates. An annealed configuration
13630 J. Phys. Chem. C, Vol. 111, No. 36, 2007
Malek et al.
Figure 2. Spontaneous microstructure formation: snapshots along the trajectory and during microstructure formation at (a) t ) 100 ps, annealed configuration, (b) t ) 40 ns, and (c) t ) 320 ns, final equilibrated configuration. The gray spheres show the carbon particles. Nafion backbones correspond to red beads, while green beads show the side chains. Small blue beads represent hydronium ions inside the pore network. Water beads are not shown for better clarity. The equilibrated structure reveals separate hydrophilic (water, hydronium ions, and side chains) and hydrophobic (carbon particles and ionomer backbones) domains. In (c), the beads of Nafion backbones are removed and Nafion side chains are represented by green transparent beads for better visualization of carbon clusters.
TABLE 1: Calculated Structural Parameters (nm) That Characterize the Phase-Segregated CL Morphology for Solvents with Different Dielectric Constants solvent dielectric constant
agglomerate size
ionomer cluster size
pore size between agglomerates
pore size inside agglomerates
ionomer pore size
2 20 80
32.788 17.067 15.036
11.230 10.163 9.744
11.534 13.864 16.590
4.795 3.444 2.397
1.971 4.320 6.068
(a) was formed after 100 ps. This configuration self-assembled within 80 ns to convert to the stabilized configuration (c). In the intermediate stages of blend formation, a rapid clustering of carbon particles into irregularly shaped aggregates was observed, while ionomer backbones assembled into a separate interconnected phase in the void spaces, Figure 2c. This structure subsequently coalesced to form agglomerated carbon particles, surrounded by ionomer and water clusters. Assuming an effective density of 1.04 nm3 for the solvent, the composition of the equilibrated structure was Nafion (32.36 wt %, 34.68 vol %), water (3.36 wt %, 4.33 vol %), carbon particles (21.42 wt %, 19.03 vol %), hydronium ions (1.99 wt %, 2.88 vol %), and solvent (40.84 wt %, 39.06 vol %). Snapshots of the final microstructures for these systems were then analyzed and evaluated. Figure 2 generally confirms the hypothesis that carbon particles form agglomerates and that the hydrated ionomer forms a separate phase that is attached to the surface of carbon aggregates. The catalyst layer thus segregates into hydrophobic and hydrophilic regions. The hydrophobic region is constituted by the carbon particles, the backbones of the ionomer, and solvent. The hydrophilic region contains water clusters, hydronium ions, solvent, and the side chains of the ionomer. The porous structure inside the ionomer evolves from the system of isolated hydrophilic clusters to a percolating three-dimensional network of irregular water-filled channels. The equilibrium structure in Figure 2c clearly shows that the side chains tend to be located at the interface between ionomer domains and waterfilled pores. Part of a pore network between agglomerates is also evident. The major finding of simulations for all three solvents is that the ionomer phase and carbon/Pt agglomerates form a phasesegregated composite. Ionomer molecules do not penetrate into
the Pt/carbon clusters. Instead they form proton-conducting volume elements that are connected to the surface of carbon agglomerates. In our equilibrated CL blend, most carbon agglomerates are closely attached to each other, surrounded by a few isolated carbon particles. The equilibrated configuration of a blend of carbon/Nafion/water that is immersed into a polar solvent (IPA, ) 20) is shown in Figure S1 (see the Supporting Information). Table 1 displays the domain sizes of the ionomer and carbon phases as well as the average pore sizes in the CL for the different values of used in the simulations. The effect of solvent on the size of the carbon agglomerates is remarkable. As increases, the average agglomerate size decreases markedly, from 33 nm for ) 2 to 15 nm for ) 80. Concomitantly, the maximum ionomer domain sizes decrease from 11 to 10 nm. The effect of the solvent on the domain sizes is, thus, more pronounced for carbon agglomerates than for the ionomer phase. Figure 3 depicts size distributions of the ionomer and carbon
Figure 3. Effect of the polarity of the solvent on the calculated domain size distributions of (a) ionomer and (b) carbon.
Self-Organization in CLs of PEFCs
J. Phys. Chem. C, Vol. 111, No. 36, 2007 13631
Figure 4. Selected structural correlation functions (or site-site radial distribution functions) between distinct components of the CL, representing the interactions between carbon particles and other components of the catalyst layer. (a-e) Pair correlation function of carbon particles (C-C), ionomer backbones (B-B), ionomer side chains (S-S), water (W-W), and hydronium ions (H-H) for the blend formed in a polar solvent ( ) 20). (f) Effect of the solvent dielectric constant on pair correlation functions for carbon particles and water (C-W) and carbon particles and ionomer side chains (C-S).
domains obtained with different values. The dielectric properties of the solvent strongly affect the size and connectivity within the ionomer and carbon clusters. Strikingly, opposite trends can be seen in Table 1 in the sizes of primary pores inside agglomerates and secondary pores between agglomerates. While the sizes of the primary pores decrease from 4.8 nm for ) 2 to 2.4 nm for )80, the sizes of the secondary pores increase from 11.5 to 16.6 nm. This indicates that in solvents with high polarity carbon particles will assemble into smaller and more compact agglomerates. This highlights the role of the solvent in controlling the size and compactness of agglomerates in CLs. It thus determines the pore size distributions and topologies of porous composite catalyst layers. As corroborated in refs 3 and 8, these fine structural details are pivotal in view of the balance between the fluxes of the reactants and water, catalyst activity, and vaporization capabilities of CLs. The findings for the size distributions of carbon agglomerates are in reasonable qualitative agreement with the experimental results (see Figure S2 in the Supporting Information). The quantitative discrepancy mostly relates to the length scale limitation of mesoscale simulations and poor experimental diagnostics. The validation of our simulation results may be further complicated due to the fact that we only considered monodispersed (initial) carbon particles. In reality, the size of
the carbon particles varies between 5 and 20 nm, leading to a polydispersed carbonaceous mixture. Moreover, the experimental method shown in Figure S2 does not distinguish between carbon and polymer clusters and therefore cannot be considered as a reliable experimental validation tool for direct comparison to our simulations. Advanced structural characterization techniques are generally needed for quantitative validations and finetuning of interaction parameters that could lead to a better correspondence. A systematically conducted experimental modeling study will provide insight into the fundamental interactions in the multicomponent CL mixture. Structural correlation functions or RDFs, gij, between distinct components of the CL provide information on phase segregation and the micro- to mesoscopic structure in the CL. Figure 4a-e shows RDFs and correlation functions for the equilibrated structure of the blend of carbon particles (C), ionomer backbones (B), ionomer side chains (S), water molecules (W), and hydronium ions (H) in a polar solvent with ) 20. We analyzed these structural functions in view of segregation of ionomer and carbon particles into distinct domains, distributions of water and hydronium ions with respect to side chains, polymer backbones and carbon particles, and ionomer and carbon network interconnections. The W-W autocorrelation function, gWW, in Figure 4a exhibits typical characteristics of a dense disordered system of
13632 J. Phys. Chem. C, Vol. 111, No. 36, 2007 interacting particles. Enhanced correlations at short range, r < 2 nm, are indicated by the sharp primary peak and pronounced secondary peaks. These features correspond to the expected short-range ordering in a fluid of hard spheres with LennardJones-type interactions. It should be noted that peak separations are determined by the sizes of the beads in the coarse-grained model. Deviations of the correlation functions from the shape expected for an ideal Lennard-Jones fluid provide information on the spontaneous structural organization of the catalyst layer. The autocorrelation functions gSS and gHH in Figure 4a exhibit a structure similar to that of gWW. This indicates a strong clustering of side chains and hydronium ions due to the aggregation and folding of polymer backbones. The primary S-S and H-H peaks are, however, suppressed compared to the primary peak in gWW due to electrostatic repulsion between these charged beads. In Figure 4b, gBB and gCC exhibit upturns toward small r, superimposed on the primary bead-bead correlations, which correspond to the characteristic dimensions of carbon particles (∼5 nm diameter) and backbone clusters (∼2-3 nm). Figure 4c reveals a strong correlation of carbon particles and polymer backbones (gCB), while gCS is seemingly suppressed. This suggests that the polymer backbones are attached to the surfaces of the carbon particles, while the side chains strive to maximize their separation from the surface of the carbon agglomerates. gCB reveals that about three layers of polymer backbones interact with the carbon surface. Parts d and e of Figure 4 illustrate the distribution of water and hydronium ions relative to distinct components of the CL. The functions gSW and gSH show that the side chains are mostly surrounded with water and hydrated protons. The height of the first peak in gSH is about twice that in gSW, due to stronger attractive interactions of the anionic side chains with hydronium ions than with water.46 Correlations between the hydrophilic species and backbones, gBH and gBW, are significantly stronger than those between these species and carbon particles, gCH and gCW. Both species W and H are expelled from the vicinity of carbon, as indicated by the reduced values of gCW and gCH for r < 2 nm. Overall, the correlation functions discussed above provide valuable structural information at the nanometer scale that allows refining the picture of the phase-segregated catalyst layer morphology. Ionomer backbones form clusters or fibers that are attached to the surface of the carbon particles. The side chains are expelled from the vicinity of carbon. Water and hydronium ions tend to maximize their separation from the carbon while trying to stay in the vicinity of the side chains. Visualization of structural conformations, e.g., those shown in Figure 2, demonstrate that the side chains, hydronium ions, and water molecules form interconnected clusters inside the ionomer domains. Figure 4f illustrates the effects of the solvent dielectric constant on the structural correlations. For the sake of comparison, we have integrated the correlation functions up to a given distance and normalized g(r) to the number of particles within that range. Overall, the magnitude of the peaks decreases from a low-polarity to a high-polarity solvent. The polar solvents ( ) 20, 80) behave similarly, while the effect of the apolar solvent is markedly different ( ) 2). A low encompasses stronger correlations between the C and B components, which indicate a stronger tendency to phase-segregate into hydrophobic and hydrophilic domains. The magnitude of short-range interactions and the carbon agglomerate size steadily decrease with increasing . Therefore, the carbon agglomerates become more separated, as the pore size between the agglomerates increases, Table 1. The peak positions for long-distance correlations are
Malek et al. shifted to large distances for the apolar solvent. A possible explanation is that a low-dielectric-constant solvent forms a separate clustered phase of carbon particles confined between thick fibrous aggregates of the ionomer. In the pair correlation function between carbon particles and water, gCW (Figure 4f), the peaks for the apolar solvent are shifted to higher distances around 4 nm. The correlation functions for the separation of carbon particles and water (C-W) and carbon particles and hydronium ions (C-H) for each of the three solvents share the same qualitative behavior. A correlation appears at around 2 nm for gCS. Figure 4f suggests that, in the presence of an apolar solvent, the structural organization into separate hydrophilic and hydrophobic domains spans to higher distances compared to that in the presence of polar solvents. We have also analyzed the pore shape and pore radius profile for typical hydrophilic pores inside ionomer domains in a CL blend. The profile of the pore radius along the pore axis shows that there are constricted zones inside the pore (Figure S3 in the Supporting Information). The pore wall is hydrophilic in all regions. The fluctuation of the ionomer molecules has a direct effect on the pore size. Moreover, the shape of the water channels in the ionomer strongly depends on the water and solvent contents and the type of solvent. The pore can be represented as consisting of big nodes (cavities) with narrow interconnections, showing an irregular structure (see the Supporting Information). For instance, the pore radius in the ionomer varies from 0.5 nm at the most constricted regions to 4 nm in the cavities. Our analysis showed that the pore radius in the ionomer decreases by using an apolar solvent ( ) 2), Table 1. Overall, explicit illustrations of optimized structural configurations (Figure 2), detailed analysis of the size distributions of the agglomerates and pores (Figure 3 and Table 1), correlation functions (Figure 4), and density profiles of different species along the simulation box (Figure S4 in the Supporting Information, contour plots, detailing the density distribution of water, polymer, and carbon particles on a typical two-dimensional cross-section) help to unravel self-organization of the carbonionomer mixture in the presence of the solvent. Fibrils of ionomer backbones are about 5 nm thick. These fibrils and clusters of side chains, water, and hydronium ions form a larger ionomer domain of size between 8 and 10 nm. “Spherical” carbon agglomerates are discernible. The size of these agglomerates increases in the presence of a low-dielectric-constant solvent. More analysis is needed to characterize the detailed structure inside ionomer clusters and aggregated carbon particles and to quantify the size distributions. The finding that the interior of the agglomerates is ionomer free has important implications for the utilization of the catalyst. Indeed, in such structures, very few Pt particles will be located directly at the three-phase interface of Pt-ionomer-carbon. This confirms an earlier conjecture that reactions in catalyst layers predominantly occur at Pt/water interfaces in water-filled pores inside agglomerates. Figures 2 and S1 (see the Supporting Information) show that at the high concentration of Nafion, the carbon agglomerates become covered by a thick (∼8 nm) ionomer layer. This observation is in good agreement with a recent study (based on atomistic MD simulations) of the structural and dynamical behavior of water and gas transport through a model polymer with a graphite-supported catalyst (Pt particles).47 The size of the ionomer domains decreases in the presence of the high-polarity solvent ( ) 80). This will affect the ability of oxygen molecules to penetrate the ionomer layer to reach catalyst sites inside the agglomerates. In the presence of an apolar solvent, the coverage of the surface of carbon
Self-Organization in CLs of PEFCs agglomerates with the ionomer layer could be complete. When a high-polarity solvent (W, ) 80) is used, the surface of the Pt/C cluster remains partially uncovered by ionomer molecules. The latter situation provides direct routes for transport of oxygen to catalyst sites inside agglomerates via gaseous diffusion. We should notice that the final microstructure also strongly relates to the choice of initial carbon particles in terms of their wetting and geometrical (size and shape) properties as well as the ionomer-carbon interactions. The transport properties of the catalyst layer are directly correlated to the porosities, pore size distributions, pore connectivity, and pore shape. Despite its simplicity and obvious limitations (e.g., we did not consider the Pt particles explicitly), this computational study provides insight into the main features of microstructure formation in catalyst layers of a PEFC. Calculated configurational properties reveal how the phase segregation into ionomer domains and carbon/Pt agglomerates occurs. The evolving structural characteristics are of vital importance for further analysis of transport of protons, electrons, reactant gas (O2), and water through catalyst layers as well as the distribution of electrocatalytic activity at Pt/water interfaces. Our study allows these properties to be related to the solvent, carbon particle choices (sizes and wettability), catalyst loading, and level of membrane hydration in the catalyst layer. Experimental results, with which we can compare our data directly, are still scarce, but we hope that this work could trigger such studies. Our findings are in qualitative agreement with experimental data on the cluster size distribution in the presence of three different solvents which show that a more polar solvent leads to smaller carbon clusters (see the Supporting Information). This principal agreement corroborates the predictive capabilities of our model. Our study highlights the importance of the applicable solvent in the structural formation process on the catalyst layer. This has important implications for the stability of the CL. Indeed, it can be expected that a large discrepancy of dielectric properties between the applicable solvent used during fabrication and water that is present in the CL during fuel cell operation will lead to the instability of the CL. Accordingly, the system specifications (e.g., type of solvent, carbon particle size, and carbon wetting properties) and parameters of the fabrication process should be considered in view of structural reorganization and irreversible performance degradation during operation of the fuel cell. Our results in combination with experimental information provide vital insights for optimization of the performance and stability of CL operation and advanced CL design. Finally, it is obviously a gross simplification to assume that interparticle interactions can be described with Lennard-Jones potentials. In the future, we plan to employ different parametrizations of pair interactions between particles in a colloidal dispersion, involving DLVO theory and ramifications thereof.48,49 DLVO-type interactions, including a screened Coulomb-repulsion term due to charging of particles (beads) in solution, will be imperative for understanding the kinetics of structure formation and for rationalizing the stability of aggregated catalyst layer structures. Despite simplifications, the present study addresses key issues such as the self-organization of carbon particles into disordered aggregates and formation of a distinct percolating phase of ionomer. We also believe that the present parametrization is sufficient for understanding the principle mechanisms of formation of stable phase-segregated morphologies in CLs.48-50 The electrostatic interactions, however, could be improved with the effective charge method developed by Gabdoulline and Wade.51 The essential parameters such as particle size and particle size distribution should be
J. Phys. Chem. C, Vol. 111, No. 36, 2007 13633 adjusted, while a more realistic parametrization refines the interaction parameters. Moreover, we plan to develop a theoretical framework for studying aggregation phenomena in multicomponent colloidal dispersions that can complement and refine our simulation approach. Conclusion We present systematic mesoscale simulations to study dynamical-structural properties of the CLs in PEFCs. The approach focuses on the structure formation of CLs. It mimics processes which occur in the fabrication process. We particularly focus on the effects of applicable solvents with varying dielectric constants and their influence in controlling the formation of microstructures with desired composite morphology. Our coarsegrained molecular dynamics simulations showed that the final microstructure moreover depends on the carbon particle choices and ionomer-carbon interactions. The important finding is that ionomer clusters do not penetrate into the carbon agglomerates. Side chains are buried inside hydrophilic domains with a weak contact to carbon domains. Ionomer domains are attached with their backbones to the surface of carbon agglomerates. They form a randomly interconnected proton-conducting network. The correlation between hydrophilic species and ionomer is significantly stronger than that between those species and carbon particles. A layer of ionomer of ∼10 nm thickness with a wellpacked morphology around Pt/C clusters was evident in the presence of polar solvents. The compactness of carbon agglomerates is affected by the polarity of the solvent. The results could be of importance in the search for advanced CL designs that could optimize the performance and stability of fuel cell operation. Acknowledgment. We gratefully acknowledge the NissanIFCI collaboration program. Supporting Information Available: Additional figures showing the final equilibrium structure of a catalyst blend (Figure S1), experimental data showing the effect of the dispersion medium on the particle size distribution (Figure S2), visualization of a hydrophilic pore (Figure S3), and twodimensional density counter plots of Nafion and carbon particles (Figure S4). This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Handbook of Fuel Cells: Fundamentals, Technology and Applications; Vielstich, W., Lamm, A., Gasteiger, H. A., Eds.; VCH-Wiley: Weinheim, Germany, 2003. (2) Weber, A. Z.; Newman, J. Chem. ReV. 2004, 104, 4679. (3) Eikerling, M. J. Electrochem. Soc. 2006, 153, E58. (4) Eikerling, M.; Kornyshev, A. A.; Kucernak, A. R. Phys. Today 2006, 59, 38. (5) Eikerling, M.; Kornyshev, A. A. J. Electroanal. Chem 1998, 453, 89. (6) Eikerling, M.; Ioselevich, A. S.; Kornyshev, A. A. Fuel Cells 2004, 4, 131. (7) Eikerling, M., Kornyshev, A.A.; Kulikovsky, A.A. Physical Modeling of Cell Components, Cells and Stacks. In Encyclopedia of Electrochemistry, Volume 5, Electrochemical Engineering; Macdonald, D. D., Ed.; in press. (8) Gode, P.; Jaouen, F.; Lindbergh, G.; Lundblad, A.; Sundholm, G. Electrochim. Acta 2003, 48, 4175. (9) Uchida, M.; Aoyama, Y.; Eda, E.; Ohta, A. J. Electrochem. Soc. 1995, 142, 463. (10) Uchida, M.; Aoyama, Y.; Eda, E.; Ohta, A. J. Electrochem. Soc. 1995, 142, 4143. (11) Uchida, M.; Fuuoka, Y.; Sugawara, Y.; Eda, N.; Ohta, A. J. Electrochem. Soc. 1996, 143, 2245.
13634 J. Phys. Chem. C, Vol. 111, No. 36, 2007 (12) Jia, N.; Martin, R. B.; Qi, Z.; Lefebvre, M. C.; Pickup, P. G. Electrochim. Acta 2001, 46, 2863. (13) Yang, T. H.; Yoon, Y. G.; Park, G. G.; Lee, W. Y.; Kim, C. S. J. Power Sources 2004, 127, 230. (14) Fischer, A.; Jindra, J.; Wendt, H. J. Appl. Electrochem. 1998, 28, 277. (15) Poltarzewski, Z.; Staiti, P.; Alderucci, V.; Wieczorek, W.; Giordano, N. J. Electrochem. Soc. 1992, 139, 761. (16) Wilson, M. S.; Valerio, J. A.; Gottesfeld, S. Electrochim. Acta 1995, 40, 355. (17) Uchida, M.; Fukuoka, Y.; Sugawara, Y.; Ohara, H.; Ohta, A. J. Elelectrochem. Soc. 1998, 145, 3708. (18) Uchida, H.; Song, J. M.; Suzuki, S.; Nakazawa, E.; Baba, N.; Watanabe, M. J. Phys. Chem. B 2006, 110, 13319. (19) Wang, Q.; Eikerling, M.; Song, D.; Liu, S. J. Electroanal. Chem. 2004, 573, 61. (20) Fernandez, R.; Ferriera-Aparicio, P.; Daza, L. J. Power Sources 2005, 151, 18. (21) Wang, C. Y. Chem. ReV. 2004, 104, 4727. (22) Springer, T.E.; Wilson, M.S.; Gottesfeld, S. J. Electrochem. Soc. 1993, 140, 3513. (23) Sahimi, M. Heterogeneous Materials II; Springer: Heidelberg, Germany, 2003. (24) van Santen, R. A. Nature 2006, 444, 46. (25) Torquato, S. Random Heterogeneous Materials: Microstructure and Macroscopic Properties; Springer: Heidelberg, Germany, 2000. (26) Wescott, J. T.; Qi, Y.; Subramanian, L.; Capehart, T. W. J. Chem. Phys. 2006, 124, 134702. (27) Galperin, D. Y.; Khokhlov, A. R. Macromol. Theory Simul. 2006, 15, 137. (28) Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423. (29) Vishnyakov, A.; Neimark, A. V. J. Phys. Chem. 2001, 105, 9586. (30) Jang, S. S.; Molinero, V.; Cagin, T.; Goddard, W. A., III. J. Phys. Chem. 2004, 108, 3149.
Malek et al. (31) Eikerling, M.; Kornyshev, A. A.; Stimming U. J. Phys. Chem. 1997, 101, 10807. (32) Vishnyakov, A; Neimark, A. V. J. Phys. Chem. 2000, 104, 4471. (33) Spohr, E.; Commer, P.; Kornyshev, A. A. J. Phys. Chem. 2002, 106, 10560. (34) Mologin, D. A.; Khalatur, P. G.; Kholhlov, A. R. Macromol. Theory Simul. 2002, 11, 587. (35) Khalatur, P. G.; Talitskikh, S. K.; Khokhlov, A. R. Macromol. Theory Simul. 2002, 11, 566. (36) Marrink, S. J.; de Vries, A. H.; Mark A. E. J. Phys. Chem. B 2004, 108, 750. (37) Izvekov, S.; Violi, A. J. Chem. Theory Comput. 2006, 2, 504. (38) Izvekov, S.; Violi, A. J. Phys. Chem. B 2005, 109, 2469. (39) Izvekov, S.; Violi, A.; Voth, G. A. J. Phys. Chem. B 2005, 109, 17019. (40) Yamamoto, S.; Hyodo, S. A. Polym. J. 2003, 35, 519. (41) Kinoshita, K. Carbon, electrochemical and physicochemical properties; John Wiley & Sons: New York, 1988; p 150. (42) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306. (43) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43. (44) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33. (45) Kisljuk, O. S.; Kachalova, G. S.; Lanina, N. P. J. Mol. Graphics 1994, 12, 305. (46) Cui, S.; Liu, J.; Selvan, M. E.; Keffer, D. J.; Edwards, B. J.; Steele, W. V. J. Phys. Chem. B 2007, 111, 2208. (47) Lamas, E. J.; Balbuena, P. B. Electrochim. Acta 2006, 51, 5904. (48) Ninham, B. W. AdV. Colloid Interface Sci. 1999, 83, 1. (49) Coniglio, A.; De Arcangelis, L.; Del Gado, E.; Fierro, A.; Sator, N. J. Phys.: Condens. Matter 2004, 16, S4831. (50) Glotzer, S. C.; Solomon, M. J.; Kotov, N. A. AIChE J. 2004, 50, 2978. (51) Gabdoulline, R. R.; Wade, R. C. J. Phys. Chem. 1996, 100, 3868.
