Simulating Startup Shear of Entangled Polymer Melts - ACS Publications

Nov 23, 2015 - Stress overshoot can be observed in entangled polymer melts during startup shear at sufficient higher shear rates. Such phenomenon can ...
1 downloads 10 Views 888KB Size
Letter pubs.acs.org/macroletters

Simulating Startup Shear of Entangled Polymer Melts Jing Cao* and Alexei E. Likhtman* Department of Mathematics and Statistics, University of Reading, Reading, U.K. RG6 6AX ABSTRACT: Start-up shear rheology is a standard experiment used for characterizing polymer flow, and to test various models of polymer dynamics. A rich phenomenology is developed for behavior of entangled monodisperse linear polymers in such tests, documenting shear stress overshoots as a function of shear rates and molecular weights. A tube theory does a reasonable qualitative job at describing these phenomena, although it involves several drastic approximations and the agreement can be fortuitous. Recently, Lu and co-workers published several papers [e.g., Lu et al. ACS Macro Lett. 2014, 3, 569−573] reporting results from molecular dynamics simulations of linear entangled polymers, which contradict both theory and experiment. On the basis of these observations, they made very serious conclusions about the tube theory, which seem to be premature. In this letter, we repeat simulations of Lu et al. and systematically show that neither their simulation results nor their comparison with theory is confirmed.

S

comparison of their stress overshoot results with the GLaMM version of the tube theory7 and found significant discrepancies, which is surprising because this theory was shown to agree well with many experimental data sets. 7,8 The controversial simulation results from Lu et al. have caused serious concern of the validation of the tube model for describing the rheological properties of entangled polymers.9−11 Later, in order to examine statements from Lu et al., Masubuchi and Watanabe10 performed primitive chain network simulation under shear flow with 1/τd < γ̇ < 1/τR, in which polymer chains are bundled by slip-links to form entangled system. They found that the stress and orientation both experience overshoot but not the chain stretching, which is in a good agreement with experiments12 and the tube theory. In addition, Watanabe et al.11 revisited the analysis of stress-optical law and demonstrated that the overshoot of chain stretching without the overshoot of orientation is contradicting to the stress-optical law, which is in accord to the results based on primitive chain network simulations. Thus, the need of examining the controversial simulations from Lu et al. was highlighted in the community, which is exactly the focus of this letter. In this work, we examine the results of Lu et al.3 based on the same molecular dynamics model. The polymer chains are represented by the conventional Kremer−Grest polymer model,6,13 which is identical to the one used by Lu et al. The equilibrium systems are first prepared by running simulations using soft nonbonded potential and harmonic bonded potential,14,15 so that the chains have the same mean square internal distance ⟨(Ri − Rj)2⟩/|i − j| as in the standard Kremer−

tress overshoot can be observed in entangled polymer melts during startup shear at sufficient higher shear rates. Such phenomenon can be well described by the tube model, introduced by de Gennes1 and Doi and Edwards,2 which accounts for topological constraints on the motion of a target chain arising from the uncrossability of surrounding chains. Under this tube picture, the motion of the chain perpendicular to the contour of the tube is suppressed and the chain follows Rouse dynamics along the tube. In the startup shear, with the applied shear rate γ̇ larger than the reciprocal of the terminal time 1/τd and smaller than the reciprocal of the Rouse time 1/τR, chain is supposed to have enough time to re-equilibrate inside its tube, and therefore, the overshoot is explained by changing tube segment orientation. In contrast, when the shear rate exceeds 1/ τR, the chain gets stretched inside the tube, which typically increases the stress level further. To highlight this, we will report our shear rates in terms of dimensionless Rouse Weissenberg number WiR = γ̇τR. From experiments, it is well-known that at WiR < 1 stress has a maximum at strain γ≈ 2, whereas as WiR increases above one, the strain at the stress maximum increases with increasing WiR. The recent simulation work from Lu et al.3−5 reported results contradicting to the tube picture mentioned above. They performed Brownian dynamics simulations of the Kremer− Grest model6 on startup shear at shear rates 1/τd < γ̇ < 1/τR in a mildly entangled system with about 10 entanglements per chain and observed shear stress overshoots. After that, they divided the chains into subchains containing Nblob beads (Nblob = 36) and calculated the chain stretching and chain orientation according to the conformation of the subchains, A clear overshoot was observed in the chain stretching instead of the orientation. Thus, Lu et al.3 argued that the chain stretching leads to the stress overshoot and not the orientation, After that, they argued that the stress-optical law is invalid as well. Lu et al. also attempted © XXXX American Chemical Society

Received: October 5, 2015 Accepted: November 19, 2015

1376

DOI: 10.1021/acsmacrolett.5b00708 ACS Macro Lett. 2015, 4, 1376−1381

Letter

ACS Macro Letters Grest model13 but are able to cross each other. After running the system with soft potentials for several Rouse time of the chains, the Lennard-Jones (LJ) and finitely extensible nonlinear elastic (FENE) potentials are switched on. The initial states of shear simulations are obtained by running the equilibrium simulation for multiple longest relaxation time of the chains. In order to simulate shear deformation, the Lees−Edwards boundary condition and the SLLOD algorithm were applied.16−18 The simulations were carried out using the MD package developed in our group, which has been proved to generate consistent results with other simulation works.6,13,18−20 We also note that the shear flow of Kremer−Grest model was studied extensively before by Kröger and co-workers.21,22 Before plunging into the investigation of the start-up shear simulations, we first analyze rheological properties of linear polymer melts in the equilibrium. This is a standard rheological procedure which the careful experimentalists use to check the consistency of the nonlinear data, and to estimate tube theory parameters of the system. According to the fluctuation− dissipation theorem, the stress response function after a small perturbation G(t) can be expressed by the correlation function of the stress in the equilibrium. An efficient method developed by Ramirez et al.20 called multiple-tau correlator is used to calculate the correlation functions of the stress. Then, the storage and loss modulus G′(ω) and G″(ω) can be easily obtained by fitting G(t) with the sum of exponential modes. These are shown by symbols in Figure 1 for flexible chains with two different chain lengths N = 350 and 512. Likhtman and McLeish23 developed a quantitative tube theory which combined self-consistent treatment for contour length fluctuations and constraint release with de Gennes−Doi− Edwards reptation theory. The simulation results from the systems with N = 512 and N = 350 are fitted with this theory simultaneously. The best-fits are shown by lines in Figure 1a with fitting tube parameters Ge = 0.0196, τe = 3290, Ne = 60, and cν = 0.1, where Ge is the entanglement modulus, τe is the Rouse relaxation time of an entangled strand, Ne is the number of monomers per entanglement and cν is the constraint release parameter, in which cν is fixed rather than fitted. We can see that the theory predicts G′ and G″ pretty well for the frequencies ω < 1/τe. Note that these parameters have large error-bars, and wider range of molecular weights is needed to make reliable estimates of tube theory parameters. As was argued in ref 19, the number of entanglements in MD chains are too small for the tube theory to work well. Thus, a better way for systematic modeling is to fit MD data with more detailed slip-springs model,24 and then compare it with the tube model at much larger number of entanglements not accessible in MD. Here however we perform a naive fit just for the purpose of roughly estimating of startup shear rheology. Note that parameters used by Lu et al. do not seem to match the data at all (Figure 1b), mainly due to the Ne being too small. Such large discrepancy of Ne used by Lu et al. and in this work is simply due to different physical observables and analysis algorithms applied in two works, which has been investigated by Pütz et al.25 that Ne values derived from monomer mean square displacement and plateau moduli have a factor of 2.3(2) difference. In addition, Everaers et al.26 carried out the primitive path analysis (PPA) on the standard Kremer−Grest model and obtained Ne value around 65 which is consistent to our result. In our simulations, the equilibrium system consists of 150 chains of length N = 512 with monomer number density ρ = 0.85σ−3 in a cubic box. As well-known, the polymer chains will be stretched and oriented under shear. Suppose the shear and

Figure 1. (a) G′ and G″ of linear monodisperse melts with different chain lengths from our simulations(symbols) and the best-fit based on LM2002 theory(lines) with tube parameters indicated in the figure. Inset contains stress relaxation function G(t). (b) Prediction of G′ and G″ based on LM2002 theory(lines) when using the tube parameters given by Lu et al. in ref 3. The symbols are the same simulation data as in part a.

velocity gradient directions are along the x- and y-axes, respectively. If the simulation box is not large enough, the size of the chains would exceed the length of the simulation box and the chains would interact with their own images through the periodic boundaries. In order to eliminate this possible defect, we first obtain configurations of the chains from an equilibrium simulation which has a cubic simulation box and then place another copy of this configuration next to the original one along the x-axis to construct the simulation box with the aspect ratio of 2:1:1. Because two parts of the box contain exactly the same conformations of the chains, we run the simulations for one terminal relaxation time τd of the component to let chains forget about the initial configurations. The Rouse time τR of the chain with N = 512 is determined from the tube parameters derived from the previous LVE theory such that τR = τe (N/Ne)2≈ 2.39× 105τ. However, in order to compare our results with Lu et al.’s data consistently, the Rouse time τR = 2.3× 105τ is used in the rest of the paper for indicating the shear rates through the Rouse−Weissenberg number WiR = γ̇τR. Note however that they used slightly different molecular weight N = 500, the difference of 2% is beyond our current accuracy. The parameters of the MD simulations used in Lu et al.’s work3 and this work are listed in Table. 1. Figure 2 shows the transient shear viscosity η+(t ) =

σxy(t ) γ̇

versus time t from the simulations in a log−log plot, in which different symbols represent different shear rates indicated by 1377

DOI: 10.1021/acsmacrolett.5b00708 ACS Macro Lett. 2015, 4, 1376−1381

Letter

ACS Macro Letters Table 1. Simulation Parameters Used in Lu et al.’s Work3 and This Work Lu et al. this work

N

Nc

τe/103τ

τR/105τ

Ne

Ge

ρσ3

box ratio

code

500 512

423 300

1.19 3.29

2.3 2.3

36 60

0.0295 0.0196

0.85 0.85

3:3:2 2:1:1

LAMMPS house made

Now, can we use standard model based on the tube picture to predict what we observed in the nonequilibrium MD simulations? The GLaMM theory7 is a tube theory that includes all main relaxation mechanisms in the tube model: reptation, contour length fluctuation, chain retraction and convective constraint release. It retains deformation information along the chain contour and so provides the most detailed stress predictions. Four input parameters τe, Ge, Ne and cν used in the GLaMM theory are obtained from the best-fit of G′ and G″ in the equilibrium simulations through LM2002 theory as indicated in Figure 1a. The predictions of the shear viscosity at different shear rates are shown by solid lines in Figure 2, as obtained from freely available software RepTate. In Figure 2, we show that the GLaMM theory predicts the overshoot of shear viscosity reasonably well for different shear rates. However, the predicted shear stresses in steady state are all slightly higher (20%) than the ones measured in the simulations. Note also small discrepancy between the dashed line and the theory predictions, which is due to additional approximations about constraint release in GLaMM theory as compared to the LM2002 theory. At this point, we do not dispute the fact that the tube model may be incomplete or requires serious modifications, but the disagreement between our simulation results and the GLaMM theory is clearly not as large as the one claimed by Lu et al.3−5 Since our system (N = 512) is similar to the one used by Lu et al. (N = 500)3 and the Rouse time τR of the chains used in both works are equal to 2.3× 105τ, it is convenient to compare our results to theirs, which is shown in Figure 3a, where σxy is plotted versus strain, instead of time, in a linear−linear plot. The blue filled symbols are extracted from Figure 1a in ref 3, and the red open symbols are the results from our simulations. The slowest

Figure 2. Shear viscosity η+(t) = σxy(t)/γ̇ from the simulations (symbols) and the predictions from the GLaMM theory (solid lines) with input parameters indicated in the figure, in which the blue dashed line is calculated using eq 1 for the linear viscoelastic regime.

Rouse−Weissenberg number WiR. The blue dashed line is obtained from the stress response function G(t) in the linear viscoelastic regime: η+(t ) =

∫0

t

G(t ′) dt ′

(1)

which is treated as an additional validation of the measurement of shear stress. We can see that, at small strains, the shear viscosity curves superimpose with the linear viscoelastic curve for all applied shear rates. If the shear rate γ̇ is faster than the reciprocal of terminal time τd of the system, the viscosity η+(t) passes through a maximum and then decreases until it reaches the steady state.

Figure 3. (a) Shear stress versus shear strain from Lu et al.(blue filled symbols) and this work(red open symbols) where N is the chain length and Nc is the number of chains in the systems. (b) Replotted shear viscosity versus time from part a. (c) Comparison of linear viscosity between our simulation (dashed line) and experiment from Auhl et al.(symbols) in mildly entangled systems. (d) Replot of linear viscosity in a well entangled system from Auhl et al. (e) Shear strain at stress maximum from experiments of different materials (open symbols) and this work (red disks). (f) Shear viscosity at steady state ηss as a function of shear rate γ̇ from Kröger et al.22 (open symbols), this work (red filled squares) and Lu et al.3 (blue filled triangles). 1378

DOI: 10.1021/acsmacrolett.5b00708 ACS Macro Lett. 2015, 4, 1376−1381

Letter

ACS Macro Letters

In ref 3, Lu et al. argued that the chain stretch induces the stress overshoot, even at the shear rate smaller than the reciprocal Rouse time 1/τR. Although we do not necessarily agree with their definition of chain stretch inside the tube through the blobs, we follow the same method for constructing the “blob primitive path” of the chain as Lu et al. The chain is divided into Z = N/ Nblob units with Nblob = 36 monomers in each unit. The “blob primitive path” of the chain is constructed by connecting the center-of-mass of the sequential blobs and the total contour length of the “blob primitive path” is calculated by the sum of the segment lengths between neighboring blobs. Figure 4 shows the ratio between the average contour length L and the equilibrium length L0 versus strain γ from our

shear rate used in our simulations is close to the fastest one used by Lu et al. and there is a factor of 6 difference between σmax xy at this specific shear rate. If we plot shear viscosity as a function of time in the same way as Figure 2 for both sets of data, the data from Lu et al. at WiR = 1/4 and 1/6 (on top of Figure 3b) somehow collapse onto each other even though the shear rates have a factor of 1.5 difference, which is usually not observed in experiment or theory. Furthermore, two sets of shear viscosity data show different slopes in the linear regime that the data of Lu et al. shows an initial slope of unity as indicated by solid line in Figure 3b, but our data has a slope considerably lower than unity. Does our data violate the basic principle of linear viscoelasticity? From eq 1, it can be seen that in order to obtain η+(t) as a linear function of time t (+1 power law), G(t) must have a clear plateau in this regime, which ranging from 105τ to 106τ according to Lu et al.’s data. However, as shown in the inset of Figure 1a, G(t) in the system with N = 512 decreases monotonically without a clear plateau. Thus, the initial slope of linear viscosity should be smaller than unity at times t > 104 τ in this mildly entangled system, just as reflected by the blue dashed line in Figure 3(b). By referring to eq 1, it is clear that the shear viscosity in linear regime can only follow a unity slope in well-entangled systems where G(t) has a clear plateau; otherwise, it has a slope less than one. In order to inspect this, it is straightforward to compare our MD shear results with experiments. Auhl et al.8 synthesized nearly monodisperse PIs with wide range of the number of entanglements and carried out start-up shear measurements for molecular weights up to 47 entanglements. Figure 3c shows the comparison of shear viscosity in linear regime between our simulation and the experiment (extracted from Figure 12d in ref 8) for the systems with similar Z values. After introducing shift factors of time (9× 104) and viscosity (3 × 10−3) to the experimental data, the two plots collapse onto each other extremely well and both experience smaller than 1 slope in the transient state. In addition, Figure 3d shows the linear viscosity of a well entangled system (Z ≈ 47) extracted from Figure 12f in ref 8, in which the LVE curve does have unity slope as expected. On the other hand, it is surprising that Lu et al. do manage to obtain a unity slope of the linear viscosity in such a mildly entangled system, which implies that the stress correlation function from their equilibrium simulations should show a clear plateau. Unfortunately, there is no G(t) data reported from Lu et al. that we can compare with. Figure 3e shows the shear strain at the stress maximum as a function of the shear rate from experiments8,12,27−29 together with our simulation data, in which all data points somehow overlap in this particular representation. Figure 3f shows the simulation results of steady-state shear viscosity ηss as a function of the shear rate for monodisperse linear chain systems represented by the Kremer−Grest model, in which open symbols are extracted from Figure 1a of ref 22, ranging from unentangled to mildly entangled systems. The model used by Kröger et al.22 is the same as the one used by Lu et al. and in this work, except that the monomer number density is 0.84σ−3 instead of 0.85σ−3. Such 1% difference of monomer number density would not affect the conclusions we give later. Our results are consistent with the ones reported from Kröger et al.22 that all data points collapse onto a master curve. However, since the simulations from Lu et al. did not reach steady states, ηss are estimated by the values of shear viscosity at the largest shear strain they reached(γ ≈ 5) and displayed by blue filled triangles. It seems that Lu et al.’s simulations overpredicted the shear stresses.

Figure 4. Renormalized contour length of the blob primitive chain as a function of strain from Lu et al. (blue filled symbols) and this work (red open symbols) under similar shear rates.

simulations (red open symbols) together with the one reported from Lu et al. (blue filled symbols), in which two sets of data have similar Rouse−Weissenberg numbers (WiR ≈ 0.5). As expected, our simulation result shows a very faint peak. If one would like to extract the maximum of the chain stretching, a horizontal line can be drawn somehow tangent to the data as a guideline which is displayed by horizontal black dashed line in Figure 4. Roughly speaking, the maximum of the chain stretching is reached around γ ≈ 3.8 with the value equal to 1.03. Under such shear rate (WiR ≈ 0.5), the “blob primitive path” of the chain is only stretched less than 3% comparing to its equilibrium value and the strain value at the stretching maximum (γ ≈ 3.8) is larger than the one at the stress maximum (γ ≈ 2) as shown in Figure 3a. Thus, the shear stress overshoot should not be due to the chain stretching. This again strongly contradicts the results of Lu et al. shown in Figure 4 by blue-filled squares. For flexible polymer chains, the stress and orientation tensors are proportional to each other in both linear and nonlinear regime, unless the chain is highly stretched and the finite extensibility becomes important, which is termed stress-optical law.11 For completeness, we examine the stress-optical law of the flexible Kremer−Grest model in both equilibrium and nonequilibrium simulations. In the equilibrium simulations, the relaxation function of total orientation tensor is supposed to be proportional to the total stress relaxation function S(t ) =

1379

Nb kBT

1 Nb 2

Nc

Nc

j=1

j=1

∑ Ojαβ(t ) ∑ Ojαβ(0) (2) DOI: 10.1021/acsmacrolett.5b00708 ACS Macro Lett. 2015, 4, 1376−1381

Letter

ACS Macro Letters G (t ) =

V ⟨σ αβ(t )σ αβ(0)⟩ kBT

Q = ρb α 2 (3)

1 G (t ) = S(t ) Q

(5)

where ρb is the number density of the bonds in the system. Using the Q value derived from the equilibrium systems, the stressoptical coefficient is found to be α = Q /ρb ≈ 0.32 for the flexible Kremer−Grest model. Figure 5b shows the comparison between shear viscosity and rescaled orientation for the systems with chain length N = 256 and different Rouse−Weissenberg number ranging from 0.58 to 57.5. We can see that the stressoptical law works extremely well for WiR ≤ 5.75, which contradicts to the results of Lu et al. In summary, we performed the molecular dynamics simulation of conventional Kremer−Grest model on startup shear to investigate the origin of the stress overshoot of linear polymer melts. The systems studied in this work were similar to the ones reported from Lu et al.;3−5 i.e., the systems are mildly entangled and the shear rates were set to be smaller than the reciprocal of the Rouse time and larger than the reciprocal of the terminal time. We found that the stress and orientation exhibited overshoot and the stress-optical law was validated consistently in both equilibrium and shear simulations, as noted experimentally. A tube-based GLaMM theory is able to predict the shear stress reasonable well comparing to the simulation results in a broad range of the shear rate. If the “primitive path” of the chain is defined in the same way as proposed by Lu et al.,3 i.e., by connecting center of mass of the subchains each containing 36 beads, no significant chain stretching is observed. Moreover, the strain at the maximum of the chain stretching is larger than the strain at the stress maximum, which demonstrates that at low rates the origin of the overshoot in shear stress should be due to the orientation of the chain instead of chain stretching. Our results are fully consistent with experimental results reported from different groups,8,12,27−29 molecular dynamics simulation results reported by Kröger et al.22 and primitive chain network simulations results reported by Masubuchi and Watanabe.10 In contrast, we found that results of Lu et al. for stress and chain stretch are not confirmed. Unfortunately, we are unable to comment on the source of problems with results of Lu et al. A more detailed study of startup shear in the systems with different chain lengths and stiffness is on the way.30

(4)

where S(t) is total orientation relaxation function, G(t) is the N−1 α β stress relaxation function, Oαβ j (t) = ∑i=1 uij (t)uij (t) is the orientation tensor of the chain j, uij is the bond vector connecting monomer i and i + 1 on chain j, σαβ(t) is the stress tensor, Nc and Nb are number of chains and bonds in the system, V is the volume of the system and Q in eq 4 is the coefficient related to the stressoptical coefficient. Figure 5a shows the stress and orientation relaxation functions of two monodisperse melts representing unentangled and mildly



AUTHOR INFORMATION

Corresponding Authors

*(A.E.L.) E-mail: [email protected]. *(J.C.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.

Figure 5. (a) Stress (lines) and orientation (symbols) relaxation functions in the equilibrium where Q = 0.0886. (b) Renormalized stress σxy/γ̇ (lines) and orientation Sxy/(γ̇α) (symbols) versus time at different Rouse−Weissenberg numbers.



ACKNOWLEDGMENTS We thank Zuowei Wang and Yuichi Masubuchi for fruitful discussions and support. This work was supported by the Engineering and Physical Sciences Research Council (EPSRC), Grant EP/K017683.

entangled systems. The coefficient is found to be Q = 0.0886 to collapse stress and orientation relaxation functions after certain characteristic time t*≈ 100 τ, which is consistent with the results reported earlier from our group.14 Under startup shear, the xy-component of the orientation 1 N tensor Sxy(t) is defined as N ∑ j =c 1 Ojxy (t ) such that if the stress-



DEDICATION J.C. dedicates this work to caring friend and esteemed adviser Prof. Alexei E. Likhtman who sadly passed away during the publication of the manuscript.

b



Sxy(t )

optical law is valid in the startup shear, σxy(t ) = α , the following relationship should be satisfied according to eqs 2, 3, and 4,

REFERENCES

(1) de Gennes, P. G. J. Chem. Phys. 1971, 55, 572.

1380

DOI: 10.1021/acsmacrolett.5b00708 ACS Macro Lett. 2015, 4, 1376−1381

Letter

ACS Macro Letters (2) Doi, M.; Edwards, S. The Theory of Polymer Dynamics; Oxford University Press: Oxford, U.K., 1988. (3) Lu, Y.; An, L.; Wang, S.-Q.; Wang, Z.-G. ACS Macro Lett. 2014, 3, 569−573. (4) Lu, Y.; An, L.; Wang, S.-Q.; Wang, Z.-G. Macromolecules 2014, 47, 5432−5435. (5) Lu, Y.; An, L.; Wang, S.-Q.; Wang, Z.-G. Macromolecules 2015, 48, 4164−4173. (6) Kremer, K.; Grest, G. J. Chem. Phys. 1990, 92, 5057−5086. (7) Graham, R. S.; Likhtman, A. E.; McLeish, T. C. B.; Milner, S. T. J. Rheol. 2003, 47, 1171. (8) Auhl, D.; Ramirez, J.; Likhtman, A. E.; Chambon, P.; Fernyhough, C. J. Rheol. 2008, 52, 801−835. (9) Graham, R. S.; Henry, E. P.; Olmsted, P. D. Macromolecules 2013, 46, 9849−9854. (10) Masubuchi, Y.; Watanabe, H. ACS Macro Lett. 2014, 3, 1183− 1186. (11) Watanabe, H.; Matsumiya, Y.; Inoue, T. NIHON REOROJI GAKKAISHI 2015, 43, 105−112. (12) Pearson, D. S.; Kiss, A. D.; Fetters, L. J.; Doi, M. J. Rheol. 1989, 33, 517−535. (13) Auhl, R.; Everaers, R.; Grest, G. S.; Kremer, K.; Plimpton, S. J. J. Chem. Phys. 2003, 119, 12718−12728. (14) Cao, J.; Likhtman, A. E. Phys. Rev. Lett. 2010, 104, 207801. (15) Likhtman, A. E. Viscoelasticity and Molecular Rheology. In Comprehensive Polymer Science, 2nd ed.; Elsevier: Amsterdam, 2011. (16) Evans, D. J.; Morriss, G. P. Comput. Phys. Rep. 1984, 1, 297−343. (17) Evans, D. J.; Morriss, G. P. Phys. Rev. A: At., Mol., Opt. Phys. 1984, 30, 1528−1530. (18) Cao, J.; Likhtman, A. E. Phys. Rev. Lett. 2012, 108, 028302. (19) Likhtman, A. E.; Sukumaran, S. K.; Ramirez, J. Macromolecules 2007, 40, 6748−6757. (20) Ramirez, J.; Sukumaran, S. K.; Vorselaars, B.; Likhtman, A. E. J. Chem. Phys. 2010, 133, 154103. (21) Kroger, M.; Loose, W.; Hess, S. J. Rheol. 1993, 37, 1057−1079. (22) Kroger, M.; Hess, S. Phys. Rev. Lett. 2000, 85, 1128−1131. (23) Likhtman, A. E.; McLeish, T. C. B. Macromolecules 2002, 35, 6332−6343. (24) Likhtman, A. E. Macromolecules 2005, 38, 6128−6139. (25) Putz, M.; Kremer, K.; Grest, G. S. Europhys. Lett. 2000, 49, 735. (26) Everaers, R.; Sukumaran, S. K.; Grest, G. S.; Svaneborg, C.; Sivasubramanian, A.; Kremer, K. Science 2004, 303, 823−826. (27) Menezes, E.; Graessley, W. J. Polym. Sci., Polym. Phys. Ed. 1982, 20, 1817−1833. (28) Kahvand, H. Ph.D. Thesis; Illinois Institute of Technology: Chicago, IL, 1995. (29) Schweizer, T.; van Meerveld, J.; Ottinger, H. J. Rheol. 2004, 48, 1345−1363. (30) Cao, J.; Likhtman, A. E. Manuscript in preparation.

1381

DOI: 10.1021/acsmacrolett.5b00708 ACS Macro Lett. 2015, 4, 1376−1381