Using coarse-grained molecular dynamics simulations, we study ionomers in equilibrium and under uniaxial tensile deformation. The spacing of ions along the chain is varied, allowing us to consider how different ionic aggregate morphologies, from percolated to discrete aggregates, impact the mechanical properties. From the equilibrium simulations, we calculate the stress-stress auto correlation function, showing a distinct deviation from the Rouse relaxation due to ionic associations that depends on ion content. We then quantify the morphology during strain, particularly the degree to which both chains and ionic aggregates tend to align. We also track the location of the ionomer peak in the anisotropic structure factor during strain. The length scale of aggregate order increases in the axial direction and decreases in the transverse direction, in qualitative agreement with prior experimental results.

Ionomers are polymers that contain a relatively small fraction of covalently bound ionic groups, widely used in a range of applications such as fuel cell membranes, packaging, and the outer covering of golf balls.1,2 Of interest here are dry ionomers (without solvent) in the melt state and above the glass transition temperature, as is often relevant during processing of these materials. Because the backbone typically has a relatively low dielectric constant, if no solvent is present, the ions tend to aggregate locally. When ions of different chains participate in the same aggregate, the chains can be held together as though they are temporarily cross-linked, leading to their favorable mechanical properties.3,4 These ionic associations are in a state of dynamic equilibrium and there is an interchange of ions between ionic aggregates at temperatures above Tg.5,6 Because of this associative nature, their rheological response changes significantly as a function of temperature in the melt regime, from fluid-like viscosity to near solid-like elastic behavior.7,8 As these properties depend on controllable features of the ionomers such as the concentration of ionic groups, their demand has increased in fields where tunable and mechanically robust polymers are needed.2,4 An understanding of molecular phenomena that lead to observable rheological properties is central to designing materials optimally.

The strong ionic aggregate formation that distinguishes ionomers from other polymers has been studied extensively as it is a key factor dictating bulk ionomer properties. Much work has focused on measurement and analysis of the “ionomer peak” observed due to ionic aggregate structures. The peak is frequently fit using the Yarusso-Cooper9 or Kinning-Thomas10 model, which assumes roughly spherical aggregates with liquid-like ordering, to estimate the size and spacing of the aggregates.11–13 To better understand material properties and their relationship with aggregate structures, numerous studies have looked into the ionomer peak as the sample is deformed by uniaxial tension using small angle X ray scattering (SAXS).13–16 Meridional (parallel to the direction of deformation) and equatorial (perpendicular to the direction of deformation) parts of the scattering intensity are analyzed separately.

Generally, it is seen that the meridional peak shifts to lower and the equatorial peak shifts to higher wavevectors up to a certain draw ratio.14,15,17 The changes in scattering peak intensity during deformation is a subject of some debate. For amorphous E-MAA ionomers, Roche et al.14 observed that the meridional peak intensity decreased and finally disappeared after a draw ratio of 2.0. The equatorial peak intensity, however, increased monotonically. Visser and Cooper performed SAXS analysis to study the morphological response of polyurethane ionomers with ionic groups regularly spaced along the polymer backbone, during uniaxial deformation.15 They observed a decrease in scattering intensity with increasing strain in both the equatorial and meridional directions and attributed this to a more disordered aggregate morphology induced by stretching. Visser and Cooper18 later observed that neither the depleted core shell model19 nor the liquid-like hard sphere model9 accurately predicted the decrease in the meridional peak intensity of the ionomer peak or the change in equatorial peak intensity with elongation. Thus, it has been difficult to quantitatively predict the deformation behavior of ionomer melts based on simple models.

Recently, Winey et al. studied a range of unneutralized poly(ethylene-co-acrylic acid) (PEAA) materials with different acid contents.20 This work can potentially shed light on the nature of aggregate alignment during deformation because the backbone architecture was precisely controlled and related samples could be directly compared. However, the unneutralized PEAA materials have no ionic groups; the acid groups tend to aggregate in a similar manner as ionic groups but much more weakly. These studies suggested that for precise, semicrystalline materials in which uncharged segments are long enough to form small crystalline regions between acid aggregates, there is a transformation from liquid-like interaggregate ordering to layered ordering of acid aggregates as increasing uniaxial tensile deformation is imposed on the material.20 Amorphous polymers (with lower backbone spacings), however, do not exhibit this transition, and their liquid-like aggregate morphology was preserved throughout deformation.20 

Here we simulate coarse-grained analogs of precisely spaced PEAA materials neutralized with sodium. These materials are relatively well understood based on a body of experimental,16,20–22 atomistic simulation,23–26 and coarse-grained simulation5,27–30 work on related materials led by Winey, Frischknecht, and others. Briefly, using a coarse-grained model similar to that used here, Hall et al. showed that short, local ion rearrangement is relatively similar for most systems, while the longer time scale rearrangements are significantly different depending on ion spacing and the resulting aggregate morphology.5 It was also shown that the structure factor obtained using this model is in good agreement with the experimental scattering results for ionomers with varying spacing between charged groups.27,28 Here, we build upon such understanding of the structure and equilibrium dynamics to study the rheological properties of the coarse-grained model ionomers in detail. We note that many prior studies have used coarse-grained MD to investigate the stress relaxation behavior of entangled polymers31–33 and their response to tensile deformation.34–36 Additionally, tensile deformation studies have been performed on polymer glasses37,38 and polymer gels39 and to understand failure at interfaces.40 However, to the best of our knowledge, stress relaxation and uniaxial tensile deformation of coarse-grained charged polymers have not been reported before. In this work, we apply uniaxial tension and track the changes in a polymer and aggregate microstructure, including ionomer scattering peak location and intensity both perpendicular and parallel to the strain axis. We also calculate bulk viscoelastic properties from equilibrium simulations.

Following Refs. 5, 27, and 28, we use a simple coarse-grained model ionomer that represents precisely spaced, fully neutralized PEAA-based ionomers. This model was shown to reproduce the general features of the ionomer peak observed in experiments on partially neutralized materials, and we use the same methodology as in Refs. 5, 27, and 28 except where noted below. Briefly, the model polymers contain an uncharged, freely jointed, linear backbone to which anions are bound at a regular spacing; for each anion, an unbound counterion of equal and opposite charge is added to the system. We map to the experimental system as shown in Fig. 1; each uncharged bead corresponds to three CH2 groups along the backbone, while the anionic pendant group represents COO and the counterion corresponds to Na+. The monomer diameter maps to approximately 0.4 nm in real units.28 This implies that the charged and uncharged monomers should be approximately the same size,41 while the counterions are half that size, corresponding to the ionic radii of Na+(approximately 0.2 nm). Note that our polymer backbones have 35-36 beads, below the entanglement length for Kremer-Grest polymers.42,43

FIG. 1.

Schematic showing (a) a PEAA-based ionomer with 9 backbone carbons between acetic acid groups, neutralized with sodium, and [(b)-(d)] the model systems with the uncharged polymer beads in gray, anionic pendant bead in red, and counterion in blue: (b) Nbb = 3 which corresponds to the system in (a), (c) Nbb = 5, and (d) Nbb = 7.

FIG. 1.

Schematic showing (a) a PEAA-based ionomer with 9 backbone carbons between acetic acid groups, neutralized with sodium, and [(b)-(d)] the model systems with the uncharged polymer beads in gray, anionic pendant bead in red, and counterion in blue: (b) Nbb = 3 which corresponds to the system in (a), (c) Nbb = 5, and (d) Nbb = 7.

Close modal

Bonded monomers interact by a nonlinear finitely extensible (FENE) spring along with the repulsive part of the Lennard-Jones (LJ) potential, given by

U F E N E r = 0.5 k R 0 2 ln 1 r R 0 2 + 4 ε L J σ L J r 12 σ L J r 6 + ε L J , r 2 1 6 σ 0.5 k R 0 2 ln 1 r R 0 2 , 2 1 6 σ < r < R 0 , r R 0 ,
(1)

where r is the distance of separation, R0 = 1.5 σ L J is the maximum possible bead-bead separation, and k = 30 ε L J / σ L J 2 is the spring constant, making the spring short and stiff to avoid backbone crossing as in the work of Kremer and Grest.42 Here, we use ε L J as the characteristic unit of energy and σ L J as the unit of length. (We have included the LJ subscript here for clarity as ε and σ are commonly used to refer to strain and stress.)

Non-bonded monomers i and j interact via the truncated LJ potential,

U L J , i j r = 4 ε L J , i j σ L J , i j r 12 σ L J , i j r 6 + S , r r c 0 , r > r c ,
(2)

where ε L J , i j is their LJ interaction strength, σ L J , i j is their LJ diameter, and S is the constant shift factor such that the potential equals 0 at the cutoff r C . We include the attractive tail of the potential out to a distance of r C = 2.5 σ L J . This is a significant difference between our model and that of Refs. 5, 27, and 28, which used a shorter cutoff to yield a fully repulsive potential. While a shorter cutoff allows for faster calculations, we expect that including the local cohesion between monomers will be important for capturing the tensile response.44,45 We run our simulations at constant pressure in the transverse directions and under deformation; in contrast, the prior work used a constant volume rather than a constant pressure and did not consider deformation.

All beads have the same interaction strength ε L J . All polymer beads are of the same size, with an LJ diameter of 1.0 σ L J , while counterions are smaller, with an LJ diameter of 0.5 σ L J . All polymer beads and counterions have unit mass. We note that to consider multiple chemical types of monomers at the coarse-grained level, one can simply adjust the interaction strengths. Our group has recently studied segmented ionenes in such a manner,46 and our future work on related PEAA-based systems will use adjusted interaction strengths to account for the chemistry of uncharged pendant groups that represent the -COOH moieties present in partially neutralized materials. In the current work, all monomer units are chemically the same except that charges are added to the pendant monomers.

Ions experience a long ranged Coulombic interaction of the form

U c r = q 1 q 2 4 π ε o ε r r ,
(3)

where εr is the dielectric constant, the two charges are q1 and q2, and εo is the vacuum permittivity. The appropriate choice of q1 and q2 relative to εoεr depends on the dielectric constant of the medium and the mapping to real charge and length units. We set ion charges to +1 and −1 to represent +1e and −1e, map σ L J to 0.4 nm, and consider a system with an effective background dielectric constant of εr = 4. This leads us to set the dielectric constant (including the unit conversion) to 0.028; i.e., the interaction between two charges at a distance of σ L J is 1/0.028 or 36εLJ. A detailed description of the dielectric constant and its mapping to reduced units can be found in Ref. 28. Note that the background dielectric constant is greater than 1 to account for the chemical nature of the medium that has not been included in our coarse-grained model; for instance, in an atomistic version of the system, local partial charges on C and H atoms that are not present in our model would act to partially screen charges. This means that if the mapping to real units works well, the effective background dielectric used here should be similar to the dielectric constant of the unneutralized (uncharged) experimental system (because the ionic interactions are explicitly included but polarizability and partial charges are not). The value of εr = 4 was chosen here, following prior work, because it allowed equilibration in a reasonable simulation time while still showing strong ionic aggregation.

The LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) molecular dynamics package is used to perform the simulations,47,48 using periodic boundary conditions and a velocity-Verlet algorithm with a time step of 0.005τ. A Nose-Hoover thermostat and barostat were applied to sample the system in the NPT ensemble; however, stress autocorrelation data were collected in the NVE ensemble, in which case the average density from the NPT simulations was used to set the volume as described below. The barostat was set to keep the average pressure at 0 using a pressure damping parameter of 100τ. A temperature damping parameter of 1.0τ was used to maintain the temperature at 1.25T*, where the reduced temperature T* = kT. This is another apparent difference between the current work and Refs. 5, 27, and 28, which used the typical value 1T*. However, the two models do not behave similarly at the same temperature—the current model’s added local cohesion (increased LJ cutoff) essentially slows down its dynamics relative to the fully repulsive model. Considering the original mapping based on room temperature, this increased temperature means the current systems are held at approximately 100 °C. Note that choosing instead to scale down both εLJ and (1/εr) by 1.25 would lead to similar effects. Scaling up T* allowed us to otherwise keep similar parameters to Refs. 5, 27, and 28 while still having a system that equilibrates in a reasonable simulation time.

Each simulated system contains 800 ionomers, with backbone spacings of Nbb = 3, 5, or 7 between pendant groups (see Fig. 1), and their counterions. A reference homopolymer with no charges and a chain length of N = 36 was also simulated to compare with ionomer systems. Initially, polymers are placed as random walks and counterions added in random positions, and a standard soft potential is briefly applied in order to push overlapping beads off each other.49 The systems were equilibrated for 175 × 103τ (in the NPT ensemble) during which we tracked the box size, mean squared displacement (MSD) (after subtracting the total system center of mass displacement), and polymer end-to-end vector autocorrelation. Specifically, the Nbb = 7 system, which has the slowest dynamics, had an average polymer center of mass MSD of 49.6σ2 at 175 × 103τ, more than 5 times its mean squared radius of gyration of Rg2 = 9.0σ2, and we chose to run for another 105τ beyond this time. By the end of the equilibration period, all systems’ average end-to-end vector autocorrelation (normalized to 1 at time 0) had decayed to less than 1/e. The end of this equilibration period is the beginning of the deformation runs described below. We also calculated the average box size from the last 105τ of the equilibration period, during which time the box size was approximately constant.

Using this average NPT box size and otherwise initializing each system the same way as above, we also ran an NVT equilibration and NVE simulation to obtain accurate equilibrium stress fluctuation data (without noise associated with the barostat and thermostat). The system was equilibrated with the same thermostat as above (but no barostat) for 23 × 104τ, at which time the thermostat was turned off. After a brief further NVE equilibration of 104τ, NVE data were collected for analysis. The constant energy data collection period was 104τ for homopolymers, Nbb = 3 system, and Nbb = 5 system and 105τ for Nbb = 7.

We compute the stress-stress autocorrelation function G(t)33 through stress fluctuations in the equilibrium simulation in the NVE ensemble (with no applied strain), given by

G t = k B T 1 V σ α β ( t ) σ α β ( 0 )
(4)

with α β. Here, V is the volume of the simulation box, α and β are Cartesian components (x, y, or z), and σαβ is the stress. We calculate stress using the virial stress tensor

σ α β = 1 V i N m i v i α v i β + i N r i α f i β ,
(5)

where v and r are the velocity and position of particle i in the α direction, f is the force on particle i in the β direction, and mi is the mass of particle i. The first term is the kinetic energy contribution and the second term accounts for all bonded and non-bonded interactions.

We computed the stress at every time step so as to accurately average bond fluctuations. We perform a running average from 0.9t to 1.1t, for each time t, which helps to smooth the data; a similar averaging window has been used in prior work31,33 to improve the signal to noise ratio. The signal is also slightly improved by averaging over the off diagonal components (σxy, σyz, and σzx) since the system is isotropic.

The viscosity was computed by the following Green-Kubo relationship:

η 0 = k B T 1 V 0 σ α β t σ α β 0 d t .
(6)

Uniaxial tensile deformation was applied along the axial (x) direction of the initially cubic simulation box with a constant engineering strain rate by increasing the x-dimension (Lx) in finite steps d L x = ε ̇ L x 0 d t , where dt is the time step, ε ̇ is the strain rate, and Lx0 is the initial dimension of the box. During deformation, the barostat is applied only in the transverse directions (y and z) to simulate constant pressure in those directions. The extension ratio (ratio of box length and original box length) is varied at every MD time step. The resulting stress, used to create stress-strain curves, was obtained by performing a block average of the axial component of the virial stress tensor ( σ x x ) with a block averaging width of 100 time steps. We thus can obtain stress-strain curves from this strain driven simulation; relative to stress-driven simulations, this approach allows for more efficient simulations and also allows us to appropriately observe the relatively constant stress plateau at higher strains.50 This type of ensemble where the dimension of the box is fixed in the axial direction and the stresses are fixed in the transverse direction (NLxσyyσzzT) mimics experimental conditions of uniaxial tension.51 

Three deformation runs for each system were performed, with strain rates of 0.326 × 10−3/τ, 0.65 × 10−3/τ, and 1.3 × 10−3/τ (each starting from the system’s equilibrated cubic box). One mapping of our systems to real units would suggest τ ≈ 1.7 × 10−12 s. However, coarse-grained simulations often require an additional factor when mapping to experimental time because removing the degrees of freedom generally speeds up the system; see discussion in the supplementary material. This naive mapping suggests that our rates are 2 × 108–8 × 108 s−1. For context, typical experimental strain rates for polymer melts undergoing uniaxial tensile deformation range from 0.01 to 1 s−1. Rates of 10-1000 s−1 are the upper limits of achievable strain rates from these experiments.52 However, strain rates in MD simulations found in the literature typically range from 106 to 1010 s−1 because lower rates require longer simulations and rates approaching those used in experiments are typically intractable.53,54 Our chosen strain rates allow us to explore the effects of varying the rate of deformation on aggregate alignment and mechanical properties, and they are within the typical limits of simulation strain rates.

Prior coarse grained ionomer work with a similar model reported the ion-ion scattering structure factor and compared it to the experimentally observed ionomer peak for precisely and randomly spaced partially neutralized PEAA materials; the comparison showed good qualitative agreement between the strength and location of the low wavevector peak and its dependence on backbone spacing. Thus, we consider the ion-ion scattering an appropriate, experimentally relevant indicator of the mesoscale order of ionic aggregates and calculate Sion-ion(k) in the same manner as the prior work (and also do not include the backbone scattering).27,28 The ion-ion partial static structure factor is

S i o n i o n k = 1 N i o n i = 1 N i o n j = 1 N i o n e i k r i r j ,
(7)

where Nion is the total number of ions in the system, ri,j are the vector positions of the ionic species (both anion and counterion), and k is the wavevector; the components of k in each Cartesian direction are discretized based on the size of the simulation box in that direction. For an isotropic system, results for wavevectors of similar lengths are typically averaged together; however, as the box is being deformed, the structure is expected to be different in the directions parallel (x) and perpendicular (y, z) to the strain axis, so we have separately calculated these structure factors. In order to improve the statistics, a thin slice of the vectors within 15° of the strain axis was considered while calculating the structure factor in the parallel direction.20 The strain points chosen for scattering analysis were 0–1 in increments of 0.1, with the addition of 0.05 and 0.15 to better show the changes at low strain. Frames within 2% strain of the chosen strain point were averaged to reduce noise in the data. The peak height and location reported are those of a parabola fit to the maximum point and the two points on either side. To further improve accuracy, we ran a second deformation simulation with its starting configuration obtained from the same equilibrated simulation as the first, but 2 × 104τ earlier in time. The peak height values were noisy, but peak locations were less noisy and were consistent across the two runs (see details in the supplementary material); peak location values reported here are an average of the two runs.

The stress-stress autocorrelation, calculated as described above from equilibrium simulations, is reported in Fig. 2. At short times, oscillations in G(t) due to bond and local fluctuations are seen,33 which are qualitatively similar for the three ionomer systems and the homopolymer reference system [Fig. 2(a), plotted on a linear scale]. G(t) is initially larger for systems with higher ion content, but at very long times the order of the ionomer curves reverses; G(t) decays slowest for ionomers with lower ion content [see Fig. 2(b), plotted on a log-log scale]. Specifically, the Nbb = 3 system shows a short and intermediate time scaling of t−1/2, similar to the Rouse relaxation for unentangled polymer melts.33,55 The Nbb = 5 and Nbb = 7 systems exhibit similar behavior at short times, though a distinct deviation from the Rouse relaxation is seen at long times due to ionic aggregation as discussed further below. Note that the stresses in the ionomer systems are considerably higher compared to the homopolymer, which is expected due to the increased overall cohesion due to the ionic interactions.

FIG. 2.

Stress-stress autocorrelation function for ionomer and homopolymer systems as labeled, plotted (a) linearly to show the behavior at early times and (b) on a log-log scale with a reference line showing t−0.5 scaling; data are connected by thin lines and every 200th point is shown with a symbol.

FIG. 2.

Stress-stress autocorrelation function for ionomer and homopolymer systems as labeled, plotted (a) linearly to show the behavior at early times and (b) on a log-log scale with a reference line showing t−0.5 scaling; data are connected by thin lines and every 200th point is shown with a symbol.

Close modal

To understand the deviation from the Rouse relaxation for ionomers with long spacers, one must consider both the morphology of the clusters and the time scales over which cluster rearrangement occurs.5,28 The ionic aggregate morphology in very closely related systems has been discussed in detail in Ref. 28; briefly, virtually all the ions in the Nbb = 3 system belong to a single percolated ionic aggregate, most of the ions in the Nbb = 5 system belong to a large percolating aggregate with some belonging to discrete clusters, and all the ions in the Nbb = 7 system belong to discrete clusters of relatively uniform size.28 There are two time scales of ion relaxation related to distinct dynamic processes.5 Local rearrangements of nearby ions within the same discrete cluster (or the same local area of a percolated cluster) occur on a relatively short time scale. A longer time, larger scale collective rearrangement of ions is prominent and easily observed in the system of discrete clusters and is also seen in the larger scale rearrangements of percolated aggregates. For percolated aggregates, ions can be transported further through multiple local relaxations, and the collective rearrangements occur more easily and on a relatively shorter time scale than those of discrete aggregates due to the high ion density. Overall, ion correlations in systems with discrete clusters (lower ion content) take longer to relax.

We surmise that the collective rearrangement required to fully relax ions within systems of discrete clusters is the reason why unlike the Nbb = 3 system, both the Nbb = 5 and 7 systems show a different scaling behavior after the initial scaling of t−0.5. The divergence from t−0.5 for Nbb = 5 sets in around 40τ, where the scaling changes to approximately t−0.3. For Nbb = 7, this divergence begins earlier, at approximately 15τ. This is because the clusters in the Nbb = 7 system are smaller and completely discrete, and the local rearrangements serve in relaxing only a small part of the stress in the system, after which only the collective rearrangement of clusters can relax the rest of the stress. Many ions in the Nbb = 5 system belong to a percolated cluster, due to which the stresses relax faster compared to the Nbb = 7 system, but complete relaxation does not occur until longer times, when the collective rearrangements can fully relax discrete clusters as well.

Because of this, the stress relaxation curves for the ionomer systems do not collapse into a master curve at early times, unlike what is expected for entangled homopolymer systems.31,56,57 The work of Leibler et al.7 discussed the features of G(t) and time scales that are important for entangled, associating polymers; associations can cause a plateau in G(t) that is qualitatively similar to that caused by entanglements, but on earlier time scales. In the current system, without entanglements, ionic associations cause a deviation from the Rouse behavior, but not a flat plateau. This is unsurprising as a completely flat plateau region requires a large, clear disparity in relaxation times and a perfectly flat region is difficult to access in MD simulations of entangled bead-spring polymer melts. In such simulations at N > 200, there is a discernible deviation from the scaling of t−0.5, which is sufficient to indicate entanglements.31,32

To further understand long time relaxation, we also calculated viscosity using the Green-Kubo relationship as outlined in the section titled Methods. The viscosities were 180, 420, and 950 for Nbb = 3, 5, and 7, respectively. This reflects the trend in both stress autocorrelation and polymer end-to-end relaxation times, reported in Fig. 3 (higher viscosities correspond to longer relaxations). We define the stress relaxation time as when the autocorrelation decays to 0.01 (crosses bottom axis of Fig. 2) and the end-to-end relaxation time as when the polymer end-to-end vector autocorrelation A C F e 2 e = R e 2 e ( t ) R e 2 e ( 0 ) / R e 2 e ( 0 ) R e 2 e ( 0 ) , which appears to roughly exponentially decay, crosses 1/e. While the trends are the same for both, with relaxation increasing with increasing spacer length, the end-to-end relaxation times are higher than the stress relaxation times. The polymer end-to-end vector decorrelates in both orientation and magnitude, but stress-stress autocorrelation is reflective of only that portion of the stress that has not relaxed. This was also observed by Sen et al.,33 who concluded that orientational decorrelation is very limiting for stress relaxation processes. It is important to note that our model does include unneutralized acid groups, which would affect dynamics particularly at long times. This may be why the trend in our viscosity results differs from the typical experimentally observed trend (for partially neutralized, lower overall ion content systems) that an increase in ion density gives rise to higher viscosities and longer stress relaxation times.6,58 We have recently created a coarse-grained model of partially neutralized systems to better match such experimental systems, which will be discussed in a future publication.

FIG. 3.

Relaxation times from the stress-stress and end-to-end autocorrelation functions as labeled, as a function of spacer length.

FIG. 3.

Relaxation times from the stress-stress and end-to-end autocorrelation functions as labeled, as a function of spacer length.

Close modal

The stress-strain curves for all systems at three strain rates are shown in Fig. 4. At each strain rate, the tensile response of the ionomers is much more pronounced than the homopolymer, which is expected due to the additional cohesion and mechanical strength imparted by the ionic groups. The three ionomer systems have a relatively similar response at the lowest strain rate but show clear differences at larger strain rates, for which the Nbb = 3 system is the toughest. This transition in the response of polymers to the applied strain rate is a phenomenon that has been studied widely, and ionomers obey the general trend of having greater yield stress at higher strain rates.59–61 Hirasawa et al.62 and Gao et al.63 also observed this trend of higher ion densities leading to an increase in the tensile strength of ionomers.

FIG. 4.

Stress-strain curves for ionomer and homopolymer systems at three different strain rates, as labeled; dashed lines are fit to the approximately linear initial region to show the apparent differences in the elastic modulus; data are connected by thin lines and only 30% of the data points are marked with a symbol for clarity.

FIG. 4.

Stress-strain curves for ionomer and homopolymer systems at three different strain rates, as labeled; dashed lines are fit to the approximately linear initial region to show the apparent differences in the elastic modulus; data are connected by thin lines and only 30% of the data points are marked with a symbol for clarity.

Close modal

Although the viscosity of the Nbb = 3 system is the lowest, these results are consistent with the fact that its stresses are the highest in the low time, high frequency regime [Fig. 2(b)]. Though the stress-strain response is similar for all systems at the low strain rate (which allows for longer time relaxation during strain), at higher strain rates, the Nbb = 3 system’s yield stress and final stress are highest (implying that the Nbb = 3 system has the highest toughness). At even lower strain rates than shown in Fig. 4, a reversal of trend can be seen, with Nbb = 3 exhibiting the least toughness of the three systems (see the supplementary material).

Overall, Fig. 4 shows that, for higher strain rates, yielding (transition from linear, elastic response to a flatter, plastic response) occurs at a higher stress and lower strain, followed by a strain-hardening upturn. Middleton et al. observed a similar change in the stress-strain response of PEAA; at low strain rates, the material began to yield at lower strain and did not show significant strain hardening, while at high strain rates, the system was tougher and exhibited strain hardening.20 

To understand the orientation of the polymer chains during deformation, we calculated a global orientation order parameter, f 2 = 3 cos 2 θ 1 / 2 , where θ is the angle between the polymer end-to-end vector and the strain axis. The results are shown for the three architectures at the same strain rate, as well as for the Nbb = 7 system at two different strain rates in Fig. 5. We see significant alignment in the direction of strain, which increases over time; this result is intuitive, and similar chain alignment has been reported previously in other simulated polymers.64 The alignment increases as the spacing between charged pendants decreases or as the clusters become percolated. It seems that as the sample is stretched, polymer backbones are able to orient more easily when their ions can rearrange within a large percolated aggregate, but this type of motion is hindered when the clusters are discrete. This trend of higher ion concentration leading to greater chain orientation was also observed experimentally by Zhao et al.71 

FIG. 5.

Orientation parameter describing backbone alignment for ionomer and homopolymer systems during strain, as labeled; data are connected by thin lines and only 20% of the data points are marked with a symbol for clarity.

FIG. 5.

Orientation parameter describing backbone alignment for ionomer and homopolymer systems during strain, as labeled; data are connected by thin lines and only 20% of the data points are marked with a symbol for clarity.

Close modal

The Nbb = 3 system consistently shows greater alignment than the other two systems at all three strain rates; however, the Nbb = 3 system only shows a significantly greater toughness at the higher strain rates. This supports the notion that overall chain alignment does not play as important a role in the mechanical properties of these systems as do the ionic associations. Also note that the orientation parameter increases quickly at low strain (almost linearly), then more slowly at higher strain, and the change in the curves occurs at a strain that is after the yield point from the stress-strain curves. This means that polymer orientation alone does not predict yielding behavior. The ionomer systems and the homopolymer have similar chain alignment until a strain of approximately 0.4 [Fig. 5(b)], after which the ionomer chains continue to align, whereas the homopolymer chains keep a fairly constant orientation. Figure 5(a) also shows that the orientation parameter increases in magnitude as the strain rate increases, implying that the chains align to a greater degree when the system is deformed more rapidly, as they do not have time to relax.

Snapshots showing backbone conformations in the box for the Nbb = 7 system as it is being deformed are given in Fig. 6. All images were created using VMD.70 Polymers initially appear coiled and are clearly seen to align in the direction of pull. Note that the system does not undergo strain induced crystallization; however, crystallization is observed in experiments for long enough spacer lengths.20 Standard Kremer-Grest chains as used here do not tend to crystallize (neither due to strain nor at low temperatures) because of their coarse-grained nature. Atomistic simulations, such as of linear polyethylene chains, naturally reproduce strain-induced crystallization behavior,65 as do specialty coarse-grained models which could be adapted for use in future ionomer work if understanding such effects is of interest.66,67

FIG. 6.

All ions (colored red and made transparent) and several selected chain backbones (opaque white) of the Nbb = 7 system are shown in their initial state (top) and after uniaxial deformation to a strain of 1.0 (bottom).

FIG. 6.

All ions (colored red and made transparent) and several selected chain backbones (opaque white) of the Nbb = 7 system are shown in their initial state (top) and after uniaxial deformation to a strain of 1.0 (bottom).

Close modal

We analyzed ionic aggregate morphological changes during deformation by performing a cluster analysis; any two ions within a distance of 0.85σ were defined to be in the same cluster. Figure 7 shows snapshots of the ionic clusters before and after uniaxial tensile deformation, up to a strain of 1.0. It is visually apparent that the clusters are somewhat aligned along the direction of pull. For Nbb = 3, the system with the highest ion density, all ion pairs belong to one percolated cluster before and after deformation. In contrast, the Nbb = 7 system contains only discrete clusters, some of which broke apart during deformation; specifically, the number of clusters increased from 157 at a strain of 0 to 169 at a strain of 1.0. The change in the radius of gyration of the clusters is quantified in the supplementary material (Fig. S3).

FIG. 7.

Ions of (a) Nbb = 3, (b) Nbb = 5, and (c) Nbb = 7 systems colored by cluster before deformation (a large percolated aggregate is shown in green for Nbb = 3, 5; for Nbb = 7, all aggregates are discrete): each system is shown in its initial state (top) and after uniaxial deformation to a strain of 1.0 (bottom). Boxes are shown to the same scale; note that different systems have slightly different sizes and different numbers of monomers and counterions.

FIG. 7.

Ions of (a) Nbb = 3, (b) Nbb = 5, and (c) Nbb = 7 systems colored by cluster before deformation (a large percolated aggregate is shown in green for Nbb = 3, 5; for Nbb = 7, all aggregates are discrete): each system is shown in its initial state (top) and after uniaxial deformation to a strain of 1.0 (bottom). Boxes are shown to the same scale; note that different systems have slightly different sizes and different numbers of monomers and counterions.

Close modal

To further quantify the structural change in the aggregates upon deformation, we show anion-cation pair correlation functions for all three systems before and during deformation at three strain rates in Fig. 8. The first sharp peak at 0.75σ corresponds to the preferred distance at contact, while the second peak is a feature of ion packing within aggregates. As the system undergoes deformation, the second peak decreases. It decreases by a greater amount when the system is deformed more quickly. Ting et al.68 observed a similar change in the pair correlation function while applying an external field to coarse grained ionomers, which indicated that the clusters were being stretched. We also note that the discrete aggregates in the Nbb = 7 system show a greater degree of stretching as compared to the percolated cluster systems, apparently due to their slower ionic aggregate relaxation time. Experimentally, ionomer chain orientation during strain was quantified by Ding et al.69 for systems with different cations; changing the cation significantly changes the aggregate morphology and degree of chain alignment.

FIG. 8.

Charged bead-counterion pair correlation functions for the three ionomer systems at different strain rates, as labeled; data are connected by thin lines and only every other data point is marked with a symbol for clarity.

FIG. 8.

Charged bead-counterion pair correlation functions for the three ionomer systems at different strain rates, as labeled; data are connected by thin lines and only every other data point is marked with a symbol for clarity.

Close modal

To show how deformation impacts the relaxation of discrete clusters, we report the cluster-cluster autocorrelation function, which is the normalized probability that two ions that were in the same cluster at time 0 are in the same cluster at time t. We see that the clusters decorrelate more rapidly when the strain rate increases (Fig. 9). The cluster relaxation time (when the autocorrelation reaches a value of 1/e) for an undeformed system (strain rate of 0) is 7133τ, which is approximately eight times longer than 920τ, the time over which the system was strained to 1.0 at a strain rate of 1.3 × 10−3/τ.

FIG. 9.

Cluster autocorrelation function for the Nbb = 7 system with discrete clusters, at different strain rates (with time 0 being the initial state just before deformation) or from the equilibrium simulation (no strain), as labeled.

FIG. 9.

Cluster autocorrelation function for the Nbb = 7 system with discrete clusters, at different strain rates (with time 0 being the initial state just before deformation) or from the equilibrium simulation (no strain), as labeled.

Close modal

The ion-ion structure factor in the coarse grained ionomer model used here was previously shown to yield the same qualitative features as the ionomer peak observed in scattering experiments; experiments and simulations also showed agreement in the change in peak location as a function of ion spacing. However, in considering the relationship to experimental scattering, one should compare peak heights only qualitatively. It is important to note that the coarse-grained model does not take into account the atomic-scale bonding and packing details, does not crystallize, and does not include unneutralized COOH groups that may be involved in aggregate formation. Here, we compute the anisotropic structure factor in the directions parallel to and perpendicular to the strain axis as described in the section titled Methods.

The main feature in the ion-ion structure factor is the first peak, corresponding to the so-called “ionomer peak” due to aggregate ordering. Representative structure factors are shown in the supplementary material. We calculated both the location and intensity of the ionomer peak as the sample was deformed. In contrast to equilibrium simulations, the data cannot be averaged over a long time interval (at least, not while considering a particular strain point), and the statistics were poor. We compared the results across two independent simulation runs; while the location of the peak changed relatively smoothly during deformation and was similar for different runs, the peak heights (see the supplementary material) were more varied and were not quantitatively reproducible across the two runs. We report only the peak location in Figs. 10 and 11, restricting strains in Fig. 10 (for the axial direction) to less than 0.5, after which the peak sometimes splits into multiple peaks with a less obvious overall location (likely due to insufficient statistics).

FIG. 10.

Position of the first peak in the ion-ion structure factor as a function of strain in the direction parallel to strain axis; (a) at a strain of 1.3 × 10−3/τ for the three ionomer systems as labeled and (b) for the Nbb = 7 system at the three different strain rates as labeled (all values are the average from two deformation runs for all systems; the difference between runs was at most 7% for any given data point).

FIG. 10.

Position of the first peak in the ion-ion structure factor as a function of strain in the direction parallel to strain axis; (a) at a strain of 1.3 × 10−3/τ for the three ionomer systems as labeled and (b) for the Nbb = 7 system at the three different strain rates as labeled (all values are the average from two deformation runs for all systems; the difference between runs was at most 7% for any given data point).

Close modal
FIG. 11.

Position of the first peak in the ion-ion structure factor as a function of strain in the direction perpendicular to strain axis; (a) at a strain of 1.3 × 10−3/τ for the three ionomer systems as labeled and (b) for the Nbb = 7 system at the three different strain rates as labeled (all values are the average from two deformation runs for all systems; the difference between runs was at most 4% for any given data point).

FIG. 11.

Position of the first peak in the ion-ion structure factor as a function of strain in the direction perpendicular to strain axis; (a) at a strain of 1.3 × 10−3/τ for the three ionomer systems as labeled and (b) for the Nbb = 7 system at the three different strain rates as labeled (all values are the average from two deformation runs for all systems; the difference between runs was at most 4% for any given data point).

Close modal

Specifically, Fig. 10 showing the peak location in the axial direction reveals that the clusters are ordering on a longer length scale (moving farther apart in the direction in which the material is pulled, as indicated by a decrease in the wavevector). This initial change in peak wavevector persists up to a certain strain, after which the peak location remains almost constant. The change in the position of clusters does not coincide with the yield point of these systems. Similarly Fig. 11 shows the peak location in the transverse direction. Clusters are initially ordering at a shorter length scale in the transverse direction (apparently moving closer together) until having a relatively constant length scale of order after a particular strain. These results are in agreement with prior experimental findings that the ionomer peak shifts to lower wavevectors in the axial direction and higher wavevectors in the transverse direction during deformation.14,15,17

We simulated simple coarse-grained model ionomers with different ion spacings in equilibrium and under uniaxial tensile deformation. We presented the stress-stress autocorrelation function, stress-strain curves, and various measures of polymer and aggregate alignment, including the anisotropic structure factor peak in directions perpendicular and parallel to strain, to understand the behavior of ionic aggregates during deformation and their relationship to mechanical properties.

Trends in the location of the ionomer peak in the scattering during deformation are in qualitative agreement with prior experimental results.14,15,17 The peak in the meridional (parallel) scattering shifts to a lower wavevector and the equatorial (perpendicular) peak shifts to a higher wavevector as the sample is deformed. This shows that the aggregate ordering shifts to longer length scales (aggregates move apart) in the direction of pull, while aggregates tend to move closer together in the perpendicular direction, as expected. These trends are enhanced at higher strain rates. At high enough strain, the S(k) peak position approaches a constant value, indicative of the fact that the ion pairs eventually are able to rearrange to keep a preferred length scale of ordering in each direction. We also quantified the degree to which both ionic aggregates and chain backbones align in the direction of pull. Chains align significantly in response to deformation, and their alignment is clearly visible in simulation snapshots, while ionic aggregates are not as clearly aligned. This is intuitive as the aggregates can break and reform during the simulation, to keep their average local structure, which is set by strong Coulombic interactions.

As has been previously understood, decreasing spacing between ions (increasing the ion density) allows for larger, percolated ionic aggregates. The ion content and resulting aggregate structure impact the mechanical properties in different ways on short and long times and length scales. While local ion rearrangements are possible and occur at relatively short times in all systems, percolated aggregates allow for easier ion rearrangements over the longer length scales of the system corresponding to full relaxation of chains. Thus, complete stress relaxation is accomplished at an earlier time in higher ion content systems, which in turn have a lower viscosity. However, increasing ion density increases the overall cohesion of the system, generally increasing strength on shorter time scales. In particular, at high strain rates (during which the system experiences a doubling of length over times shorter than the overall chain relaxation time), the higher ion content systems show greater toughness. This work considered fully neutralized materials only, though typical experimental materials are partially neutralized. In future work, we plan to consider the impact of uncharged associating groups (e.g., the COOH groups present in PEAA materials).

See supplementary material for additional data sets (stress-strain curve at a low strain rate, orientation parameter, radius of gyration, and anisotropic structure facture) and time mapping between LJ and real units for the model studied in this paper.

We thank Amalie Frischknecht, Karen Winey, and L. Robert Middleton for helpful discussions during the course of this work; Jonathan Brown for providing the cluster analysis script; Connor Barber for his help in updating that script; and Ting Ge for advice regarding calculating the stress autocorrelation function. This material is based on work supported by the National Science Foundation under Grant No. 1463103. This work benefited from an allocation in computing time from the Ohio Supercomputer Center. We also acknowledge partial support of this work from the H.C. “Slip” Slider Professorship in Chemical and Biomolecular Engineering.

1.
M. R.
Tant
,
K. A.
Mauritz
, and
G. L.
Wilkes
,
Ionomers: Synthesis, Structure, Properties and Application
(
Blackie Academic & Professional
,
1997
).
2.
A.
Eisenberg
and
J. S.
Kim
,
Introduction to Ionomers
(
Wiley Interscience
,
New York
,
1998
).
3.
R. A.
Weiss
and
H.
Zhao
,
J. Rheol.
53
,
191
(
2009
).
4.
L.
Zhang
,
N. R.
Brostowitz
,
K. A.
Cavicchi
, and
R. A.
Weiss
,
Macromol. React. Eng.
8
,
81
(
2014
).
5.
L. M.
Hall
,
M. J.
Stevens
, and
A. L.
Frischknecht
,
Macromolecules
45
,
8097
(
2012
).
6.
N. K.
Tierney
and
R. A.
Register
,
Macromolecules
35
,
6284
(
2002
).
7.
L.
Leibler
,
M.
Rubinstein
, and
R. H.
Colby
,
Macromolecules
24
,
4701
(
1991
).
8.
J.
Le Meins
and
J.
Tassin
,
Colloid Polym. Sci.
281
,
283
(
2003
).
9.
D. J.
Yarusso
and
S. L.
Cooper
,
Macromolecules
16
,
1871
(
1983
).
10.
D. J.
Kinning
and
E. L.
Thomas
,
Macromolecules
17
,
1712
(
1984
).
11.
S.
Kutsumizu
,
H.
Tagawa
,
Y.
Muroga
, and
S.
Yano
,
Macromolecules
33
,
3818
(
2000
).
12.
Y.
Tsujita
,
M.
Yasuda
,
M.
Takei
,
T.
Kinoshita
,
A.
Takizawa
, and
H.
Yoshimizu
,
Macromolecules
34
,
2220
(
2001
).
13.
T.
Hashimoto
,
S.
Sakurai
,
M.
Morimoto
,
S.
Nomura
,
S.
Kohjiya
, and
T.
Kodaira
,
Polymer
35
,
2672
(
1994
).
14.
E. J.
Roche
,
R. S.
Stein
,
T. P.
Russell
, and
W. J.
Macknight
,
J. Polym. Sci., Polym. Phys. Ed.
18
,
1497
(
1980
).
15.
S. A.
Visser
and
S. L.
Cooper
,
Polymer
33
,
4705
(
1992
).
16.
M. E.
Seitz
,
C. D.
Chan
,
K. L.
Opper
,
T. W.
Baughman
,
K. B.
Wagener
, and
K. I.
Winey
,
J. Am. Chem. Soc.
132
,
8165
(
2010
).
17.
K. A.
Page
,
F. A.
Landis
,
A. K.
Phillips
, and
R. B.
Moore
,
Macromolecules
39
,
3939
(
2006
).
18.
S. A.
Visser
and
S. L.
Cooper
,
Macromolecules
25
,
2230
(
1992
).
19.
M.
Fujimara
,
H.
Takeji
, and
H.
Kawai
,
Macromolecules
15
,
136
(
1982
).
20.
L. R.
Middleton
,
S.
Szewczyk
,
J.
Azoulay
,
D.
Murtagh
,
G.
Rojas
,
K. B.
Wagener
,
J.
Cordaro
, and
K. I.
Winey
,
Macromolecules
48
,
3713
(
2015
).
21.
A.
Taubert
and
K. I.
Winey
,
Macromolecules
35
,
7419
(
2002
).
22.
E.
Boz
,
A. J.
Nemeth
,
I.
Ghiviriga
,
K.
Jeon
,
R. G.
Alamo
, and
K. B.
Wagener
,
Macromolecules
40
,
6545
(
2007
).
23.
C. F.
Buitrago
,
D. S.
Bolintineanu
,
M. E.
Seitz
,
K. L.
Opper
,
K. B.
Wagener
,
M. J.
Stevens
,
A. L.
Frischknecht
, and
K. I.
Winey
,
Macromolecules
48
,
1210
(
2015
).
24.
D. S.
Bolintineanu
,
M. J.
Stevens
, and
A. L.
Frischknecht
,
ACS Macro Lett.
2
,
206
(
2013
).
25.
D. S.
Bolintineanu
,
M. J.
Stevens
, and
A. L.
Frischknecht
,
Macromolecules
46
,
5381
(
2013
).
26.
S. M.
Waters
,
J. D.
McCoy
,
A. L.
Frischknecht
, and
J. R.
Brown
,
J. Chem. Phys.
140
,
014902
(
2014
).
27.
L. M.
Hall
,
M. J.
Stevens
, and
A. L.
Frischknecht
,
Phys. Rev. Lett.
106
,
127801
(
2011
).
28.
L. M.
Hall
,
M. E.
Seitz
,
K. I.
Winey
,
K. L.
Opper
,
K. B.
Wagener
,
M. J.
Stevens
, and
A. L.
Frischknecht
,
J. Am. Chem. Soc.
134
,
574
(
2012
).
29.
D.
Ruan
and
D. S.
Simmons
,
Macromolecules
48
,
2313
(
2015
).
30.
D.
Ruan
and
D. S.
Simmons
,
J. Polym. Sci., Part B: Polym. Phys.
53
,
1458
(
2015
).
31.
A. E.
Likhtman
,
S. K.
Sukumaran
, and
J.
Ramirez
,
Macromolecules
40
,
6748
(
2007
).
32.
Q.
Zhou
and
R. G.
Larson
,
Macromolecules
39
,
6737
(
2006
).
33.
S.
Sen
,
S. K.
Kumar
, and
P.
Keblinski
,
Macromolecules
38
,
650
(
2005
).
34.
J. D.
Davidson
and
N. C.
Goulbourne
,
Modell. Simul. Mater. Sci. Eng.
24
,
065002
(
2016
).
35.
C.
Svaneborg
,
R.
Everaers
,
G. S.
Grest
, and
J. G.
Curro
,
Macromolecules
41
,
4920
(
2008
).
36.
F.
Leonforte
,
Phys. Rev. E
82
,
041802
(
2010
).
37.
R. A.
Riggleman
,
G. N.
Toepperwein
,
G. J.
Papakonstantopoulos
, and
J. J.
De Pablo
,
Macromolecules
42
,
3632
(
2009
).
38.
R. A.
Riggleman
,
K.
Yoshimoto
,
J. F.
Douglas
, and
J. J.
De Pablo
,
Phys. Rev. Lett.
97
,
045502
(
2006
).
39.
Y. R.
Sliozberg
,
R. A.
Mrozek
,
J. D.
Schieber
,
M.
Kröger
,
J. L.
Lenhart
, and
J. W.
Andzelm
,
Polymer
54
,
2555
(
2013
).
40.
T.
Ge
,
G. S.
Grest
, and
M. O.
Robbins
,
Macromolecules
47
,
6982
(
2014
).
41.
G.
Kamath
,
F.
Cao
, and
J. J.
Potoff
,
J. Phys. Chem. B
108
,
14130
(
2004
).
42.
K.
Kremer
and
G. S.
Grest
,
J. Chem. Phys.
92
,
5057
(
1990
).
43.
R. S.
Hoy
,
K.
Foteinopoulou
, and
M.
Kröger
,
Phys. Rev. E
80
,
031803
(
2009
).
44.
G. J.
Papakonstantopoulos
,
M.
Doxastakis
,
P. F.
Nealey
,
J. L.
Barrat
, and
J. J.
De Pablo
,
Phys. Rev. E
75
,
031803
(
2007
).
45.
H.
Chao
and
R. A.
Riggleman
,
Polymer
54
,
5222
(
2013
).
46.
P.
Vijayaraghavan
,
J. R.
Brown
, and
L. M.
Hall
,
Macromol. Chem. Phys.
217
,
930
(
2016
).
47.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
48.
S. J.
Plimpton
,
R.
Pollock
, and
M. J.
Stevens
, in
Proceedings of Eighth SIAM Conference on Parallel Processing for Scientific Computing
,
Minneapolis
,
1997
.
49.
R.
Auhl
,
R.
Everaers
,
G. S.
Grest
,
K.
Kremer
, and
S. J.
Plimpton
,
J. Chem. Phys.
119
,
12718
(
2003
).
50.
B. M.
Aguilera-Mercado
,
C.
Cohen
, and
F. A.
Escobedo
,
Macromolecules
47
,
840
(
2014
).
51.
S.
Lee
and
G. C.
Rutledge
,
Macromolecules
44
,
3096
(
2011
).
52.
J.
Richeton
,
S.
Ahzi
,
K. S.
Vecchio
,
F. C.
Jiang
, and
R. R.
Adharapurapu
,
Int. J. Solids Struct.
43
,
2318
(
2006
).
53.
L.
Yang
,
D. J.
Srolovitz
, and
A. F.
Yee
,
J. Chem. Phys.
107
,
4396
(
1997
).
54.
F. M.
Capaldi
,
M. C.
Boyce
, and
G. C.
Rutledge
,
Polymer
45
,
1391
(
2004
).
55.
J. D.
Halverson
,
W. B.
Lee
,
G. S.
Grest
,
A. Y.
Grosberg
, and
K.
Kremer
,
J. Chem. Phys.
134
,
204905
(
2011
).
56.
J.
Hou
,
C.
Svaneborg
,
R.
Everaers
, and
G. S.
Grest
,
Phys. Rev. Lett.
105
,
068301
(
2010
).
57.
M.
Rubinstein
and
R. H.
Colby
,
Polymer Physics
(
Oxford
,
2003
).
58.
Q.
Chen
,
S.
Liang
,
H.
Shiau
, and
R. H.
Colby
,
Macromolecules
2
,
970
(
2013
).
59.
Y.
Shui-sheng
,
L.
Yu-bin
, and
C.
Yong
,
Lat. Am. J. Solids Struct.
10
,
833
(
2013
).
60.
S.
Deschanel
,
B. P.
Greviskes
,
K.
Bertoldi
,
S. S.
Sarva
,
W.
Chen
,
S. L.
Samuels
,
R. E.
Cohen
, and
M. C.
Boyce
,
Polymer
50
,
227
(
2009
).
61.
R.
Duckett
,
S.
Rabinowitz
, and
I.
Ward
,
J. Mater. Sci.
5
,
909
(
1970
).
62.
E.
Hirasawa
,
Y.
Yamamoto
,
K.
Tadano
, and
S.
Yano
,
J. Appl. Polym. Sci.
42
,
351
(
1991
).
63.
Y.
Gao
,
N. R.
Choudhury
, and
N. K.
Dutta
,
J. Appl. Polym. Sci.
124
,
2908
(
2011
).
64.
N. M.
Lacevic
,
R. S.
Maxwell
,
A.
Saab
, and
R. H.
Gee
,
J. Phys. Chem. B
110
,
3588
(
2006
).
65.
M. S.
Lavine
,
N.
Waheed
, and
G. C.
Rutledge
,
Polymer
44
,
1771
(
2003
).
66.
C.
Luo
and
J. U.
Sommer
,
Comput. Phys. Commun.
180
,
1382
(
2009
).
67.
C.
Luo
and
J. U.
Sommer
,
Phys. Rev. Lett.
112
,
195702
(
2014
).
68.
C. L.
Ting
,
M. J.
Stevens
, and
A. L.
Frischknecht
,
Macromolecules
48
,
809
(
2015
).
69.
Y. S.
Ding
,
R. A.
Register
,
C.-Z.
Yang
, and
S. L.
Cooper
,
Polymer
30
,
1204
(
1989
).
70.
W.
Humphery
,
A.
Dalke
, and
K.
Schlten
,
J. Mol. Graph.
14
,
33
(
1996
).
71.
Y.
Zhao
,
C. G.
Bazuin
, and
R. E.
Prud’homme
,
Macromolecules
22
,
3788
(
1989
).

Supplementary Material