Gert
Dorenbos
*ab and
Kei
Morohoshi
*b
aKnowledgenet Co., Lofty Chuo Bldg. (9F), 1-17-24, Shinkawa, Chuo-ku, Tokyo 104-0033, Japan. E-mail: gert@nano.tec.toyota.co.jp
bToyota Motor Co., Future Project Division, 1200, Mishuku, Susono, Shizuoka 410-1193, Japan. E-mail: kei@morohoshi.tec.toyota.co.jp
First published on 4th August 2010
Using dissipative particle dynamics we model phase separation within block and grafted polymers composed of hydrophobic (A) and hydrophilic, acid-containing (C) beads. The grafted polymers have their hydrophilic beads located at the end of the grafted side chains. Pore morphologies are calculated at a hydration level λ of 4 H2O molecules/C bead. Monte Carlo tracer diffusion calculations are used to model the restricted movement of water within the pore networks. For the block polymers we find that at fixed C bead fractions, or ion exchange capacity (IEC), an increase in C block length results in larger pores and increased water diffusion. For grafted polymers of equal IEC, increasing the side chain length results in a better connected pore network and increased long-range water mobility.
Broader contextBecause of their environmental friendliness, polymer electrolyte fuel cells are intensively studied nowadays, as a potential clean energy source. Widespread commercial use requires a drastic cost reduction, which has focused researchers on developing a (simplified) thermal and water management system that enables high temperature operation. However, high temperatures result in dehydrated polyelectrolytes with increased proton resistivity, leading to poor fuel cell performance (low power output). To overcome this problem requires the development of a new membrane that can be kept hydrated by facilitating diffusion of water generated at the cathode towards the anode. Since in practice experimental production and characterization of membranes takes a long time, in this work we study the correlation between polymer design and water diffusion by means of molecular simulation. The correlations are evaluated by a newly developed method that combines dissipative particle dynamics with Monte Carlo (tracer) diffusion calculations. This method allows a fast screening in which, by systematic variation of the polymer design, interesting relations are obtained between polymer sequence, pore structure and diffusive properties. |
Despite their use in fuel cells, PFSA membranes, of which Nafion® is the most popular, face several problems in technical and commercial sense. As soon as protons flow through the membrane, electro-osmotic drag8–10 can result in a non-uniform water profile,11,12 which may cause the membrane to dry out at the anode side or cause cathode flooding, especially at increased current densities. Since the proton conductivity reduces rapidly with decreasing water content2 this prevents fuel cell designers from achieving high current densities. Another drawback is that the cathode reactions should preferably take place at a high temperature, which also can cause membrane dehydration. Moreover, membrane lifetime decreases with increasing temperature due to destructive oxidation reactions. From a commercial point of view, the cost of Nafion® is too high in order for fuel cells to become widely used in every day life.
The need for high performance and low cost membranes capable to operate in fuel cells has focused research groups onto alternatives such as sulfonated poly(ether–ether–ketone) (SPEEK),13,14 sulfonated styrene/ethylene–butylene/styrene (SSEBS),15 polybenzimidazole (PBI),16 and sulfonated polyaryleneethersulfone or BPSH.17,18 These membranes also exhibit phase separation to some extent. Since a good connected water cluster network provides the essential diffusive pathways for (long-range) proton and water transport it is worth predicting pore morphologies within candidate membranes.
Preferably such morphologies should be obtained from full-scale molecular dynamics (MD) simulations. However, for mainly two reasons MD is not yet capable to probe the (mesoscale) connectivity between pores. First, the time scale at which equilibrium morphologies are established is too slow: in real membranes the polymer lengths reach up to several hundreds to thousands monomer fragments, which means that polymer backbone motions are governed by (slow) Rouse dynamics which decrease inversely with polymer length. Therefore prediction of mesoscale equilibrium pore morphologies using conventional MD requires too much CPU time. A second restriction is the large system size that is needed to probe the pore connectivity (or tortuosity): the MD simulation box should be at least larger than the characteristic length scale (i.e. the distance between pores) of the simulated morphology. This requires a system volume ∼10,000 nm3 involving O (105) atoms.
Despite these shortcomings, MD has been applied in modeling PFSA membranes within a rather small simulation box.19–25 Early work of Elliot et al.19 demonstrated phase separation for the Nafion®–water system. That study involved several hundred water molecules and very short polymers constructed of two CF3 fragments connected with a single side chain end linked with a SO3 group (i.e. very low equivalent weight (EW), cf. high ion exchange capacity (IEC)), thus ignoring any connectivity between the side chains. A MD study by Vishnyakov and Neimark20 also predicted water clustering using a more realistic presentation of the Nafion® polymers. In that study the system was limited to ∼125 nm3 and involved ∼103 water molecules. More recent MD studies also predict water clustering which can be considered as the initial onset for the formation of long-range pore networks in large systems.
MD simulations aiming to predict the sequence design dependence of structural properties of membranes have also been reported. For example, Cui et al.21 compared calculated water channel networks obtained for a short side chain PFSA polymer (Dow (EW977)) with those obtained for the longer side chain Nafion® (EW1143) polymer. The network within the Dow membrane was found to be more dispersed and better connected than within Nafion®. Jang et al.22 applied coarse grained MD (CGMD) and predicted differences between pore networks for Nafion® oligomers that had their side chains distributed uniformly and highly non-uniformly (blocky) along their backbones. They found that the phase separation involving the blocky polymer was most pronounced.
The above MD studies also report on water and/or proton diffusion coefficients that were obtained by sampling the mean square displacements (MSD) during a sampling period of, at maximum, a few ns. However, with regards to diffusion, on this timescale only local environments can be probed by MD and calculated diffusion constants cannot be representative for long-range diffusion which is governed by pore connectivity and bottleneck phenomena. A powerful demonstration of this was obtained from quasi-elastic neutron scattering (QENS) studies26,27 of diffusion in Nafion®1100. Pivovar et al.26 revealed on a picosecond level that the mobility of water molecules within the pores are very high and resembles that of pure water. They showed that the local diffusion constant of water molecules increases from 1.22 × 10−5 cm2 s−1 to 2.14 × 10−5 cm2 s−1 when the hydration increases from λ = 1 (∼3 vol%) to λ = 16 (∼33 vol%). These values are remarkably close to the pure water value of 2.3 × 10−5 cm2 s−1. At a hydration level of λ = 5 (∼14 vol%) the local water mobility is ∼1.6 × 10−5 cm2 s−1, i.e. 0.7 times that of pure water. Since the long range (μ-scale) water diffusion constant at this water content is near 2.3 × 10−6 cm2 s−1, a further decrease of a factor ∼7 must therefore be caused by the pore network morphology.
The above findings were confirmed by QENS experiments performed by Perrin et al.27 who could distinguish between 2 types of protons. For λ > 2, the major type consists of protons that diffuse very fast within the pores and their diffusion constant is nearly equal to that of protons in pure water. Independent of hydration level, a remaining small fraction corresponding to 3 protons per sulfuric acid site moves much slower (probably related to H3O+). Perrin et al.27 also showed that long-range diffusion is significantly slower than short-range diffusion. From these QENS studies it is evident that the pore connectivity determines the overall long-range diffusion constant. Modeling long-range diffusion within membranes therefore requires simulation of large-scale pore networks.
Since large size (104–106 nm3) phase separated structures cannot be obtained by conventional MD methods, coarse grained methods have been frequently applied. Mologin et al.28 simulated the morphology of Nafion® using the bond fluctuation model. Khalatur et al.29 applied the reference interaction site model. Both studies concluded that for increasing water content a bi-continuous phase appears in which elongated water–sulfuric clusters are connected to each other. Wescott et al.30 implemented a self-consistent mean field theory in the MESODYN simulation code. They found their pore morphologies were consistent overall with experimental data, although the hydrophilic phases did not form a continuous network. Yamamoto et al.31 applied dissipative particle dynamics (DPD) and revealed a bi-continuous phase, with water and sulfuric groups together forming a connected network within the Teflon® phase. The above CG studies did not consider water transport within the pore networks, such as recently challenged by Malek et al.32 They compared their results obtained from CGMD and DPD for system size up to ∼103 nm3 . For 2 < λ < 9 the diffusion constants obtained from the CGMD were in reasonable agreement with NMR data. For λ = 9, DPD derived diffusion coefficients were significantly higher as compared with CGMD.
We also used DPD to predict pore networks within Nafion® for a system size of 5.8 × 103 nm3.33,34 The pore connectivity was probed by Monte Carlo (MC) trajectory (tracer) calculations. From that study we predicted that, without knowledge of the QENS work (at that time unpublished), water mobility within the pores should be about the same as that of pure water. This result was confirmed later for a larger system size of 2.4 × 104 nm3.35Interestingly the pore networks were found to depend strongly on the EW (or chain architecture) of the Nafion® polymer. Experimental Bragg spacings obtained from literature at several water contents and for various EW agreed very well with those calculated from the DPD morphologies (within ∼7%, see Fig. 23 and Fig. 24 in ref. 35). Since, with regards to PFSA polymers, our approach turned out to be suitable in predicting trends, we extent it here much more generally to block and grafted polymers.
Our approach divides the simulation of diffusion in membranes into two parts by calculating pore morphologies and water diffusion separately. Morphologies are mapped onto a high-density grid on which a pore network is defined. Each pore network is subjected to MC trajectory (or tracer) calculations in which the particle movement is confined to the hydrophilic water containing pores. Since frozen lattices are used as input for the MC calculations, no computational effort is spent in calculating rearrangements of the polymer during the course of water diffusion. The partitioning essentially assumes that the pore networks do not change morphology while water is moving through it. This is justified for Nafion®1100 since the water diffusion constants vary between 5% and 25% compared to that of pure water when the membrane contains 5 to 30 vol% water, respectively. Slow polymer backbone motions cannot account for these high values.
A similar freezing of phase separated pore networks combined with MC tracer diffusion has also been performed recently by Zhdanov36 who used a lattice gas approach involving hydrophobic and hydrophilic interactions to create pore structures in a 3D lattice (8 × 106 nodes). The pore network was calculated for the Nafion®–water system and the percolation threshold for long-range diffusion within Nafion® was found to be ∼23 vol%. This is appreciably higher than the experimental values (around 5–10 vol% for Nafion®1100). This overestimation is probably due to the lack of connectivity between hydrophobic particles. As demonstrated earlier,35 polymer connectivity is an essential point that significantly affects long-range diffusion and percolation thresholds (see e.g. Fig. 11(d) in ref. 35). Therefore, for equal porosity (or water volume fraction, Fr(W)) different membranes are expected to show very different diffusion constants.
Our objective here is to predict how sequence design is expected to affect pore morphologies and transport properties within membranes. For the block polymers we study the dependence of diffusion and pore morphology on hydrophilic block length. For the grafted polymers these properties are studied as function of side chain separation and side chain length.
![]() | ||
| Fig. 1 Bead representations of the membrane polymers considered. A beads are hydrophobic, C beads are hydrophilic. (a) Block polymer. (b) Grafted polymer. | ||
The block polymers (Fig. 1(a)) are composed of hydrophobic and hydrophilic blocks with each repeat unit containing x A beads and y C beads. The number of repeat units is given by z. Since we are searching for trends, we do not give an explicit chemical formula for A and C beads. When one wants to model a polymer with know chemical formula this requires specification of each bead type and an adequate parameterization of the DPD interactions between them.
Water is represented by W beads, each of them contains Nw water molecules. The number of W beads is thus (λ/Nw) times the number of C beads. At fixed hydration level λ, the W bead fractions Fr(W), defined as the number of W beads divided by the total number of beads, can be calculated from the polymer sequence using eqn(1).
![]() | (1) |
For the grafted polymers the side chains are grafted regularly along the backbones (Fig. 1(b)). One repeat unit contains x consecutive A beads along the backbone. Each side chain is composed of y A beads terminated with a hydrophilic C bead. The distance between branching points is therefore equal to x and the length of the side chain is equal to y. W bead fractions are calculated according to eqn (2).
![]() | (2) |
Previously for grafted PFSA polymers we predicted the dependence of the pore morphology (and transport properties) on the value of x, while keeping the length of the side chain fixed and the number of repeat units equal to 5 (z = 5).35 The interaction parameters and bead sizes were chosen such as to model hydrated Nafion® and Dow membranes. At equal water contents (Fr(W) = constant) we found that for increasing values of x (i.e. increase in EW or decrease in IEC) water clusters become separated further apart and larger in size. In this work we will also vary the value of y (side chain length) and z (number of repeat units) systematically.
To allow a comparison between different membranes, λ is kept fixed at a value of 4, and each W bead contains Nw = 2 water molecules. Selecting a much higher value for Nw would a priori exclude the formation of certain pore networks such as the hypothetical morphology in which each hydrophilic C bead is associated with only 1 W bead dispersed within a hydrophobic A matrix. This means that the number of W beads in a simulation should in principle be larger than the number of C beads, which implies Nw < λ.
The volume of a single water molecule is 0.03 nm3 (mass density = 1 g cm−3). Therefore, the volume of each W bead equals V = 0.06 nm3. This is the same volume which we adapt for each polymer bead and is roughly equal to half the volume occupied by an aromatic ring or a CF2CF2 fragment in Nafion® (estimated from Teflon® mass density ∼2.2 g cm−3). The bead density is ρ = 3, therefore one DPD length-unit corresponds to rc = (3 × 0.06)1/3 nm3 = 0.565 nm.
Selecting a value of Nw = 1 would require a simulation that contains twice as many beads for the same system volume. Furthermore, since the time step Δt scales proportionally to Nw5/3,42 the simulation time needed to reach equilibrium-like morphologies also increases, this would diminish the advantages of using DPD. Since DPD is a coarse graining technique, hydrogen bond networks cannot be taken into account.
The time evolution is given by Newton's equations of motion,
, where ri, vi, mi, are position, velocity and mass of bead i, respectively. For non-connected beads, the interaction between two beads i and j consists of the sum of three forces: a conservative, a dissipative (drag) and a random (noise) force. The conservative force is soft repulsive and linear with the bead–bead separation:
![]() | (3a) |
rij = ri − rj,rij = |rij|, ij = rij/|rij|. | (3b) |
The dissipative force is proportional to the relative velocity between beads, while the (pair wise) random force corresponds to thermal noise that pumps energy into the system. Together they act as a thermostat. These forces are given by:
FDij = −γωD(rij)( ij·vij) ij | (4) |
FRij = σωR(rij)θij(Δt)−0.5 ij | (5) |
![]() | (6) |
The noise (σ) and friction parameter (γ) are related according to σ 2 = 2γkBT. Here σ = 3, γ = 4.5, kBT = 1. All three forces act along the line of centers conserving linear and angular momentum.
To keep adjacent beads that belong to the same polymer joined together, an additional spring force applies:
| FSij = −Crij | (7) |
The spring constant C is set at the default value of 4, at which the distance between polymer beads that are connected with each other is about the same as the first peak in the pair correlation function of non-connected DPD beads.37 The appearance of Δt−0.5 in eqn (5) ensures consistent diffusion.37 Since the total force on each bead depends on the velocity viaeqn (4), a modified Verlet integration algorithm is used31,35,37 with the empirical factor κ = 0.6535 and time step Δt = 0.05. In this way the temperature control was kept constant near kBT = 1.0. The mass and length (i.e. interaction range) of the DPD beads are scaled to a value of 1.
The repulsions aij are chosen such that the compressibility of water is reproduced for a system that contains only W beads. Since Nw = 2 and the bead density ρ = 3, we obtain for interacting W beads a value of aij(WW) ∼ 51 according to eqn (14) and eqn (16) in ref. 37. Within a DPD simulation the interactions are always repulsive. This implies that hydrophilic interactions are those for which the repulsions are lowest, and repulsions between incompatible (hydrophobic) beads will be highest. Here we adapt the following values for the repulsions:
| a(WW) = a(CC) = a(AA) = 51 |
| a(CW) = 41 (hydrophilic; Δa = −10) |
| a(AW) = a(AC) = 60 (hydrophobic; Δa = 9), or |
| a(AW) = a(AC) = 70 (strongly hydrophobic; Δa = 19) |
In earlier studies by Groot and Warren37 and Groot and Rabone42 obtained several linear relationships between the excess repulsions (Δa) and the corresponding Flory–Huggins χ-parameters:
| χij = const × Δaij | (8) |
The values for the constant prefactor depend on Nw and whether each bead represents a single monomer or a connected polymer bead. We did not intend to determine the exact value of the prefactor for Nw = 2. Based on the values given in ref. 37 (eqn (24) & (27)) and ref. 42 (eqn (18)) a reasonable guess for the prefactor value is ∼0.28 ± 0.02, for which the adapted repulsions can be converted to the following χ-parameters:
| χ(CW)∼ −2.8 (hydrophilic) |
| χ(AW) = χ(AC)∼ +2.5 (hydrophobic), or |
| χ(AW) = χ(AC)∼ +5.3 (strongly hydrophobic) |
The system size is cubic with cell dimension of L = 30 or 35 DPD length units (case dependent). Each simulation involves the dynamics of 81
000 or 128
625 beads ( = ρ × L3) and contains several hundreds to a few thousand polymer chains. The volumes are 4860 nm3 (L = 30) and 7718 nm3 (L = 35). The lengths of the box edges are 17.0 nm (L = 30) and 19.8 nm (L = 35). Periodic boundaries are applied in each orthogonal direction. Morphologies sampled after 20
000 time steps (1000 DPD time units) are used as input for the MC tracer diffusion calculations. This period is the same as in other studies31,33–35,38 and based on the evolution of the W bead pair correlation functions assumed to be sufficient to establish equilibrium-like morphologies.31,35 Rather than sampling our results over an ensemble of morphologies per sequence we have chosen to present the results that were obtained for a single morphology per sequence. In Section III of the results and analysis this will be discussed in more detail.
In order to construct a pore network from the DPD bead positions, a cubic lattice that contains 2.7 × 107 (3003) or 4.3 × 107 (3503) nodes is defined for system size L = 30 or L = 35, respectively. The resolution of the grid is therefore 0.1 × rc = 0.0565 nm and is considered high enough to cover all morphological details that can be deduced from the DPD bead positions. The large number of beads (O (105)) and huge number of nodes (O (107)) requires some programming effort. During the pore network construction procedure we calculate for each node the distance to the nearest polymer (A or C) bead and the distance to the nearest W bead. Periodic boundaries are applied when calculating these distances. Nodes for which the nearest DPD bead is a W bead belong to the water pore network. The remaining nodes are part of the polymer phase. For all morphologies the pore channel volume fractions (fraction of total number of nodes that are part of the pore network) are only very slightly (<0.005) larger than the nominal water contents Fr(W). This also implies that locally within both water and polymer phase the bead densities are equal to ρ = 3.
In the MC calculations N tracer particles (L = 30: N = 1000; L = 35: N = 2000) are put randomly on nodes that belong to the pore network. At each time step, a jump trial towards a nearest node in one of the 6 orthogonal directions is randomly selected for each particle independently. Move trials are accepted for those particles for which the nearest node in the selected direction is also part of the pore network. The diffusion constant is obtained from the mean square displacement (MSD) of N trajectories:
![]() | (9) |
Fig. 2 shows displacement curves obtained for 3 block polymers of equal length but differing in the way hydrophilic C beads are distributed along the backbone. The curve for pure water (obtained by accepting each jump trial) is also shown. The horizontal line drawn at 1.225 × 105 (=3502) is indicative for the system size. For all three membranes the whole pore network has been probed indicating that within each membrane a percolating pore network exists.
![]() | ||
| Fig. 2 MSD sampled over N = 2000 particle trajectories versus MC time obtained for the block polymers sequences A16C4, (A8C2)2, and (A4C)4. System size L = 35. Water contents Fr(W) = 0.286. The horizontal line drawn at 1.225 × 105 is an indication of the system size. Inset: short time behavior of MSD in membranes relative to pure water case. | ||
The calculated tracer diffusion constants are expressed relative to their values in pure water as follows:
![]() | (10) |
D′(H2O)rel. can thus be obtained from the ratio of the slopes of the MSD curves such as those displayed in Fig. 2. The prime in eqn(10) is introduced to emphasize that local water mobility within the pores is assumed to be equal to that of pure water and that, for the reasons just mentioned, extra care should be taken when calculated water diffusion constants are compared with experimental data. Nevertheless, we expect that predicted trends are valid when comparing values of D′(H2O)rel. obtained for membranes that belong to the same “family” (e.g. Teflon®-like or hydrocarbon-like) and that carry the same acid moieties. The inset in Fig. 2 demonstrates that the short time behavior for all three membranes indeed approaches that of pure water.
At first sight, one might be tempted to estimate water transport directly from the diffusion of W beads during a DPD simulation. However, we selected the MC approach mainly for two reasons. First, DPD allows some extent of bond crossing of polymer beads that lift up entanglements. DPD polymers are therefore diffusing faster than real polymers (or polymers in a realistic MD calculation), resulting in rearrangements of the pore structures while water is diffusing through it. In our approach the pore network is fixed in time. A second reason is that even if during a DPD calculation the pore network does not change significantly, any estimation of long-range water diffusion sampled from the movements of DPD W beads would require very long DPD runs. In order to estimate diffusion constants, the whole system should be probed. Since the system size considered involves a volume of ∼8000 nm3, this would require the water beads to travel ∼20 nm.
We calculated for various values of Nw(= 1,2,3,4,5,6) and for ρ = 3 the diffusion of DPD beads using the repulsions that approximately reproduce the water compressibility by setting aij = 25 × Nw.37 Using the definition given in eqn (9) (see also ref. 42), we obtained D∼0.3 × Nw−0.59. For Nw = 2 the diffusions are: D(Nw=2)∼0.2. For these beads to probe a system size of L = 35 would require 352 × 6 × 0.2 ∼O(103) DPD time units, or ∼20
000 additional time steps. However, within the pore networks presented in the next section, the overall diffusion is one or more orders of magnitude smaller than that of pure water, which would require a simulation of O(105) to O(106) additional time steps for the W beads to probe the whole pore network.
![]() | ||
| Fig. 3 (a) D′(H2O)rel. plotted against water volume fraction, Fr(W). The results were obtained for (AxCy)z block polymers. Symbols: filled circles: a(AC) = a(AW) = 60, triangles: a(AC) = a(AW) = 70. System size L = 30. (b) Comparison of data in (a) (filled circles) with results obtained for system size L = 35 (crosses). Repulsions a(AC) = a(AW) = 60. | ||
We verified and confirmed the results that were obtained for a(AC) = a(AW) = 60 by performing additional simulations for a larger system size (L = 35) and varying number of repeat units (z). This was done for sequences for which x/y = 4, 5, 6, and 7. For each value of x/y (or Fr(W)), the total length ((x + y) × z) of the polymers was kept fixed but the distribution of the hydrophilic C beads along the backbones was systematically varied with the hydrophilic blocks containing respectively 1, 2, or 4 consecutive C beads. The results are shown in Fig. 3(b) (crosses) and compared with the data from Fig. 3(a) for which z = 5 (filled circles in Fig. 3(b)). It is confirmed in Fig. 3(b) that an increase in hydrophilic block length from 1 to 2 consecutive C beads results in increased diffusion. Note that for the polymers with hydrophilic block lengths that consist of 4 consecutive C beads (sequences A28C4, A25C4, A20C4, and A16C4 in Fig. 3(b)) the water mobility increases significantly.
The results shown in Fig. 3 can be explained by inspecting the pore morphologies. DPD structures that were obtained for the sequences A16C4, (A8C2)2 and (A4C)4 are displayed in Fig. 4. For each of these sequences the total number of A beads equals 16 and the number of C beads is 4 (x = 4y). According to eqn(1), the water volume fraction is Fr(W) = 0.2857. Since ρ = 3 and L = 35, the number of W beads equals 0.2857 × ρ × L3 = 36
749. Although we believe that the morphologies are equilibrated, we cannot classify them in a particular structure (lamellar, cylindrical-like, spherical, perforated lamellar etc.). From Fig. 4 it is evident that the sequence with the longest hydrophilic C-block (A16C4) shows the most pronounced, well connected, pore channel network. The morphology appears bi-continuous with both polymer and water phases forming a percolating and connected network. For the sequence (A8C2)2, for which the lengths of the C blocks are reduced by a factor 2, the pores still appear quite well connected, but remarkably smaller in size and more irregular in character. When reducing the hydrophilic block lengths further, ((A4C)4 polymer) the pore structure is no longer well defined and W beads seem to be nearly uniformly distributed, but water clustering is still apparent for the (A4C)4 polymer. Since long-range diffusion also occurs for this structure, this indicates the presence of a percolating and continuous water network. Some isolated water clusters might be present that cannot participate in long range diffusion, and strictly speaking this structure might not be bi-continuous.
![]() | ||
| Fig. 4 DPD structures obtained for the sequences A16C4, (A8C2)2 and (A4C)4. A beads (red), C beads (yellow), W beads (blue). Iso-surfaces of the W bead density are shown at a density of 1.5 (A16C4) and 1.0 ((A8C2)2 and (A4C)4). The repulsions are a(CW) = a(AW) = 60. Fr(W) = 0.286. System size L = 35. | ||
The sequences for which x = 5y ((A20C4, (A10C2)2 and (A5C)4), x = 6y (A24C4, …), and x = 7y (A28C4, …) essential show a similar trend when decreasing the hydrophilic C block length from 4 to 1 consecutive C beads and are not displayed here.
To estimate the size of the water clusters and the average separation between them we calculated the pair correlation function g(r) of the water beads. Fig. 5 displays g(r) plots that were calculated for the structures shown in Fig. 4. The pair correlation functions obtained for the sequences with x = 7y (A28C4 etc.) are also displayed in Fig. 5. Points located at distances from a water bead where the values of g(r) are larger than 1 are more likely to be part of the pore network than not. Distances for which g(r) drops below the value of 1, correspond to anticorrelation, and are an indication of the cluster boundary, i.e. pore radius Rpore. The average separation between clusters, Dcl-cl, can be derived from the position of the 2nd maximum in the g(r) plots. The distance at which g(r) starts to drop below the value 1 increases with increase of C block length. Also the 2nd maximum shifts towards larger distances with increase of the length of the C blocks. This indicates that for these membranes the cluster (or pore) radius and the distance between clusters increase with the length of the hydrophilic blocks, in agreement with the observations in Fig. 4. Note that the repulsions between beads (viaeqn 3(a)) cause the 1st maximum in g(r) to also persist for a highly dispersed distribution of W beads. This is demonstrated by the two curves in Fig. 5 that were calculated for (A7C)4 sequence and Fr(W) = 0.286 (upper curve) and Fr(W) = 0.2 (lower curve) by setting all repulsions at a value of 51.
![]() | ||
| Fig. 5 Upper symbols: g(r) of the W beads for the morphologies shown in Fig. 4. Lower symbols: g(r) obtained for (A7C)4, (A14C2)2, and A28C4. For A28C4 sequence the approximate positions for Rpore and Dcl-cl are indicated. Solid curves are g(r) (W bead) calculated for water dispersed in the (A7C)4 polymer matrix for Fr(W) = 0.2 (lower curve) and Fr(W) = 0.286 (upper curve) by setting all repulsions at aij = 51. System size L = 35. | ||
Rpore or Dcl-cl are plotted versus D′(H2O)rel. in Fig. 6(a) and Fig. 6(b), respectively. The graphs clearly demonstrate various relations among the transport properties (water mobility), morphology (pore size and distance between pores) and sequence design (hydrophilic C block-length).
![]() | ||
| Fig. 6 (a) Pore radius and (b) inter cluster distance derived from the pair correlation function plotted against D′(H2O)rel. Dashed lines serve as eye guide and connect data points obtained for membranes that contain the same amount of water (in volume fraction Fr(W)). Error bars for Rpore and Dcl-cl are estimated to be about ±0.3 nm. | ||
Cluster size distributions of the DPD water beads were calculated for a range of fixed values of Lbond (Fig. 7). Here Lbond is defined as the largest distance for which two water beads are considered to be connected with each other and belong to the same cluster. For a large value of Lbond (e.g.Lbond = 0.6 nm), nearly all water beads are part of the same percolating cluster. When selecting a very small value for Lbond, then the possibility of two nearby water beads being connected with each other will become very small (due to repulsive interactions between DPD beads, eqn (3a)). The number of W beads contained within the largest cluster increases therefore as function of Lbond. This is shown in Fig. 7(a–b) where for several block polymers the largest cluster size is plotted against Lbond. The size of the largest cluster is expressed relative to the total number of W beads that is contained within each membrane, i.e. (number of W beads within largest cluster)/(total number of W beads). For membranes that contain an equal fraction of C beads over A beads (i.e. y/x = constant), an increase in hydrophilic C block length results in a shift of the sigmoidally shaped curves towards the left. This means that for a fixed value of Lbond the size of the largest water cluster increases with increasing block length. This trend is observed for each of the four polymer series for which x = 4y, x = 5y, x = 6y, or x = 7y. Fig. 7(c) expresses D′(H2O)rel. against those values of Lbond for which 80% of the total number of water beads are located within the same percolating cluster. The relations between Lbond (80%) (from Fig.7(c)) and Rpore (from Fig. 6(a)) or Dcl-cl (from Fig. 6(b)) are compiled in Fig.8(a) and Fig.8(b). We conclude that for the block polymer sequences an increase in block length results in:
![]() | ||
| Fig. 7 (a–b) Fraction of W beads contained within largest water cluster versusLbond calculated for several polymers. (c) D′(H2O)rel.versus distance Lbond (80%) that is needed such that 80% of the W beads take part in the largest (percolating) cluster. | ||
![]() | ||
| Fig. 8 (a) Rpore and (b) Dcl-clversus the value of Lbond that is necessary in order for the largest cluster to contain a fraction of 0.8 of the total number of water beads. | ||
i. Larger pores (Rpore) (Fig. 8(a)) and larger distance between them (Dcl-cl) (Fig. 8(b)).
ii. Better connected water clusters leading to smaller Lbond (80%) values (Fig. 7(c), Fig. 8(a–b)), and therefore less isolated clusters.
iii. Enhanced water diffusion (Fig. 6(a–b), Fig. 7(c)).
Fig. 9(a) displays the calculated values of D′(H2O)rel. against water content, Fr(W). These results were obtained for strongly hydrophobic A beads (a(AW) = a(AC) = 70). An interesting trend can be observed: for polymers with a fixed number of A beads per repeat unit (i.e. x + y = constant), the water mobility depends strongly on the distribution of A beads within the polymer sequence. For these polymers, the ones with long side chains and a short inter branching distance show faster water diffusion than those with short side chains and a long inter branching distance, as illustrated for two cases in the inset of Fig. 9(a). This trend is observed independent of the number of A beads that are contained within one repeat unit. When the data are re-plotted against Q = y/x (i.e. the ratio of the length of the side chains (y) and the inter-branching point distance (x)) logarithmic fits can be obtained for a fixed number (x + y) of A beads per repeat unit within the polymers. These fits are given as solid lines in Fig. 9(b).
![]() | ||
| Fig. 9 (a) D′(H2O)rel.versus water volume fraction. Inset: DPD bead representations of the sequences (A[A6C])5 and (A6[AC])5. (b) D′(H2O)rel. plotted against Q( = y/x). Lines are fits obtained for data points with same fraction of C beads (i.e.x + y = constant) within the polymer sequence. | ||
In order to explain the trend shown in Fig. 9, we took a closer look at the pore morphologies. Morphologies that were obtained for the sequences (A7[AC])4 and (A[A7C])28 are shown in Fig. 10(a).
![]() | ||
| Fig. 10 (a) Pore morphologies obtained for the sequences (A7[AC])4 and A[A7C])28. DPD bead representations of two repeat units are included. The water volume fraction Fr(W) = 0.182. Iso-surfaces are shown for a W bead density of 1.0. System size L = 35. (b) Pore area against enclosed volume fractions that were calculated at several iso-values. Dashed vertical line is drawn at Fr(W) = 0.182. Inset: specific area against enclosed volume. | ||
The curved surfaces drawn in Fig. 10(a) enclose volumes where the (local) W bead density is larger than 1 (i.e. iso-value is 1.0) and are illustrative for the pore network. The sum of the enclosed volumes are a fraction, 0.2257 and 0.2157, of the total system volume for sequences (A[A7C])28 and (A7[AC])4, respectively. The networks do not differ much in appearance, although the W bead density iso-surfaces seem to be better connected and less spherical for the sequence with the longer side chains, (A[A7C])28. We calculated for iso-values ranging from 0 to 3 the areas of iso-surfaces and the summed volume fractions that these surfaces enclose. From Fig. 10(b) we observe that the area of iso-surfaces for the (A[A7C])28 sequence is less than for the (A7[AC])4 sequence when the enclosed volume fractions are larger than ∼0.15. Since specific pore areas (inset Fig. 10(b)), defined as (total pore area)/(total pore volume), tend to be smaller for the longer side chain sequence we would expect them to be more spherical in character. Furthermore, for pores with a more spherical character, we would expect a lowering in diffusion, as was predicted for the Nafion® membranes.35 However, from Fig. 9 it is in general the sequences with the longer side chains that show enhanced diffusion as compared to the shorter side chain sequences. In Fig. 11(a) this is once again demonstrated where for the sequence (A[A7C])28 diffusion is much faster than for (A7[AC])4.
Apparently an increase in side chain length does not result in more spherically shaped pores and consequently lower diffusion constants. Therefore it is likely that the decrease in specific areas is due to an overall increase in pore size (possibly accompanied with less isolated pores). Indeed, from the pair correlation functions (Fig.11(b)) some clear differences can be drawn. For the polymer with the longer side chains, (A[A7C])28, a larger Dcl-cl and also, although less pronounced, a larger value for Rpore is observed. This explains the trends with regards to water diffusion that we predict in Fig.9, and the more specific cases displayed in Fig.11(a).
From the above we conclude that for the grafted polymers an increase in pore size (or Dcl-cl) correlates with enhanced water diffusion. Since this trend is essentially the same as was obtained for the block polymers we might also anticipate that the size of the largest water clusters for fixed Lbond is larger for those sequences that exhibit higher water diffusion. This is confirmed in Fig. 11(c), where for the polymer sequence with the longer side chains, (A[A7C])28, the number of water beads within the largest cluster is indeed higher than for the shorter side chain sequence. This is observed for each value of Lbond.
The consistency of these 2 trends (Fig. 11(b)–(c)) was found to be generally valid for the grafted sequences (for which x + y = constant) that we considered. As a final illustration in Fig. 11(d) the size of the largest water cluster is plotted versusLbond for the sequences (A6[AC])5 and (A[A6C])5 (system size L = 30). Again the sequence with the longest side chains, i.e. (A[A6C])5, contains the highest number of water beads within the largest W cluster for fixed values of Lbond.
![]() | ||
| Fig. 11 (a) MSD versus MC time. (b) Pair correlation function of the water beads. (c) W bead fraction contained within the largest water cluster versusLbond. Data in (a–c) were calculated for the sequences (A7[AC])4 and A[A7C])28 and system size L = 35. (d) W bead fraction contained within largest water cluster versusLbond obtained for (A6[AC])5 and (A[A6C])5 and system size L = 30. | ||
In the figures each datum point was obtained for only one pore morphology. The system size was slightly smaller (L = 30 or 35) than in previous work35 (L = 40) where we found that convergence of calculated diffusion constants within acceptable error bars was reached for system size L∼30 and that it was sufficient to calculate only one morphology. We estimate that the relative error for the obtained (tracer) diffusion coefficients for most morphologies to be less than ∼10%.
To illustrate this, we generated for each of the sequences (A[A6C])5 and (A6[AC])5 four independent DPD structures which served as an input for the MC calculations. MSD curves sampled for these eight structures are shown in Fig. 12(a). The average deviation of the slopes from their mean value is 7% and 5% for sequence (A[A6C])5 and (A6[AC])5, respectively. Fig. 12(b) displays the MSD at t = 2 × 106 MCS as function of the number of trajectories. For N>∼500 the MSD does not vary substantially, which indicates that probing over 1000 trajectories is sufficient to reach convergence for each single morphology. We note that when one wants to make accurate estimates for percolation thresholds then an averaging should be performed over results obtained for several morphologies at various values of Fr(W), which was not the intention of this study.
![]() | ||
Fig. 12 MSD curves calculated for sequence (A[A6C])5 and (A6[AC])5 containing 20 vol% water. Data points represented by symbols were obtained for two separate DPD run sampled at 80 000 time steps. (b) MSD at 2 million MCS versus the number of trajectories. | ||
All results presented in this work were obtained for DPD structures sampled at t = 20.000 time steps which seems to be long enough in order to obtain equilibrium-like structures. To demonstrate this we sampled pore morphologies up to 80.000 time steps for the sequences (A[A6C])5 and (A6[AC])5. MSD curves obtained for these 2 sequences are included in Fig. 12(a) (symbols). Water diffusion is not drastically different for DPD morphologies sampled over these long equilibration periods.
The trend that for the longer side chain sequences water clusters are larger for fixed values of Lbond was verified also for these 8 morphologies. Fig. 13 shows the W bead fractions contained in the largest cluster versusLbond. The 8 curves split up in two groups, each group containing either long side chain or short side chain sequences.
![]() | ||
| Fig. 13 W bead fraction contained within the largest water cluster versusLbond obtained for four independent morphologies for sequences (A6[AC])5 (empty symbols) and A[A6C])5 (solid symbols), respectively. Drawn lines are averages obtained for the two sequences. | ||
For the block polymers the pore connectivity, or D′(H2O)rel., increases with hydrophilic block length. This result was obtained for both type of interactions a(AC) = a(AW) = 60 and becomes more pronounced for a(AC) = a(AW) = 70 (Fig. 3(a)). Our prediction that for the grafted polymers the short side chain sequences reveal the slowest diffusion were also obtained, although again less pronounced, when the morphologies were generated for the repulsions a(AC) = a(AW) set at a value of 60 (results not presented here).
In an earlier study35 pore morphologies of the short side chain Dow1000 membrane and long side chain Nafion®1000 membrane were calculated. The polymers were modeled by hydrophobic backbone A beads, representing (CF2)4 fragments, with side chains composed of a single hydrophilic C bead (CF2CF2SO3H) (Dow) or a hydrophobic B bead (OCF2CFCF3O) with a pending hydrophilic C bead (Nafion®). Each W bead contained 4 H2O molecules (Nw = 4), with the repulsions more than twice higher than assumed here. Despite the differences in bead volumes and adapted repulsions, also in that study the diffusion of water within the short side chain (Dow1000) membrane was predicted to be less than within the longer side chain (Nafion®1000) membrane for water contents Fr(W) = 0.1, 0.2 and 0.3.
Apparently, hydrated grafted polymers composed of hydrophobic backbone and hydrophobic side chains with hydrophilic side chain terminus, generally tend to result in larger pores for sequences with the longer side chains. This can be rationalized by realizing that phase separation is governed mainly by two processes; (1) minimization of the hydrophilic-hydrophobic interface and (2) maximization of the number of hydrophilic-hydrophilic (C–W) contacts. Clearly process (1) stimulates phase separation by formation of large water clusters which are located within a hydrophobic, A bead containing, phase. However, the size of these clusters is restricted due to the presence of the hydrophilic C beads which due to process (2) tend to arrange themselves at the water–polymer interface. For two polymers of the same EW (or IEC) but different side chain length it is therefore those containing the longer side chains which have the highest flexibility (or entropic freedom) to accommodate to a morphology with larger water clusters.
The results obtained for the grafted polymers are more salient in nature and significantly more interesting. Of particular interest is the trend shown in Fig. 9(b) that for polymers with equal IEC (i.e. x + y = constant), those with the longer side chains are showing faster water diffusion. Since the protons diffuse through the same pore network as water does, this means that increased conductivity is expected for the longer side chain sequences. Note that this is only valid when comparing membranes that have the same EW or IEC and thus equal H+ concentration. Because λ = 4, this is valid for those data points in Fig. 9(a) obtained at equal Fr(W), or likewise for those data points that are connected with solid lines (i.e. x + y = constant) in Fig. 9(b).
We tried to find confirmation of our predictions in literature. Sumner et al.45 and Creager et al.46 compared the proton conductivities of bis[(perfluoroalkyl) sulfonyl] imide membranes for several hydration levels with those measured for Nafion®1100. Both membranes have the same IEC but differ in side chain length (being shorter for the Nafion®1100 membrane) and distance between the branching points (being on average ∼30% longer for the Nafion® membrane). At a RH of 31% both membranes contain the same amount of water, (λ∼4) but the conductivity of the sulfonyl ionomer was found to be twice that of Nafion®1100. The aim of their study might have been to increase the acidity of the ionic group, thus increasing the amount of solvated (mobile) protons within the pore network leading to higher proton diffusion and higher conductivities as well. But these results can also be explained by our prediction, and as suggested by Sumner et al.,45 that the pore networks within the membrane containing the longer side chains are better connected.
Creager et al.46 estimated the error bars for λ to be about equal to 1, which therefore does not allow us to make strong conclusions. Interestingly, the same research group reconfirmed their findings and measured at a hydration level λ = 3 that the Imide1075 (λ = 3 at RH = 58%) membrane has a 3.6 times higher conductivity than the Nafion®1100 membrane (λ = 3 at RH = 31%), see Table 1 and Table 2 in ref. 47.
It would be more convincing to make a comparison between diffusion constants in grafted membranes that have exactly the same end grafting acid moieties. Such membranes are available for the PFSA membranes. However, with exception of Nafion® and Dow membranes their (average) inter-branching point distance and side chain lengths are not precisely known, which makes validation of our predictions difficult.
The fact that Dow membranes with smaller side chains perform better in fuel cells than Nafion® membranes which are composed of longer side chains, at first sight seems to be in conflict with our prediction for grafted polymers, but in fact is not. The reason why Dow membranes perform better has been thoroughly discussed by Kreuer et al.48 One reason is that a higher crystallinity and glass transition temperature makes Dow membranes mechanically more stable than Nafion® membranes of the same IEC. For fuel cell applications this means that a Dow membrane with higher IEC, or lower EW, than Nafion® can be used at higher temperatures. Furthermore, these low EW membranes have two advantages: (1) a higher proton density and (2) because they absorb (at fixed RH) about the same amount of water/sulfuric acid site, they contain a higher water volume fraction, Fr(W). Overall, this means that for low EW membranes more protons can diffuse in better connected pores, which makes Dow membranes more attractive than Nafion® membranes.
In Fig. 14 the solid lines are logarithmic fits to the data points in Fig. 9 that are obtained for membranes with same EW (x + y = constant). The thick arrow illustrates that reducing the side chain length from y = 3 to y = 2 while keeping the inter-branching distance x fixed at a value of 5, results in increased diffusion. This is observed for all possible transitions for which Δy = −1. This means in general that lowering the membranes' EW by solely reducing the side chain length results in faster diffusion at equal λ, and better connected pores.
![]() | ||
| Fig. 14 General trends obtained for grafted polymers. Thick arrow: decreasing the side chain length while keeping x fixed results in enhanced diffusion. Thin arrow: Increasing the side chain length while decreasing x leads to enhanced diffusion. | ||
However, when comparing membranes with the same EW (i.e. x + y = constant) then the ones with the longer side chains show the highest diffusion coefficients, as illustrated by the thin arrow in Fig. 14. Also this trend can be deduced from experiments when comparing water diffusion or proton conduction of Dow membranes with those obtained for Nafion® membranes. Kreuer et al.48 made such a comparison when they plotted water diffusion and proton conductivities within Dow1084 and Nafion®1100 membranes in Fig. 4 and Fig. 5 of ref. 48, respectively. Up to water contents of Fr(W)∼0.25, water diffusion within Dow1084 is systematically lower than within Nafion®1100. At Fr(W)∼0.08 the water diffusion constant in Dow1084 is less than half of that within Nafion®1100, while at Fr(W) = 0.2 this difference is much less pronounced but still several tens of percent. Also the proton diffusion coefficient and conductivity for these two membranes is lower for the Dow membrane, as can be observed from Fig. 5 and Fig. 6 (for λ<∼7.5) in ref. 48. Since these two membranes have nearly the same EW or IEC but with the shorter side chain Dow membrane revealing the slowest diffusion, this experimental observation is in line with the general trend displayed in Fig. 14 (thin arrow). Although the DPD parameterization leading to the trends in Fig. 14 was not chosen to model Dow or Nafion® membranes, it is interesting that when such a parameterization is done also slower diffusion for Dow1000 compared to Nafion®1000 was predicted in recent DPD-MC calculations.35 For Fr(W) = 0.1, 0.2 and 0.3 diffusion within Dow1000 was found to be systematically less than within Nafion®1000 membranes. 35
DPD morphologies were obtained by assuming very flexible polymer (main and/or side) chains. This allowed us to focus on one “driving force” that is involved with pore network formation. The use of bending potentials may affect the pore networks and consequently the estimated values for D′(H2O)rel. For very stiff polymers the bending forces will compete with the forces that govern phase separation and may play a more decisive role in establishing the morphologies. From design point of view it will be interesting to see to what extent the trends such as displayed in Fig. 3 and Fig. 14 are affected when bending forces are included. The introduction of different bending forces for hydrophilic and hydrophobic parts in block polymers, or for side chain and main chain fragments in grafted polymers, are expected to be interesting design parameters in order to enhance water or proton diffusion within PEM's.
One disadvantage of DPD is the incapability to simulate crystallization. This occurs in membranes and might to some extent obstruct phase separation. Several studies on PFSA membranes reveal the presence of a certain crystalline fraction.6,47–50 X-ray diffraction measurements have also revealed an appreciable amount of crystallinity in the sulfonyl imide membranes,46,47 and even for SPEEK membranes a small crystalline fraction has been suggested.51 Realizing that polymer crystallization is a non-equilibrium process involving various extents of crystalline order,52 means that the amount and size of the crystallites may, besides the processing history of the membrane, also depend on chain architecture. Current DPD cannot deal with processing history and crystallization phenomena, so this should be challenged by other methods.
Summarizing, for block and grafted type polymer membranes composed of hydrophobic and hydrophilic fragments we find that upon hydration the pore morphologies strongly depend on chain architecture. For the block polymers longer hydrophilic blocks lead to larger pores that are better connected resulting in increased water diffusion. For membranes (that belong to the same family) that are composed of grafted polymers with end-grafted hydrophilic acidic sites, at equal IEC (or EW), the ones with the longest side chains reveal the highest water diffusion constants.
| This journal is © The Royal Society of Chemistry 2010 |