Modeling and Simulation for Fuel Cell Polymer Electrolyte Membrane

We have established methods to evaluate key properties that are needed to commercialize polyelectrolyte membranes for fuel cell electric vehicles such as water diffusion, gas permeability, and mechanical strength. These methods are based on coarse-graining models. For calculating water diffusion and gas permeability through the membranes, the dissipative particle dynamics–Monte Carlo approach was applied, while mechanical strength of the hydrated membrane was simulated by coarse-grained molecular dynamics. As a result of our systematic search and analysis, we can now grasp the direction necessary to improve water diffusion, gas permeability, and mechanical strength. For water diffusion, a map that reveals the relationship between many kinds of molecular structures and diffusion constants was obtained, in which the direction to enhance the diffusivity by improving membrane structure can be clearly seen. In order to achieve high mechanical strength, the molecular structure should be such that the hydrated membrane contains narrow water channels, but these might decrease the proton conductivity. Therefore, an optimal design of the polymer structure is needed, and the developed models reviewed here make it possible to optimize these molecular structures.


Introduction
In the field of industrial polymer materials, high functional polymer membranes that require high physical strength and good gas separation properties have been investigated.In order to develop and commercialize new polymeric materials, a number of demanding properties must be satisfied simultaneously.In this development process, new polymer synthesis is needed repeatedly, which makes it difficult to commercialize new functional polymers within a short period of time.To avoid such a time delay, molecular simulation is frequently used, thus enabling the synthesis and evaluation of membranes virtually with a computer.However, since there are not so many physical properties that can be evaluated by simulation, its utilization is often limited to the analysis of physical phenomena.Thus, we have considered that if certain properties could be evaluated by molecular simulation, it would improve the polymer development process.
The polyelectrolyte membrane, which is used for fuel cell vehicles, satisfies various properties (Figure 1), and as such it was used as a study subject for improving the development process.Among the required properties, it is already known that the deterioration factor of chemical stability exists in the catalyst layer [1] and that dimensional stability against swelling can be evaluated by simulation [2].Therefore, we developed new methods to evaluate physical properties that are associated with proton conductivity [3], gas permeability [4], and mechanical strength [5].These methods are reviewed in this work.Polyelectrolyte membranes for fuel cell applications play a decisive role in separating hydrogen from oxygen, which are fuel gases, while it also functions as a proton conductor.Nafion, which is a perfluorosulfonic acid membrane and consists of hydrophobic Teflon backbone and side chains that contain a hydrophilic sulfonic group, is well known as a representative membrane [6].The sulfonic groups absorb water under a moist environment, resulting in water that is phase separated from Teflon, and eventually results in a water cluster network with clusters up to several nanometers in diameter [7,8].The phase-separated morphology contains both the gas barrier for oxygen diffusion by Teflon phase and the proton conductivity pathway through the water phase.
Already the earliest molecular dynamics simulations of hydrated Nafion in which all atoms are treated directly indeed revealed phase-separated structures [9,10].However, in praxis, it is difficult to obtain an appropriate water cluster network, due to following two reasons.At first, the relaxation of phase separation to its equilibrium structure is slow and requires very long computation time.This is because, in the case of polymer, the relaxation time increases with chain length.Second, it is necessary to consider a system the size of dozens of nanometers to evaluate the network of water clusters, meaning that more than 100,000 atoms are needed, which requires a drastic increase in computation time.An advantage of molecular dynamics simulation is that proton conductivity can be calculated directly from the mean square displacement of the protons, but due to the limited sampling time of only a few ns, it can evaluate only their motion within certain parts of the water cluster.As a result, the diffusion constants derived from simulations would be different by an order of magnitude [11] when comparing experimental long-range data with the short-range diffusion obtained from simulation studies.Therefore, several methods have been developed by which some aggregates of atoms are coarse grained, and larger time and space scales can be treated.Water cluster size derived from experiments can be reproduced well by mean field methods, though well-resolved network structures are not obtained [12].Water cluster networks that correspond well with those obtained from experiments can be acquired by using dissipative particle dynamics (DPD) [13].Also, by combining DPD simulations with a statistical analysis of small angle X-ray and neutron-scattering experiments [14], good agreement with the experiment was obtained.However, since DPD, in principle, allows for the crossing of polymer chains, serious care should be taken when directly applied to study dynamic phenomena.Molecular simulation along with ab initio calculation was applied to proton transport phenomena [15].This multi-scale treatment is a straightforward approach, but evaluating long-range motion through the water cluster network for various membranes at several hydration levels needs a lot of computation time.
In this paper, we review the coarse-graining DPD and Monte Carlo (MC) methods [16] that were used for the evaluation of proton conductivity and gas permeability by considering the motion of tracer particles by MC through the phase-separated structures obtained from DPD.Another coarse-graining strategy, coarse-grained molecular dynamics (CGMD), was applied to understand the microstructure within the membrane [17] and catalyst layer [18].These studies compared water cluster size with small angle X-ray scattering data.In addition to equilibrium structures also polymer chain dynamics can be analyzed by CGMD [19].For this reason, we also used CGMD to evaluate mechanical strength which is the resultant property of the motion of the phase -separated structure itself.

Calculation Method
In coarse-graining Nafion, shown in Figure 2, several aggregates of atoms were represented by beads.Here the fragments -CF 2 -CF 2 -CF 2 -CF 2 -, -O-CF 2 -CF[CF 3 ]-O-and -CF 2 -CF 2 -SO 3 H were defined as beads A, B and C, respectively, and H 2 O were represented by W beads, with each W bead containing four water molecules [13].The volume of the beads was around 0.12 nm 3 .

DPD and MC
In DPD [20], forces that act on each particle consist of the conservative force C , the dissipative force D and the random force R , each of which is described in Equations (1)(2)(3).
Here, the interaction distance is expressed by r c , the relative position and velocity between particle i and j is given by r ij = r i -r j and v ij = v i -v j , ω D and ω R are weight functions, and the randomly fluctuating variable θ ij introduces randomness into the system with Gaussian statistics.For particles bonded to each other within the same polymer chain, Equation ( 4) is added as the bonding force S .

S (4)
In Equation ( 4), the spring constant C was set to 4 so that the mean distance between the bound particles was about the same as that between the unbound particles [20].The bead density was set at the default value 3. Since for Nafion the bead volumes are about 0.12 nm 3 , this results in the interaction distance in of r c = 0.71 nm.For Nafion, the values in Table 1 were used for the interaction parameter a ij in Equation (1).They were derived by a ij = a ii + 3.27 χ ij .By setting the diagonal terms at a value of a ii = 104, the water compressibility is reproduced for a system solely composed of W beads. Table 1 also shows the χ parameters in parenthesis which were calculated from atomistic simulations of mixing energies [13].The simulation cells were cubic with box edge lengths of more than 35 times the interaction distance r c .In all DPD calculations, equilibrium phase-separated structures that served as input for the MC calculations were obtained at more than 20,000 time steps, using a time step of 0.05.
In the DPD-MC strategy, the pore morphologies at various hydration levels were calculated by DPD and diffusion constants within the membrane was calculated separately by mapping static DPD morphologies onto a cubic grid on which lattice MC calculations of tracer particles through the pore networks occurs.In the MC trajectory calculations, the particle distribution obtained by DPD was converted onto a three-dimensional grid on which the motion of MC particles through the phase separated structure was simulated.The grid resolution was one-tenth of the interaction distance r c in DPD so that the phase-separated structure of DPD could be rebuilt.The type of each node was defined as that type of bead which is most close to that node.As a result, every node belongs to the polymer phase or the water phase.When investigating proton conductivity, the pore network is defined as a summation of nodes with a nearest bead being W or C. For examining errors by conversion, it was confirmed that the difference between the pore bead fraction within the DPD and pore node fraction of the MC grid was within 0.5%.The diffusion of water and permeation of hydrogen and oxygen gas through the membranes was mimicked by following the motion of tracer particles through the mapped morphologies.This was done by initially placing 1000 or 2000 virtual water particles at nodes that belong to the pore networks, or in case of gas permeation by placing 4000 virtual gas particles at randomly selected nodes in both pore and polymer phase.During each subsequent step, each particle performs one jump trial towards a randomly selected nearest node.In case of water diffusion, the jump trial is only accepted when this node also belongs to the pore space.For gas diffusion, the probability to accept the jump trial depends on both the property of the node on which the particle is located and the property of the node to which it tries to jump.The diffusion coefficient of water within membrane is obtained through Equation (5).
Here, N is the number of tracer particles and R i (t) is the position of tracer particle i on the lattice at time t.According to Equation ( 5), the slopes obtained by plotting the mean square displacement (MSD) against MC step (MSC), are proportional to the diffusion constant.For pure water, each jump trial is accepted resulting in slopes equal to 1.For diffusion through the hydrated membranes the slopes will be <1, since not each jump trial is accepted.Therefore, for the model membranes, the calculated diffusion constants are expressed relative to the pure water case and are given by D' rel.= d(MSD/d(MCS).A comparison with experimental values for diffusion in Nafion was achieved by multiplication of D' rel. with the pure water diffusion constant of 2.3 × 10 −5 cm 5 /s.

CGMD
Each particle motion is described by the Langevin equation in CGMD, and the forces that act on the particles are obtained from the potentials given in Equations (6)(7)(8).
U bond = k r ( 7) The Lennard-Jones type potential in Equation ( 6) was used to describe the non-bond interaction between particles, while the harmonic and linear oscillator type potentials in Equations ( 7) and ( 8) were used for describing the bond length and bond angle interactions, respectively.Since setting of potential parameters in Equations (6-8) is important for an analysis of actual polymer, these were determined from atomistic modeling in which all atoms were considered.
In order to determine the non-bond potential, reversible work [21] is used to map the explicit atoms onto coarse grained particles.This is obtained, using Equation ( 9), from the potential energy U between coarse-graining units sampled over lots of configurations [21].
The parenthesis indicates that potential energy between two coarse graining units separated by a distance r ij is averaged over all orientations.
For coarse-grained particles that are bonded to each other, the bond and angle potentials are determined by the bond and angle average of the Boltzmann's factor [21].Equations ( 10) and ( 11), based on the length and angle distribution P obtained by sampling a large number of conformations, are used [22] to obtain these potentials.
In obtaining the potential energies U LJ and distributions P, energy calculations more than 10,000 times for potential energy and molecular dynamics simulations longer than 10 ns for calculating the distributions were performed for sampling, where COMPASS [23] force field was applied to all atoms.From obtained U and P, the potentials derived from Equations (9)(10)(11) were fitted to the function forms of Equations (6)(7)(8) to determine the parameters.Thus, obtained fitted parameters are listed in Table 2.These values are normalized by the energy of 0.533 kcal/mol and the length of 0.544 nm and the angle is measured as external angle (the angle of linearly connected atoms is 0 degree).
Table 2. Potential parameters of coarse-grained particles.

Proton Conductivity
Since protons are conducted through the hydrophilic water-containing phase, good correlation is recognized between proton conductivity and water diffusion.Therefore, water cluster shapes in membranes were calculated by the DPD module contained in Materials Studio commercialized by Accelrys Inc., and water diffusion though these clusters was evaluated by MC tracer diffusion simulations.A comparison between calculated diffusion constants with those obtained from experiment in Nafion117 are shown in Figure 3.In the calculations, water was allowed to diffuse within the hydrophilic pore network.The calculations agree well with the experimental values over a wide range of moisture absorption This indicates that the DPD-MC approach can be utilized for quantitative estimation of diffusion constants.As for the required CPU time, one type of membrane can be evaluated within about one day and therefore a systematic search can be performed for membranes with various molecular structures generated in a computer virtually.This allowed us to deduce persistent trends for diffusion when a certain parameter is systematically varied, such as membrane EW [16], side chain length, block length, parameters between different molecular fragments [3], etc.For this reason, polyelectrolyte membranes shown in Figure 4 were investigated.The effect of the block chain length of hydrophobic A particles and hydrophilic C particles for the block type (a), and the effect of main chain length and side chain length for the grafted type (b), were estimated by changing values of x and y in Figure 4.
Figure 5 shows the water cluster networks that were obtained for three types of block polymer architectures that contain 28.6 vol % water.Since the C bead fraction within each type is the same, the equivalent weight (EW), defined as the average molecular weight per acid molecule, was also the same for all membranes.For {A 16 C 4 } 1 , which had the longest hydrophilic and hydrophobic block chain lengths, well-connected water clusters with a large diameter were seen.For{A 8 C 2 } 2 , for which the block lengths were two times smaller, water clusters were connected, but their diameter was small.On the other hand, for {A 4 C 1 } 4 , extremely small clusters that seem to be badly connected were seen and some of them were even isolated.The calculated water diffusion coefficient increases with an increase of hydrophilic block length.In order to explore the effect of water cluster shapes on water diffusion, the cluster radius and the separation between clusters were estimated from the pair correlation function of the W beads.The inter-cluster distance is defined as the position where the second peak in the pair correlation function occurs, while the pore radius is defined as that distance for which the pair correlation function of the W beads starts to drop below the value 1.The relationship between cluster shape and water diffusion is presented in Figure 6 for block types with various molecular structures.For each block type, the level of moisture absorption expressed as the number of water molecules per acidic site was at the same level of λ = 4.In Figure 6, dotted lines are drawn as guide lines that connect the results obtained for different block chain lengths, but with the same EW.  Figure 6 indicates that water diffusion is enhanced with the increase of water cluster radius and also with cluster separation.Moreover, since in the graph water diffusion is enhanced by moving to the different dotted line on the right hand side, rather than by moving up to the upper right along the same dotted line, it is understood that effective molecular design to enhance water diffusion is achieved by reducing EW at the first step and lengthening the block length at the next step.
Figure 7 shows water cluster networks for two grafted type membranes that differ in side chain length without changing EW.The moisture absorption was 18.2 vol % for both.Although the water cluster shapes are almost the same, the cluster radius and separation distance increase were found to increase with side chain length and the water diffusion is also enhanced.
In order to figure out why the water diffusion greatly differs in spite of similar morphologies, a percolation analysis was performed for the network structures (Figure 8).The fraction of the W beads that is contained within the largest cluster is plotted against the maximum distance for which two W particles are considered to be connected and join the same cluster.It enables the grasping of the critical distance for which almost all W particles are contained in one network from the inflection points of the graph, and the capturing of the quantity of the isolated water clusters which are not included in the network from the convergent values for which the graph becomes flat at large-defined length of W bead connection.Figure 8 indicates that this critical distance is large and the W bead fraction at which the graphs converge is small in the membrane with short side chains.In other words, water clusters are difficult to be connected and 30% of the water exist as isolated clusters, so this is the cause of the drastically low water diffusion.As a result of a systematical search, in which diffusion constants were derived for several tens of polymer architectures, interesting trends were obtained.Figure 9 shows the relationship between grafted type molecular structures and water diffusion.The line of constant EW is indicated with a solid line, that of constant main chain length is indicated with a broken line, and that of constant side chain length is indicated with a dotted line.
As shown by the arrows in the figure, there are three molecular design indices to enhance water diffusion.Among these indices, lowering EW by shortening the main chain (reducing x) has the most significant effect on enhancing water diffusion, lowering EW by shortening the side chain (reducing y) has the middle effect, and the effect of lengthening the side chain while keeping the EW fixed is lowest.Moreover, if the molecular structure of the membrane is drastically changed, water diffusion can be improved greatly, even by means of lengthening the side chain with keeping EW at constant.

Gas Permeability
First, the solubility S and diffusion coefficient D of gas species within the hydrated membrane were estimated, before calculating the gas permeability K as K =S × D. Contrary to water and protons, gas dissolves not only within the water phase, but also within the polymer phase, and gas diffusion occurs within both phases.Furthermore, in the diffusion process, gas species can enter one phase while leaving the other.In the hydrated Nafion, the solubility within these two components were obtained based on the assumption that the gas solubility is determined as the volume averaged value based on the solubility and volume fraction of each component.The volume averaged value of solubility is expressed in Equation ( 12) where f Nafion is the volume fraction of Nafion.
Equations ( 13a) and (13b) were used for solubility of H 2 and O 2 in dry Nafion, respectively.These equations were derived from the experimental data [20].The H 2 and O 2 solubility in water were deduced from the temperature dependencies of the molar fraction solubility [24,25] and water mass density [25].Thus obtained temperature dependent solubilities are given by Equations (14a) and (14b).
Equations (15)(16)(17)(18) with T in Kelvin were used for diffusion in the Nafion phase [26] and water phase [27].Ratio of side and main chain length As an example, for 40 degrees Celsius, the solubilities of oxygen are 1.02 mmol/L in water and 6.11 mmol/L in Nafion, while the diffusions of oxygen are 3.04 × 10 −5 cm 2 /s in water and 0.014 × 10 −5 cm 2 /s in Nafion.
The calculation of diffusion coefficients was performed by MC.Since the diffusion of both gas species is always highest in Nafion, a move trial of a gas particle to jump from water node towards a nearby selected water node is always accepted, while within the Nafion phase this probability is reduced by D Nafion /D water .The gas solubility within the Nafion phase is always larger than within the water phase.Therefore, for a gas particle that tries to leave a water node and enter a nearest selected Nafion node, the trial is always accepted, while a jump trial in the reverse direction is reduced and becomes S water /S Nafion .These jumping rules ensured that the diffusion of the gas particles within both phases was well described, and that during the simulations, the gas particle concentration within each phase matched exactly the experimental pure phase component solubility.The gas diffusion constants within the membrane were obtained by multiplication of the slopes derived from the MSD curves with the diffusion constant of oxygen and hydrogen in the pure water phase given by Equations ( 17) and (18), respectively.The permeability K was calculated by multiplying of the thus-obtained diffusion constant with the solubility values given by Equation (12).
Figure 10 shows the comparison of hydrogen and oxygen permeability through Nafion 125 with experimental values.Three different hydration levels at three temperature conditions were considered, volume fractions of Nafion are 0.5, 0.6 and 0.7 and temperature conditions are 40, 55 and 70 degrees Celsius.Hydrogen Oxygen The results indicate that gas permeability can be predicted quantitatively by the developed method for various gaseous species, component fractions and temperatures.There may be several reasons why the predicted values are less than in the experiment.For instance, due to the plasticizing effects of the diffusion coefficient, the Nafion component might increase with the increase of moisture sorption.However, when considering the gas transport through hydrated membranes, it is difficult to distinguish between diffusion constants in Nafion phase and water phase from an experiment.Realizing this would likely make the model prediction better.A similar argumentation can be made for the actual gas solubility within the Nafion phase with the water uptake.The interesting finding presented in Figure 10 is that, without invoking any fitting parameters, gas permeability values can be approached by the developed coarse-graining DPD-MC strategy.
In order to know the direction to raise gas barrier properties of oxygen and hydrogen, which are fuel gases, gas transport phenomena were analyzed.The Maxwell model for gas transport in heterogeneous media was used for the analysis [28] of the MC permeation results.The Maxwell model enables the association of the permeability with inner structures of the two-component system by the parameter G in Equation ( 19), as shown in Figure 11.In Equation ( 19), K water and K Nafion are the permeability of the gas species in the pure water and Nafion components, respectively.
Based on the calculated hydrogen and oxygen permeability, the parameter G was calculated by Equation ( 19) for various Nafion fractions (Figure 12).If mobility between different components would be the same, one G value can be associated with a certain structure.However, for the same phase-separated structures at the same Nafion volume fraction, values of G greatly differ depending on gaseous species, so mobility between two components are different for hydrogen and oxygen.Moreover, as the volume fraction of Nafion decreases, the water cluster networks are better connected leading to an increased G value.This is indeed clearly seen for oxygen, but not for hydrogen.This is because the migration or transfer between two components occurs more easily for hydrogen, while the migration of oxygen mainly occurs within the single water component.That is, the migration-location differs depending on gaseous species.In order to confirm the migration-location directly, we investigated the component where gas particles mainly stay during 10,000 steps.The results are presented in Figure 13.• Hydrogen ▲ Oxygen

Nafion volume fraction f nafion Structural parameter G
The horizontal and vertical axes indicate time and components in which the particles stay, respectively, and the line across the components indicates a jump between two components.Figure 13 also shows that, although the jump events occur for both hydrogen and oxygen, the frequency of which they occur is higher for hydrogen.Also, the oxygen gas particles stay longer in the Nafion component than the hydrogen particles.
Since the extent to which gas migration occurs within one phase as compared to the other is different for hydrogen and oxygen, this will affect the direction of polymer design in order to control permeation of fuel gasses.Therefore, it has been clarified that gas barriers can be enhanced by control of percolation and interface structures of water clusters for hydrogen and also be enhanced by lowering the diffusivity in the polymer phase for oxygen.

Mechanical Strength
Mechanical strength was calculated for membranes at various moisture levels.Each Nafion 117 chain consisted of 25 repeat units and the water volume fraction within the Nafion membrane was varied as listed in Table 3.After generating coarse-graining particles of Nafion and water randomly at the initial state, energy was converged within the NPT ensemble at 25 degrees Celsius for more than 10 6 time steps to obtain equilibrium structures of hydrated membranes.The time step was set at 0.006 so that relative errors of the temperature were less than 0.3%.When evaluating elastic modulus from a stress-strain curve by extending the hydrated membranes, the cells were deformed by 5% with a speed of 0.0176 m/s along one direction while keeping the cell volume constant under the deformation.To save the calculation cost the deformation rate used here is almost twenty times higher as the deformation rate in the experiment [29].The general coarse graining molecular dynamics simulator OCTA/COGNAC [30] was used for the calculation.At first, before evaluating elastic modulus, we confirmed whether the experiment could be reproduced with regards to the bulk density and water cluster network structures for hydrated Nafion membranes.The experimental Nafion mass density in dry state was reported as 2.05 g/cm 3 [31] while our calculated value was 1.94 g/cm 3 .Under the water absorption conditions at λ = 3, λ = 4, λ = 5 and λ = 14, the mass densities shown in Figure 14 were obtained, this indicates that the calculated bulk density agreed well with the experimental values.On the other hand, absorbed water forms cluster network structures inside the membranes and the separation distance between the clusters have already been estimated from small-angle X-ray scattering experiments [7,8].We estimated the cluster separations from the radial distribution function of the coarse-graining water particles, and the results are compared with the experimental results in Figure 15.Overall, the inner structures were reproduced well quantitatively for a wide range of hydration levels, which reveals that the calculation method used in this study is adequate.At high water volume fraction the obtained cluster spacing is slightly lower than the experimental values.The reason for this might be that the atomistic force field used for the water molecules to obtain the coarse-grained potentials is the 3-site model, which might somewhat underestimate the interaction between water, resulting in more dense water and therefore smaller cluster spacings at high hydration levels.Second, the membranes with moisture absorption of λ = 0, λ = 4, λ = 8 and λ = 20 were stretched and elastic modulus was evaluated.Since deformation volume at the longest extension is as small as 5%, while the yield point in the experiments appears above 5% deformation, we calculated the elastic modulus from the gradient of the stress-strain curve under the assumption that stress distortion occurred within the elastic domain.The moisture absorption dependence of the elastic modulus is shown in Figure 16.While from the simulations the tendency appears in which elastic modulus has a maximum when plotted against moisture absorption as observed in experiment, the absolute values of the elastic modulus cannot be reproduced quantitatively.As mentioned above, this is because the extension speed was twenty times faster than in the experiments.The cause of the apparent maximum in the elastic modulus in Figure 16 was investigated by analyzing the internal structure of the hydrated membranes.Since higher-order structure such as appearance of crystals was not observed during our calculations for any hydrated membrane, we tried to analyze from the viewpoint of the molecular theory of rubber elasticity.
In membranes, polymer chain configurations are complicated and intertwined with each other, and this intertwinement acts as a crosslink under deformation.Since we presumed that the extent in which intertwinements occur depends on the amount of absorbed water, we analyzed them.In the analysis, chains were pulled by shortening the equilibrium bond length while keeping both ends of all polymer chains fixed in space, and we evaluated the number of entanglement points per chain by counting how many of them are intertwined with each other [32].The results are shown in Figure 17.

Water content λ Number of entanglements per chain
There are 1.5-3.6 entanglements per polymer chain and overall they decrease monotonously with the amount of absorbed moisture.This change in the number of entanglements with water moister uptake observed in Figure 17 seems not to be related with the elastic modulus data in which a maximum was recognized.
The effect of water content on crosslink formation was examined.Although Nafion does not have an intrinsic crosslinked molecular structure, it is possible that the sulfonic groups in the structure, along with water molecules, generate tightly bound hydrogen bonds, and that these hydrogen bonds are involved in quasi-crosslinked structures between two sulfonic groups that are bridged by water molecules.Presuming that these quasi-crosslinks are related to the enhancement of mechanical strength, we estimated the number of quasi-crosslinked C-W-C structures that are generated between the coarse-grained C particles and the coarse-grained water W particles.A quasi-crosslink is defined as a C-W-C bridge for which both C-W and W-C distances are shorter than the equilibration distance 1.1σ obtained from the interaction potential.Figure 18 shows moisture absorption dependence of the number of quasi crosslinks per polymer chain.The number of quasi crosslinks has a maximal value for moisture absorption near λ = 5 which is near the λ value for which a maximum in the elastic modulus occurs.Moreover, near this maximum, more than a half of the sulfonic groups in one polymer chain are involved in the formation of quasi-crosslinked structures.Since the number is about 1.5 times larger than in the case of large amount of water, it is likely that quasi-crosslinked structures consisting of sulfonic groups and water molecules are the source of the maximum value observed for the elastic modulus, which affects mechanical strength of membranes.If the quasi-crosslinks affect the elasticity significantly, the quasi-crosslink points that are once formed are supposed to persist during elongation.Thus, we investigated how much of the initial crosslinked structures are maintained when they are stretched.Since the result showed that about 70% of the structures were maintained after extension by 5%, it is understood that quasi-crosslinks cause the generation of elasticity.Further, the reason why the number of quasi crosslinks decreases for large moisture absorption is that the quantity of water contained between sulfonic groups increases, which makes it difficult to form quasi-crosslinks. From the above results, an increase of quasi-crosslinked structures is one of the methods to enhance the mechanical strength of membranes and to improve its durability.In order to increase quasi-crosslinks also under high moisture absorption conditions, it is important to generate narrow water channels as cluster networks, where a large quantity of water cannot easily be located between the acidic sites.It is possible to perform molecular design for membranes that realize the above from the relationship between cluster sizes and polymer structures [3].However, water cluster networks are the migration path of protons.Therefore, proton conductivity simultaneously could fall for narrow channels.Hence, it is important to clarify the cluster size that realizes both higher mechanical strength and better proton conductivity.Utilization of polymer simulation enables the deriving of molecular structures which can have optimal values of cluster sizes within a short simulation time.For the mechanical strength, further study in regards to the variation of internal polymer conformations and elastic property after swelling and shrinking the membrane repeatedly will be needed.

Conclusions
For polyelectrolyte membranes for fuel cell electric vehicles, we have established the method to evaluate key properties to commercialize the membrane, such as water diffusion, gas permeability, and mechanical strength.Based on coarse-graining models, for water diffusion and gas permeability through the membrane, the DPD-MC approach was applied, and for mechanical strength of hydrated membrane, CGMD was developed.As for the water diffusion, the molecular structure of the membrane with long side chain presents high diffusivity.The map that reveals the relationship between many kinds of molecular structures and water diffusion was investigated, and the direction to enhance the diffusivity by improving membrane structure can be understood clearly.With regards to gas permeability, hydrogen diffuses through both polymer and water phase, and, on the other hand, oxygen moves mainly through the polymer phase.Therefore, when improving gas permeability, the direction of molecular design is different for hydrogen and oxygen.The CGMD model clarified that the source of elasticity of hydrated membranes are hydrogen bonds which act as quasi-crosslinks between two sulfonic acids bridged by water.In order to increase mechanical strength, we should design a molecular structure with narrow water channels, but this also decreases proton conductivity.Thus, the optimal design of a polymer structure is needed, and the developed model would be useful to optimize the molecular structure in a shorter period than experiments.There are a number of industrial polymer materials other than fuel cell polyelectrolyte membranes for which good gas barrier characteristics and high mechanical strength are required.It is considered that the molecular simulation would be widely used and also the effective polymer development would be realized by adding the specified models in the future.

Figure 1 .
Figure 1.Characteristics required for fuel cell polyelectrolyte membranes.

Figure 2 .
Figure 2. Chemical formula of Nafion monomer and coarse-graining strategy employed.

Figure 5 .
Figure 5. DPD structures obtained for the sequences {A 16 C 4 } 1 , {A 8 C 2 } 2 and {A 4 C 1 } 4 .A bead is in gray, C bead in white and W bead in black.Reproduced with permission from [3].Copyright 2010 Royal Society of Chemistry.

4 Figure 6 .
Figure 6.(a) Pore radius and (b) inter-cluster distance derived from pair correlation function plotted against D' rel. .Reproduced with permission from [3].Copyright 2010 Royal Society of Chemistry.

Figure 7 .
Figure 7. Pore morphologies obtained for the sequences {A 7 [A 1 C]} 4 and {A 1 [A 7 C]} 28 .DPD bead representations of two repeat units are included.Reproduced with permission from [3].Copyright 2010 Royal Society of Chemistry.

Figure 8 .
Figure 8. W bead fraction contained within the largest water cluster vs. length which defines connection between W beads. Reproduced with permission from [3].Copyright 2010 Royal Society of Chemistry.

Figure 9 .
Figure 9. D' rel.plotted against y/x (the ratio of side chain length and branching point distance).Reproduced with permission from [3].Copyright 2010 Royal Society of Chemistry.

Figure 11 .
Figure 11.Permeation through heterogeneous media.Relation between permeability and internal heterogeneous structure.

Figure 13 .
Figure 13.Tragectory of gas particles during diffusion process through hydrated Nafion at 40 degrees Celsius and 0.7 volume fraction of Nafion.

Figure 17 .
Figure 17.Number of entanglements per chain in hydrated Nafion as function of water content.

Figure 18 .
Figure 18.Number of quasi-crosslinks per chain in hydrated Nafion as a function of water content.

Table 1 .
Dissipative particle dynamics (DPD) repulsions used for Nafion and corresponding χ parameters in parenthesis.

Table 3 .
Simulation conditions of hydrated Nafion membrane at various water content values.