Dependence of pore morphology and diffusion on hydrophilic site distribution within hydrated amphiphilic multi block co-polymer membranes

Polymer 54 (2013) 5024e5034

Contents lists available at SciVerse ScienceDirect

Polymer journal homepage: www.elsevier.com/locate/polymer

Dependence of pore morphology and diffusion on hydrophilic site distribution within hydrated amphiphilic multi block co-polymer membranes G. Dorenbos Knowledgenet Co., Lofty Chuo Bldg. (9F), 1-17-24, Shinkawa, Chuo-ku, Tokyo 104-0033, Japan

a r t i c l e i n f o

a b s t r a c t

Article history: Received 28 March 2013 Received in revised form 26 June 2013 Accepted 2 July 2013 Available online 9 July 2013

Dissipative particle dynamics (DPD) combined with Monte Carlo (MC) tracer diffusion calculations are used to study phase separation and diffusion within hydrated (amphiphilic) alternating multi block copolymer membranes. The co-polymers are composed of hydrophobic (A) and hydrophilic (C) fragments. The hydrophobic A block length is alternating short (x A fragments) and long (y A fragments). One repeat unit is represented as AxCAyC, with y
x. The phase separated morphologies that were generated for 18 architectures by DPD at a water content of 25 percent by volume reveal that water is contained within a pore network with hydrophilic C fragments located near the pore boundary. The morphologies are mapped onto a cubic grid on which MC (tracer) trajectory calculations are performed in which particle movement is restricted to the water containing pore networks. For architectures for which the hydrophilic fragment fractions are the same (same ion exchange capacity), an increase in difference of hydrophobic block lengths (or y  x) result in a linear increase in inter pore distance and a significant increase in long-range diffusion through the pore networks. Ó 2013 Elsevier Ltd. All rights reserved.

Keywords: Diffusion Morphology Ionomer

1. Introduction Currently NafionÒ is mostly applied as a proton exchange membrane (PEM) in fuel cells. Its purpose is to separate the cathode from the anode while simultaneously acting as a conductor for protons. NafionÒ is a per-fluorinated sulfonic acid (PFSA) polymer composed of a hydrophobic Teflon backbone to which side chains are grafted with a pendant hydrophilic SO3H group (Fig. 1(a)). The water uptake, expressed as the number of water molecules l per sulfonic site, increases with relative humidity [1e8]. Phase separation results in a water containing pore network with pores a few nm in diameter surrounded by polymer phase [6e8]. Since these pores are the diffusive pathways for water and protons this makes their degree of connectivity essential for obtaining high proton conductivities. Because at low humidity the proton conductivity decreases rapidly [1e4,6,9,10], this requires a fuel cell design strategy that provides sufficient hydration of the PFSA membranes during operating conditions. The development of membranes with good connected pore networks would be an important step toward wide scale application of fuel cells.

E-mail address: [email protected]. 0032-3861/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.polymer.2013.07.007

In order to find alternative polymers that can provide sufficient proton conductivity and mechanical strength at working conditions several alternative grafted or block polymer membranes have been synthesized. As for the PFSA membranes, these alternatives have acidic, hydrophilic, sites located within the backbones of the block polymers or within the side chains of the grafted polymer. Examples are the block polymers poly-benzimidazole (PBI) [11], sulfonated poly-ether-ether-ketone (SPEEK) [12,13], sulfonated styrene/ ethyleneebutylene/styrene (SSEBS) [9] sulfonated poly aryl ether ketones (SPAES) [14] and BPSH [15], sulfonated polyimide (SPI) [16], and the grafted [bis] (perfluoroalkyl) sulfonyl imide perfluorinated ionomers [17]. Due to the presence of hydrophilic (acidic) sites within the molecular architectures, membranes composed of these ionomers swell and phase separate under humid conditions, resulting in water containing pore networks. The extent to which the pores are connected may, besides the hydration level, depend on several parameters like the molecular architecture, acidic site density, side chain length, side chain distribution, backbone and/or side chain stiffness, etc. Theoretical approaches such as molecular dynamics (MD) [18e30], coarse grained MD [31,32] and meso scale methods such as dissipative particle dynamics (DPD) [33e45] may be helpful to gain insight in how the pore size and their connectivity is affected by these parameters. With respect to fuel cell membranes MD has frequently been applied

G. Dorenbos / Polymer 54 (2013) 5024e5034

Fig. 1. (a) Repeat unit of the grafted NafionÒ (EW ¼ 1143) polymer and bead parameterization determined in previous work [39] A0 : CF2eCF2eCF2eCF2, B0 : OeCF2e CF(CF3)eO, C0 : CF2eCF2eSO3H. (b) DPD representation of a repeat unit of the model multi-block co-polymers considered here which are composed of two hydrophobic blocks that contain x and y A beads, respectively, covalently bonded by a hydrophilic C bead.

to model phase separation within the grafted (PSFA) [18e30] membranes and predicts indeed water clustering in agreement with experiment. From MD diffusion constants can be obtained by sampling mean square displacements of water and/or protons during a period of at most several ns. During this time only the local environment of the diffusing species can be probed and obtained diffusion constants do not represent experimentally obtained long range diffusion constants. For this reason simulation of long range diffusion requires a meso scale probing of the pore connectivity. This essentially means that the system size should be significantly larger than the characteristic distance of the pore network. As an example, for the PFSA membranes the characteristic distance as estimated from the Bragg spacing is around 5 nm [6e8,10], and therefore the obstructions (bottlenecks) in the diffusion pathways can only be probed for system volumes of the order of 104 nm or more. For limited system size, water and proton diffusion through hydrated Nafion1143 (1143 is the equivalent weight (EW) and refers to 1143 g polymer/mol SO3H) was calculated by Devanathan et al. [25]. They observed an appreciable increase in diffusion of these species near l ¼ 5. The size of the simulation box was at most w60 nm3 (at l ¼ 20) containing 800 water molecules and 40 Nafion side chains. For around 3 times larger system size Cui et al. [26] obtained pore morphologies and found an increase of water and proton diffusion through Dow (EW ¼ 977) and Nafion (EW ¼ 1143) membranes with water uptake, in accordance with experiment. Karo et al. [27] simulated for a system size of w400 nm3 water diffusion within Dow and Nafion membranes at l ¼ 15. They noted that the pore morphology depends strongly on system size, and therefore questioned results deduced from earlier MD simulations [28] obtained for smaller simulation box volumes. The largest system size, comprising w2  106 atoms (volume w27,000 nm3), for which MD modeling has been performed on PFSA membranes was reported by Knox and Voth [29]. They pre-assumed 5 different morphologies and a random morphology. With exception of the random morphology, each of them evolved such that the characteristic scattering peak could be reproduced. Perhaps the most interesting observation was the formation of wires containing hydronium ions, water and sulfonate side chain groups. These wires connect the hydrophilic clusters that were originally imposed on the pre-assumed morphologies. Since water and proton diffusion takes place through the thus established hydrophilic network the extent to which the water clusters are interconnected by such wires will affect the long range transport properties. Park et al. [30] performed MD for Nafion1143 and several SPI block ionomers which differ in the length of the hydrophilic (sulfonated) blocks. Their work revealed that near a water content of 25 vol% (fw ¼ 0.25) for Nafion a well phase separated morphology occurs. The distribution of water within SPI membranes turned out

5025

to be much more dispersed. This was attributed to the stiff backbone and the possibility of water molecules to form also hydrogen bonds with the ether oxygen in the hydrophobic blocks. Since such an association of water with the Teflon backbone in Nafion is absent this resulted in a much more phase separated morphology within Nafion. From experiments of Park et al. [30] it was found that above a hydrophilic block length fraction of f ¼ 0.6 a sharp increase in water uptake (up to 45 vol% at f ¼ 0.8) and proton conductivity occurs, suggesting that water becomes better connected. Also their MD study predicted for f > 0.6 a sharp increase in proton diffusivity at water contents comparable with experiment. Jang et al. [31] applied coarse grained MD (CGMD) to study the effect of side chain distribution along the polymer backbones on morphology and diffusion. They found that at fw ¼ 0.2 for Nafion polymers with uniform side chain attachments phase separation is less pronounced and the pores are smaller than for polymers with all side chains located at one end of the backbone, but water and proton diffusion was not drastically affected by these differences in sequence design. Allahyarov et al. [32] concluded from CGMD simulations that for PFSA like polymers with similar EW those with the longer side chains form larger clusters and higher proton conductivities, in line with DPD studies [35,38] in which side chain length and density were varied. Dorenbos et al. [33e38] combined DPD with Monte Carlo (MC) trajectory calculations and simulated the diffusion of neutral species (water O2, N2, H2) through hydrated Nafion. The static pore networks obtained from DPD (using the parameterizations proposed by Yamamoto and Hyodo [39]) were mapped on a cubic grid on which MC trajectory calculations through the pore networks were performed. Using this strategy experimental long range water diffusion [34] and O2 and H2 gas permeation rates [35] could be approached. The interesting finding was the capability of predicting trends regarding diffusive transport when the EW or side chain length of the polymers is varied. Water diffusion and O2 and N2 permeation within hydrated Nafion membranes was expected at similar fw to decrease with increase in EW or decrease of ion exchange capacity (IEC), where IEC is the inverse of EW or mole equivalent of acidic sites per gram polymer. The explanation for these findings was that while the pore size increases with decrease of IEC, the connections between them become narrower, which act as bottlenecks for diffusion [33e35]. The same findings were deduced from theoretical studies [36,37] of model polymers similar in architecture as PFSA but in which the side chain lengths were varied. For architectures of the same IEC the ones that contain the longer side chains revealed at fixed fw the largest water clusters, largest inter cluster spacing and fastest diffusion. For architectures with same side chain lengths a decrease in IEC, due to increase in inter branching distance along the backbones, also resulted in increased cluster size and distance between them, but slower diffusion (see Fig. 8 in Ref. [37]). DPD-MC modeling aimed to predict percolation thresholds for diffusion, revealed that, at similar IEC, for the longer side chain architectures also much lower percolation thresholds are expected than those containing short side chains [37]. This was explained due to differences in topological distance between hydrophilic sites within the architectures: An increase in topological distance allows the formation of larger and better connected pores that result in higher diffusion constants at same water content and thus a lower percolation threshold. All of the above studies assumed idealized architectures with branching points (grafted polymers) or hydrophilic sites uniformly distributed along the backbone (block polymers). For realistic polymers for which these distributions might be non-uniform or statistical, this might affect the pore morphology. Indeed, a recent DPD-MC study [38] predicted that for grafted model ionomers,

5026

G. Dorenbos / Polymer 54 (2013) 5024e5034

similar in architecture as PFSA, a less uniform side chain distribution results in a better connected pore network and increased diffusion. Interestingly, also the agreement between experimental and calculated O2 permeability and water diffusion constants through hydrated Nafion was better for architectures for which side chains are statistically distributed as compared with Nafion for which the side chains are uniformly distributed. The purpose of the current study is to predict how diffusion through the pore networks of membranes composed of flexible multi-block copolymers might be affected by the distribution of the hydrophilic fragments along the polymer chain. To keep this study systematic, consecutive hydrophobic block lengths are chosen to be alternating long and short separated by a hydrophilic fragment (i.e. bimodal hydrophobic block length distribution, see Fig. 1(b)). DPD is used to generate the pore morphologies for chain architectures that differ in hydrophilic site density (IEC) and for each IEC all possible locations of these sites along the polymer backbones are considered. The morphologies are generated at a fixed hydration level in which 25 vol% of the system is occupied by water. As in previously published work [33e38], the diffusion through the generated pore networks is modeled separately by employing Monte Carlo (MC) tracer calculations. In Section 2 the DPD parameterization of the multi block copolymers is outlined. In Section 3 the obtained pore morphologies are analyzed and clear relations between morphology, chain architecture and obtained diffusion constants are presented. Diffusion constants are found to increase significantly when the difference in length between the consecutive hydrophobic blocks increases. From cluster size distributions it is deduced that the increase in diffusion is not due to more water being located within the percolating clusters, but rather should be explained by the presence of narrow connections (or bottlenecks) within the connected (percolating) pore network. In Section 4 the results are discussed and several items regarding convergence are addressed. 2. Computational details 2.1. Polymer chain architecture In DPD molecular fragments or several atoms are coarse grained by grouping them together into beads. A single repeat unit of the block co-polymers is shown in Fig. 1(b). The polymers are composed of hydrophobic blocks that contain x and y consecutively connected hydrophobic A beads separated by a hydrophilic C bead. The number of repeat units is 5, each architecture contains thus 10 hydrophilic C beads. The values of x and y for all architectures are listed in Table 1 and can be split into 3 groups for which x þ y ¼ constant (¼10, 12 or 14) and for which the C bead fractions and thus IECs are the same. The fractional hydrophilic bead density is thus 61 for the x þ y ¼ 10 and 81 for the x þ y ¼ 14 architectures. Within each group all possible combinations of x and y for which y
x are considered which results in a total of 18 architectures. The synthesis of these kinds of alternating hydrophobic block length distribution might be a challenge for experimentalists. In the present computer experiments these were selected because they are (1) among the simplest architectures with varying block lengths and (2) it could be verified whether similar trends are obtained for various IEC membranes. Statistical distributions of block lengths are not considered here. Water is modeled by W beads. The volume of each (A, C and W) bead is assumed to be 0.12 nm3 and corresponds to the volume occupied by 4 water molecules. The (dry) acidic site densities are thus 2.31 mmol/cm3 for the x þ y ¼ 10 and 1.73 mmol/cm3 for the x þ y ¼ 14 architectures. In principle the IEC expressed in mmol/

Table 1 Inter cluster distance DCleCl (4th column), Monte Carlo diffusion constants DMC at fw ¼ 0.25 (5th column), and amount of W beads contained within largest cluster (6th column). (x,y)

xþy

yx

DCleCl (nm)

DMC

Ncl

(1,9) (2,8) (3,7) (4,6) (5,5) (1,11) (2,10) (3,9) (4,8) (5,7) (6,6) (1,13) (2,12) (3,11) (4,10) (5,9) (6.8) (7,7)

10 10 10 10 10 12 12 12 12 12 12 14 14 14 14 14 14 14

8 6 4 2 0 10 8 6 4 2 0 12 10 8 6 4 2 0

5.04 4.75 4.61 4.46 4.46 5.69 5.54 5.33 5.26 5.11 5.04 6.19 5.98 6.05 5.90 5.62 5.54 5.54

0.181 0.118 0.082 0.047 0.044 0.150 0.103 0.031 0.019 0.021 0.009 0.106 0.095 0.02 0.006 0.004 0.002 0.0

31,907 31,746 31,533 31,073 31,191 31,907 30,969 29,454 28,900 30,650 29,504 31,887 31,828 30,849 21,069 14,641 12,390 5639

g can be obtained when the mass densities of the molecular fragments are known: For a polymer with a mass density of 1 g/cm3 (dry weight) the IEC for a x þ y ¼ 10 architecture becomes 2.31 mmol/g, which corresponds with an EW of 433 g/mol. For the purpose of keeping this work general the chemical formulas of the A and C beads are not explicitly given. 2.2. Dissipative particle dynamics Two decades ago Hoogerbrugge and Koelman [40] introduced DPD in order to study hydrodynamics. Since then DPD has been refined and is used in areas such as phase separation [41], membrane rupture [42], vesicles [43,44], etc. Here the mathematical framework described by Groot and Warren is followed [45] as outlined below. The beads evolve according to Newton’s equations of motion: dri/dt ¼ vi, midvi/dt ¼ fi, with mi, ri and vi, being the mass, position and velocity of bead i, respectively. The interaction between non-connected beads i and j is given by the sum of a R conservative force FCij, a dissipative force FD ij , and random force Fij. The total force fi acting on bead i is given by Eq. (1).

fi ¼

X

FijC þ FijR þ FijD

(1)

jsi

The sum is over all particles j located within a cutoff distance rc. FCij is soft repulsive and decreases linearly with the distance between beads and becomes zero beyond rc. This distance also defines the length scale:

( FijC

¼

  r r ij aij 1  rijc b

rij ¼ ri  rj ;

0

 

rij < rc rij
rc

  rij ¼ rij ;

 

  b r ij ¼ rij =rij 

(2a)

(2b)

In Eq. (2) b r ij is the unit vector in the direction of bead j toward bead i: FRij corresponds to thermal noise (Eq. (3)), and FD ij is proportional to the beads’ relative velocity (Eq. (4)).

  r ij FijR ¼ suR rij zij ðDtÞ0:5 ðkB TÞ1 b    D D Fij ¼ gu rij b r ij $vij b r ij

(3) (4)

G. Dorenbos / Polymer 54 (2013) 5024e5034

vij ¼ vi  vj, uD and uR are weight functions, zij introduces randomness into the system involving a randomly fluctuating variable with Gaussian statistics with zero mean and unit variance: ¼ 0, and ¼ (dikdjl þ dildjk)d(t  t0 ). For each pair of beads there is an independent random function. uD and uR are given by Eq. (5).

8 2 < 1  rij h  i2   r c uD rij ¼ uR rij ¼ : 0

  rij < rc   rij
rc

(5)

Noise (s) and friction (g) parameters are related as s2 ¼ 2gkBT [46], with s ¼ 3, g ¼ 4.5, kB the Boltzmann constant and T the temperature of the system. All three forces act along the line of centers conserving linear and angular momentum. Adjacent beads that belong to the same polymer are joined by a spring force:

FijS ¼ Crij

(6)

The constant C is 4. As for the interaction range rc, the mass of all beads and the unit of energy (kBT) are scaled to 1, as a result also the unit of time s ¼ rc(m/kBT)0.5 is equal to 1. A modified Verlet integration algorithm is used [45] to solve the equations of motion with empirical factor 0.65 and time step Dt ¼ 0.05. In this way the temperature is kept constant near kBT ¼ 1.0. The bead density r is 3. The adapted repulsions aij between beads are given in Table 2. Hydrophilic interactions are those for which the repulsions are lowest (e.g. CeW), while repulsions between incompatible (hydrophobic vs hydrophilic) beads are highest (AeW and AeC). The repulsions between similar beads are set at 104 which reproduces the water compressibility for a system composed solely of water beads [38,39]. The non-diagonal repulsions are calculated from the assumed FloryeHuggins c parameters. These were chosen to be comparable to the c-values that were calculated for the molecular fragments of Nafion by Yamamoto and Hyodo [39]. In their DPD parameterization a repeat unit of Nafion was modeled as A0 4[B0 C0 ] (Fig. 1(a)) with assumed bead volumes of 0.12 nm3 (similar as in this work) and hydrophobic backbone A0 beads composed of (CF2eCF2e CF2eCF2) and side chain B0 beads (OeCF2eCF(CF3)eO) and hydrophilic (side chain) C0 beads composed of CF2eCF2eSO3H fragments. They calculated from atomistic simulations of the binding energies the c parameters between dissimilar beads, and obtained cA0W ¼ 5.79, cB0 W ¼ 4.9, cC0 W ¼ 2.79, cA0 B0 ¼ 0.02, cA0 C0 ¼ 3.11, and cB0 C0 ¼ 1.37. In order to reduce in this work the parameter space, each model block-copolymer is composed of only one type of hydrophobic and one type of hydrophilic bead. Furthermore, the c parameters between hydrophobic (A) beads and hydrophilic (C and W) beads are chosen to be the same with cAC ¼ cAW ¼ 4.9. These values are in between the values cA0W ¼ 5.79 and cA0 C0 ¼ 3.11 and equal to the cB0 W ¼ 4.9 estimate [39], and were also adapted in Ref. [37]. Also the selected cCW ¼ 2.6 is close to the estimated cC0 W ¼ 2.79 in Ref. [39]. The motivation to select these c-values was that for Nafion well phase separated morphologies could be obtained [34,35,39]. Therefore it can be anticipated upon that the

5027

c-values in Table 2 will also result in well phase separated morphologies for the block polymer architectures in Fig. 1(b). In order to predict general trends and for the sake of generality, atomistic calculations to obtain c parameters for given chemical formulas of the beads were omitted. Here, the most important assumption is that A and C beads are incompatible with each other, with A fragments unfavorably interacting with water, while C fragments have a favorable interaction with water. Another motivation to adapt the c parameters listed in Table 2 is that they are the same as in systematical studies on phase separation and diffusion within model block and grafted polymers [36,37], thus expanding a database that might guide fuel cell membrane development. Fig. 1(b) might suggest that A and C beads are the same as in Nafion, resulting in a block polymer with all carbon atoms covalently bonded along the backbone. But these beads might as well contain ether linkages, (fluoro) aromatic and/or hydro or fluorocarbon molecular fragments, different acid groups, etc. Certainly, for all these synthesized polymers the exact c-values will deviate to various extent from those listed in Table 2. The DPD repulsions that correspond with the c-values were calculated from Eq. (7) and are given in Table 2.

Daij ¼ ð4:16  0:15Þ  cij

In Eq. (7) Daij ¼ aij  aii is the excess repulsion between beads i and j. Eq. (7) was obtained by using a similar procedure as described by Groot et al. [42,45] by calculating the mutual solubility of N-mers (N ¼ 2, 3, 4) composed of A beads and N-mers composed of C beads for various repulsions aAC > 104. Using a simulation box of 30  8  8 the left hand side was initially filled with 5760/2N Atype N-mers, and the right hand side was filled with the same amount of C-type N-mers. With the diagonal repulsions aAA and aCC set at 104, the minority A (C) bead concentration rA (rC) within the majority C (A) phase was averaged over the time window 1.3  105Dt  1.6  105Dt for several repulsions aAC > 104. For each set of repulsions the effective c parameters were calculated according to Eq. (8) [45,47].

cN ¼

  ln rA=r C

(8)

1  2rA

In Fig. 2 c is plotted against the excess repulsion DaAC. Since bead volumes are the same, the water volume fractions fw are given by Eq. (9) (NA, NW and NC are the total number of A, W and C beads).

Table 2 DPD repulsions aij between beads. The corresponding FloryeHuggins cij-parameters are given in parenthesis.

A C W

A

C

W

104 (0) 124.4 (4.9) 124.4 (4.9)

104 (0) 93.2 (2.6)

104 (0)

(7)

Fig. 2. c plotted against DaAC.

5028

G. Dorenbos / Polymer 54 (2013) 5024e5034

fw ¼

NW

NW þ NA þ NC

(9)

In this work fw ¼ 0.25. Since each W bead represents 4 water molecules, the water content expressed in terms of l (number of water molecules/C beads) is given by Eq. (10).

l ¼

NW NC

4

(10)

For the architectures in Fig. 1(b) the relation between fw and l is then given by Eq. (11).

l ¼

2  fw  ðx þ y þ 2Þ ð1  fw Þ

(11)

For fw ¼ 0.25 the corresponding l values are l ¼ 8 (x þ y ¼ 10), l ¼ 9.33 (x þ y ¼ 12) and l ¼ 10.67 (x þ y ¼ 14). The system size is cubic with cell dimension of L ¼ 35. One DPD length unit corresponds to rc ¼ (3  0.12)1/3 nm3 ¼ 0.71 nm, therefore box edge lengths are L  rc w25 nm and system volumes are w15,000 nm3. The DPD simulations involve the dynamics of 128,625 beads with 1608 (x þ y ¼ 10), 1379 (x þ y ¼ 12) or 1206 (x þ y ¼ 14) polymer chains. Periodic boundaries are applied in each orthogonal direction. Morphologies generated at 2  104Dt (or t ¼ 1000s) were used as input for the MC calculations presented in Section 3. An estimate of the physical time to which s corresponds can be obtained by matching the DPD diffusion constant of a system solely composed of W beads with the experimental diffusion constant of pure water [42]. For aww ¼ 104 s corresponds to w130 ps. The generated morphologies were thus obtained at a corresponding physical time of O(0.1 ms). It was verified that the presented results are not drastically affected when performing longer DPD runs up to 4  104Dt. 3. Results and analysis

Figs. 3 and 12 additional sequences listed in Table 1, are plotted against y  x in Fig. 5. For membranes of similar IEC (x þ y ¼ constant), DCleCl increases linearly with y  x. The DCleCl values for all 18 chain architectures are listed in Table 1. 3.2. Monte Carlo tracer diffusion The pore connectivity can be probed by MC trajectory calculations by restricting the tracer particle movement to the pore phase [34]. By following the random movement of the particles, the tracer diffusion coefficients can be determined. The particle movement takes place on a cubic grid of size 3503 thus containing almost 43 million nodes. On this grid the W pore network is constructed in such a way that nodes for which the nearest bead is a W bead belong to the pore phase. After construction of the pore network 2000 particles (Ntracer ¼ 2000) are initially put at randomly selected nodes that belong to the pore phase. At each Monte Carlo step (MCS), a jump trial toward a nearest node in one of the six orthogonal directions is randomly selected for each particle independently. A trial is successful when the aimed site belongs to the W (or pore) phase. Diffusion constants are usually derived from the mean square displacement (MSD) of Ntracer trajectories (Eq. (12)):

D ¼

1 6Ntracer

lim

t/N

N D E d X 2 ½Ri ðtÞ  Ri ðt 0 Þ dt i ¼ 1

(12)

MSD curves obtained for the lhs morphologies in Fig. 3 and additional sequences (2,8) and (4,6) are displayed in Fig. 6(a). In Fig. 6(a), the pure water case, for which every jump trial is accepted, is also included. For pure water the slopes are per definition equal to 1. The diffusion constants within the membranes, DMC, can therefore be expressed relative to the diffusion in pure water value according to Eq. (13) [34,37].

DðMSDÞ DðMCSÞ

3.1. Pore morphology

DMC ¼

In Fig. 3 morphologies obtained at fw ¼ 0.25 for x þ y ¼ 10 ((5,5), (3,7), (1,9)) and x þ y ¼ 14 ((7,7), (4,10), (1,13)) architectures are displayed. For each morphology the water beads are located within the polymer A matrix while C beads are located near the boundaries of the water clusters (W/A interface). All pore structures appear to be the same, in the sense that they are sponge-like. But it seems that the pore size and (since fw ¼ constant) the distances between the pores are largest for the architectures for which the y values are highest (i.e. the bottom (1,9) and (1,13) morphologies in Fig. 3). This suggests that an increase in difference (y  x) between consecutive hydrophobic block lengths along the polymer backbone favors the formation of larger pores. Also when a comparison is made between the rhs x þ y ¼ 14 morphologies and the lhs x þ y ¼ 10 morphologies, the pores within the high IEC x þ y ¼ 10 membranes are smallest. An indication of the average size of the pores can be deduced from the W-bead pair correlation functions, G(r). Those calculated for the morphologies in Fig. 3 are shown in Fig. 4. The distance for which G(r) drops below 1, is indeed smallest for the x þ y ¼ 10 architectures, and ranges between w1.9 nm (for (5,5) architecture) and w2.1 nm ((1,9) architecture). For the x þ y ¼ 14 architectures these values range between w2.4 nm ((7,7) architecture) and w2.7 nm ((1,13) architecture). An estimate of the distance between water clusters DCleCl, is derived from the location of the second maximum in G(r), located near 5 nm (see inset Fig. 4). DCleCl increases with decrease in IEC (or increase of x þ y). DCleCl obtained for the morphologies shown in

The D(MSD)/D(MCS) slopes were determined by linear fitting over the time window 4  106e6  106 MCS [37]. From Fig. 6(a) it is observed that diffusion through the pore networks is fastest for the (1,9) architecture for which the hydrophilic C beads are highly non-uniformly distributed along the backbones. When the difference in consecutive hydrophobic block lengths is decreased, diffusion decreases also systematically in the order (1,9) > (2,8) > (3,7) > (4,6) > (5,5). From the deduced DMC values listed in Table 2 a relative increase of a factor 4.1 in diffusion is observed for the (1,9) architectures as compared with the (5,5) architecture. This trend is getting more pronounced when the IEC is decreased (or EW is increased): For the x þ y ¼ 14 series (Fig. 6(b)) long range diffusion is absent for the (7,7) architecture but increases gradually when y is increased. In Fig. 7 the obtained DMC values are plotted against y  x. Two main trends are observed. First all architectures for which y  x is the same diffusion decreases with increase of y þ x. This means that a decrease of IEC (increase of EW) results in slower diffusion. Secondly, it is found that for membranes of the same IEC (x þ y ¼ constant) an increase in difference in length between the consecutive hydrophobic blocks (given by y  x) enhances diffusion. In Fig. 8 DMC is plotted against DCleCl. For sequences of similar IEC (x þ y ¼ constant) DMC increases with inter cluster distance (DCleCl). The reversed trend, in which increased diffusion correlates with a decrease in DCleCl, occurs for the uniform (y ¼ x) and most non-uniform ((1,9), (1,11), (1,13)) architectures.

(13)

G. Dorenbos / Polymer 54 (2013) 5024e5034

5029

Fig. 3. Pore morphologies at fw ¼ 0.25 obtained for (AxCAyC)5 architectures with x þ y ¼ 10 (lhs figures) and x þ y ¼ 14 (rhs figures). A beads are colored red (gray), C beads yellow (light gray) and W beads blue (dark gray). Iso-surfaces of the W bead densities are drawn at an iso-value of 0.8.

For membranes of the same IEC an increase in (y  x) results in higher DCleCl values (Fig. 5) and a significant increase in diffusion (Fig. 7). Apparently, for membranes of the same IEC an increase of the difference (y  x) between the consecutive hydrophobic block lengths allows for (1) significantly more water to participate within a connected network and/or (2) less severe bottlenecks for diffusion through the percolating network. However, as also was deduced from Fig. 8 an increase of DCleCl does not necessarily result in increased diffusion: When comparing membranes whose IEC is different such as the sets ((1,9), (1,11), and (1,13)) or ((5,5), (6,6), and (7,7)) diffusion systematically decreases with DCleCl. Since this trend is reversed to the one obtained for membranes of the same IEC, this implies that larger DCleCl, and thus larger pores (since fw ¼ constant), not always result in better connected pore

networks that favor fast diffusion. This might be due to less water participating in a percolating network or that within the existing percolating network there are more dead ends, or bottlenecks that slow down diffusion. 3.3. Cluster size distribution The distribution of water beads within the DPD morphologies was investigated by calculating the cluster size distributions (CSD) of the W beads. Due to the repulsive interactions between DPD beads (Eq. (2)) W beads are not randomly distributed within the water phase. For very small Lbond values W beads can thus not be contained within large clusters and each W bead will be isolated. But for Lbond values near or well above the first peak in the pair

5030

G. Dorenbos / Polymer 54 (2013) 5024e5034

Fig. 4. Pair correlation functions of the W beads for the morphologies shown in Fig. 3. (5,5) (Crosses); (3,7) (filled triangles); (1,9) (filled circles); (7,7) (); (4,10) (open triangles); (1,13) (open circles).

correlation function (Fig. 4) water clustering can be expected. Here, a water cluster is defined to be of size Ncl when it contains Ncl W beads and the distance of each bead within that cluster to at least one of the remaining Ncl  1 beads within that cluster is less than Lbond ¼ 1 DPD length unit. Within a DPD simulation this distance is equal to the interaction range of the DPD beads (see, Eq. (2)), slightly larger than the distance of w0.6 nm (w0.8 DPD length unit) for which the first peak in the W bead pair correlation function occurs. Note that more than one W bead can be within 1 DPD length unit of a certain W bead. Such W beads will automatically belong to the same cluster. Fig. 9(a) displays the CSD for the uniform (y ¼ x) and most non-uniform architectures. The total number of W beads within each morphology was 32,157 (¼fwrL3). For each architecture a large majority of the clusters contain less than 10 W beads (varying between w86% for the (7,7) architecture and 99.5% for the (1,13), (1,11) and (1,9) architectures), but the accumulated number of W beads located within these clusters turns out to be at most 1.2% of the total amount of W beads. The number of W beads contained within the largest cluster for all 18 morphologies are listed in Table 1. For most of the architectures the values are close to the total number of water beads

contained within the DPD simulation box (32,157). Notable exceptions are the (4,10), (5,9), (6,8) and (7,7) architectures. For these four architectures also the lowest MC tracer diffusion constants are obtained. Fig. 9(b) displays the fraction (fisol) of W beads that is not contained within the largest cluster. From Fig. 9(b) it is observed that for architectures for which the differences in block lengths (y  x) are the same fisol increases with increase of y þ x. With exception of the 4 sequences just mentioned, more than 90% of the W beads are contained within the largest (percolating) cluster. This means that, besides the exceptions mentioned above, the differences in diffusion constants as observed in Figs. 7 and 8 are not due to a large amount of water being located within isolated clusters, but rather should be explained by the presence of narrow connections (bottleneck) within a single percolating cluster. The fraction of isolated water clusters can in principle also be estimated from long MC trajectory calculations by inspecting the maximum displacement of each tracer particle during its trajectory. Then the fraction of tracer trajectories not capable of performing maximum displacements larger than the size of the simulation box will converge to the isolated fraction. In a recent study it was confirmed that the isolated fraction derived from long MC calculations approach those derived from the DPD W bead cluster size distributions obtained for Lbond near 1 DPD length [38].

4. Discussion The MC tracer diffusion constants, DMC, are found to depend strongly on the size of the hydrophobic block lengths x and y. These results could be obtained because the selected coarse graining strategy allowed a minimal computational effort with respect to (1) evolvement toward equilibrium-like (DPD) morphologies and (2) the calculation of diffusion through the mapped morphologies. This combination makes the DPD-MC strategy very efficient. In principle diffusion coefficients can also be calculated from the , of the displacements of W beads ensemble average, DMembrane Wbead during the DPD simulations. The molecular water diffusion coeffican then be obtained by cient through the membrane DMembrane H2 O comparing DMembrane with the W bead diffusion coefficient Wbead DMembrane for a system solely composed of W beads according to Eq. Wbead (14). In Eq. (14) DPurewater is the experimental pure water diffusion H2 O coefficient of 2.3  105 cm2/s.

DMembrane ¼ H2 O

Fig. 5. (a) DCleCl at fw ¼ 0.25 versus y  x. Lines are linear fits through data for which x þ y are the same.

DMembrane Wbead water DPure Wbead

water  DPure H2 O

(14)

Instead of deriving diffusion coefficients from the DPD W bead movements the MC approach was applied in Section 3 for two reasons. First, DPD allows some extent of bond crossing of polymer beads that can lift up entanglements. This causes DPD polymers to diffuse faster than real polymers, resulting in fast evolvement to equilibrium like morphologies. But consequently the DPD pore morphology might also alter while the water beads are diffusing through it. For this reason, frozen (static) DPD morphologies were taken as input for the MC trajectory calculations and the pore networks were thus fixed in time. A second reason is that even when during a DPD calculation the pore networks would not change at all, or less (eg. by forbidding of bond-crossing or adapting high masses to polymer beads), any estimation of long range water diffusion sampled from the movements of W beads would require very long simulation runs. In order to estimate (long-range) diffusion constants, the whole system should be probed. This means that the W beads must be capable of traveling distances longer than the box edge length. For

G. Dorenbos / Polymer 54 (2013) 5024e5034

5031

Fig. 6. MSD curves obtained for fw ¼ 0.25. (a) All 5 MSD curves obtained for architectures for which x þ y ¼ 10. The pure water case (fw ¼ 1.0), for which the slopes are equal to 1, is also included (filled circles). (b) All 7 MSD curves obtained for architectures for which x þ y ¼ 14.

the present system size the edge lengths are L ¼ 35 DPD length units. The diffusion of W beads for a system solely composed of W beads (ie. fw ¼ 1.0) for r ¼ 3 and aww ¼ 104 was calculated according to Eq. (12) to be D ¼ 0.13. The average time needed for these beads to travel through the whole system is then 352/6D w 1600s, which requires a simulation of duration of 32,000 time steps. However, the time for W beads to probe the pore networks at fw ¼ 0.25 is expected to be much longer: The diffusion constants (Table 1) through the static pore networks are up to O(102) less than that of pure water. This requires an exceptional long DPD run of a duration of O(106) time steps for the hypothetical case that the pore morphology does not rearrange during the DPD simulations. But since the polymers can rearrange continuously, it is expected that , are higher than the MC diffusivities W bead diffusivities, DMembrane Wbead DMC listed in Table 1. The above was verified by restarting all 18 DPD calculations for another 20,000 time steps starting from the morphologies that were generated at 20,000Dt. As mentioned, during this period W beads can only probe a portion of the simulation cell. In Fig. 10 obtained W bead diffusivities for the x þ y ¼ 10 (Fig. 10(a)) and x þ y ¼ 14 (Fig. 10(b)) architectures, are compared with the MC diffusivities. Interestingly, both MC diffusion constants and W bead diffusion constants exhibit the same trends: (1) diffusion constants within x þ y ¼ 10 membranes are faster than within the x þ y ¼ 14 membranes (for similar y  x); (2) diffusion increases with y  x. Indeed as expected, the diffusivities of the DPD W beads are larger than the long range diffusivities obtained by MC. Besides the

Fig. 7. DMC plotted against y  x. fw ¼ 0.25.

argumentation in the previous paragraph, this is also partly caused by the fact that during the simulations the probed length scale is of pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the order 6D1000 DPD length units (rc). From Fig. 10 water ranges between 0.13 and 0.3 and the abso=DPure DMembrane Wbead Wbead

lute W bead diffusivity ranges between DMembrane ¼ 0:017 and Wbead DMembrane ¼ 0:045, the average distances traveled by the W beads Wbead are thus w10rc and w16rc, respectively. Therefore at most half of the simulation box size (L ¼ 35rc) could be probed during the period of 20,000Dt. The DMembrane/DPure water values derived from the W bead motions in Fig. 10 should therefore be considered as upper bounds. The mass of each bead was the same. When supposing that the parameterization of an A bead would be similar to the A0 beads in Nafion [34,35,39], its weight would correspond to 200 atomic mass units (amu) (mass of 4 CF2 fragments), while the mass of each W bead (72 amu, corresponding to 4H2O) is nearly 3 times less. While this will not affect the eventual pore morphology this may lead to further decrease of the W bead diffusivities and approach the MC values. This was not verified in the present study. Within the presented figures each datum (DCleCl, DMC, cluster size) was obtained for only one pore morphology rather than from ensemble averages obtained from various morphologies. Since very persistent trends were found the conclusions are not expected to be affected when DCleCl, DMC, etc., are calculated by performing a sampling over the results obtained from several separate DPD runs. As an illustration, for the (5,5), (2,8) and (1,9) architectures DMC (Fig. 11(a), DCleCl Fig. 11(b)), and the water bead fraction (fisol) contained within isolated clusters (Fig. 11(c)) are plotted against

Fig. 8. Tracer diffusion coefficients DMC obtained for the architectures listed in Table 1 plotted against DCleCl. fw ¼ 0.25.

5032

G. Dorenbos / Polymer 54 (2013) 5024e5034

Fig. 9. (a) Cluster size distributions obtained for uniform (y ¼ x) (filled symbols) and non-uniform (open symbols) architectures. (b) Fraction of isolated W beads that are not contained within the largest W cluster plotted against y  x.

DPD time. Beyond t ¼ 1000s, DMC and DCleCl do not vary significantly and also fisol for these cases is always below 0.03, meaning that more than 97% of W beads is always located within one W cluster. For all times (500, 1000, 1500, 2000) DCleCl and DMC increases while fisol decreases when the hydrophobic block distribution becomes less uniform. The selected value of 2000 for the number of tracer particles, Ntracer, in obtaining the DMC is sufficient to reach convergence. This is illustrated in Fig. 11(d), where the MSD at t ¼ 6  106 MCS is plotted against Ntracer. The polymers were flexible with no bending constraints that could impair the process that governs the formation of the micro phase separated morphologies. These are (1) the maximization of the number of hydrophilicehydrophilic (WeW and CeW) and hydrophobicehydrophobic (AeA) contacts, and (2) minimization of hydrophobicehydrophilic (WeA) interface. As illustrated in Fig. 12, when all bonds in the architectures would be removed but with the A, C, and W bead fractions kept the same, these processes result in a morphology with a mixture of C and W beads fully separated from the A bead phase. The present work allowed a systematic study of the impact of the constraints of bead connectivity and the location of the hydrophilic C beads within the architectures on the phase separation process. For the present block polymers both processes result in different pore networks for membranes of the same IEC. When the hydrophobic block length distributions are highly non-uniform (asymmetric), the diffusion through the pore networks is more facile than when the blocks are uniformly distributed. A plausible explanation for this is that C beads which are separated by the longer hydrophobic blocks within the asymmetric architectures are more likely to join a different pore than those C beads that are separated by the shorter hydrophobic block. This can be most easily understood for a (1,13) architecture: C beads separated by a short hydrophilic block are separated by only 1 A bead (see Fig. 1(b)) and will join the same

pore. C beads that are separated by 13 A beads are expected to have a much higher probability to join a different pore. With increase of asymmetry (going from (7,7) to (1,13)) this leads to an increased distance between the pores as observed in Fig. 5. Since fw was fixed at 0.25 this resulted thus also in larger pores when asymmetry (y  x) is increased. For each IEC membrane this increase in DCleCl and pore size due to higher asymmetry results in easier diffusion (Figs. 5, 7 and 8). Diffusion can be increased even further when the IEC is increased as observed for the (1,13), (1,11), and (1,9) architectures (Fig. 8) despite the fact that DCleCl decreases with increase of IEC. Clearly, an increase in C bead fraction (or decrease in hydrophobic A bead fraction) also plays an essential role in the pore network formation. For fuel cell applications sufficient membrane rigidity is needed. This can be obtained if a small crystalline fraction is present. It is well possible that for similar IEC membranes, the crystalline fraction depends on the length of the blocks. Since for the non-uniform (y > x) architectures the distance between water clusters tend to significantly increase (see Fig. 5), this might also favor the formation of micro crystals within the polymer phase. Especially, architectures for which the hydrophilic fragments are statistically distributed along the backbones (not studied here), will contain relatively long hydrophobic fragments that might group together and crystallize in the hydrophobic phase regions. These issues were beyond the scope of this work. Since crystallization processes cannot be taken into account with DPD these should be challenged by other (CG) molecular dynamics methods. This study predicts systematic trends for hypothetical model multi block co-polymer membranes. The impact of amount of incompatibility/compatibility between beads was not addressed here, but a systematic study of the dependency of the diffusive property on the c-values might be worthwhile, thus opening the

Fig. 10. (a and b) Monte Carlo tracer diffusion DMC, (diamonds) within the static DPD morphologies sampled at 20,000Dt compared with (short range) DPD W bead diffusion constants (triangles) sampled over time window 20,000Dte40,000Dt for (a) x þ y ¼ 10 and (b) x þ y ¼ 14 architectures.

G. Dorenbos / Polymer 54 (2013) 5024e5034

5033

Fig. 11. Examples of obtained time dependence of DMC (a), DCleCl (b), the fraction of isolated W beads (c) for the (5,5) (diamonds), (2.8) (squares) and (1,9) (triangles) architectures. (d) MSD at t ¼ 6  106 MCS plotted against number of tracer particles for architectures for which x þ y ¼ 10 within the DPD morphologies generated at 20,000Dt.

way to verify if there are c-values that for fixed polymer architectures form the best connected pore network. Molecular fragments that correspond to such optimal c-values could then be chosen as building blocks in the synthesis of new membranes. In order to increase the confidence in such an approach it would be helpful to know whether the trends reported in this work can also be obtained from full atomistic MD simulations, or from experimental data on synthesized polymers in future.

5. Conclusions Phase separated morphologies of hydrated multi-block copolymer membranes were obtained by DPD while diffusion through the pore networks was probed by MC tracer diffusion calculations. The polymers are composed of alternating short (x consecutively connected hydrophobic beads) and long (y consecutively connected hydrophobic beads) hydrophobic blocks that are separated by (or linked together with) a hydrophilic bead. For fixed water volume fraction fw ¼ 0.25 the majority of the water containing pores are for most architectures contained within a percolating cluster that allows for long range diffusion. For architecture of similar ion exchange capacity the inter cluster distances increase approximately linearly with difference in hydrophobic block lengths (y  x). Even more interesting is that an increase in (y  x) also results in increased long range diffusion through the pore networks. Since most water is located within one large (percolating) cluster this means that the extents to which bottlenecks in the (connected) pore networks hinder diffusion can be decreased by increase of y  x. These results might be of relevance with respect to optimizing the pore structures in future block type PEM for fuel cell applications. Acknowledgment I gratefully thank Toyota Motor Corporation for stimulating discussions and providing assistance in publishing this work. References

Fig. 12. Distribution of A (red; dark gray), C (yellow; light gray), and W beads (blue; dark) at t ¼ 500s. System size is cubic with L ¼ 15. The fractional bead occupancies are similar as for the simulations involving the x þ y ¼ 10 architectures, ie. fW-bead ¼ 0.25, fA-bead ¼ 0.67, fC-bead ¼ 0.08.

[1] Saito M, Arimura N, Hayamizu K, Okada T. J Phys Chem B 2004;108:16064e70. [2] Zawodzinski Jr TA, Derouin C, Radzinski S, Sherman RJ, Smith VT, Springer TE, et al. J Electrochem Soc 1993;140:1041e7. [3] Hinatsu JT, Mizuhata M, Takenaka H. J Electrochem Soc 1994;141:1493e8. [4] Zhao Q, Carro N, Ho YR, Benziger J. Polymer 2012;53:1267e76. [5] Bass M, Freger V. Polymer 2008;49:497e506. [6] Gierke TD, Munn GE, Wilson FC. J Polym Sci Polym Phys 1981;19:1687e704.

5034 [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]

G. Dorenbos / Polymer 54 (2013) 5024e5034 Mauritz KA, Moore RB. Chem Rev 2004;104:4535e85. Schmidt-Rohr K, Chen Q. Nat Mater 2008;7:75e83. Edmonson CA, Fontanella JJ. Solid State Ionics 2002;152e153:355e61. Kreuer KD, Schuster M, Obliers B, Diat O, Traub U, Fuchs A, et al. J Power Sources 2008;178:499e509. Li Q, He R, Jensen JO, Bjerrum NJ. Fuel Cells 2004;4:147e59. Xing PX, Robertson GP, Guiver MD, Mikhailenko SD, Wang KP, Kaliaguine S. J Membr Sci 2004;229:95e106. Li X, Zhao C, Lu H, Wang Z, Na H. Polymer 2005;46:5820e7. Li N, Shin DW, Hwang DS, Lee JM, Guiver MD. Macromolecules 2010;43: 9810e20. Lee M, Park J, Lee H-S, Lane O, Moore RB, McGrath JE, et al. Polymer 2009;50: 6129e38. Miyatake K, Zhou HZ, Matsuo T, Uchida H, Watanabe M. Macromolecules 2004;37:4961e6. Atkins JR, Sides CR, Creager SE, Harris JL, Pennington WT, Thomas BH, et al. J New Mater Electrochem Syst 2003;6:9e15. Elliot JA, Hanna S, Elliot AMS, Cooley GE. Phys Chem Chem Phys 1999;1: 4855e63. Vishnyakov A, Neimark AV. J Phys Chem B 2001;105:9586e694. Jinnouchi R, Okazaki K. J Electrochem Soc 2003;150:E66e73. Seeliger D, Hartnig C, Spohr E. Electrochim Acta 2005;50:4234e40. Urata S, Irisawa J, Takada A, Shinoda W, Tsuzuki S, Mikami M. J Phys Chem B 2005;109:4269e78. Zhou X, Chen Z, Delgado F, Brenner D, Srivastava R. J Electrochem Soc 2007;154:B82e7. Cui S, Liu J, Selvan ME, Paddison SJ, Keffer DJ, Edwards BJ, et al. J Phys Chem B 2007;111:2208e18.

[25] Devanathan R, Venkatnathan A, Rousseau R, Dupuis M, Frigato T, Gu W, et al. J Phys Chem B 2010;114:13681e90. [26] Cui S, Liu J, Selvan ME, Paddison SJ, Keffer DJ, Edwards BJ. J Phys Chem B 2008;112:13273e84. [27] Karo J, Aabloo A, Thomas JO, Brandell D. J Phys Chem B 2010;114:6056e64. [28] Brandell D, Karo J, Liivat A, Thomas JO. J Mol Model 2007;13:1039e46. [29] Knox CK, Voth GA. J Phys Chem B 2010;114:3205e18. [30] Park CH, Lee CH, Sohn J-H, Park HB, Guiver MD, Lee JM. J Phys Chem B 2010;114:12036e45. [31] Jang SS, Molinero V, Cagin T, Goddard III WA. J Phys Chem B 2004;108: 3149e57. [32] Allahyarov E, Taylor P. J Polym Sci Part B Polym Phys 2011;49:368e76. [33] Dorenbos G, Suga Y, Kume H. JSAE Trans 2006;37:99e104. [34] Dorenbos G, Suga Y. J Membr Sci 2009;330:5e20. [35] Dorenbos G, Morohoshi K. J Chem Phys 2011;134:044133. [36] Dorenbos G, Morohoshi K. Energy Environ Sci 2010;3:1326e38. [37] Dorenbos G, Morohoshi K. J Mater Chem 2011;21:13503e15. [38] Dorenbos G, Morohoshi K. J Chem Phys 2013;138:064902. [39] Yamamoto S, Hyodo S-A. Polym J 2003;35:519e27. [40] Hoogerbrugge PJ, Koelman JMVA. Europhys Lett 1992;19:155e60. [41] Groot RD, Madden TJ. J Chem Phys 1998;108:8713e24. [42] Groot RD, Rabone KL. Biophys J 2001;81:725e36. [43] Yamamoto S, Murayama Y, Hyodo S. J Chem Phys 2002;116:5842. [44] Wang H, Liu Y-T, Qian H-J, Lu Z-Y. Polymer 2011;52:2094e101. [45] Groot RD, Warren PB. J Chem Phys 1997;107:4423e35. [46] Espanol P, Warren P. Europhys Lett 1995;230:191e6. [47] Strobl G. The physics of polymers: concepts for understanding their structure and behaviour. Heidelberg, Berlin, Germany: Springer-Verlag; 1997.