Mesoscale Simulation of Proton Transport in Proton Exchange

Article pubs.acs.org/JPCC

Mesoscale Simulation of Proton Transport in Proton Exchange Membranes Ryan Jorn and Gregory A. Voth* Department of Chemistry, James Franck Institute and Computation Institute, University of Chicago, Chicago, Illinois 60637, United States, and Computing, Environment, and Life Sciences Directorate, Argonne National Laboratory, Argonne, Illinois 60439, United States

ABSTRACT: Previous efforts to model proton transport through fuel cell membranes have largely focused on disparate length scales: molecular dynamics at the atomistic level and fuel cell stack engineering approaches at the macroscale. A new multiscale approach to bridge these extremes is proposed in this work which combines concepts from coarse-grained (CG) modeling with smoothed particle hydrodynamics (SPH) to capture the qualitative morphology and transport behavior of a proton exchange membrane at the length scale of tens of nanometers. This method allows for connection to atomistic simulations via the inclusion of transport coefficients from molecular dynamics and coarse-grained forces derived for the polymer backbone, side chain, proton, and water interactions. Information pertaining to macroscopic conductivity is obtained by volume averaging based on the flux and chemical potential fields within the membrane. Proton transport is effectively coarse-grained via introduction of a composition variable associated with each interacting site which carries the field information. By combining this technique with local electrostatics and coordinate dependent diffusion constants, the effects of double layer formation within the water pores and the influence of proximity to sulfonate groups on transport is recovered. The combined CG-SPH method is validated and subsequently applied to an equilibrated hydrated Nafion structure with a box length of 40 nm. The resulting conductivities calculated for the material agree very well with trends from experiment and provide insight into the complex interplay of morphology, proton distribution, and diffusion coefficients at a length scale that can be expanded beyond feasible atomistic molecular dynamics simulations to capture the effects of mesoscopic morphology on proton conduction.

I. INTRODUCTION The need for alternatives to fossil fuels is becoming increasingly compelling as the impact of their consumption on air quality,1 global climate change,2 and energy security3 continues to be manifested worldwide. While the need for new energy technologies is apparent, optimizing the materials used to build “greener” systems remains a significant challenge. The incorporation of electrochemical systems to address modern energy demands has yet to achieve energy densities and power output comparable to gasoline combustion at a competitive cost.4 Fuel cells in particular have been promoted for replacing the internal combustion engine in automobiles;5 however, their successful implementation requires engineering a complex array of catalytic and transport processes.6−9 Fuel cells are powered by the redox chemistry of hydrogen and oxygen gases which are introduced at either end of a layered structure consisting of gas diffusion electrodes, catalyst layers that facilitate the redox © XXXX American Chemical Society

reactions, and an ion-conducting electrolyte that separates the cathode and anode assemblies.4 The continued improvement of fuel cell technology thus depends on optimizing the properties of each component in this stack, as well as their interfaces.6,10 Improvement of the ion-conducting electrolyte requires advancement of materials with high ion conductivity, low permeability to reactant gases, negligible electronic conductance, and chemical and mechanical stability.6,10−12 Proton exchange membranes (PEMs) have been a popular choice for realizing these features, of which Nafion is considered the standard for its reasonably high proton conductivity and mechanical stability.13 Nafion consists of a polytetrafluoroethylene (PTFE) backbone augmented with sulfonic acid side Received: January 2, 2012 Revised: March 2, 2012

A

dx.doi.org/10.1021/jp300040w | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

domains can be connected by narrow water wires;30 however, fluctuations in the polymer backbone may cause these connections to rapidly vary in time, providing a dynamic conduction pathway.33 In general, the structure of the water network is also influenced by the length of the sulfonated side chains, the equivalent weight of the polymer, and the dispersity of side chains along the PTFE backbone.28,31,32,37 Clearly a great deal of structural information has been gleaned from atomistic-level simulations. However, to correctly model proton transport at larger length scales, the computational limitations of molecular dynamics, i.e., system size and length of the simulation, become problematic.12 Pore sizes are estimated to be on the order of nanometers; thus, using simulation boxes of a similar scale will not allow for an accurate representation of the tortuosity and porosity of the hydrated morphology, both of which can influence the effective proton diffusion.54 Furthermore, for large simulations with polymer strands representative of the actual molecular weights observed in experiment, e.g., 106 g/mol,55 it is not possible to fully equilibrate the structure on a computationally feasible time scale. Thus, it is necessary to combine atomistic information and data from experiment with more CG approaches to attain realistic membrane structures which allow for accurate predictions of conductivity. Several CG approaches have been implemented to overcome the limitations of atomistic methods and capture pore size distributions.38−48 These methods provide qualitative agreement with experiment and the results from MD simulations; i.e., nanophase segregation is observed, and the cluster sizes agree with reported scattering profiles. Dissipative particle dynamics (DPD) is of particular relevance to the work presented in this article and has been applied to the phase separation in diblock copolymers56,57 as well as hydrated Nafion.40,43−47 DPD is based on a Brownian-type equation of motion in which stochastic and dissipative terms act between particle pairs to drive the system toward thermal equilibrium and thus produce an NVT ensemble.57−60 While the pairwise force for DPD is typically chosen to be of a simple linear repulsion term,57 it has recently been compared with a more detailed CG force field.39 Efforts have also been made to connect DPD with kinetic Monte Carlo to study water and gaseous species transport within equilibrated mesoscale structures;45,61,62 however, extension of this work to the study of proton transport on the mesoscale has not yet been addressed. The use of theoretical and computational methods to improve fuel cell PEMs requires not only an accurate method for treating the local and global structure of the polymer membrane but also proper treatment of the proton transport within the context of a reasonable morphology. Proton transport in water is often discussed in terms of vehicular conduction in addition to structural diffusion, also known as Grotthus shuttling.63 Vehicular diffusion can be captured by standard MD within linear response theory64 and corresponds to the translation of the entire excess proton−water complex in a fashion common to all mobile species in solution.21,23,24,36 Paddison et al. have extended these studies to account for higher order nonequilibrium effects and have reported transport coefficients which vary with distance from the pore walls in the case of simple lamellar geometries.65−67 While insightful at low hydration levels, these methods are questionable with increasing water content, since experiments suggest that the Grotthus shuttling mechanism becomes

chains attached randomly along the polymer. The membrane thus consists of a moiety with hydrophilic and hydrophobic groups, the sulfonic acid groups and the Teflon backbone, respectively, which undergo microphase separation when exposed to water.13 Decades of research have been invested into understanding the morphology of Nafion and its impact on the proton conduction mechanism with the motivation of driving future PEM design. While agreement has been reached regarding general features of the percolation behavior with increasing water content, an outstanding issue remains the precise morphology of the hydrated Nafion domains during proton conduction. At low water contents (λ < 5), where λ is the number of water molecules per sulfonate group, the sulfonated side chains cluster together with the added water to form isolated hydrophilic regions. As the water content increases, the hydrophilic domains expand to form spanning water channels which are capable of efficiently transmitting protons. Measurements from a variety of neutron and X-ray scattering experiments have been used to propose a number of models for the shape of the water domains and their expansion behavior with increasing hydration, including lamellar layers of polymer and water,14,15 cylinders of water surrounded by polymer chains as well as cylinders of polymer surrounded by water domains,16,17 and spherical water clusters connected by narrow channels,18 among several other structures.13,19 Given the ambiguity of the experimental data, several simulation efforts have attempted to provide a more detailed understanding of the morphology of the water network in perfluorosulfonic acid (PFSA) membranes in general and Nafion in particular.20 These efforts have relied on different levels of resolution from fully atomistic molecular dynamics (MD),21−32 in which all the atoms in the polymer chain and solvent are treated explicitly, to united atom approaches in which the CF2 groups in the polymer are replaced with single interaction sites,33−37 and finally more coarse-grained (CG) models in which multiple carbons are grouped into collective degrees of freedom along with their substituents.38−48 From the former two categories, detailed information has been obtained for the water structure surrounding the sulfonate groups and the behavior of the hydrophilic domains with increasing water content. Not surprisingly, it has been shown that water favors coordination with the headgroup of the side chain and the side chains tend to be oriented perpendicular to the hydrophobic polymer backbone.34−36 At low hydration levels, sulfonate groups share water and hydronium ions through bridging structures which allow for the formation of small hydrophilic clusters in spite of the repulsion between the side chains.22−24,34−36 Ab initio molecular dynamics simulations49−53 have explored the threshold for dissociation of the proton from the sulfonic acid headgroup and demonstrated that for water contents around λ = 3 the chemical bond is replaced by a sulfonate−proton contact ion pair. With increased water content, the sulfur groups move farther apart and the protons become more loosely associated with the conjugate base side chain. For water contents above λ = 9, a substantial fraction of protons in the system are found on average beyond the first hydration shell of the sulfonate group.22,23 In regard to the water cluster sizes and distribution, MD has shown that at low hydration the morphology consists of isolated clusters of tens of water molecules, which swell to hundreds of water molecules with increasing water content and yield estimates for the percolation threshold between five and six waters per sulfonate.26,36 At intermediate hydration, the hydrophilic B



Article

accomplished here within the context of our previous work on the Nafion PFSA membrane30,71,72 along with work from our group on the composition-dependent mechanics of biological membranes.84,85 We make use of a CG model to describe the water network of Nafion and combine it with nonequilibrium thermodynamics to describe proton conduction through the resulting complex heterogeneous structure. Our method is distinct from previous attempts in that we do not restrict the system to ideal geometries or the simulation size to the scale of individual pores. This method is therefore a CG model of the PEM which can also describe the proton flow and dynamically influence transport. By connecting to results from previous MSEVB calculations, we also include the fundamental physics of transport at the microscale. In section II, we introduce our combined CG approach with smoothed particle hydrodynamics (SPH) and provide its validation for a simple pore geometry. In section III, we discuss the morphologies derived from CG DPD simulations and the resulting conductivities with comparison to the previous literature. In section IV, we conclude by describing several future extensions of this work and discuss areas for improvement of the CG-SPH method.

increasingly relevant as the protons become more loosely associated with the sulfonate groups.68 Accounting for structural diffusion requires a method capable of describing changes in bonding topology between the water and mobile protons and hence lies outside the scope of standard molecular dynamics.12,20,63 Several ab initio studies have investigated the effect of λ and chain length on the dissociation of hydrogen from individual sulfonate groups in Nafion;49−53,69,70 however, scaling these simulations up to the size of pores with flexible polymer backbones is computationally intractable. Alternatives which have been applied successfully at the pore level to account for Grotthus shuttling of protons include empirical valence bond theory (EVB),71−75 reactive force fields,29 and stochastic hopping algorithms within standard MD.26,76,77 Multistate EVB (MS-EVB) is a particularly attractive method, since the model parameters can be derived directly from quantum simulations and the method has already been successfully applied to model transport through disordered Nafion−water networks with highly concentrated proton environments.71,72 Unfortunately, the computational cost associated with such simulations remains prohibitive to application on mesoscopic simulation scales. Modeling proton transport at length scales beyond those feasible with atomistic methods has largely relied on continuum mechanics and nonequilibrium thermodynamics. Macroscopic conservation equations have been integrated to study water transport within PEM fuel cells and to evaluate cell performance via polarization curves.8 The proton dynamics giving rise to the observed currents can be studied by deriving simple models for the transport coefficients78 and incorporating them into the Poisson−Nernst−Planck theory for electromigration of ions, which accounts for convective flow from electrophoretic effects as well as diffusive motion.8,79,80 The Poisson− Boltzmann method has been applied to describe equilibrium proton distributions within cylindrical pores and has been discussed in the limit of rapid diffusive motion with respect to the established flow.81−83 The disadvantage of these approaches is that they require as input transport coefficients and assume idealized membrane structures. In addition to these challenges, there is usually not any clear connection with atomistic information to form a multiscale bridge for studying charge transport. Fuel cell engineering approaches, which incorporate the entire membrane electrode assembly, also neglect the detailed structure of the membrane itself and replace it with material parameters which are often either fit from experiment or derived from simple arguments. There is no clear connection to the detailed structural information provided by the atomistic methods discussed previously, and attempts to reintroduce membrane structure often rely on highly averaged quantities. Hence, the ability for these continuum approaches to optimize materials for future PEM design and predict performance is rather limited. Clearly there is a need for new mesoscopic methods that can connect realistic membrane morphologies with the correct physics guiding the proton transport process in order to help inform PEM design and advance fuel cell research. Such methods must be sufficiently coarse-grained in their description of the system such that mesoscale structural changes can be included at reasonable computational cost. Our goal in this work is to lay the foundation for such a multiscale method that can connect information from atomistic MS-EVB simulations with mesoscopic descriptions of both the membrane morphology and the proton conduction process. This effort will be

II. CG-SPH METHOD The model PEM considered in this work is the hydrated Nafion 117 membrane in which the PTFE polymer has an equivalent weight of 1100 g/(mole SO3). For simplicity, the sulfonate side chains are assumed to be equally spaced which allows the implementation of the large-scale coarse-graining strategy of Wescott et al.38 The repeat unit for the polymer strand is shown in Figure 1. In the CG model, the backbone is divided

Figure 1. Schematic representation of the CG model of Wescott et al.38 in which the atomistic repeat unit for Nafion 117 is replaced with a simple bead representation including a side chain bead and two polymer backbone beads. In the atomistic representation, the turquoise atoms are carbons, the white correspond to fluorines, the red are oxygens, and the yellow is a sulfur atom.

into two [CF2CF2]4 groups which are replaced by single interaction sites, and the entire sulfonate side chain is replaced by a single CG particle. As discussed previously,38 the CG beads have a volume of 0.315 nm3 and a diameter of 0.844 nm. In order to account for the presence of water in the hydrated polymer, water CG particles of the same volume are also introduced which contain on average 10.5 water molecules. DPD simulations were performed using the parameters in C



Article

Table 1 in the LAMMPS molecular dynamics package86 to model random water networks for subsequent conduction

variables and has been shown to establish thermal equilibrium at sufficiently long times:58,60 σ=

Table 1. Parameters Used in the DPD Simulations of Hydrated Nafion in Accord with the Model of Wescott et al.38 parameter Δt (fs) T (K) rc (nm) C γ [(kJ fs)/(mole nm2)] aWW [(kJ fs)/(mole nm2)] aWP [(kJ fs)/(mole nm2)] aWS [(kJ fs)/(mole nm2)] aPS [(kJ fs)/(mole nm2)]

× 102 × 10−1 × × × × ×

Fijspring = C

104 102 102 102 102

∑ (Fijc + Fijr + Fijd + Fijspring ) (1)

j

where the total force on particle i is written as a sum over interactions with the remaining j particles in terms of a conservative linear repulsion term, Fcij, a random force, Frij, a dissipative frictional force, Fdij, and a bonding term for the CG particles which comprise the polymer chains, Fspring . For ij distances less than the defined cutoff, rc, the expression for the first three terms is given by57

aij = aii + 3.096

(4)

RT χ rc ij

(5)

Since there are three beads in the model, there are three unique Flory parameters which can be evaluated on the basis of cohesive energies for the individual materials. The interaction coefficients chosen in this study are listed in Table 1 and were selected on the basis of atomistic MD calculations of the solubility parameters.38 The relative magnitudes of these interactions account for the hydrophobic interaction between the fluoronated backbone and water, the favorable interaction between the sulfonate side chain and water, and the repulsive interaction between the backbone and the side chain. It is interesting to note that in previous DPD studies interaction parameters have been used which reflect two different possibilities for the relative magnitudes of the side chain− backbone interaction, aSP, and the water−backbone interaction, aWP. The selection of a stronger side chain−PTFE repulsion as compared to the water−PTFE interaction39,43,44 (aWP < aSP) and the opposite case40,45−47,61,62 (aWP > aSP) both produced plausible structures emphasizing the need for additional connections with atomistic simulation and the need for continual refinement of the CG models. As such, the structures generated by the model presented here are not intended to be a

Fijc + Fijr + Fijd = aijw(rij)riĵ + σw(rij)ξij(Δt )−1/2 riĵ − γw(rij)2 (riĵ ·vij)riĵ

RT rij rc 2

where C is a constant that is adjusted to give an average bond length that matches the first peak in the radial distribution function from simulations of nonbonded particles.57 The values for the self-interaction coefficients in eq 2, aii, are determined by matching CG simulation results for the pressure as a function of bead density to the compressibility of the material. Simulations with water beads were performed, and the pressure was found to depend on the square of the bead density, in agreement with the previously cited DPD studies. The aWW parameter was extracted from a quadratic fit of the pressure versus density data for water particles and subsequently adjusted to give the correct compressibility of water at the bead density specified above, i.e., ρ = 3/rc3. Care must be taken in this process to account for how many water molecules are contained in the CG particles,87 as failure to do so results in matching the compressibility of the water beads to experimental values rather than the water itself. The resulting repulsion term will then be too soft and yield incorrect density fluctuations. For simplicity, the value of the self-interaction coefficient obtained for water (aWW) was used for each of the other particle types, as in previous polymer simulations.40,43−47,61,62 When the CG particles correspond to different chemical entities, water and the Teflon backbone, for example, Flory parameters, denoted χij, are introduced to provide the cross interaction coefficients according to57

studies. Polymer chains with 15 of the repeat units from Figure 1 were created by a random walk algorithm, and water beads were randomly inserted into the simulation box to achieve the desired water content. In each simulation, 198 000 coarsegrained beads were included and the structure was equilibrated for 2.5 million time steps using the Intrepid Blue Gene P supercomputer at the Argonne Leadership Computing Facility. A brief overview of the DPD methodology is provided subsequently, and the interested reader is referred to the cited literature for further details.57−60 There are three types of CG beads in our model system: water, side chain, and polymer backbone, which will be referred to in the following notation as W, S, and P particles, respectively. The CG particles interact via the following force equation: Fi =

(3)

where R is the gas constant and T is the desired temperature of the system. For values of rij greater than the cutoff, the weighting function in eq 2 is equal to zero and the first three terms in eq 1 do not contribute to the motion of the particle. The bonding term in eq 1 is chosen as a simple linear restoring force akin to a Gaussian chain:

value 10.0 3.00 9.81 3.00 9.95 7.05 7.41 7.05 7.61

2RT γ

(2)

where w(rij) = (1 − rij)/rc is a weighting function, rij is the distance between the two particles, r̂ij is the unit vector pointing from particle j to particle i, and vij is their relative velocity. The cutoff distance is chosen following previous studies by Groot and co-workers56,57,87 and Wu et al.43,44 to fulfill the relationship ρrc3 = 3, which corresponds to a density of 3 beads per rc3. This density regime was shown to be convenient for reproducing bulk liquid properties with a minimal cutoff range for the DPD force terms and results in a value of roughly a nanometer for the DPD cutoff in the current simulations. There are three additional parameters in eq 2 which must be included to control the thermostat and the interaction energies between beads i and j, i.e., aij, σ, and γ. The following relationship has been used to eliminate one of the thermostat D



Article

dC+ = −C+∇⃗ ·vwater − ∇⃗ ·J+diff dt

definitive identification of the morphology of hydrated Nafion but rather provide plausible structures that demonstrate the behavior of the method to describe transport through complex networks at the mesoscale with future opportunities to connect to more accurate representations of the morphology. As mentioned at the outset, having a structure for the water network in the membrane system from DPD simulations is only one piece of developing a mesoscopic model for PEM proton conduction; one must also properly account for the migration of protons through the system.88 It has been argued previously that continuum type approaches to diffusion are justified in molecular systems as long as the pore size is greater than 2 D screening lengths,80,89 a criteria that is likely satisfied for protons in the highly acidic Nafion environment as a result of its nanometer-sized pores and highly ionic environment.13,90 While transport has recently been studied at the nanometer length scale using finite element methods for idealized cylindrical pores,79,80 a second approach which is more amendable to connection with molecular and CG dynamics is smoothed particle hydrodynamics (SPH).91−94 While SPH has its origins in astrophysics and star modeling,95,96 it has since been successfully applied to a variety of problems such as shock wave propagation,97 diffusion and flow in porous media,98−101 heat conduction, 102 and a variety of other fluid flow problems.92,93 Beyond providing an additional technique for solving the differential equations of fluid dynamics, the fundamental role of the quasi-particles used in this method has also been discussed in connection with DPD103 and as a conceptual bridge between the microscopic reversible molecular dynamics and the macroscopic irreversible mechanics of liquids.91 Some of the advantages of SPH as a numerical method include the similarity in propagation algorithm to standard MD, which allows for efficient parallelization based on linked-cell methods,86 ease of incorporating complex boundary conditions which are also time dependent, and description of the fluid motion in the Lagrange frame of reference. As with any numerical method, there are also challenges that have to be overcome including the proper treatment of field values at the boundary, identification of sufficient particle support for the weighting functions, and optimizing the weighting functions used in the calculations to maximize numerical accuracy.92 Below, it shall be shown that DPD morphologies can be used in synergy with SPH to provide a mesoscopic model for proton transport. The key connection is made by mapping the location of water DPD beads to the interpolation quasi-particles used in the SPH algorithm. The description here of proton transport is based on nonequilibrium thermodynamics within SPH and begins with considering the mass conservation equations for the three entities in our simulation: polymer, water, and protons. To model transport, we consider their respective fluxes denoted as Jpoly, Jwater, and J+. In the following discussion, the flow of polymer will be neglected because of the temperature range considered such that the polymer strands can be treated as fixed during proton transport. Proton flux can be partitioned into a diffusive component arising from dispersion of the species and a convective flux arising from fluid flow with an averaged velocity denoted vfluid. The fluid flow arises from water transport within the membrane which is an additional mass transport problem that has attracted a great deal of interest.8,104 Changes in proton concentration within regions of the water network can then be written as

(6)

where a material derivative is implied in eq 6 on the left-hand side and the diffusive flux is denoted Jdiff + . In principle, the solution of eq 6 requires integration of the coupled equation for water flow under the influence of proton motion and a pressure gradient generated by water creation at the cathode and depletion at the anode. If the effects of electroosmosis and pressure-driven flow are neglected, the proton transport equation is reduced to diffusive motion and the superscript on the flux in eq 6 will be dropped in the subsequent discussion. Within linear response theory, the diffusive flux can be written in terms of the motion of individual protons and the local gradient of the chemical potential in the membrane as105,106 J+ =

1 (∑ 3VRTNA k , k ′

∫0

∞

⟨vk ·vk ′(τ)⟩ dτ)∇⃗μ+ (7)

where the quantity preceding the gradient is known as the Onsager coefficient where V is the volume of the system, NA is Avogadro’s number, vk (vk′) denotes the velocity of individual protons k and k′, μ+ is the chemical potential for the proton at a given point, and the brackets ⟨···⟩ correspond to an ensemble average. The evaluation of the correlation function in eq 7 requires atomistic information which can be obtained from MD simulations in conjunction with MS-EVB methods to account for vehicular diffusion and Grotthus hopping.71,72 For the moment, only the dilute limit for this expression is considered: J+ = −

C+D+ ∇⃗μ+ RT

(8)

where the proton tracer diffusion constant, also known as the self-diffusion constant, is given as D+. This quantity has been studied in detail by previous MD simulations,21,23,24,36 and recently with the inclusion of MS-EVB methods,71,72 where good agreement was found with previous NMR, conductivity, and scattering experiments.107−110 Rather than consider the diffusion coefficient as a global “constant” for the system, there is increasing evidence from atomistic simulations that these transport coefficients should be treated as dependent on distance from the pore walls in Nafion and similar porous materials.21,65,67,111 The physical origins for including coordinate-dependent diffusion coefficients arise from the formation of a Stern layer at the pore wall composed of sluggish water molecules as well as the influence of the local water structure and sulfonate side chains on the Grotthus hopping mechanism within the channel. Recently, MS-EVB simulations have demonstrated that proximity to the sulfonate groups in proton diffusion results in a complex and anticorrelated interplay of vehicular transport and Grotthus shuttling.71,72 To recapture these effects on proton dynamics, a simple model for the diffusion constant is implemented which is heuristically related to previous atomistic MS-EVB simulations. A linear relationship between distance to the nearest sulfonate group and the effective diffusion coefficient was assumed and parametrized by comparing averaged diffusion coefficients from lamellar geometries to results from similar water contents in MS-EVB simulations.72 This approach results in the following expression for the proton diffusion coefficient: D+(r ) = αr + D+0 E

(9)



Article

where α = 2.69 nm/ns, D0+ = 2.69 × 10−2 nm/ns, r is the distance to the nearest sulfonate group in units of nanometers, and D is given in units of nm2/fs. In this work, the diffusion constants of ref 29b were also used when the coordinate dependence of the transport coefficient was neglected for comparison. The highly charged environment of the Nafion pores not only influences the local transport properties of the water channels but also the local chemical potential experienced by the protons as they traverse the membrane. While electroneutrality can be enforced in a global sense, at the level of individual pores, double layers of protons distributed around the sulfonate groups will create anisotropic proton disbritubions. In order to incorporate the local differences in electrostatic potential experienced by each proton as it traverses the membrane, the so-called quasi-electrochemical potential is used which expresses chemical potential as a combination of the infinite dilution entropic term as well as an interaction with the local charge densities:112 μ+ = RT ln C+ + F Φ

common problem of incorporating electrostatics in MD methods must also be addressed in these simulations. The potential-shifting method of Wolf et al. has been described as an effective replacement to Ewald summation in certain cases and can also be treated as an effective means to describe electrostatics in a screening environment.115 This potentialshifting approach was implemented in our code with a cutoff between 4 and 6 nm as determined by numerical convergence for the lamellar pore structures. Given the general equations for proton flux and expressions for the quantities entering eq 8, the SPH method can now be applied to generate the appropriate equations of motion for the field variables. The general concept in SPH is to discretize the velocity, density, and energy fields from the conservation equations of continuum mechanics using a set of interpolation points, or quasi-particles, at which the field values are known in conjunction with smoothing functions to provide the field values continuously throughout the problem domain. As an example, consider an arbitrary field variable evaluated at position qj, A(qj):

(10)

A(q j) =

where the first term corresponds to the noninteracting diffusion limit, F is Faraday’s constant, and Φ is the electrostatic potential. Calculating the electric potential, Φ, requires knowledge of the charge distribution in the Nafion pore and presents a unique challenge to mesoscopic calculations as compared with traditional MD, since the exact ion positions have been subsumed into CG entities. Recent efforts to add electrostatics to DPD simulations have addressed CG potentials by assuming a model distribution for the charges in each particle. Two examples that have been considered are a linear concentration drop from the center of the particle113 and an exponential decay of charge density from the center of the particle.114 Both forms provide an analytic result for the bead−bead interaction energy which is amendable to Ewald summation and particle mesh methods. For the moment, a simpler case for the bead concentration profile is considered: the sulfonate side-chains are given a unit negative charge uniformly distributed across the volume of the particle and similarly the water beads contain a uniform positive charge proportional to the local proton concentration. While temporal coarse-graining and the charge delocalization of the protons contribute to an effective smearing out of these charges across the coarse-grained particles, assuming a uniform distribution is clearly an idealization that will improve only with decreasing particle size. Given that our present motivation is to understand the qualitative role of anisotropic proton distributions in the Nafion pore, this assumption can be treated as a first approximation subject to future refinement. The electrostatic potential for a sphere of uniform charge is identical to a point-charge located at its center for points outside of the sphere, and hence, calculation of the local electrostatic potential at each SPH particle position requires a simple superposition of the potential arising from each neighboring bead in addition to a self-interaction term originating from the charge inside the bead itself. Within this approach, each proton interacts with the average field of the other protons in the system and hence some of the activity of the protons in the acidic environment is retained in a meanfield sense, but the explicitly correlated atomistic motion is not. Since a form for the electrostatic potential is recovered which still contains the long-range nature of Coulomb’s law, the

≈

∫ A(q′)δ(q j − q′) dq′ ∫ A(q′)W (q j − q′) dq′/∫ W (q j − q′) dq′ (11)

where the weighting or smoothing function W(q − q′) behaves just as the entity in eq 2, but in practice has a different form, and a correction term has been included in the denominator to normalize the result. With the introduction of the smoothing function, the approximation in the second line becomes exact in the limit that W(q − q′) approaches a delta function. In practice, the weighting function has a finite support length, or cutoff distance, which means that it approaches zero for sufficiently long particle separations and is strictly zero beyond the cutoff, similar to the weighting function in DPD. The smoothing function will be chosen to depend only on the distance between points, and hence, W(q − q′) will be replaced by W(|q − q′|). In what follows, a quintic spline will be used for the smoothing function due to its stability for random particle distributions.93 The cutoff for this function was varied to confirm numerical convergence of the resulting concentration profiles. The integral in eq 11 can be evaluated numerically by a Riemann sum in which the q space is discretized on a set of grid points: A(q j) ≈

∑ A(q i)W (|q j − q i|)Δq i/∑ W (|q j − q i|)Δq i i

i

(12)

where the volume element Δqi can vary at each point and can be rewritten in terms of the density and mass of the material: m m A(q j) ≈ ∑ i A(q i)W (|q j − q i|)/∑ i W (|q j − q i|) ρ ρi i i i (13)

General schemes have been proposed for calculating gradients and higher order derivatives of the field variable using Taylor expansions which correct the SPH equations in the presence of boundaries to within first order accuracy.92,116 These so-called corrected SPH approaches will be used in this work to calculate the proton fluxes, i.e., gradients of chemical potential, given by eq 8. While this method has been extended to higher order derivatives, Monaghan has pointed out a F



Article

thickness of the water layer and the hydration level of the membrane, resulting in linear expansion of the water layer with increasing water content. For the lowest water content considered, λ = 5.25, the water layer was 0.844 nm wide, corresponding to one SPH particle diameter between the polymer sheets. At the highest water content considered, the pore width was 3.376 nm, or four SPH particle diameters. An example structure is shown in the inset of Figure 2a with the

simpler algorithm based on the integral approximant of Brookshaw117 which can be used to evaluate second order derivative expressions relevant to transport with comparable accuracy and less computational effort: ∇⃗ ·(A(q)∇⃗B(q)) ≈

∫ (A(q) + A(q′))(B(q) − B(q′)) F(|q − q′|) dq′

(14)

where A and B are arbitrary functions and F(|q − q′|) is related to the weighting function of eq 13 via ∇⃗ W(|q − q′|) = (q − q′) F(|q − q′|). As previously discussed by Chen et al.,116 the accuracy of the SPH approximation for integrals is improved by inclusion of a correction term accounting for particle inconsistencies which to lowest order is simply the normalization factor given in eq 11 for W. This concept can be extended to eq 14 by including a division by the integral of the gradient of the smoothing function to provide increased numerical accuracy in the presence of boundaries. Taken together with eqs 6, 8−10, and 14, the final transport equation for the evolution of concentration with respect to time is given by dC+i = dt

∑ j

mj RT ρj

(C+j D+j + C+i D+i )(μ+i − μ+j )Fij/∑ j

mj ρj

Fij (15)

Transport equations based on eq 15 were found to be quite stable and accurate, even in the presence of significant particle inconsistency at the boundaries of the water domain within model pores as demonstrated below. Equation 15 describes the concentration dynamics of the system as composition exchanges between SPH particles i and j based on their differences in proton chemical potential and the average of their diffusion coefficients. The connection between the SPH treatment and the DPD model discussed above arises from first examining the morphologies from the DPD trajectory and then substituting the water CG beads with SPH particles of the same size. In the studies involving the complete electrostatics of the system, the equilibrated proton distributions are obtained by initializing the simulation with local proton density from sulfonate groups donating their proton to the nearest water. Starting from these initial concentration profiles, eq 15 was then integrated for the water particles using a linked-cell method within a self-contained code developed in our group while keeping the polymer structure static. Ongoing work in our group is continuing to develop the CG-SPH method to describe the effect of fluctuating polymer backbones. Proton conductivity was calculated using established volume averaging techniques to explore the effective transport properties.118−120 The membrane is broken up into smaller pieces at steady state, and the superficial volume average of the flux is divided by the derivative of the intrinsic volume average of the chemical potential to obtain the proton conductivity. It is anticipated in future work that these averaged quantities can then be bridged into engineering models for the entire fuel cell stack to complete the multiscale connection from atomistic simulation to macroscale modeling.

Figure 2. Comparison of results for the proton concentration (a) and proton flux (b) for a lamellar model (inset) in which the concentration is fixed at either end of the structure. The analytical (lines) and SPH results (symbols) are compared at 5 ns (black solid lines), 15 ns (red dashed lines), and 50 ns (blue dotted lines) for SPH smoothing lengths of 2 (triangles) and 3 (solid circles) times the average nearestneighbor water particle distance.

colors corresponding to the model in Figure 1 and the water layer represented by regularly spaced SPH particles on a lattice of points. The smoothing length chosen for these simulations, and the subsequent investigation of transport through DPD Nafion structures, was a multiple of the average nearestneighbor water particle distance obtained from the radial distribution function. In the lamellar simulations, the smoothing length was readily selected, since the particles are evenly spaced on a lattice of sites separated by 0.844 nm. Initially, the concentration on either side of the water layer was fixed and the electrostatics within the pore were neglected in addition to the coordinate dependence of the diffusion coefficient, resulting in simple Fick’s law diffusion which can be compared with analytical results.121 The fixed concentration boundary conditions were incorporated by creating a layer of SPH particles with fixed concentration at the edges of the simulation box which function as a source and a sink for protons. Fixing the concentration at either end is intended to resemble, to a rough approximation, the operation of a PEM fuel cell in which protons are being continuously produced at the anode by dissociating hydrogen gas and being likewise continually consumed at the cathode by water formation to generate regions of proton rich and proton depleted water within the membrane. In future work, the effect of varying the flux conditions at either boundary to account for more realistic reaction rates can be considered as well as proper account for the electric fields near the electrode assembly. One can see from Figure 2a that the agreement between the SPH simulations and the analytical results is excellent for the

III. RESULTS AND DISCUSSION To validate the present method and demonstrate the basic physics, we first considered the simple lamellar model for the hydrated Nafion morphology proposed by Litt.14 Within this model, a linear relationship was proposed between the G



Article

concentration profile along the direction of the applied concentration gradient and is not extremely sensitive to the choice of SPH smoothing length. The results for a cutoff of 2 times the average particle neighbor distance, i.e., 1.688 nm, show deviations of 3% or less from the analytical result, while a cutoff of 3 times the neighbor particle distance produces profiles which differ by less than 1% from the analytical result. These errors increase substantially for smoothing lengths below 1.688 nm, and as one continues to decrease the SPH cutoff to 1 nm and below, the particles cease to exchange composition entirely and no proton flux is seen through the channel. The result in Figure 2a is an important validation of our approach because there is substantial particle inconsistency experienced by each SPH quasi-particle in the simulation region; however, proximity to the boundaries does not affect the accuracy of the results. Interestingly, the proton flux is more sensitive to the choice of SPH cutoff, as seen in Figure 2b. The deviation from the analytical results for a smoothing length of two particle distances is as high as 25% before the establishment of a steady state. Fortunately, since the propagation of the concentration profiles is independent of the evaluation of fluxes, errors in the calculated flux do not contaminate the observed concentration profiles. Given that the proton distribution is correct at steady state, the flux eventually returns to excellent agreement with the analytical solutions and validates our use of the shorter smoothing length to save on computational time. The one exception to this observation remains directly at the interface between the constant concentration boundary and the water region where the calculated gradient does not improve in accuracy as a function of simulation time or as a function of smoothing length. This effect arises from the presence of a corner in the concentration profile, and likewise a discontinuity in the first derivative, at the boundary where the edge particles have a fixed concentration and the internal particles are evolving in time. It is not surprising that the presence of a corner would contribute to significant error in the calculated derivative and could be minimized in future work by exploring different implementations of the boundary conditions.116 In order to avoid edge effects in the present simulations, the flux was only considered one smoothing length in from the boundaries in all the conductivity calculations and the calculated flux was verified to converge with respect to distance from the fixed surface. The conductivities calculated for the lamellar pores are shown in Figure 3 and demonstrate the importance of incorporating both diffusion and morphological effects.88 The value for the conductivity varies by a factor of about 20 between water contents of 5.25 and 15.75, in agreement with expectations based on the changes in diffusion coefficient and water volume fraction. The diffusion constant increases by an order of magnitude in going from the lower to the higher water content,72 and the volume fraction of the water in the channel increases by about a factor of 2 with the increased water content. The blue curve in Figure 3 shows the results for addition of the electrostatics to the same structure in Figure 2, as well as accounting for coordinate dependence of the diffusion coefficient from eq 9. The increase in conductivity is not as pronounced in this case as it was for the lamellar pores without electrostatics. When the charges in the system are included as well as the coordinate dependent diffusivity, the conductivity only increases by a factor of 2 in going from a water content of

Figure 3. Calculated conductivities as a function of water content for the lamellar model of Figure 2. Results are shown using the diffusion coefficients of ref 29b without charged side chains (red) and with charged side chains and eq 9 for the diffusion coefficient (blue).

5.25 to 15.75, in contrast to the previous scenario in which there was an increase by a factor of 20. This more gradual increase can be attributed to two factors. First, the majority of the proton concentration is now pushed against the walls of the pore by the electrostatic interaction between the protons and the sulfonate groups. Second, the variation in the diffusion constant is now based on distance from the nearest sulfonate which, combined with the double layer formation at the pore walls, ensures that the majority of protons in both the low and high water content simulations experience similar transport coefficients within the resolution of the present coarse-grained model. At higher water content, the expanded water regions contain higher proton mobility in the center of the pore; however, these regions are not substantially occupied and therefore do not alter the observed conductivity as profoundly as seen in the red curve of Figure 3. These findings are further supported by Figure 4a,b which shows the resulting concentration profiles at steady state for the lamellar structure with a uniform negative charge density along the walls of a slightly larger 4.22 nm pore. Parts a and c of Figure 4 show the concentration profiles for “equilibrated” water layers in which protons were added to the pore to enforce electroneutrality in the extended water plus polymer

Figure 4. Concentration profiles for lamellar structures with uniform charge distributions along the pore walls (a, b) and with corregated surfaces formed by charged side chains (c, d). The profiles in (a) and (c) correspond to equilibrated structures without an applied gradient, while (b) and (d) include a proton gradient. All concentrations are given in units of protons/nm3, and each frame shares the orientation given by the inset to Figure 2. H



Article

Figure 4c,d also shows the impact of accounting for roughened surfaces in which the side chain groups are treated as protruding from the polymer matrix into the water layer rather than being embedded in the Teflon-rich regions, a behavior which has been noted in previous MD simulations.34−36,71 In comparing with Figure 4a,b, one can see that the protons now favor the bridging positions between the sulfonate groups and the concentration fields are deformed around the side chains to incorporate the corregation of charge density at the surface. Cross sections of the concentration profile are plotted in Figure 5 with blue lines corresponding to two slices through the pore at the extremes of the roughened pore walls: the dotted line connecting bridging positions between sulfonate groups on either side of the pore and the solid line corresponding to a similar slice running from the top of one side chain to the other along the direction of the pore height. The dashed line shows a similar order of magnitude difference between the impermeable walls of the pore and the center of the channel, while the solid line shows a smaller ratio. The factor of 2 difference in magnitude between the two curves reflects the difference in the number of effective sulfonate interactions available to the protons at the two different positions along the pore wall. The concentration profile under the application of a proton gradient in Figure 4d also reveals a qualitatively different mechanism for proton transport from the smooth wall case in Figure 4b. In the latter case, the protons were allowed to slide along the walls of the pore, while in the former case the ions effectively hop between the sulfonate sites along the walls. Up to this point, we have assumed a highly idealized geometry for the pores of the membrane and showed that, within these well-defined systems, our method provides physically reasonable results. However, the true strength of the CG-SPH method lies in its ability to interface with complex water networks derived from CG simulations within the polymer matrix to describe transport through tortuous materials. Before discussing our results for random water networks based on CG DPD simulations, we first provide a validation that our method is not restricted to evenly spaced grids of SPH particles. Figure 6 shows a snapshot of the concentration profile across a water layer of randomly distributed SPH sites at a density comparable to the bead density used in the subsequent DPD simulations. As one can see, the deviations from the analytical result in Figure 6a are minimal and on the order of a few percent relative error. It is particularly noteworthy that the contours for the concentration are nearly straight lines across the width of the pore, which implies that the concentration field is uniform in both directions perpendicular to the applied gradient as anticipated by the unidirectionality of the proton migration. As noted in the case of the uniformly spaced lattice points, Figue 6b shows larger errors in the calculated flux prior to steady state. In agreement with Figure 2b, these errors are appreciably high at the boundaries of the simulation region, i.e., 20% relative to the analytical result. Not only are the relative errors higher for the flux calculation but also the variation in flux across the pore width. Near the center of the simulation region, the flux contours bounce around the average value with a larger amplitude than seen in the concentration profiles, and the amplitude increases near the boundaries. Also, in agreement with the uniform grid spacing results, the accuracy of the flux and the uniformity of its values across the pore width improve as the system approaches steady state.

system, and their density was allowed to evolve in time subject to periodic boundary conditions along the pore length. Clearly, the electrostatic attraction between the protons and the sulfonate groups causes the formation of a nanometer wide diffuse layer along the top and bottom polymer planes in the absence of an applied concentration gradient (see Figure 4a). Suprisingly, even when the proton gradient is applied across the pore in Figure 4b, the majority of additional protons are still pushed against the walls of the pore, suggesting that the interaction with the sulfonate groups overcomes additional repulsion from protons already present near the walls. Figure 5 provides a slice in the concentration profile from Figure 4a along the height of the pore and shows a factor of 10

Figure 5. Concentration in the water layers of Figure 4a (red solid line) and Figure 4c (blue solid and dashed lines) plotted perpendicular to the direction of proton flux. The solid blue line corresponds to a slice taken that hits a sulfur particle on both sides of the lamellar layer, while the dotted line connects two bridging positions on both sides of the same structure.

increase in the proton concentration from the center of the water channel to the edge of the polymer layer. The ratio of concentrations in the double layers and the center of the pore is on the lower end of reports from previous Poisson−Boltzmann (PB) calculations79,80,122 which typically show at least an order of magnitude difference, or more, in ion concentration between the center of the channel and the edge. It is interesting to note that due to the finite resolution of our current model, which is physically motivated by the size of the water clusters used to represent the hydrophilic domains, the behavior of the proton concentration at the edges of the pore is quite different from the exponentially sharp features seen in PB simulations. In Figure 4a,b, the concentration plateaus before the polymer interface, in better agreement with the shifted maximum in ion concentration seen in atomistic simulations.123,124 However, the details of the shifted peak in ion distribution are clearly not seen due to the finite resolution. Recently, the deficiencies of the Poisson−Boltzmann method have been investigated by direct comparison with atomistic simulations111,123,124 and the effects of ion−wall interactions and finite ion size have been extensively discussed in connection with improving PB calculations. Our method might provide an alternative to overcoming these limitations by connecting a more accurate CG model that accounts for such interactions within our combined CG-SPH approach. Furthermore, by implementing coordinate-dependent transport coefficients, including viscosity, that are derived from atomistic simulations, the effects of Stern layers on water and proton transport will be possible in future work. I



Article

merge to form larger pores which are not well described by any of the standard models proposed in the literature, in agreement with previous DPD simulations.40,43,44 The full percolation of the water network is not captured in these plots, since they only show two-dimensional projections of the three-dimensional network, but the threshold behavior is apparent when transport is considered through these materials in the following discussion. In order to explore the size of the resulting water pores formed within these morphologies, radial distribution functions (RDFs) were also considered; see Figure 8. In previous work,

Figure 6. Contour plots of the proton concentration (upper panel) and flux (lower panel) for a lamellar water layer 4.22 nm thick comprised of randomly distributed SPH points at 15 ns. The units for concentration and flux are the same as those used in Figure 2, and the contour values are given on the right-hand side of the plots.

The DPD equilibrated structures for the hydrated Nafion morphology are shown in Figure 7 after 25 ns of equilibration time subject to the parameters of Table 1. One can see that the correct qualitative features are recovered: clustering of the sulfonate groups with the water to form isolated hydrophilic domains at the lowest hydration level considered, λ = 5, followed by expansion of the water clusters to form percolating and complex networks at the higher water contents of λ = 10 and 15. To further highlight the differences between these structures, slices were also taken through each structure to display the water distribution; see the lower panels of Figure 7b. These plots reinforce the observations from the equilibrium snapshots that at low water contents the clusters appear to be around 1−1.5 nm in diameter and well isolated from each other. At higher hydration levels, these clusters clearly begin to

Figure 8. Radial distribution functions for the water beads in the DPD structures averaged over a 25 ns production run after equilibration at varying hydration levels: λ = 5 (red), 10 (blue), and 15 (green). The inset shows an expanded view of the indicated region to highlight the oscillatory structure arising from microphase separation in the polymer.

water RDFs were obtained which showed a clear oscillatory structure corresponding to the underlying water domain formation and spacing.39,44 A similar structure can be seen in Figure 8; however, the locations of extrema for the lowerfrequency phase separation component are masked by a higher frequency oscillation corresponding to the excluded volume repulsion between the DPD beads. In previous work, the size

Figure 7. Snapshots from the equilibrated DPD structures for hydrated Nafion with water contents of λ = 5 (left-most panels), 10 (middle panels), and 15 (right-most panels). The upper panels correspond to the equilibrated DPD structure, color-coded in agreement with Figure 1, while the lower panels correspond to water densities taken from a slice through the middle of each structure using the SPH summation of eq 13. J



Article

and separation of water regions were estimated from the first point at which the RDF dropped below a value of 1 and the location of the first subsequent maximum, respectively. It is clear from Figure 8 that the location of the second maximum shifts with increasing water content, which implies that the centers of the water rich regions are moving farther apart as the clusters expand in size. In addition to the shift of the second maxima, the width of this feature also changes. If one measures the distance between the first time the RDF crosses 1 after 2 nm and the point at which it falls below 1 again, one can estimate the changes in hydrophilic domain sizes to increase from 1.7 nm at λ = 5 to 3.0 nm at λ = 15. An increase in the size of water channels on the order of nanometers is in good agreement with previous studies of the Nafion morphology. Proton transport through the DPD structures shown in Figure 7 was studied by using the same method applied to the lamellar morphology; i.e., a layer of SPH particles was placed on either side with fixed concentration. Figure 9 shows

Figure 10. Conductivity of the DPD equilibrated structures from volume averaging across segments of the structure both without electrostatics (red) and including electrostatics with eq 9 (blue). Experimental results from ref 108 (empty triangles) and ref 109 (solid triangles) are included for comparison.

a steady state with concentration dependent largely on distance from the boundaries at λ = 15, implying minimal tortuosity effects. It is anticipated that the tortuosity is of greater relevance at water contents closer to the percolation threshold. In going from λ = 10 to 15, the water volume fraction increases but there is no significant change to the concentration profiles. Figure 10 shows the resulting conductivities obtained by averaging over four sample volumes in the membrane along the direction of applied proton gradient without accounting for the internal electric fields and demonstrates a factor of about 20 between the lowest and highest water contents considered. This result is surprisingly close to the results from the simple lamellar structure studied initially. In agreement with the lamellar structure, the change in diffusion coefficient as a function of water content accounts for an order of magnitude and the remaining factor arises from the increased porosity at higher water content. The inclusion of electrostatics in the simulation produces a similar effect on the conductivity as seen in the case of the lamellar geometry, namely, the conductivity with respect to hydration level becomes more linear. Both Figures 3 and 10 show that, without the inclusion of electrostatics, the conductivity has a noticeable curvature which is not seen in experiment. Regarding the quantitative agreement with experiment in Figure 10, all three of our results are within an order of magnitude of the measurements, with the two higher water contents being within a factor of 2 of experiment. The trend of higher accuracy at higher values of λ is most likely valid and is related to the level of resolution used in the CG model. By restricting the water to occupy spheres that are tens of molecules in size, the formation of smaller water channels for smaller values of λ is precluded and the effects of porosity may be underestimated. Beyond comparison with conductivity measurements, the distribution of protons within the Nafion pores has also been investigated spectroscopically and can likewise be examined within the DPD structures presented here. At high hydration levels, it was found that the proton concentration within each pore can vary by a factor of 2−3 between the center and the pore wall, while at lower water contents the distributions may be more uniform.90 The results in Figure 11 agree with this trend, and specifically, the boxed regions demonstrate regions in which the proton concentration varies by a factor of 2 between the middle of the pore and the edges. The magnitudes of the concentrations are also close to the range reported by

Figure 9. Concentration profiles for slices through the middle of the equilibrated DPD structures at λ = 5 (left panel) and 15 (right panel) parallel to the applied concentration gradient.

representative results for the proton concentration in the hydrated DPD model using the diffusion coefficients of ref 29b at each water content. The percolation behavior for Nafion is clearly captured as seen by comparing the profile for λ = 5 with the results for a higher hydration level. The darkened middle region reflects minimal interaction with the boundaries and significant proton depletion, while the edges of the box have equilibrated with the boundary conditions on either side. The water network at λ = 5 does not seem to possess substantial continuous channels to allow for significant proton flux through the entire system, but the proximity to the percolation threshold is indicated by the small degree of connectivity that allows for proton penetration into the first 6 nm from the lefthand side and the nonzero conductivity calculated for the system; see Figure 10. A more definitive identification of the threshold will require further investigation of additional water contents both below and above the value of 5 waters/SO3. Once spanning channels form and expand at higher hydration levels, the protons migrate through the entire structure to form K



Article

Figure 11. Proton concentration in DPD structures with local electrostatics and eq 9 for the diffusion coefficient. Boxed areas highlight pores which demonstrate protons favoring wall positions over the middle of the domain in analogy to double-layer formation.

Argonne Leadership Computing Facility on the Intrepid Blue Gene P system.

experiment with the highest proton concentrations corresponding to about 6 M and the lower values closer to 1.5 M. Within the boxed region shown for λ = 15, the value at the pore walls corresponds to about 3 M, while the concentration away from the walls approaches 1.7 M.

■

IV. CONCLUSIONS In this work, we have introduced an approach to address modeling proton conduction through the Nafion PEM at the mesoscale with a unique combination of coarse-grained molecular dynamics as well as continuum mechanics using dissipative particle dynamics and smoothed particle hydrodynamics, respectively. The method is capable of interfacing with detailed information from atomistic MD simulations and offers clear directions for refinement in development of more reliable CG potentials and more accurate descriptions of the coordinate dependence of proton transport coefficients. This method is also highly parallelizable, since the main equations carry the same physical picture used to describe traditional MD, namely, the pairwise interaction of particles. We have applied the present CG-SPH method to two candidate morphologies for Nafion 117 and found that in both cases the inclusion of local electrostatics within the membrane plays a substantial role in affecting the conductivity as a function of hydration behavior of the membrane. Comparison with experiment revealed that, by including double-layer effects and coordinate-dependent proton transport coefficients, the linear conductance profile previously observed could be replicated at a semiquantitative level. As part of ongoing efforts to understand proton conduction in Nafion and related PEMs, the interplay between proton transport and polymer fluctuations will be examined in the future, as well as the effect of using higher resolution in the CG model.

■

REFERENCES

(1) Wang, K.; Dickinson, R. E.; Liang, S. Science 2009, 323, 1468− 1470. (2) “Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, 2007”, Cambridge University Press, 2007. (3) Annual Energy Review 2009; U.S. Energy Information Admistration, 2010. (4) Winter, M.; Brodd, R. J. Chem. Rev. 2004, 104, 4245−4269. (5) Stone, C.; Morrison, A. E. Solid State Ionics 2002, 152−153, 1− 13. (6) Steele, B. C. H.; Heinzel, A. Nature 2001, 414, 345−352. (7) Djilali, N.; Sui, P. C. Int. J. Comput. Fluid Dyn. 2008, 22, 115− 133. (8) Weber, A. Z.; Newman, J. Chem. Rev. 2004, 104, 4679−4728. (9) Goddard, W.; Merinov, B.; Duin, A. V.; Jacob, T.; Blanco, M.; Molinero, V.; Jang, S. S.; Jang, Y. H. Mol. Simul. 2006, 32, 251−268. (10) Paddison, S. J. Ann. Rev. Mater. Res. 2003, 33, 289−319. (11) Kreuer, K. D. J. Membr. Sci. 2001, 185, 29−39. (12) Kreuer, K. D.; Paddison, S. J.; Spohr, E.; Schuster, M. Chem. Rev. 2004, 104, 4637−4678. (13) Mauritz, K. A.; Moore, R. B. Chem. Rev. 2004, 104, 4535−4585. (14) Litt, M. H. Polym. Prepr. 1997, 38, 80. (15) Haubold, H.-G.; Vad, T.; Jungbluth, H.; Hiller, P. Electrochim. Acta 2001, 46, 1559−1563. (16) Rubatat, L.; Gebel, G.; Diat, O. Macromolecules 2004, 37, 7772− 7783. (17) Schmidt-Rohr, K.; Chen, Q. Nat. Mater. 2008, 7, 75−83. (18) Hsu, W. Y.; Gierke, T. D. J. Membr. Sci. 1983, 13, 307−326. (19) Hwang, G. S.; Kaviany, M.; Gostick, J. T.; Kientiz, B.; Weber, A. Z.; Kim, M. H. Polymer 2011, 52, 2584−2593. (20) Elliott, J. A.; Paddison, S. J. Phys. Chem. Chem. Phys. 2007, 9, 2602−2618. (21) Karo, J.; Aabloo, A.; Thomas, J. O.; Brandell, D. J. Phys. Chem. B 2010, 114, 6056−6064. (22) Devanathan, R.; Venkatnatham, A.; Dupuis, M. J. Phys. Chem. B 2007, 111, 8069−8079. (23) Devanathan, R.; Venkatnathan, A.; Dupuis, M. J. Phys. Chem. B 2007, 111, 13006−13013. (24) Venkatnathan, A.; Devanathan, R.; Dupuis, M. J. Phys. Chem. B 2007, 111, 7234−7244. (25) Idupulapati, N.; Devanathan, R.; Dupuis, M. J. Phys. Chem. B 2011, 115, 2959−2969. (26) Devanathan, R.; Venkatnathan, A.; Rousseau, R.; Dupuis, M.; Frigato, T.; Gu, W.; Helms, V. J. Phys. Chem. B 2010, 114, 13681− 13690.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Funding for this project was provided by the U.S. Department of Energy, Office of Basic Energy Sciences (DOE-BES grant number DE-FG02-10ER16171), and under Contract DEAC02-06CH11357. Computer time was provided by the L



Article

(27) Brandell, D.; Karo, J.; Thomas, J. O. J. Power Sources 2010, 195, 5962−5965. (28) Brandell, D.; Karo, J.; Liivat, A.; Thomas, J. O. J. Mol. Model. 2007, 13, 1039−1046. (29) Hofmann, D. W. M.; Kuleshova, L.; D’Aguanno, B. J. Phys. Chem. B 2009, 113, 632−639. (30) Knox, C. K.; Voth, G. A. J. Phys. Chem. B 2010, 114, 3205− 3218. (31) Jang, S. S.; Molinero, V.; Ç aǧin, T.; Goddard, W. A. J. Phys. Chem. B 2004, 108, 3149−3157. (32) Hristov, I. H.; Paddison, S. J.; Paul, R. J. Phys. Chem. B 2008, 112, 2937−2949. (33) Vishnyakov, A.; Neimark, A. V. J. Phys. Chem. B 2001, 105, 9586−9594. (34) Urata, S.; Irisawa, J.; Takada, A.; Shinoda, W.; Tsuzuki, S.; Mikami, M. J. Phys. Chem. B 2005, 109, 4269−4278. (35) Cui, S.; Liu, J.; Selvan, M. E.; Keffer, D. J.; Edwards, B. J.; Steele, W. V. J. Phys. Chem. B 2007, 111, 2208−2218. (36) Cui, S.; Liu, J.; Selvan, M. E.; Paddison, S. J.; Keffer, D. J.; Edwards, B. J. J. Phys. Chem. B 2008, 112, 13273−13284. (37) Liu, J.; Suraweera, N.; Keffer, D. J.; Cui, S.; Paddison, S. J. J. Phys. Chem. B 2010, 114, 11279−11292. (38) Wescott, J. T.; Qi, Y.; Subramanian, L.; Capehart, T. W. J. Chem. Phys. 2006, 124, 134702. (39) Malek, K.; Eikerling, M.; Wang, Q.; Liu, Z.; Otsuka, S.; Akizuki, K.; Abe, M. J. Chem. Phys. 2008, 129, 204702. (40) Yamamoto, S.; Hyodo, S.-a. Polym. J. 2003, 35, 519−527. (41) Galperin, D. Y.; Khokhlov, A. R. Macromol. Theory Simul. 2006, 15, 137−146. (42) Komarov, P. V.; Veselov, I. N.; Chu, P. P.; Khalatur, P. G.; Khokhlov, A. R. Chem. Phys. Lett. 2010, 487, 291−296. (43) Wu, D.; Paddison, S. J.; Elltiot, J. A.; Hamrock, S. J. Langmuir 2010, 26, 14308−14315. (44) Wu, D.; Paddison, S. J.; Eliott, J. A. Energy Environ. Sci. 2008, 1, 284−293. (45) Dorenbos, G.; Morohoshi, K. Energy Environ. Sci. 2010, 3, 1326−1338. (46) Dorenbos, G.; Suga, Y. J. Membr. Sci. 2009, 330, 5−20. (47) Dorenbos, G.; Morohoshi, K. J. Mater. Chem. 2011, 21, 13503− 13515. (48) Qi, Y.; Lai, Y.-H. Polymer 2011, 52, 201−210. (49) Ilhan, M. A.; Spohr, E. J. Phys.: Condens. Matter 2011, 23, 234104. (50) Glezakou, V.-A.; Dupuis, M.; Mundy, C. J. Phys. Chem. Chem. Phys. 2007, 9, 5752−5570. (51) Paddison, S. J.; Zawodzinski, T. A. Solid State Ionics 1998, 113− 115, 333−340. (52) Paddison, S. J.; Elliott, J. A. J. Phys. Chem. A 2005, 109, 7583− 7593. (53) Paddison, S. J.; Elliott, J. A. Phys. Chem. Chem. Phys. 2006, 8, 2193−2203. (54) Zhao, Q.; Majsztrik, P.; Benziger, J. J. Phys. Chem. B 2011, 115, 2717−2727. (55) Curtin, D. E.; Lousenberg, R. D.; Henry, T. J.; Tangeman, P. C.; Tisack, M. E. J. Power Sources 2004, 131, 41−48. (56) Groot, R. D.; Madden, T. J. J. Chem. Phys. 1998, 108, 8713− 8724. (57) Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423− 4435. (58) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Europhys. Lett. 1992, 19, 155−160. (59) Schlijper, A. G.; Hoogerbrugge, P. J.; Manke, C. W. J. Rheol. 1995, 39, 567−579. (60) Español, P.; Warren, P. Europhys. Lett. 1995, 30, 191−196. (61) Dorenbos, G.; Morohosi, K. J. Chem. Phys. 2011, 134, 044133. (62) Dorenbos, G.; Pomogaev, V. A.; Takigawa, M.; Morohoshi, K. Electrochem. Commun. 2010, 12, 125−128. (63) Knight, C.; Voth, G. A. Acc. Chem. Res. 2011, in press.

(64) Frenkel, D.; Smit, B. Understanding Molecular Simulation, 2nd ed.; Academic Press: San Diego, CA, 2002. (65) Paul, R.; Paddison, S. J. J. Chem. Phys. 2005, 123, 224704. (66) Paddison, S. J.; Paul, R.; Zawodzinski, T. A. J. Chem. Phys. 2001, 115, 7753−7761. (67) Paddison, S. J.; Paul, R. Phys. Chem. Chem. Phys. 2002, 4, 1158− 1163. (68) Zawodzinski, T. A.; Neeman, M.; Sillerud, L. O.; Gottesfeld, S. J. Phys. Chem. 1991, 95, 6040−6044. (69) Choe, Y.-K.; Tsuchida, E.; Ikeshoji, T.; Yamakawa, S.; Hyodo, S.-a. Phys. Chem. Chem. Phys. 2009, 11, 3892−3899. (70) Roudgar, A.; Narasimachary, S. P.; Eikerling, M. Chem. Phys. Lett. 2008, 457, 337−341. (71) Petersen, M. K.; Voth, G. A. J. Phys. Chem. B 2006, 110, 18594− 18600. (72) Feng, S.; Voth, G. A. J. Phys. Chem. B 2011, 115, 5903−5912. (73) Seeliger, D.; Hartnig, C.; Spohr, E. Electrochim. Acta 2005, 50, 4234−4240. (74) Spohr, E. Mol. Simul. 2004, 30, 107−115. (75) Spohr, E.; Commer, P.; Kornyshev, A. A. J. Phys. Chem. B 2002, 106, 10560−10569. (76) Selvan, M. E.; Keffer, D. J.; Cui, S. J. Phys. Chem. C 2011, 115, 18835−18846. (77) Selvan, M. E.; Keffer, D. J.; Cui, S.; Paddison, S. J. J. Phys. Chem. C 2010, 114, 11965−11976. (78) Choi, P.; Jalani, N. H.; Datta, R. J. Electrochem. Soc. 2005, 152, E123−E130. (79) Berg, P.; Ladipo, K. Proc. R. Soc. London, Ser. A 2009, 465, 2663−2679. (80) Ladipo, K. O.; Berg, P.; Kimmerle, S.-J.; Novruzi, A. J. Chem. Phys. 2011, 134, 074103. (81) Guzman-Garcia, A. G.; Pintauro, P. N.; Verbrugge, M. W.; Hill, R. F. AIChE J. 1990, 36, 1061−1074. (82) Pintauro, P. N.; Tandon, R.; Chao, L.; Xu, W.; Evilia, R. J. Phys. Chem. 1995, 99, 12915−12924. (83) Yang, Y.; Pintauro, P. N. Ind. Eng. Chem. Res. 2004, 43, 2957− 2965. (84) Ayton, G. S.; McWhirter, J. L.; McMurty, P.; Voth, G. A. Biophys. J. 2005, 88, 3855−3869. (85) Ayton, G. S.; McWhirter, J. L.; Voth, G. A. J. Chem. Phys. 2006, 124, 064906. (86) Plimpton, S. J. Comput. Phys. 1995, 117, 1−19. (87) Groot, R. D.; Rabone, K. L. Biophys. J. 2001, 81, 725−736. (88) Selvan, M. E.; Calvo-Muñoz, E.; Keffer, D. J. J. Phys. Chem. B 2011, 115, 3052−3061. (89) Corry, B.; Kuyucak, S.; Chung, S.-H. Biophys. J. 2000, 78, 2364− 2381. (90) Spry, D. B.; Fayer, M. D. J. Phys. Chem. B 2009, 113, 10210− 10221. (91) Hoover, W. G.; Hoover, C. G. Mol. Phys. 2003, 101, 1559− 1573. (92) Liu, M. B.; Liu, G. R. Arch. Comput. Methods Eng. 2010, 17, 25− 76. (93) Monaghan, J. J. Rep. Prog. Phys. 2005, 68, 1703−1759. (94) Monaghan, J. J. Annu. Rev. Astron. Astrophys. 1992, 30, 543−574. (95) Lucy, L. B. Astron. J. 1977, 82, 1013−1024. (96) Gingold, R. A.; Monaghan, J. J. Mon. Not. R. Astron. Soc. 1977, 181, 375−389. (97) Monaghan, J. J. J. Comput. Phys. 1983, 52, 374−389. (98) Y, Z.; Fox, P. J. Transp. Porous Media 2001, 43, 441−471. (99) Y, Z.; Fox, P. J. J. Comput. Phys. 2002, 182, 622−645. (100) Tartakovsky, A. M.; Meakin, P. Adv. Water Resour. 2006, 29, 1464−1478. (101) Tartakovsky, A. M.; Ward, A. L.; Meakin, P. Phys. Fluids 2007, 19, 103301. (102) Cleary, P.; Monaghan, J. J. J. Comput. Phys. 1999, 148, 227− 264. (103) Español, P.; Revenga, M. Phys. Rev. E 2003, 67, 026705. (104) Yan, L.; Ji, X.; Lu, W. J. Phys. Chem. B 2008, 112, 5602−5610. M



Article

(105) Groot, S. R. d.; Mazur, P. Non-equilibrium thermodynamics; Interscience Publishers: New York, 1962. (106) Steele, W. A. Time-Correlation Functions. In Transport Phenomena in Fluids; Hanley, H. J. M., Ed.; Marcel Dekker: New York, 1969. (107) Ochi, S.; Kamishima, O.; Mizusaki, J.; Kawamura, J. Solid State Ionics 2009, 180, 580−584. (108) Zawodzinski, T. A.; Neeman, M.; Sillerud, L. O.; Gottesfeld, S. J. Phys. Chem. 1991, 95, 6040−6044. (109) Kreuer, K. D.; Schuster, M.; Obliers, B.; Diat, O.; Traub, U.; Fuchs, A.; Klock, U.; Paddison, S. J.; Maier, J. J. Power Sources 2008, 178, 499−509. (110) Perrin, J.-C.; Lyonnard, S.; Volino, F. J. Phys. Chem. C 2007, 111, 3393−3404. (111) Ho, C.; Qiao, R.; Heng, J. B.; Chatterjee, A.; Timp, R. J.; Aluru, N. R.; Timp, G. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 10445−10450. (112) Newman, J. S. Electrochemical Systems; Prentice-Hall: Englewood Cliffs, NJ, 1973. (113) Groot, R. D. J. Chem. Phys. 2003, 118, 11265. (114) González-Melchor, M.; Mayoral, E.; Velázquez, M. E.; Alejandre, J. J. Chem. Phys. 2006, 125, 224107. (115) Wolf, D.; Keblinski, P.; Phillpot, S. R.; Eggebrecht, J. J. Chem. Phys. 1999, 110, 8254−8282. (116) Chen, J. K.; Beraun, J. E.; Carney, T. C. Int. J. Numer. Methods Eng. 1999, 46, 231−252. (117) Brookshaw, L. Proc. Astron. Soc. Aust. 1985, 6, 207−210. (118) Whitaker, S. AIChE J. 1967, 13, 420−427. (119) Vidts, P. D.; White, R. E. J. Electrochem. Soc. 1997, 144, 1343− 1353. (120) Gupta, A.; Seo, J. H.; Zhang, X.; Du, W.; Sastry, A. M. J. Electrochem. Soc. 2011, 158, A487−A497. (121) Crank, J. The mathematics of diffusion; Clarendon Press: Oxford, U.K., 1979. (122) Eikerling, M.; Kornyshev, A. A. J. Electroanal. Chem. 2001, 502, 1−14. (123) Kim, D.; Darve, E. Phys. Rev. E 2006, 73, 51203. (124) Qiao, R.; Aluru, N. R. J. Chem. Phys. 2003, 118, 4692−4701.

N


Mesoscale Simulation of Proton Transport in Proton Exchange

Recommend Documents