Slip-Spring Model for the Linear and Nonlinear Viscoelastic Properties

Article pubs.acs.org/Macromolecules

Slip-Spring Model for the Linear and Nonlinear Viscoelastic Properties of Molten Polyethylene Derived from Atomistic Simulations A. P. Sgouros, G. Megariotis, and D. N. Theodorou* School of Chemical Engineering, National Technical University of Athens (NTUA),GR-15780 Athens, Greece S Supporting Information *

ABSTRACT: Atomistic simulations have been very useful for predicting the viscoelastic properties of polymers but face great difficulties in accessing the dynamics of dense, well entangled longchain melts with relaxation times longer than μs due to the high computational cost required. A plethora of coarse-grained models have been developed to address longer time scales. In this article we present a multiscale simulation strategy that bridges detailed molecular dynamics (MD) simulations to slip-spring based Brownian dynamics/kinetic Monte Carlo (BD/kMC) simulations of long-chain polymer melts. The BD/kMC simulations are based on a mesoscopic Helmholtz energy function incorporating bonded, slip-spring, and nonbonded interaction contributions (Macromolecules 2017, 50, 3004). Bonded contributions are expressed as sums of stretching and bending potentials of mean force derived from detailed MD simulations of shorter-chain melts, while nonbonded interaction contributions in the absence of slipsprings are derived from an equation of state that is consistent with thermodynamic properties predicted by detailed MD and measured experimentally. Monodisperse linear polyethylene melts of chain lengths C260 to C2080 are used as a test case. Estimates of the chain self-diffusivity, the longest relaxation time, the stress relaxation modulus, and the zero-shear viscosity from ms-long equilibrium BD/kMC simulations are in excellent agreement with MD results for the shorter-chain melts and with experiment. The BD/kMC scheme is extended to simulate Couette flow using Lees−Edwards periodic boundary conditions over a range of Weissenberg numbers (Wi) from 10−2 to 105. Predictions for the shear viscosity as a function of shear rate, the first and second normal stress difference coefficients, the startup shear stress, as well as for changes in chain conformation and entangled structure with increasing Wi are in favorable agreement with experimental and atomistic simulation evidence.

1. INTRODUCTION The study of rheological properties of polymer melts is of major importance for the scientific community and for industry.1 Simulations have enabled significant advances in understanding various heretofore elusive microscopic phenomena regarding structure−property relations, thermodynamics, and dynamics in polymer melts. MD simulations are employed routinely in studying polymer melts under quiescent2−8 and flow conditions.3,9,10 Usually, due to the enormous computational cost involved in performing such simulations, the examined systems are small and the applied shear rates are relatively large, while simulation times are on the order of hundreds of nanoseconds. The longest relaxation times of high molar mass polymers, however, are on the order of milliseconds up to several seconds and thus not easily amenable to MD simulations. This has motivated the development of coarsegrained (or mesoscopic) simulations that offer access to mesoscopic time and length scales in reasonable simulation time. The role of topological constraints (TCs) in long-chain polymer rheology has been recognized early. One of the most significant models in describing TCs is the tube model first developed by Edwards11 and later studied by de Gennes.12 This © XXXX American Chemical Society

model envisions the polymer chains as reptating within an artificial tube formed by entanglements with their neighboring chains; lateral motion of the chain against the tube is suppressed. The tube model has undergone several refinements over the years through the incorporation of important mechanisms such as constraint release and contour length fluctuations.13−15 An excellent review of tube models has been written by McLeish.16 In MD simulations the polymeric chains are unable to cross each other due to excluded volume effects, hence giving rise to TCs that form the artificial tube.17−20 The modeling of TCs in mesoscopic descriptions of coarse-grained model systems, on the other hand, is more complicated since in many cases much softer effective potentials are involved which do not prevent such chain crossings; thus, the predicted viscoelastic response of the melt becomes inaccurate. References 21 and 22 report some pioneering multibead− spring MD simulations which have been shown to predict very accurately several linear and nonlinear effects regarding the Received: April 3, 2017 Revised: May 11, 2017

A

DOI: 10.1021/acs.macromol.7b00694 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

approaching bead dimensions. The Helmholtz energy of chain strands is derived as a potential of mean force from detailed united-atom MD simulations of the polymer. In comparison to the early multibead−spring simulations,21,22 slip-spring approaches such as the one adopted here can be made more economical computationally by aggressive coarse-graining; that is, they enable the simulation of very long times for given computer time and system size. For example, in simulating a system with 520,000 carbon atoms (10,000 beads) with a time step of 50 ps, our code can cover ∼0.4 ms per day on five processors (Intel Xeon CPU E5645 @ 2.40 GHz); faster codes have been developed by other groups. The present work differs from earlier slip-spring simulation approaches in that (a) it develops the slip-spring model in parallel with detailed molecular simulations of a united atom (hereafter called “atomistic”) model of the polymer, thus introducing a strategy for mapping the latter model onto the former and establishing a hierarchical approach relying on both models; (b) it validates predictions through comparisons with both detailed atomistic simulations and experiments on molten PE; (c) it employs the mesoscopic model within nonequilibrium BD/kMC simulations of shear flow, employing Lees−Edwards periodic boundary conditions over a wide range of Weissenberg numbers, in order to predict shear thinning, shear stress evolution upon startup of shear flow, and the first and second normal stress difference coefficients. This article is organized as follows. Section 2 describes the atomistic and mesoscopic models invoked and explains the mapping between the two models and the treatment of entanglements at the mesoscopic level. Section 3 provides some details of the atomistic MD and mesoscopic BD/kMC simulations. Section 4 presents predictions for structure and rheological properties from equilibrium simulations of the quiescent melts, while section 5 presents predictions from nonequilibrium simulations of the sheared melts. Finally, section 6 summarizes the main findings of the work.

structure and the viscoelastic properties of polymer melts. In the models employed, the uncrossability effect arises naturally due to the incorporation of repulsive nonbonded interactions among the coarse grained particles, thus making chain crossings very improbable. However, the mapping to an uncrossable bead−spring model must certainly have a limit; upon making the coarse-grained beads bigger and bigger (i.e., increasing the number of actual monomers being represented by a bead), the potential of mean force describing their interactions becomes softer and softer and the uncrossability is lost, unless it is added by hand. In more recent coarse-grained simulations the effect of TCs has been described through several extensions of the tube model which are designed to capture the main mechanisms that are necessary to follow the evolution of the chains within the field of the tube. In the innovative work by Padding and Briels23,24 leading to the “Twentanglement” method, the coarse-grained particles constituting a chain are connected with interacting elastic bands that can give rise to TCs. The incorporation of such TCs was shown to restore the proper viscoelastic response in polyethylene melts. Another successful class of coarse-grained models are the single- and multichain slip-link models, in which entanglements are modeled as zero dimensional TCs and form a tube along the contour path of the chain.11,25−33 In such formulations, chains are crossable, and the effect of uncrossability is added “by hand” through the introduction of slip-links. Those sliplinks can be either stationary or mobile, while their creation/ destruction processes constitute dynamical phenomena that evolve though a reptation-like mechanism. In later years, another class of such models emerged, the “slip-spring models”.34 Slip-springs are entropic springs connecting the contours of neighboring chains, which are created and destroyed at chain ends and can slide along the contours of their chains, hence mimicking the mechanism of reptation. Both single- and multichain slip-spring and slip-link models have been used extensively in the literature to model the viscoelastic properties of polymeric systems under equilibrium and nonequilibrium conditions.35−45 Excellent reviews of coarse-grained simulations of entangled polymer melts and of slip-link and slip-spring models in particular have been written by Masubuchi41 and by Schieber,33 while more recent work is reviewed in ref 35. In this work a mixed particle- and field-based mesoscopic BD/kMC simulation methodology35 utilizing slip-springs is further developed and applied to monodisperse linear polyethylene (PE) melts of chain lengths C260 to C2080 to predict linear and nonlinear rheological properties. As in ref 35, the dynamics is derived from a mesoscopic Helmholtz energy function which incorporates contributions from the elasticity of chain strands, slip-springs, and nonbonded interactions. In ref 35 the nonbonded interactions were expressed as a spatial integral of a Helmholtz energy density functional depending on the local density and derived from an equation of state that is consistent with the equilibrium volumetric properties of the melt. In the current work, however, the nonbonded interactions are described with the excess Helmholtz energy density functional, the “communal” entropy associated with exchange of beads from one cell to another being taken care of by the simulation. This formulation involving the excess Helmholtz energy is more correct theoretically; practically, it has negligible consequences for cells such as the ones used here that occupy large regions in space, but would be substantial for cells

2. MODEL SYSTEMS 2.1. Atomistic Model. The atomistic simulations were conducted using the TraPPE force field,46,47 wherein the hydrogen atoms are lumped onto the carbon atoms to which they are connected to form “united atoms”. CH2 and CH3 united atoms are treated as single interaction sites. The intramolecular interactions between united atoms that are separated by more than three bonds, as well as intermolecular interactions are described with the Lennard-Jones potential: ⎡⎛ σ ⎞12 ⎛ σ ⎞6 ⎤ VLJ(r ) = 4ε⎢⎜ ⎟ − ⎜ ⎟ ⎥ ⎝r⎠ ⎦ ⎣⎝ r ⎠

(1)

where ε = 0.0914 kcal/mol, σ = 3.95 Å. The cutoff distance of the potential was set to 9.085 Å, while the attractive tail contributions were computed by direct integration assuming isotropy.48 Even in the presence of strong chain orientation, such schemes provide accurate corrections to the potential energy and pressure if applied in conjunction with large cutoffs, where the radial atomic pair distribution function in directions both parallel and perpendicular to the oriented chains has reached unity. The interactions between neighboring atoms were described with a harmonic potential with spring constant kb = 191.7652 kcal/(mol Å2) and equilibrium length req = 1.54 Å: B


Article

Macromolecules

intramolecular springs) which model the mainly entropic elasticity of the polymer. Effective bending potentials depending on the angles formed by consecutive strands may be introduced to modulate single chain stiffness. There are two types of beads: internal nodal points and end points. The effect of the entanglements is modeled with slip-springs.35 More details regarding the slip-springs are provided in section 2.4. Our mesoscopic model systems are described thermodynamically by a Helmholtz energy function: A total = A strands + A angles + A nonbonded + Aentanglements

Astrands is the free energy of the strands that connect consecutive beads along the same chain and Aangles is a bond angle bending free energy. The functional expressions and parameters for those free energy effective potentials (intramolecular springs and angles) were extracted directly from the atomistic simulations using the mapping strategy discussed in the next section 2.3. Anonbonded is the excess free energy potential, relative to an ideal gas of beads, and models the nonbonded interactions. It is calculated from an equation of state (EoS); in this work, the Sanchez−Lacombe EoS50 is used for this purpose, with parameters for linear PE from ref 51. An advantage of this scheme is that it reproduces several thermodynamic properties of the system, such as the experimental isothermal compressibility.35,52 Details regarding the mesoscopic representation of nonbonded interactions can be found elsewhere.35 Aentanglements is the contribution of the entanglements to the total free energy. Ideally, Aentanglements should be very small, otherwise substantial perturbations to the chain conformations may occur. Recently, Chappa and collaborators proposed a solution to this problem wherein the effect of the slip-springs (SS) on chain conformations can be fully negated by using a compensating potential such that38

Figure 1. Mesoscopic representation of entangled polyethylene. Beads represent segments with a length of several Kuhn steps. Beads that belong to the same chain are connected by permanent strands (blue springs). The mesoscopic system presented herein cannot form “real” entanglements due to the softness of the effective nonbonded potential, so the entanglements are introduced via slip-springs (thin red springs). Each end of a slip-spring resides on a bead and can hop between adjacent beads along a chain. In the figure, beads with the same color belong to the same chain, while beads with darker shades denote the chain ends.

Zslip ‐ springs + compensating = Zno slip ‐ springs

⎛ A ss(rij) ⎞ ⎟⎟ Acompensating = zactivkBT ∑ exp⎜⎜ − ⎝ kBT ⎠ i>j

where zactiv is the activity (which is equal to e with μ being chemical potential of a slip-spring and tunes the number of slipsprings); i, j are indices that run over all possible pairs of beads that are able to form a slip-spring; and Ass(rij) is the free energy of a slip-spring with vector rij. Indeed, several test runs have been performed under equilibrium conditions with and without nonbonded interactions; chain conformations in all studied systems in the presence of the compensating potential remained unperturbed regardless of the activity and the strength of the slip-spring interactions. Using the compensating potential, a zero contribution to the pressure from the entanglements was observed, while the chain end-to-end distance and radius of gyration were unaffected. Unfortunately, a practical problem with using a compensating potential is that it requires the computation of a large number of pairwise interactions. Even with the use of neighbor lists, the cost of this computation increases as r3 with r being the maximum extension of the slip-springs. When relatively long slip-springs are allowed with lengths that are commensurate to the tube diameter d of real polymeric materials (d ≥ 40 Å) these computations become impractical.

(2)

(3)

Each dihedral angle φ was described by the Toxvaerd49 torsional potential: 8

Vtorsion(φ) =

∑ cn cosn(φ) n=0

(7)

μ/kBT

The bond angle bending potential of eq 3 was used with ka = 124.2003 kcal/(mol rad2) and θeq = 114°. Vangle(θ ) = 0.5ka(θ − θeq)2

(6)

with Z being the configurational partition function of the system. The compensating potential is given by the following expression:

Figure 2. From the atomistic to the mesoscopic representation. (a) Mapping of a real polymer chain with 260 carbon atoms onto a sequence of five coarse-grained beads (n = 52). (b) Two approaches for mapping atomistic polymer segments onto beads. A bead can be assigned either to the center of mass of a segment (top, rcm,n coordinates) or to the coordinates of the central atom(s) of a segment (bottom, rmid,n coordinates).

Vb(r ) = 0.5k b(r − req)2

(5)

(4)

where c0 = 1.9892, c1 = 4.2328, c2 = −0.6021, c3 = −7.1778, c4 = 4.4255, c5 = 3.9068, c6 = −8.9206, c7 = −3.4498, and c8 = 5.5980 in kcal/mol. 2.2. Mesoscopic Model. In the mesoscopic representation of polyethylene (Figure 1), segments of n carbon atoms are lumped into coarse-grained beads. Consecutive beads in the same chain are connected permanently by strands (i.e., C


Article

Macromolecules

Figure 3. (a, b) Probability density ρangle(cos θ) of cosines of the angles formed by three consecutive beads along the same chain and (c, d) the probability density of the distance between neighboring beads ρstrand(r) for segments with various lengths n, in carbon atoms, as computed directly from atomistic MD simulations of C260H522 melts. In parts a and c, the coordinates of the beads follow the center of mass, while in parts b and d, they follow the atom(s) in the middle of the represented segment. The markers represent beads with n = 10 (circles), n = 26 (squares), and n = 52 (diamonds) carbon atoms. The dashed line in parts a and b represents the freely jointed chain distribution ρangle(cos θn→∞) = 0.5.

Table 1. Mean Strand Length ⟨rstrand⟩, End-to-End Distance ⟨rchain,ete⟩ and the Corresponding Mean Square Values from Atomistic MD Simulations, Mesoscopic BD and BD/kMC Simulations, and from the Analytical Expressions of the Mesoscopic Model Obtained from the Mapping Procedurea MD, rcm,52 MD, rmid,52 BD, freely jointed chains BD, semiflexible chains BD/kMC, semiflexible chains with slip-springs theory, inverse Langevin

⟨rstrand⟩

⟨rstrand2⟩

⟨rchain, ete⟩

⟨rchain, ete2⟩

24.5 28.7 28.6 28.7 28.7

670 921 923 925 925

58.0 58.9 56.3 59.1 59.1

3904 4040 3695 4038 4044

28.6

923

−

point of view, however, with our choice of parameters the computed viscoelastic properties were the same, whether a compensating potential was used or not. This is because in this work, as well as in a previous study of cis-1,4 polyisoprene,35 the slip-springs employed were soft and long (r ≈ d), so the perturbations of chain conformations due to slip-springs (which were partially compensated by the nonbonded free energy potential) were minimal. For this reason, and since we checked that the viscoelastic properties remain unaffected whether a compensating potential was used or not, a compensating potential was not employed in generating the results presented in this work. However, it should be noted that the use of a compensating potential is mandatory in systems with high densities of strong slip-springs. Some results and discussion regarding the insensitivity of the viscosity to the compensating potential are included in section I of the Supporting Information to the present paper. 2.3. Mapping between the Atomistic and the Mesoscopic Representations. The quantities extracted from MD simulations were used as parameters in the mesoscopic BD/kMC simulations cast exclusively in terms of the beads, each bead representing a segment consisting of n consecutive monomers (methylene or methyl groups). Upon tracking the trajectories of the beads in the course of atomistic equilibrium MD simulations, we extracted the expressions and the parameters for the free energy effective potentials that describe the interactions among the coarse-grained particles in BD, thus allowing the replication of the conformations between

−

In all cases, each chain consists of five beads with size n = 52 methylene or methyl groups.

a

An empirical method to speed up the compensating potential calculations is to compute the pairwise interactions in eq 7 up to a length that is computationally feasible, and then account for the rest by multiplying the compensating potential with a scale factor such that the contribution of the entanglements to the total pressure is zero. Indeed, we found that by using this approach in some test cases the chains remained unperturbed, even when the density of the entanglements was set really high, and in the absence of the nonbonded interaction scheme employed here. The implementation of this approximation and some test results are discussed in Appendix A. From a practical D


Article

Macromolecules

Figure 5. Slip-spring length distribution from eq 22 (red line) and BD simulations (circles). For comparison, a length distribution from Gaussian springs is displayed (dashed line) with a mean squared slipspring length equal to 1600 Å2.

atom(s) of the represented segments instead of their center of mass (see Figure 2b, bottom). From now on, the coordinates/ representations of the n-sized segment’s center of mass will be denoted as rcm,n and those of the central atom(s) as rmid,n. Regarding the size of the coarse-grained particles, many different approaches have been followed in the literature. In relevant works that deal with coarse-graining,3,20,24,53,54 the beads were chosen to represent segments with n = 10 up to n = 50 backbone carbon atoms. In the production phase of this work, we chose to map chain segments onto mesoscopic beads with size equal to n = 52 carbon atoms (∼5 Kuhn segments for PE) from C260H522 melts at T = 450 K. This is done in a uniform fashion (the first n segments form one bead, the next n segments a second bead, etc., as shown in Figure 2a) and requires that the polymerization degree be a multiple of n. The chosen value of n is large enough to allow for large integration steps in the mesoscopic dynamics, while the size of the represented segments is less than the typical tube diameter of polyethylene (35−45 Å)4 in the “reptation picture”.24 Most recently, the distribution of the cosine of the angles (defined between three consecutive beads as shown in Figure 2a) has been assumed to be uniform (represented by the dashed line in Figure 3, parts a and b), corresponding to a random relative direction of successive strands; thus, no effective bending potentials have been employed at the mesoscopic level.20,24,53,54 Interestingly, analysis of the trajectories obtained by direct mapping of the atomistic MD simulations to the coarse-grained level shows otherwise. In particular, considerable skewness was observed in the angle distribution for a wide range of bead sizes. Figure 3a displays the angle distribution ρangles(cos θcm,n) from rcm,n coordinates. In cases where the bead size is small (e.g., n = 10) the distribution is skewed toward large angles due to the rigidity of segments with lengths commensurate with the Kuhn length of polyethylene, which is about 14.5 Å using the experimental characteristic ratio from neutron diffraction C∞ = 7.8.55,56 However, even when the bead size is considerably larger, the effect is still apparent. On the other hand, it was observed that the angle distribution extracted from rmid,n coordinates is much less skewed, as shown in Figure 3b. It should be noted that similar distributions were obtained from MD simulations on

Figure 4. (a) Strand length distribution ρstrand, (b) the free energy effective potential of the strands Astrand, and (c) the force Fstrand from rmid,52 coordinates. Results from atomistic simulations are depicted with circles. Theoretical calculations from Gaussian springs with variance ⟨r2⟩ = 921 Å2 are shown with dashed lines. Fittings with the inverse Langevin function on the MD data are depicted with dotted red lines, while fittings with the “soft” inverse Langevin function with a tolerance length rtol = 53.9 Å are shown with continuous red lines. In part a, the strand length distribution from equilibrium BD simulations with the bead−spring model is shown with squares, while results from BD/kMC simulations with the semiflexible bead−spring model in the presence of slip-springs are shown with blue ×-markers.

the two representations. This was done by performing Boltzmann inversion of the atomistically calculated distribution functions ρ(r), thus obtaining potentials of mean force A(r) that could adequately describe these effective interactions up to an additive constant: A(r) = −kBT ln ρ(r)

(8) 20,24,53,54

In many relevant works, the coarse-graining of atomistic trajectories is performed by tracking the center of mass of the represented segments (Figure 2b, top). In this work, another coarse-graining option was examined as well, wherein the bead coordinates are set to track the central E


Article

Macromolecules both C260H522 and C520H1042 melts considering angles that are located near the middle or the ends of the polymer chains (see Figure S2 of the Supporting Information). Parts c and d of Figure 3 present the strand length distributions ρstrand(rn) (i.e., of the distance between two consecutive beads) for beads with various sizes n, placed at the rcm,n and rmid,n coordinates, respectively. Even though the distributions from those two representations look qualitatively similar, they present significant differences, especially for larger values of n. For example, if we examine the mean and the mean squared value for bead sizes equal to n = 52 we find that for rcm,52 coordinates ⟨rcm,52⟩ = 24.5 Å and ⟨rcm,522⟩ = 670 Å2, while for rmid,52 coordinates ⟨rmid,52⟩ = 28.6 Å and ⟨rcm,522⟩ = 921 Å2 (see Table 1). In the following, the mapping will be performed from the rmid representation, since the extracted strand lengths are in close match with the average end-to-end distance of the fine-grained chain segments and the angle distribution shows little bias compared to the one obtained from the rcm representation. The circles in Figure 4a depict the strand length distribution ρstrand(rmid,52) from MD simulations. From this distribution we extracted the free energy, up to an additive constant, and forces (see circles in Figure 4, parts b and c, respectively) by using the following equations: A strand (r) = −kBT ln

L−1(x) ≈ L301−1(x) x

=

4πr 2

3 r2 kBT 2 2 ⟨r ⟩

U301(x) =

(13)

∫ dx L301−1(x)

= −ln(1 − x 2) +

x2 x4 x6 − − , 2 20 15 (14)

Δ = 0.11%

were Δ is the maximum relative error. From eqs 13 and 14 we have derived analytical expressions for the free energy and the forces of the strands: −1 (xext) λ L Fstrand = − r 301 xext 3

xext

A strand = −

(10)

(15)

⎧ r , r = |r| ≤ rtol ⎪ ⎪ rmax =⎨ ⎪ rtol , r = |r| ≥ r tol ⎪r ⎩ max

(9)

∫

(16)

⎧ λ 2 rmaxU301(xext), r ≤ rtol ⎪ ⎪ 3 Fstrand dr=⎨ −1 ⎪ λ 2 L301(xext) + Bs , r ≥ rtol r ⎪6 xext ⎩ (17)

Bs =

⎡ ⎤ λ 1 rmax 2⎢U301(xext) − xextL301−1(xext)⎥ ⎣ ⎦ 3 2

(18)

where rmax is the maximum length of the strand, xext is the relative extension ratio of the strand, and λ is the slope of the force with respect to distance in the vicinity of small strand extensions. The free parameters rmax and λ were fitted to 55 Å and 7.2 × 10−3 kJ mol−1 Å−2, respectively, so that the mean and the variance of the strand length distribution from MD simulations are reproduced in the mesoscopic description. Brownian dynamics is driven by systematic and stochastic forces, so there is always a chance (albeit tiny if a proper time step is chosen) that the strand may overextend and generate very large forces. For this reason we introduced a tolerance length rtol which is used to soften the potential near the singularity at xext = 1, while Bs is a constant that ensures continuity of the force at |r| = rtol. A value of rtol equal to 53.9 Å (corresponding to rtol/rmax = 98%) was adequate for our purposes, although nearby values provided identical distributions. The red lines in Figure 4 depict the analytical calculations using eqs 9, 15 and 17 while the square markers in Figure 4a depict the strand length distributions obtained from BD simulations with the bead−spring model discussed above. The mean and the mean squared values of the strand length distribution from the mesoscopic model expressions and from the BD simulations coincide with those from MD simulations (see Table 1). Finally, a weak free energy bending potential for our semiflexible bead spring model was obtained by performing Boltzmann inversion on the angle distributions from the molecular dynamics simulations. The free energy bending

(11)

In eq 11, the only free parameter is the variance of the strand length distribution, a quantity that is readily available from MD simulations. A disadvantage of Gaussian springs, however, is that they fail to describe properly the behavior of strands at large extensions (as shown with dashed lines in Figure 4, parts a−c). This is especially important for the nonequilibrium BD/ kMC simulations employed here to study steady shear flow; hence finitely extensible strands were used in this work. The forces due to strands were described by an inverse Langevin function. The Langevin function (see eq 12) does not have an analytical inverse, so one has to perform either numerical inversion or use analytical approximations.57,58 1 x

,

Δ = 0.28%

In the definitions of eqs 9 and 10 the interbead distance vector r between two beads i and j is defined as r = ri − rj, while Fstrand(r) is the force on bead i due to bead j. As seen in parts b and c of Figure 4, the extracted free energy and magnitude of the force from the rmid,52 coordinates ascend abruptly at low distances. This is attributed to excluded volume effects between the midsections of the segments; this is not observed when mapping from rcm,52 coordinates. In our earlier work35 the interactions between consecutive beads in a chain were modeled with Gaussian entropic springs, where the free energy was computed by the Gaussian expression:

L(x) = coth(x) −

1 − x2

ρstrand (r )

Fstrand(r) = −∇r A strand

A strand (r) =

3x − 5 (6x 2 + x 4 − 2x 6)

(12)

Here we chose to work with an accurate and efficient approximant L301−1 and its integral U301 developed recently by M. Kröger:58 F


Article

Macromolecules

distribution is replicated with very good accuracy in BD simulations. For lengths larger than 80 Å the probability density is very low, hence the chosen value for the length tolerance does not have any appreciable impact on the length distribution. CReTA analysis of atomistically sampled melt configurations has uncovered a rather complex picture regarding the distribution, the mobility, and the lifetimes of TCs.17−20 In particular, successive strands between topological constraints were not uncorrelated directionally, while TCs of various strengths were detected: strong constraints which persist for a long time, and weak binary contacts which “blink” over times much shorter than are needed for the primitive path to change appreciably. In this work, we opted to adopt a simplified slipspring model as used by other researchers, where TCs can pass through each other and are only created and destroyed at chain ends. The hopping frequency vhop governs the mobility and the lifetimes of the slip-springs, hence it dictates the long time dynamics of the system. Here it was treated as an adjustable parameter, whereby the complex behavior of the entanglements and the effects related to the decreased number of degrees of freedom due to the coarse graining are described collectively. In particular, vhop was estimated through trial and error to vhop = 5 × 105 s−1 upon following the procedure below: (1) An upper bound for vhop was obtained through the rough argument presented in Appendix D of ref 35. (2) The self-diffusion coefficient D was calculated using several values of vhop for at least three systems with different molar masses (e.g., chains with 5, 10, and 20 beads). (3) A thorough comparison was performed of the dependence of D on M as obtained from BD/kMC and as obtained from MD simulations and the best value, or an interpolated intermediate value, was picked for vhop. The final parameter of our model is the number of slipsprings which is tuned through the activity. The number of entangled strands along each chain, zE̅ S, is estimated through the following formula:

potential that was extracted from the rmid,52 coordinates by using the expression Aangle(θ) = −kBT ln(ρangle(cos θ)), was fitted with a harmonic potential of the form: A angle = ka(θ − π )2

(19)

where ka was fitted to 0.2 kJ mol−1 rad−2. This parameter reproduces the effective angle distribution from atomistic MD (see Figure 3b, red diamonds) very well. The angle distributions of MD and BD are depicted in Figure S3 of the Supporting Information to the present paper. 2.4. Modeling of the Entanglements. In the current work, the effect of entanglements is modeled with slip-springs using a modified version of the scheme developed in ref 35. This model is governed by three fundamental parameters: the length distribution of the springs, the pre-exponential factor vhop which dictates the dynamics of the slip-springs, and the number of slip-springs or, in a grand canonical ensemble formulation, their activity zactiv = exp(βμ). Since in this work we dealt with strong shear flows, the free energy of the slip-springs was modeled using a soft FENE59 force field instead of a Gaussian one35 following the footsteps of ref 38: Fss = −r

λss Css(r )

⎧ r2 ⎪1 − , r ≤ rtol,ss rmax ,ss 2 ⎪ Css(r ) = ⎨ ⎪ rtol,ss 2 − 1 , r ≥ rtol,ss ⎪ rmax ,ss 2 ⎩

A ss = −

∫

(20)

(21)

2 ⎧−0.5λ r ss max ,ss ln[Css(r )], r ≤ rtol,ss ⎪ ⎪ Fss dr=⎨ λssr 2 ⎪ 0.5 + Bss , r ≥ rtol,ss ⎪ C (r ) ss tol,ss ⎩

(22)

⎡ 1 − Css(rtol,ss) ⎤ Bss = −0.5λssrmax ,ss 2⎢ + ln Css(rtol,ss)⎥ ⎢⎣ Css(rtol,ss) ⎥⎦

z ES ̅ = (23)

Nmonomers/chain N̅ES

(24)

where N̅ ES is the “the ensemble average of the number of monomers” in a strand between topological constraints.18 Interestingly, CReTA predicts up to ∼2.5 more topological constraints than estimated via the experimental molar mass between entanglements Me, which equals 1150 g/mol at T = 443 K for polyethylene.60 In particular, according to analysis with CReTA, N̅ ES ≈ 0.4Ne where Ne is the number of monomers in a primitive path Kuhn segment. Ne from CReTA is estimated as ∼75, in very good agreement with the experimental Ne = 82 at 443 K.17,61 For comparison the ensemble average of the entanglement spacing through analysis of PE melts with the Z-code61 was computed as 61 ± 6, for chain lengths larger than the characteristic value which is about 200 skeletal carbon atoms.62 Brownian dynamics simulations conducted in the past reinforce the prediction that N̅ ES < Ne in order to reproduce the experimental plateau modulus of polybutadiene solutions, Masubuchi et al. used half the experimental value of Ne.63 In the current study, we wanted to test this prediction with the current mesoscopic model, so we performed some initial

where rmax,ss is the maximum extension of the spring, λss is the slope of the force versus extension curve at small extensions and Css is the argument of the logarithm in the FENE potential. rmax,ss and λss are the free parameters of our model; they were set to 100 Å and 5.15 × 10−3 kJ mol−1 Å−2, respectively, in order to achieve a slip-spring length distribution with second moment equal to ⟨d2⟩ = 1600 Å2, where d is the tube diameter of high molar mass polyethylene (see Figure 5a).4 Furthermore, rmax,ss and λss were chosen so that the shape of the FENE length distribution is close to the shape of the length distribution of Gaussian slip-springs with the same variance (see Figure 5a, dashed line). In the same spirit as with the inverse Langevin springs, we introduced a tolerance length rtol,ss in order to avoid the generation of high forces due to overextension of the springs, and a constant Bss that ensures continuity of the forces at r = |r| = rtol,ss. The value of rtol,ss was set to 95 Å (so that rtol,ss/rmax,ss = 95%), since at these lengths the magnitude of the force ascends steeply. As seen in Figure 5, the theoretical FENE length G


Article

Macromolecules Table 2. Parameters and Statistics of the Mesoscopic Simulations system

Nbeads

zentang/chain

zactiv

τ1 (ns)

⟨Rg2⟩ (Å2)

⟨Rete2⟩ (Å)

⟨Rete2⟩/⟨Rg2⟩

C∞

C260H522 C520H1042 C1040H2082 C2080H4162

5 10 20 40

4.76 9.51 19.03 38.05

0.004 87 0.005 52 0.005 93 0.006 01

191 1163 12525 171421

785 1672 3433 7248

4050 9396 20443 42932

5.16 5.62 5.96 5.92

7.91 8.40 8.71 8.93

simulations with the experimental number of Me. Interestingly, the plateau modulus from the stress relaxation function turned out to be about half the experimental one when using the experimental Me. Subsequently, we followed the predictions of CReTA and experimented with a larger number of entangle2 ments, using NES ̅ ≈ 3 Ne (×1.5 more entanglements), hence, obtaining a better prediction of the plateau modulus and stronger scaling of the viscosity with molar mass. During our test experiments, we observed that the plateau modulus is roughly proportional to the number of slip-springs in the regime of large molar masses. We chose to work with N̅ ES = 54.7 carbon atoms, although an even lower value could have been used. If the number of entangled strands per chain is known, then the average number of slip-springs is obtained from the following relation:

z̅ − 1 ⟨nss⟩ = Nchains ES 2

Figure 6. Schematic representation of an elementary hopping event. The ends of the slip-spring (red) can only hop to the first neighbors ai‑1, ai+1, and bi‑1 since the distance between beads ai and bi+1 is larger than the attempt radius.

(25)

the box were adjusted by trial and error so that the time averaged pressure of the system was close to zero (±10 atm). Initially, the melts were thoroughly equilibrated via connectivity-altering Monte Carlo simulations,70−72 which allow vigorous sampling of the long length scale conformational properties of chains in the dense melt. The equilibrated PE melts were then subjected to molecular dynamics (MD) simulations for up to 375 ns using the LAMMPS package.73,74 The MD simulations were performed at constant temperature and volume, employing a Nosé−Hoover thermostat with an effective relaxation time of 0.1 ps to sample a canonical (NVT) ensemble.75,76 The equations of motion were integrated with the velocity-Verlet algorithm and the RESPA multitime scale integrator with two hierarchical levels.74,77,78 The forces due to bonds, bond angle bending and torsional force fields were computed every 1 fs (inner level), while the forces due to nonbonded interactions were computed every 4 fs (outer level). 3.2. Brownian Dynamics/Kinetic Monte Carlo. In our mesoscopic simulations we studied monodisperse polyethylene melts at T = 450 K under equilibrium (EBD/kMC) and steady shear flow (NEBD/kMC). Four systems were examined with different chain lengths (see Table 2) in order to probe the dependence of the structural and dynamical properties on the molar mass. The systems with the lowest molar masses were the same ones as those examined with MD simulations (C260H522 and C520H1042) and were used to perform the mapping between the atomistic and mesoscopic simulations. The initial configurations were extracted directly from atomistic representations generated with the MAPS package69 that were coarse-grained with the procedure described in section 2.3. In all cases, the dimensions of the box were set to be at least four times as large as the root-mean-square radius of gyration of the examined polyethylene chains. In BD simulations79,80 the evolution of the trajectory of each bead i is governed by the following stochastic differential equation:

In our simulations the number of slip-springs fluctuates as a result of entanglement destruction and creation events at chain ends. In other words, the system is grand canonical with respect to slip-springs, and the average density of slip-springs is controlled through the activity zactiv. Appendix B presents two methodologies for the thermodynamic and dynamical estimation of zactiv, for a system with a constant number of slipsprings. The activity has been calculated by both approaches and was found to be practically identical. Finally, section II of the Supporting Information to this paper contains extensive results regarding the distribution of the entanglements, which are compared thoroughly with the atomistic simulation-based predictions of CReTA.17

3. SIMULATION DETAILS 3.1. Molecular Dynamics. In our MD simulations three systems of linear monodisperse polyethylene melts were considered with different molar mass per chain, at a constant temperature T = 450 K and zero average pressure.64 First we examined two systems in the unentangled regime with short chains (72 × C16H34 and 32 × C52H106, where the first number in the multiplication indicates the number of parent chains present in the model system) in order to check the validity of our model by comparing our results to the literature.2,4−6,65 In the subsequent simulations, systems with larger molar masses were examined (53 × C260H522 and 37 × C520H1042) to probe their behavior in the entangled polymer regime and obtain useful information that was later used for our mesoscopic simulations. Initial configurations for the polyethylene melts were generated with the Amorphous Builder program66,67 based on an extension of the method of Theodorou and Suter68 available in the MAPS software.69 The dimensions of the simulation boxes were chosen to be at least twice as large as the chain’s root-mean-square radius of gyration in order to minimize system size effects. Furthermore, the dimensions of H


Article

Macromolecules ζi

dri(t ) = Fi({ri(t )}) + ξi(t ) dt

Monte Carlo (kMC) simulation running in parallel with the Brownian dynamics, using a modified version of the fluctuating slip-spring scheme developed in ref 35 (see Appendix C for more details). According to this scheme, the probability of an end of the ith slip-spring to hop to a new candidate site equals:

(26)

The systematic forces Fi were obtained as the negative of the gradient of the free energy in eq 5 with respect to ri, while the stochastic forces ξi are random forces with zero mean which obey the fluctuation−dissipation theorem. The friction coefficient ζi of bead i with mass mbead is calculated using the following relation: mbead ζi = ζmonomer mmonomer (27)

⎛ A (r ) ⎞ Phop, i = v hop exp⎜ ss i ⎟ΔtkMC ⎝ kBT ⎠

where Ass(ri) is the free energy of the ith slip-spring with length ri and ΔtkMC is the elapsed time between hopping events. Even though the value ΔtkMC does not affect the dynamics of the slipsprings, ΔtkMC should be large enough to allow for efficient simulations and low enough so that Phop,i does not exceed unity. If Phop,i becomes greater than 1, then some elementary events are missed and the dynamics of slip-springs is retarded. Here, ΔtkMC was set to 10 time steps of the BD; values which were by a factor of 10 smaller provided identical results. The process of hopping of slip-spring ends along their respective chains requires some care to preserve microscopic reversibility. In this work, since we dealt with finitely extensible springs, we consider as hopping candidates the first neighbors of the slip-spring ends under the condition that the length of the candidate slip-spring does not surpass an attempt radius rattempt. In theory, rattempt should be equal to the maximum extension of FENE slip-springs, since slip-springs with larger lengths constitute impossible states with infinite energy. In practice, in order to avoid the formation of extended slipsprings that generate huge forces, the attempt radius was set equal to the tolerance length rtol,ss discussed above, since at these lengths the probability distribution tends to zero and the free energy rises steeply. For clarity, Figure 6 illustrates an example of an elementary hopping event. The end of the slipspring that resides on bead ai can hop either to bead ai‑1 or ai+1, while the end located at bead bi can only hop to bead bi‑1. Bead bi+1 is excluded from the candidates list since the length of the resulting spring that connects beads ai with bi+1 would be longer than the hopping attempt radius. The BD/kMC simulations were conducted using the EMSIPON code for Brownian dynamics simulations.35,52 The atomistic representations and the figures presented herein were generated using relevant software.84,85

where mmonomer is the mass of a methylene group and ζmonomer is the monomeric friction coefficient per methylene or methyl segment. ζmonomer is a quantity that is readily available from MD simulations and its calculation constitutes a natural step of the proposed multiscale strategy.2,4 Harmandaris et al. estimated the friction factor from equilibrium MD simulations from both the self-diffusion (ζD ≈ 4.7 × 10−13 kg/s) and the Rouse time (ζτ ≈ 4.0 × 10−13 kg/s) in the limit of large molar masses below the entanglement molar mass (C90−C150), where it reached an asymptotic chain length-independent value that was in excellent agreement with experiment.2 In this study, we chose to work with the experimental value of the monomeric friction coefficient ζexp ≈ 4.15 × 10−13 kg/s, since it lies within the predictions of the atomistic simulations.81 Assuming that the velocities of the particles decorrelate much faster than the duration of a time step Δt of the BD simulation, the equations of motion can be integrated using the algorithm of van Gunsteren and Berendsen82 shown below: r′i , α (tn + Δt ) = ri , α(tn) +

⎤ 1⎡ 1 Fi , α(tn)Δt + Fi ,̇ α(tn)(Δt )2 ⎥ ⎦ ζi ⎢⎣ 2

+ R i , α(Δt )

(30)

(28)

where i is an index that runs over all particles, α is the dimension (x, y, z), Ḟ is the time derivative of the systematic force and Ri,α is a random variable obtained from a Gaussian 2k T distribution with zero mean and variance ⟨R i2, α(Δt )⟩ = ζB Δt . i

The time step of the equilibrium simulations was set to 50 ps. Simulations with much shorter time steps (1 and 5 ps) were performed and yielded identical structural and dynamical properties in small systems. The nonequilibrium properties of the studied systems were examined by performing BD/kMC simulations under steady shear flow (NEBD/kMC). The Lees − Edwards83 boundary conditions were implemented along the direction of the gradient (y-axis), while a suitable flow field was imposed along the x-axis. For a given shear rate γ̇ the positions of the particles are updated as follows:

4. PROPERTIES OF THE QUIESCENT MELTS 4.1. Structural Characteristics. The local structure of the chains is governed by the free energy effective potentials in eq 5 extracted from atomistic simulations and experimental data. The mesoscopic simulations conducted here reproduce the theoretical length distributions of the strands (Figure 4a), slipsprings (Figure 5a) and bending angles (see Figure S3 of the Supporting Information) accurately and so, in principle, our model should be able to predict the long-length scale conformations of the chains accurately. Figure 7a displays the end-to-end distance distribution from MD, EBD, and EBD/ kMC simulations for a short chain made of five beads, and Figure 7b shows the computed ⟨Rete2⟩ for various molar masses. In the absence of bending potentials, the end-to-end distance distribution from EBD deviates from the one extracted from MD simulations (see Figure 7). The discrepancy between EBD and the extrapolated MD data increases with the molar mass since the proportionality constant is slightly different. A similar behavior has been observed in the work of Padding and Briels (Table 1 in ref 24). In that work,24 based on coarse-graining PE

ri , α(tn + Δt ) ⎧ r′i , α (tn + Δt ), α = y , z ⎪ =⎨ ⎪ ̇ t, α = x ⎩ r′i , α (tn + Δt ) + 0.5[ri , y(tn) + r′i , y (tn + Δt )]γ Δ

(29)

The mobility of the beads at short times, where the chains have not yet experienced the confinements due to the “tube” is governed by the friction coefficient. The long time behavior of the chains, on the other hand, is governed by the mobility and the lifetime of the slip-springs which act as TCs. The connectivity of the slip-springs is updated through a kinetic I


Article

Macromolecules

Figure 7. (a) End-to-end distance distribution of coarse-grained chains with five beads (C260H522) in a melt and (b) the dependence of the mean squared end-to-end distance ⟨Rete2⟩ on the molar mass as computed from MD, EBD and EBD/kMC simulations. Results from this work have been obtained from atomistic systems (circles) and from mesoscopic systems with strand (triangles), strand−angle (×-markers) and strand−angle−slip-spring (squares) contributions to the effective potential. Purple markers depict MD results from the literature (crosses,2 dashes,3 asterisks4), while diamonds show results from the BD study of Padding and Briels.24 The data points from BD simulations in (b) have been modified using eq 31. The dashed line is a linear fit based on the MD results from the current work.

Figure 8. Dependence of (a) the self-diffusion coefficient and (b) the longest relaxation time τ1 (obtained from the autocorrelation function of the end-to-end vector) on molar mass, from this work (MD, circles; EBD/kMC, red squares), as well as from various simulation works in the literature (green triangles,2 blue crosses,3 orange ×-markers,4 black dashes,5 purple diamonds6). The dotted line in part a shows a fitting on experimental results performed by Lodge,87 while the dashed lines in parts a and b are guides to the eye.

chains using the center of mass approach and using only strand effective free energy potentials, even though the scaling of Rete vs M is in excellent agreement with MD simulations, the obtained Rete is slightly underestimated (see Figure 7b, blue diamonds). Please note that, in order to compare Rete between MD and EBD one has to account for the fact that the mesoscopic chains miss half a chain segment on each chain end due to coarse-graining. An approximate way to correct Rete from EBD and make comparisons with MD is to perform the following estimation: nbeads per chain + 1 R′ete,BD = R ete,BD nbeads per chain (31)

Ideally, the introduction of slip-springs should only affect the viscoelastic properties of the system and not the conformational properties. At low molar masses (up to ∼8000 g/mol) the effect of slip-springs on the chain conformations is very weak (see Figure 7 and Table 1) since the slip-spring density is quite low and the chains are short. At higher molar masses, however, the effect of slip-springs is stronger and combined with the nonbonded interactions (Sanchez−Lacombe EoS) leads to a minor swelling of chains. By accounting for excluded volume effects, nonbonded interactions maintain a uniform density in the system. Similar behavior was observed in our previous work.35 When the nonbonded free energy is on, without any compensating potential, the perturbation of chain conformations caused by slip-springs is not strong and the scaling between ⟨Rete2⟩ and molar mass remains linear. Moreover, our model follows the theoretical relation between end-to-end distance and radius of gyration86 ⟨Rete2⟩ ≈ 6⟨Rg2⟩, while the characteristic ratio (obtained from Rete,BD ′ in eq 31) ranges from

In the presence of the effective bending potentials the results from EBD coincide with those from MD, since the slight skewness of the angle distribution toward 180° (see Figure 3b, blue diamonds) stiffens the coarse-grained chain, and Rete increases. J


Article

Macromolecules

Figure 9. (a) Stress relaxation function from equilibrium EBD/kMC simulations. (b) Zero shear viscosity from MD6 (triangles), NEBD/kMC (squares), EBD/kMC (diamonds), and experimental data (circles,6 asterisk88 and cross89) versus molar mass. The dotted line in part a and the dashed line in part b are guides to the eye and illustrate the t−0.5 and M3.4 scalings, respectively. The red error bars in part b display the error estimate of the extrapolated viscosity from NEBD/kMC calculations.

Another quantity that characterizes the dynamics of polymer melts is the longest relaxation time, usually denoted as τ1. The longest relaxation time was computed from the autocorrelation function of the end-to-end vector by performing a numerical fit with the expression y = A exp(−t/τ1) at time regimes where the function shows a clear exponential behavior. Results from this work and from the literature are displayed in Figure 8b. The longest relaxation time from EBD/kMC is in good match with MD simulations and scales as M3.6 at large molar masses, similarly with ref 40. Figure 9a displays the time evolution of the shear relaxation modulus G(t) from EBD/kMC simulations calculated from the stress tensor using the multiple tau correlator algorithm.35,90 G(t) was computed from the off-diagonal elements of the stress tensor which contains contributions from strands, angles, slipsprings, and nonbonded interactions, the latter affecting only the diagonal elements.35 Because of the fact that the calculation of G(t) from EBD/kMC is a slowly converging process− especially when studying polymer chains with large molar masses−the data was fitted with Maxwell modes at long time regimes.35 In the quasi-entangled regime (M ≈ 3500 g/mol), G(t) decreases steeply; at these molar masses the slip-spring density is low, giving rise to a Rouse-like behavior at short times (G(t) ∼ t−0.5). At larger molar masses (M > 15000 g/mol), however, G(t) features an onset of a rubbery plateau that broadens with rising molar mass. The rubbery plateau GN from EBD/kMC simulations is located at about 2 MPa, while the experimental value for HDPE is GN = 2.6 MPa at 413 K.91 Figure 9b depicts the zero shear viscosity η0 from MD,6 EBD/kMC, NEBD/kMC, and experimental data.6,88,89 In EBD/kMC simulations, η0 was calculated directly from the shear relaxation modulus using a Green−Kubo expression92 (see Figure 9b, diamonds). At low molar masses, even though there is a reasonable qualitative agreement between EBD/kMC and experiment, the scaling of η0 vs M is rather weak. This can be attributed to the fact that these systems are weakly entangled and the rubbery plateau on G(t) has not been formed yet. However, as the molar mass increases, the plateau modulus starts to stabilize and the scaling gets stronger.

8 (crossover regime) to 9 (well entangled regime), in favorable agreement with atomistic simulations.2−4,6 4.2. Rheology. The model incorporated in this work captures smoothly the scaling transitions of the mean square displacement over a vast range of time scales.35 Figure 8a presents the dependence of the self-diffusion coefficient on molar mass from this work (MD, black circles; EBD/kMC, red squares) and from the literature.2−6 The dotted curve in Figure 8a is a fit to experimental data reported by Lodge87 from about 2500 up to 40000 g/mol with a scaling exponent equal to −2.3. In this study the self-diffusion coefficient was computed from the mean square displacement as D = lim

g1(t )

t →∞ 6t

. The self-

diffusion coefficient obtained from MD (Figure 8a, circles) is in good match with the data points from the literature and captures the behavior at the crossover region between the Rouse and the entangled regime6 at approximately 1500 g/mol. Additionally, MD simulations are in good qualitative agreement with the experimental findings of Lodge,87 although the values of the predicted diffusion coefficient are slightly larger. The results from EBD/kMC are in reasonable agreement with MD simulations, although the scaling of D versus molar mass is slightly stronger than the one obtained from experiments and MD simulations. Similar behavior has been observed in many relevant works in the literature. In ref 40 the self-diffusion coefficient from algorithms A1 and A3 (that can be calculated from Table 2 in ref 40 as Dk = R02/τk) shows a scaling stronger than −2.3. References 36 and 41 report a scaling close to −2 for chains with 10−70 beads, while the longer chains in ref 36 display a scaling slightly stronger than −2.3. The discrepancy/ consistency of MD, BD and experimental data with the predictions of the tube model can be further intensified through additional plots that present the dependence of DN2 on N, where N is the length of the chain.34,37 A relevant plot (DM2 versus M) is displayed in Figure S9 of the Supporting Information to the present paper. Generally, the results from EBD/kMC simulations are in good quantitative agreement with experimental data, hence the quiescent properties computed here are relevant to reality. K


Article

Macromolecules

5. PROPERTIES OF THE SHEARED MELTS Besides EBD/kMC simulations, the linear polyethylene melts were studied under steady shear deformation by employing the

Figure 11. (a) Ψ1 and (b) − Ψ2 versus Wi from NEBD simulations with chains of molar mass 3649 (squares), 7296 (triangles), 14590 (diamonds), and 29177 (circles) g/mol. Inset: the extrapolated Ψ1 for Wi → 0 versus molar mass from NEBD (circles) and experiment89 (cross). The experimental scaling of Ψ1 is shown with the broken line. The error bars represent error estimates obtained through the block averaging method.95

Figure 10. (a) Shear viscosity and (b) the normalized number slipsprings per chain (z̃) versus Wi. Markers depict results from chains with 5 (squares), 10 (triangles), 20 (diamonds), and 40 (circles) beads. The asterisks in (a) denote results from chains with 5 beads in the absence of slip-springs. The crosses and the x×-markers in (b) were obtained from refs 42 and 38, respectively. The error bars represent error estimates obtained through the block averaging method.95

simulations over a large range of molar masses. Asterisks in Figure 10a depict results from a system with unentangled, fivebead chains. In this case, the calculated viscosity is lower by an order of magnitude as compared with the similar entangled system. Thus, the contribution of slip-springs to the viscoelastic properties of the material is very strong, even for the shorter PE chains examined here. The employed model can reproduce shear thinning phenomena that occur at large shear rates. As illustrated in Figure 10b, for Wi ≤ 1, the normalized number of slip-springs per chain z̃ which is defined here as z ≐ ̃ Nsls/chain(Wi)/ Nsls/chain(Wi → 0), remains relatively unaffected. For the longer chains examined here, for intermediate values of Wi (10−1 − 10−2), z̃ scales as Wi−0.2, in good agreement with the Brownian dynamics simulations from ref 42 (see Figure 10b, crosses) and with the data from molecular dynamics simulations10 presented therein. Qualitatively similar behavior was observed in ref 38, although shear thinning effects took place at lower Wi values of about Wi ≈ 10−2 in that reference. For even larger values of Wi, z̃ decreases with a larger rate and scales as Wi−0.6. Finally, a similar response is observed for the lower molar mass melts, although the shear thinning phenomena there take place at

Lees−Edwards boundary conditions and the integration scheme discussed in section 3.2. The melt systems were simulated over a large range of Weissenberg values (Wi), where Wi = γ̇τ1 and τ1 is the relaxation time extracted from the autocorrelation of the end-to-end vector (displayed in Table 2). On the basis of these NEBD/kMC simulations, the nonlinearity of several structural and dynamical properties was examined. The shear viscosity was computed for various shear rates using the xy component of the stress tensor: η = τxy/γ̇. As shown in Figure 10a, for Wi > 1, the shear viscosity can be described by a power law of the form: η ∼ Wi−α, where α ≈ 0.7 ± 0.5; this result compares well with experiments89,93 and relevant slip-spring,38 slip-link94 and classical multibead−spring approaches.21,22 For values of Wi less than one, η forms a plateau around the value of η0. The extrapolated value of η0 in the limit Wi → 0 is presented in Figure 9b (squares) and is in excellent agreement with the calculations from EBD/kMC L


Article

Macromolecules

Figure 12. Shear stress τxy versus simulation time for (a) M = 7296 g/mol and (b) M = 29177 g/mol for various values of Wi.

Figure 13. (a, b) Distribution of the end-to-end distance and (c, d) the second rank order parameter of the end-to-end vector versus Wi. Parts a and c present results from systems with 10 beads per chain while parts b and d show results from chains that have 40 beads.

molar masses, cohere with the experimental findings discussed above. As Wi → 0, Ψ1 and −Ψ2 gradually form a plateau while the scaling of Ψ1 at large shear rates is close to the experimental Wi−1.37 and numerical simulations.96,97 The extrapolated value of Ψ1 as Wi → 0, is depicted in the inset of Figure 11a, and seems to be in reasonable agreement with experiment,93 although the extrapolated value for larger molar masses seems underestimated. The ratio − Ψ2/Ψ1 was calculated as 0.17 ± 0.03 at low shear rates, a value that is within the prediction of atomistic simulations (according to ref 10, Ψ2/Ψ1 = 0.18 ± 0.10). Moreover, the material constant θ can be obtained from the ratio of the normal stress differences at low shear rates as

lower Wi numbers. This effect can be attributed to the fact that chains with lower molar mass have a slightly decreased slipspring activity (see Table 2), and hence, the replenishment rate of the destroyed slip-springs is lower. In this work, the first (Ψ1) and the second (Ψ2) normal σzz − σyy σxx − σyy stress coefficients are defined as Ψ1 = and Ψ2 = 2 2 γ̇

γ̇

in view of the geometry of the current system. Experiments on commercial HDPE (Chemplex 6109)93 at T = 453 K show that in the regime of large shear rates, Ψ1 scales as Wi−1.37. Furthermore, experimental works89,93 also show that both Ψ1 and Ψ2 form a plateau as Wi → 0, while Ψ2 is strictly negative due to conformational asymmetries between −y and −z directions. The results in Figure 11, which shows the dependence of Ψ1 and −Ψ2 on Wi over a large range of

Ψ2 Ψ1

M

=

θ . 1−θ

The material constant was calculated as θ = −0.20 ± DOI: 10.1021/acs.macromol.7b00694 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules 0.04. For comparison Mitsoulis89 adopted a value of θ for HDPE equal with θ = −0.11. Figure 12 presents the evolution of τxy from melts with chains comprised of 10 (7296 g/mol) and 40 (29177 g/mol) beads over the time span of a NEBD/kMC simulation under different shear rates. For large shear rates (Wi > 1), τxy displays a strong stress overshoot with a duration that scales as Wi−0.50 and Wi−0.65 for chains with 10 and 40 beads, respectively. During this transition period the density of entanglements decreases and the chains align in parallel with the direction of the flow. In contrast, for Wi ≤ 1 the stress overshoot disappears and τ xy increases gradually over a timespan that is commensurate with the longest relaxation time τ1. Very similar response has been observed in several numerical calculations97,98 and in experimental data from entangled polystyrene melts.99 The application of strong shear deformations is expected to have strong effects on the long-range conformations of the chains. Figure 13 shows the distribution of the end-to-end distance ρ(Rete) and the second rank order parameter (P2,α) of the end-to-end vector, from two systems with 10 (7296 g/mol) and 40 (29177 g/mol) beads, for various values of Wi. In this study, the second rank order parameter is defined as P2, α =

3 1 ⟨cos2(θα)⟩ − , α = x , y , z 2 2

were restored through the use of slip-springs, which are artificial topological constrains that mimic the effect of the entanglement. The nonlinear response of polyethylene was studied through steady shear flow simulations by using the Lees− Edwards boundary conditions and by imposing a suitable velocity field along the direction of the flow. Our results were compared with relevant works from the literature and experimental data when available. The current study is a stepping stone toward future studies in systems with polymer/ solid interfaces such as confined polymer melts and nanocomposite materials.

■

APPENDIX A: A LONG-RANGE APPROXIMATION FOR THE COMPENSATING POTENTIAL The compensating potential developed in ref 38 is a repulsive potential between all possible pairs that can be connected by slip-springs. The free energy, force, and resulting virial between a pair ij that is able to form a slip-spring are given by ⎛ A ss(rij) ⎞ ⎟⎟ Acomp(rij) = zactivkBT exp⎜⎜ − ⎝ kBT ⎠

(33)

⎛ A ss(rij) ⎞ ⎟⎟Fss(rij) Fcomp(rij) = −∇rij Acomp(rij) = −zactiv exp⎜⎜ − ⎝ kBT ⎠

(32)

where θα is the angle between the end-to-end vector and a direction α. A value of P2,α = 1 means that Rete is perfectly aligned along direction α, P2,α = 0 shows that Rete is randomly oriented, while, when P2,a= −0.5, Rete is oriented perpendicular to direction α. For Wi ≤ 0.1, R ete distributions remain relatively unperturbed; thus, they are similar to those obtained from EBD/kMC simulations. For such low shear rates P2,α ≃ 0, hence the chains do not present any preference in their orientation. Upon increasing the shear rate, ρ(Rete) broadens significantly due to the elongation of the chains along the direction of the flow. As shown in Figure 13, parts c and d, P2,x gradually increases with Wi, hence the chains tend to orient in parallel with the flow direction. Interestingly, even though both P2,y and P2,z take negative values with increasing Wi, the effect is stronger along the y-direction. This phenomenon correlates strongly with the fact that Ψ2 is strictly negative since, on average, chains experience larger stresses along the y- than along the z- direction, σyy > σzz, due to the fact that the velocity gradient is imposed along the y-axis (see eq 29), hence generating increased stress along this direction. For very large Wi values (Wi ≥ 100) the distribution becomes bimodal. The bimodality of the distributions is attributed to the fact that the strands employed in this study are finitely extensible, and thus the chains can only be extended up to a maximum finite value. The obtained behavior of the end-to-end distributions cohere with relevant MD and coarse-grained BD simulations.3,10

(34)

Vcomp(rij) = Fcomp(rij) ·rij

(35)

where Ass(rij) is the free energy and Fss(rij) is the force from a slip-spring with length rij. The contribution of the compensating potential to the total free energy and virial is ⎛ A ss(rij) ⎞ ⎟⎟ Acomp,total = zactivkBT ∑ exp⎜⎜ − ⎝ kBT ⎠ i>j

(36)

⎡ ⎛ A (r ) ⎞ ⎤ ss ij ⎟⎟Fss(rij) ·rij⎥ Vcomp,total = −zactiv ∑ ⎢exp⎜⎜ − ⎢ ⎝ kBT ⎠ ⎥⎦ i>j ⎣

(37)

Using eq 37, one can compute the contribution of the compensating potential to the pressure (either analytically or numerically) which is equal to the negative of the contribution of slip-springs Vss,total to the pressure. However, the computation of Acomp,total in systems with long slip-springs is computationally expensive, since the number of operations scales as r3. A way to accelerate the calculation is to compute the interactions up to a cutoff distance rc and then account for the omitted interactions by multiplying Fcomp(rij) by a scale factor λ so that V ccomp,total(λ) = −λzactiv

6. CONCLUSIONS A multiscale strategy involving atomistic (MD) and mesoscopic (BD/kMC) simulations was developed in order to study monodisperse polyethylene melts over a broad range of time and length scales. The bridging between those different levels of description was performed through a rigorous mapping procedure that led to an accurate replication of the local and long-range conformations obtained from atomistic simulations. The viscoelastic properties of the melt at the mesoscopic level

⎡ ⎛ A (r ) ⎞ ⎤ ⎢exp⎜⎜ − ss ij ⎟⎟Fss(rij) ·rij⎥ ⎢ ⎝ kBT ⎠ ⎥⎦ i > j , rij < rc ⎣

∑

(38)

Given that the contribution of Vcomp,total to the pressure from eq 37 equals the negative of the contribution of Vss,total, it would be reasonable to choose λ such that the contribution of the approximate compensating potential Vccomp,total(λ) equals the contribution of − Vss,total as well, hence negating the effect of slip-springs. Since Vccomp,total(λ) is proportional to λ, the scale factor can be computed as N


Macromolecules

■

⎤ ⎡ A ss(rij) ∑i > j ⎢exp − k T Fss(rij) ·rij⎥ ⎦ ⎣ B λ= c = ⎤ ⎡ A ss(rij) V comp,total(1) ∑i > j , r < r ⎢exp − k T Fss(rij) ·rij⎥ ij c ⎣ ⎦ B

(

Vcomp,total

APPENDIX C: OPTIMIZATION OF THE FLUCTUATING SLIP-SPRING SCHEME In the fluctuating slip-spring scheme developed in a previous study,35 the algorithm responsible for the creation of slipsprings at the Nends chain ends of a system is designed as follows: (1) Iteration over all the chain ends (2) Iteration over all possible candidates within an attempt radius rattempt (3) Creation of a new slip-spring with probability:

)

(

)

(39)

Finally, eq 38 reduces to eq 37 upon substituting λ from eq 39.

■

APPENDIX B: THERMODYNAMIC AND MEAN-FIELD CALCULATION OF THE ACTIVITY THROUGH SIMULATIONS WITH CONSTANT NUMBER OF SLIP-SPRINGS The configurational partition function of the system of polymers and slip-springs is written as38 ∞

Z=

∫ +[{r}]e−A({r})/(k T) ∏ ∑ B

i > j nij = 0

Article

Pform, i = v hopzactivncands, iΔtkMC

(44)

where ΔtkMC is the time span between successive iterations of kMC and ncands,i is the number of candidate slip-springs that can be formed from the ith end. In homogeneous systems, the average number of candidates equals:

n

ij zactiv e−nijAss(|ri − rj|)/(kBT ) nij!

(40)

with A symbolizing Astrands+Aangles+Anonbonded and Ass being the Helmholtz energy of a single slip-spring. ri is the position vector of bead i, while nij is the number of slip-springs connecting beads i and j. The total number of slip-springs in a configuration would be ∑i > j nij , so the mean number of slipsprings equals 1 ⟨nss⟩ = ⟨∑ nij⟩ = Z i>j ∞

∫ +[{r}]e

⟨ncands⟩ =

i > j nij = 0

⟨Ncreate⟩ = Nends⟨Pform⟩

−A({r})/(kBT )

(∑ nij)

ij zactiv e−nijAss(|ri − rj|)/(kBT ) nij!

= zactiv⟨∑ e−Ass(|ri − rj|)/(kBT )⟩ (41)

hence, zactiv =

⟨∑i > j e

⟨nss⟩ −A ss(|ri − rj|)/(kBT )

⟩

P′form, i = Pform, i /Pselect

(42)

′ ⟩ = NendsPselect⟨Pform ′ ⟩ ⟨Ncreate

(48)

which exactly reduces to eq 46 upon substituting eq 47. Thus, the number of newly created slip-springs is the same. Some care is needed in setting the value of Pselect, since, if P′form surpasses ′ ⟩ unity, then some elementary events are missed and ⟨Ncreate decreases (our code halts and issues a warning). In this study, Pselect was set to 0.005, thus the algorithm presents a speedup of ×200 times.

■

v hop⟨ncands⟩ΔtkMCNends ⟨Ncreate⟩

(47)

After completion of the algorithm, the average number of the slip-springs equals

A Widom-like expression has been obtained which relates the mean number of slip-springs to the activity in the considered grand canonical model with respect to slip-springs. This allows the estimation of zactiv for a desired ⟨Nss⟩ from a simulation with a constant number of slip-springs (i.e., in which the destruction events are immediately followed by creation events).35 Another way to compute zactiv is through the calculation of the average number of creation or destruction events at each kMC step, from an equilibrium simulation with constant number of slip-springs of a homogeneous system. With ⟨Ncreate⟩ known, the activity can be computed from eqs 44-46 as zactiv =

(46)

This algorithm is computationally consuming, since a total of Nends iterations have to be performed in order to create new slip-springs within a given instance. In this work, a modified version of the algorithm was used, which is the following: (1) Iteration over all the chain ends (2) A chain is selected to try to form a new slip-spring with a probability Pselect. (3) If the selection is rejected, then the algorithm returns to step (1). (4) If the selection is accepted then a new slip-spring is formed with a probability:

i>j

i>j

(45)

with ρ being the overall number density of beads. After completion of the algorithm, the average number of new slip-springs is

n

∏∑

4 πrattempt 3ρ 3

ASSOCIATED CONTENT

* Supporting Information S

(43)

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.7b00694. . (Figure S1) Viscosity with and without the compensating potential; (Figure S2) bond-angle distributions from MD; (Figure S3) bond-angle distributions from MD and BD; (Figure S4) distributions of entanglements per chain

Even though the dynamical estimate of zactiv from eq 43 converges at a much slower rate than its thermodynamic estimate by eq 42, the activities obtained through both methodologies are practically identical. Activities employed in the fluctuating slip-spring scheme simulations reported in this article are listed in Table 2. O


Article

Macromolecules

■

Polyethylene Melt in Steady Shear. Macromolecules 2010, 43, 6886− 6902. (11) Edwards, S. F. The Statistical Mechanics of Polymerized Material. Proc. Phys. Soc., London 1967, 92, 9−16. (12) de Gennes, P. G. Reptation of a Polymer Chain in the Presence of Fixed Obstacles. J. Chem. Phys. 1971, 55, 572. (13) Doi, M. Explanation for the 3.4-Power Law for Viscosity of Polymeric Liquids on the Basis of the Tube Model. J. Polym. Sci., Polym. Phys. Ed. 1983, 21, 667−684. (14) Milner, S. T.; McLeish, T. C. B. Reptation and Contour-Length Fluctuations in Melts of Linear Polymers. Phys. Rev. Lett. 1998, 81, 725−728. (15) Likhtman, A. E.; McLeish, T. C. B. Quantitative Theory for Linear Dynamics of Linear Entangled Polymers. Macromolecules 2002, 35, 6332−6343. (16) McLeish, T. C. B. Tube Theory of Entangled Polymer Dynamics. Adv. Phys. 2002, 51, 1379−1527. (17) Tzoumanekas, C.; Theodorou, D. N. Topological Analysis of Linear Polymer Melts. Macromolecules 2006, 39, 4592−4604. (18) Tzoumanekas, C.; Theodorou, D. N. From Atomistic Simulations to Slip-Link Models of Entangled Polymer Melts: Hierarchical Strategies for the Prediction of Rheological Properties. Curr. Opin. Solid State Mater. Sci. 2006, 10, 61−72. (19) Anogiannakis, S. D.; Tzoumanekas, C.; Theodorou, D. N. Microscopic Description of Entanglements in Polyethylene Networks and Melts: Strong, Weak, Pairwise, and Collective Attributes. Macromolecules 2012, 45, 9475−9492. (20) Tzoumanekas, C.; Lahmar, F.; Rousseau, B.; Theodorou, D. Onset of Entanglements Revisited. Topological Analysis. Macromolecules 2009, 42, 7474−7484. (21) Kröger, M.; Hess, S. Rheological Evidence for a Dynamical Crossover in Polymer Melts via Nonequilibrium Molecular Dynamics. Phys. Rev. Lett. 2000, 85, 1128−1131. (22) Kröger, M.; Loose, W.; Hess, S. Rheology and Structural Changes of Polymer Melts via Nonequilibrium Molecular Dynamics. J. Rheol. 1993, 37, 1057. (23) Padding, J. T.; Briels, W. J. Zero-Shear Stress Relaxation and Long Time Dynamics of a Linear Polyethylene Melt: A Test of Rouse Theory. J. Chem. Phys. 2001, 114, 8685−8693. (24) Padding, J. T.; Briels, W. J. Time and Length Scales of Polymer Melts Studied by Coarse-Grained Molecular Dynamics Simulations. J. Chem. Phys. 2002, 117, 925−943. (25) Doi, M.; Edwards, S. F. Dynamics of Concentrated Polymer Systems. Part 2.-Molecular Motion under Flow. J. Chem. Soc., Faraday Trans. 2 1978, 74, 1802−1817. (26) Han, C. C.; Akcasu, A. Z. On the Dynamic Structure Factor, S (Q, Ro), of Polymers in Dilute Solution. Polymer 1981, 22, 1019− 1025. (27) Edwards, S. F.; Vilgis, T. The Effect of Entanglements in Rubber Elasticity. Polymer 1986, 27, 483−492. (28) Rubinstein, M.; Panyukov, S. Elasticity of Polymer Networks. Macromolecules 2002, 35, 6670−6686. (29) Schieber, J. D.; Neergaard, J.; Gupta, S. A Full-Chain, Temporary Network Model with Sliplinks, Chain-Length Fluctuations, Chain Connectivity and Chain Stretching. J. Rheol. 2003, 47, 213. (30) Schieber, J. D. Fluctuations in Entanglements of Polymer Liquids. J. Chem. Phys. 2003, 118, 5162−5166. (31) Nair, D. M.; Schieber, J. D. Linear Viscoelastic Predictions of a Consistently Unconstrained Brownian Slip-Link Model. Macromolecules 2006, 39, 3386−3397. (32) Khaliullin, R. N.; Schieber, J. D. Analytic Expressions for the Statistics of the Primitive-Path Length in Entangled Polymers. Phys. Rev. Lett. 2008, 100, 188302. (33) Schieber, J. D.; Andreev, M. Entangled Polymer Dynamics in Equilibrium and Flow Modeled Through Slip Links. Annu. Rev. Chem. Biomol. Eng. 2014, 5, 367−381. (34) Likhtman, A. E. Single-Chain Slip-Link Model of Entangled Polymers: Simultaneous Description of Neutron Spin-Echo, Rheology, and Diffusion. Macromolecules 2005, 38, 6128−6139.

and per reduced chain contour; distribution of (Figure S5) the number of beads and (Figure S6) the distance between entanglements; RDF of (Figure S7) the entangled beads and (Figure S8) the center of mass of slip-springs (RDF0/RDF1/RDF2); (Figure S9) DM2 versus M (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

D. N. Theodorou: 0000-0002-4763-9739 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The valuable contribution of Dr. Georgios G. Vogiatzis in all phases of this work is deeply appreciated. The authors gratefully acknowledge financial support through the project “Multiscale Simulations of Complex Polymer Systems” (MuSiComPS) by the Limmat Foundation, Zurich, Switzerland. They would also like to thank Prof. Marcus Müller (Georg-August-Universität Gö ttingen) and Dr. Christos Tzoumanekas for fruitful discussions. Part of this work was supported by computational time granted from the Greek Research & Technology Network (GRNET) in the National HPC facilityARISunder Project ID pr003013 (MSMS).

■

REFERENCES

(1) Denn, M. M. Extrusion Instabilities and Wall Slip. Annu. Rev. Fluid Mech. 2001, 33, 265−87. (2) Harmandaris, V. A.; Mavrantzas, V. G.; Theodorou, D. N. Atomistic Molecular Dynamics Simulation of Polydisperse Linear Polyethylene Melts. Macromolecules 1998, 31, 7934−7943. (3) Nafar Sefiddashti, M. H.; Edwards, B. J.; Khomami, B. Individual Chain Dynamics of a Polyethylene Melt Undergoing Steady Shear Flow. J. Rheol. 2015, 59, 119−153. (4) Harmandaris, V. A.; Mavrantzas, V. G.; Terzis, A. F.; Paspalakis, E. Segmental Dynamics in Polyethylene Melts through Atomistic Molecular Dynamics Simulations. Recent Research Topics and Developments in Chemical Physics: From Quantum Scale to Macroscale 2008, 661, 179−196. (5) Hur, K.; Jeong, C.; Winkler, R. G.; Lacevic, N.; Gee, R. H.; Yoon, D. Y. Chain Dynamics of Ring and Linear Polyethylene Melts from Molecular Dynamics Simulations. Macromolecules 2011, 44, 2311− 2315. (6) Harmandaris, V. A.; Mavrantzas, V. G.; Theodorou, D. N.; Kröger, M.; Ramírez, J.; Ö ttinger, H. C.; Vlassopoulos, D. Crossover from the Rouse to the Entangled Polymer Melt Regime: Signals from Long, Detailed Atomistic Molecular Dynamics Simulations, Supported by Rheological Experiments. Macromolecules 2003, 36, 1376−1387. (7) Theodorou, D. N.; Vogiatzis, G. G.; Kritikos, G. Self-ConsistentField Study of Adsorption and Desorption Kinetics of Polyethylene Melts on Graphite and Comparison with Atomistic Simulations. Macromolecules 2014, 47, 6964−6981. (8) Kritikos, G.; Sgouros, A.; Vogiatzis, G. G.; Theodorou, D. N. Molecular Dynamics Study of Polyethylene under Extreme Confinement. J. Phys.: Conf. Ser. 2016, 738, 012012. (9) Jeong, S.; Cho, S.; Kim, J. M.; Baig, C. Molecular Mechanisms of Interfacial Slip for Polymer Melts under Shear Flow. J. Rheol. 2017, 61, 253−264. (10) Baig, C.; Mavrantzas, V. G.; Kröger, M. Flow Effects on Melt Structure and Entanglement Network of Linear Polymers: Results from a Nonequilibrium Molecular Dynamics Simulation Study of a P


Article

Macromolecules

Small Angle Neutron Scattering. Macromolecules 1997, 30, 4973− 4977. (57) Cohen, A. A Padé Approximant to the Inverse Langevin Function. Rheol. Acta 1991, 30, 270−273. (58) Kröger, M. Simple, Admissible, and Accurate Approximants of the Inverse Langevin and Brillouin Functions, Relevant for Strong Polymer Deformations and Flows. J. Non-Newtonian Fluid Mech. 2015, 223, 77−87. (59) Kremer, K.; Grest, G. S. Dynamics of Entangled Linear Polymer Melts: A Molecular Dynamics Simulation. J. Chem. Phys. 1990, 92, 5057. (60) Fetters, L. J.; Lohse, D. J.; Colby, R. H. CHAPTER 25 Chain Dimensions and Entanglement Spacings. Phys. Prop. Polym. Handb. 2006, 445−452. (61) Fetters, L. J.; Lohse, D. J.; Milner, S. T.; Graessley, W. W. Packing Length Influence in Linear Polymer Melts on the Entanglement, Critical, and Reptation Molecular Weights. Macromolecules 1999, 32, 6847−6851. (62) Foteinopoulou, K.; Karayiannis, N. C.; Mavrantzas, V. G.; Kröger, M. Primitive Path Identification and Entanglement Statistics in Polymer Melts: Results from Direct Topological Analysis on Atomistic Polyethylene Models. Macromolecules 2006, 39, 4207−4216. (63) Masubuchi, Y.; Ianniruberto, G.; Greco, F.; Marrucci, G. Entanglement Molecular Weight and Frequency Response of Sliplink Networks. J. Chem. Phys. 2003, 119, 6925−6930. (64) Paul, D. R.; DiBenedetto, A. T. Diffusion in Amorphous Polymers. J. Polym. Sci., Part C: Polym. Symp. 1965, 10, 17−44. (65) Hur, K.; Winkler, R. G.; Yoon, D. Y. Comparison of Ring and Linear Polyethylene from Molecular Dynamics Simulations. Macromolecules 2006, 39, 3975−3977. (66) Ramos, J.; Peristeras, L. D.; Theodorou, D. N. Monte Carlo Simulation of Short Chain Branched Polyolefins in the Molten State. Macromolecules 2007, 40, 9640−9650. (67) Vogiatzis, G. G.; Theodorou, D. N. Local Segmental Dynamics and Stresses in Polystyrene-C60 Mixtures. Macromolecules 2014, 47, 387−404. (68) Theodorou, D. N.; Suter, U. W. Detailed Molecular-Structure of a Vinyl Polymer Glass. Macromolecules 1985, 18, 1467−1478. (69) MAPS; Scienomics: Paris, 2015. (70) Pant, P. V. K.; Theodorou, D. N. Variable Connectivity Method for the Atomistic Monte Carlo Simulation of Polydisperse Polymer Melts. Macromolecules 1995, 28, 7224−7234. (71) Mavrantzas, V. G.; Boone, T. D.; Zervopoulou, E.; Theodorou, D. N. End-Bridging Monte Carlo: A Fast Algorithm for Atomistic Simulation of Condensed Phases of Long Polymer Chains. Macromolecules 1999, 32, 5072−5096. (72) Karayiannis, N. C.; Mavrantzas, V. G.; Theodorou, D. N. A Novel Monte Carlo Scheme for the Rapid Equilibration of Atomistic Model Polymer Systems of Precisely Defined Molecular Architecture. Phys. Rev. Lett. 2002, 88, 105503. (73) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (74) LAMMPS Molecular Dynamics Simulator. Available fromhttp://lammps.sandia.gov/, 2017. (75) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (76) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (77) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible Multiple Time Scale Molecular Dynamics. J. Chem. Phys. 1992, 97, 1990−2001. (78) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87, 1117−1157. (79) Ermak, D. L.; Yeh, Y. Equilibrium Electrostatic Effects on the Behavior of Polyions in Solution: Polyion-Mobile Ion Interaction. Chem. Phys. Lett. 1974, 24, 243−248. (80) Ermak, D. L.; McCammon, J. A. Brownian Dynamics with Hydrodynamic Interactions. J. Chem. Phys. 1978, 69, 1352.

(35) Vogiatzis, G. G.; Megariotis, G.; Theodorou, D. N. Equation of State-Based Slip-Spring Model for Entangled Polymer Dynamics. Macromolecules 2017, 50, 3004−3029. (36) Langeloth, M.; Masubuchi, Y.; Böhm, M. C.; Müller-plathe, F. Recovering the Reptation Dynamics of Polymer Melts in Dissipative Particle Dynamics Simulations via Slip-Springs. J. Chem. Phys. 2013, 138, 104907. (37) Uneyama, T.; Masubuchi, Y. Multi-Chain Slip-Spring Model for Entangled Polymer Dynamics. J. Chem. Phys. 2012, 137, 154902. (38) Chappa, V. C.; Morse, D. C.; Zippelius, A.; Müller, M. Translationally Invariant Slip-Spring Model for Entangled Polymer Dynamics. Phys. Rev. Lett. 2012, 109, 148302. (39) Ramírez-Hernández, A.; Peters, B. L.; Schneider, L.; Andreev, M.; Schieber, J. D.; Müller, M.; de Pablo, J. J. A Multi-Chain Polymer Slip-Spring Model with Fluctuating Number of Entanglements: Density Fluctuations, Confinement, and Phase Separation. J. Chem. Phys. 2017, 146, 014903. (40) Ramírez-Hernández, A.; Detcheverry, F. A.; Peters, B. L.; Chappa, V. C.; Schweizer, K. S.; Müller, M.; de Pablo, J. J. Dynamical Simulations of Coarse Grain Polymeric Systems: Rouse and Entangled Dynamics. Macromolecules 2013, 46, 6287−6299. (41) Masubuchi, Y. Simulating the Flow of Entangled Polymers. Annu. Rev. Chem. Biomol. Eng. 2014, 5, 11−33. (42) Masubuchi, Y. Effects of Degree of Freedom below Entanglement Segment on Relaxation of Polymer Configuration under Fast Shear in Multi-Chain Slip-Spring Simulations. J. Chem. Phys. 2015, 143, 224905. (43) Masubuchi, Y.; Takimoto, J. I.; Koyama, K.; Ianniruberto, G.; Marrucci, G.; Greco, F. Brownian Simulations of a Network of Reptating Primitive Chains. J. Chem. Phys. 2001, 115, 4387−4394. (44) Uneyama, T.; Masubuchi, Y. Detailed Balance Condition and Effective Free Energy in the Primitive Chain Network Model. J. Chem. Phys. 2011, 135, 184904. (45) Ramírez-Hernández, A.; Peters, B. L.; Andreev, M.; Schieber, J. D.; de Pablo, J. J. A Multichain Polymer Slip-Spring Model with Fluctuating Number of Entanglements for Linear and Nonlinear Rheology. J. Chem. Phys. 2015, 143, 243147. (46) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of N -Alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (47) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. On the Simulation of Vapor−liquid Equilibria for Alkanes. J. Chem. Phys. 1998, 108, 9905. (48) Dodd, L. R.; Theodorou, D. N. Atomistic Monte Carlo Simulation and Continuum Mean Field Theory of the Structure and Equation of State Properties of Alkane and Polymer Melts. Adv. Polym. Sci. 1994, 116, 249−281. (49) Toxvaerd, S. Equation of State of Alkanes II. J. Chem. Phys. 1997, 107, 5197. (50) Sanchez, I. C.; Lacombe, R. H. An Elementary Molecular Theory of Classical Fluids. Pure Fluids. J. Phys. Chem. 1976, 80, 2352− 2362. (51) Sanchez, I. C.; Lacombe, R. H. Statistical Thermodynamics of Polymer Solutions. Macromolecules 1978, 11, 1145−1156. (52) Megariotis, G.; Vogiatzis, G. G.; Schneider, L.; Müller, M.; Theodorou, D. N. Mesoscopic Simulations of Crosslinked Polymer Networks. J. Phys.: Conf. Ser. 2016, 738, 012063. (53) Chen, L. J.; Qian, H. J.; Lu, Z. Y.; Li, Z. S.; Sun, C. C. An Automatic Coarse-Graining and Fine-Graining Simulation Method: Application on Polyethylene. J. Phys. Chem. B 2006, 110, 24093− 24100. (54) Lahmar, F.; Tzoumanekas, C.; Theodorou, D. N.; Rousseau, B. Onset of Entanglements Revisited. Dynamical Analysis. Macromolecules 2009, 42, 7485−7494. (55) Horton, J. C.; Squires, G. L.; Boothroyd, A. T.; Fetters, L. J.; Rennie, A. R.; Glinka, C. J.; Robinson, R. A. Small-Angle Neutron Scattering from Star-Branched Polymers in the Molten State. Macromolecules 1989, 22, 681−686. (56) Fetters, L. J.; Graessley, W. W.; Krishnamoorti, R.; Lohse, D. J. Melt Chain Dimensions of Poly(ethylene-1-Butene) Copolymers via Q


Article

Macromolecules (81) van Meerveld, J. A Method to Extract the Monomer Friction Coefficient from the Linear Viscoelastic Behavior of Linear, Entangled Polymer Melts. Rheol. Acta 2004, 43, 615−623. (82) van Gunsteren, W. F.; Berendsen, H. J. C. Algorithms for Brownian Dynamics. Mol. Phys. 1982, 45, 637−647. (83) Lees, A. W.; Edwards, S. F. The Computer Study of Transport Processes under Extreme Conditions. J. Phys. C: Solid State Phys. 1972, 5, 1921−1928. (84) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (85) Kraus, D. Consolidated Data Analysis and Presentation Using an Open-Source Add-in for the Microsoft Excel® Spreadsheet Software. Med. Writ. 2014, 23, 25−28. (86) Theodorou, D. N.; Dünweg, B.; Landau, D. P.; Milchev, A. I. Computer Simulations of Surfaces and Interfaces. Polymers at Surfaces and Interfaces 2003, 329−419. (87) Lodge, T. P. Reconciliation of the Molecular Weight Dependence of Diffusion and Viscosity in Entangled Polymers. Phys. Rev. Lett. 1999, 83, 3218−3221. (88) Ansari, M.; Alabbas, A.; Hatzikiriakos, S. G.; Mitsoulis, E. Entry Flow of Polyethylene Melts in Tapered Dies. Int. Polym. Process. 2010, 25, 287−296. (89) Mitsoulis, E. Effect of Viscoelasticity in Fountain Flow of Polyethylene Melts. Int. Polym. Process. 2009, 24, 439−451. (90) Ramírez, J.; Sukumaran, S. K.; Vorselaars, B.; Likhtman, A. E. Efficient on the Fly Calculation of Time Correlation Functions in Computer Simulations. J. Chem. Phys. 2010, 133, 154103. (91) Fetters, L. J.; Lohse, D. J.; Colby, R. H. Physical Properties of Polymers Handbook; Mark, J. E., Ed.; Woodbury: New York, 1996. (92) Hartkamp, R.; Daivis, P. J.; Todd, B. D. Density Dependence of the Stress Relaxation Function of a Simple Fluid. Phys. Rev. E 2013, 87, 032155. (93) Min, K.; White, J. L.; Fellers, J. F. High Density Polyethylene/ polystyrene Blends: Phase Distribution Morphology, Rheological Measurements, Extrusion, and Melt Spinning Behavior. J. Appl. Polym. Sci. 1984, 29, 2117−2142. (94) Ramírez Hernández, A.; Müller, M.; de Pablo, J. J. Theoretically Informed Entangled Polymer Simulations: Linear and Non-Linear Rheology of Melts. Soft Matter 2013, 9, 2030−2036. (95) Flyvbjerg, H.; Petersen, H. G. Error Estimates on Averages of Correlated Data. J. Chem. Phys. 1989, 91, 461. (96) Stephanou, P. S.; Kröger, M. Solution of the Complete CurtissBird Model for Polymeric Liquids Subjected to Simple Shear Flow. J. Chem. Phys. 2016, 144, 124905. (97) Ilg, P.; Kröger, M. Molecularly Derived Constitutive Equation for Low-Molecular Polymer Melts from Thermodynamically Guided Simulation. J. Rheol. 2011, 55, 69. (98) Stephanou, P. S.; Tsimouri, I. C.; Mavrantzas, V. G. FlowInduced Orientation and Stretching of Entangled Polymers in the Framework of Nonequilibrium Thermodynamics. Macromolecules 2016, 49, 3161−3173. (99) Schweizer, T.; van Meerveld, J.; Ö ttinger, H. C. Nonlinear Shear Rheology of Polystyrene Melt with Narrow Molecular Weight distributionExperiment and Theory. J. Rheol. 2004, 48, 1345.

R


Slip-Spring Model for the Linear and Nonlinear Viscoelastic Properties

Recommend Documents