Viscoelastic Properties of Polymer Melts from Equilibrium Molecular Dynamics Simulations

Introduction. The dynamics of polymer chains are now generally very well understood in the pioneering frameworks of the Rouse model for short chains, and the reptation model for longer chains.1-5 While these theories yield well-known expressions for the variation of the diffusivity D and viscosity η with the chain length N, an open question is the chain length at which one crosses over from Rouse-like behavior to reptation. This length, termed the entanglement length, Ne, is relatively easy to estimate experimentally, but has been hard to access in a simulation due to the fact that, up until recently, there was no satisfactory definition of an entanglement.6,7 In numerical studies, most estimates of the entanglement length involve conducting simulations and examining dynamic quantities, such as the frequency dependent storage modulus. Extensive work has been done by Kremer and Grest8 to identify the onset of chain entanglement, with somewhat mixed results. Examination of the chain length dependence of chain diffusion yields a cross over from the Rouse scaling (N-1) to reptation scaling (N-2) behavior in the vicinity of chains of length 35. In contrast, estimates of storage modulus from nonequilibrium molecular dynamics simulations yield an entanglement chain length in the vicinity of 80. We approach this problem with the tools of equilibrium molecular dynamics (MD).9 Equilibrium stress fluctuations which result naturally from these simulations directly provide estimates of the viscosity through the Green-Kubo equation, as suggested by Smith et al.10 The storage and loss modulus are also computed from the stress autocorrelation function, and the hints of a plateau are seen in the storage modulus. The plateau is the strongest indication of onset of reptation dynamics, and is seen for chain lengths of 80 and higher. Estimates of the entanglement length Ne, provide a number closer to 30, consistent with the first estimate provided by Kremer and Grest. Simulation Model and Methods. The MD simulation employs a standard chain model. Interaction between nonbonded monomers are described by a shifted, purely repulsive Lennard-Jones (LJ) potential: U(r) ) 4 [(σ/r)12 (σ/r)6] + for r 21/6σ. Adjacent bonded monomers interact via a stiff FENE potential, in the form of VFENE ) k(R0/2) ln(1 (r/R0)), which constrains the distance between adjacent monomers to about 1σ (we use the same parameters for the FENE potential as Grest and Kremer in refs 5 and 8). We performed constant volume simulations of monodisperse polymer melts of varying chain lengths, and will focus our attention on the crossover regime, specifically N ) 20, 40, 80 and 120, respectively. In a typical simulation we use a total of 2400 monomers embedded in the periodic simulation box, of size 14.1σ in each direction, which corresponds to the reduced segment density of F* ) 0.85. In a few simulations we have doubled the number of monomers, and for the longest chain system, N ) 120, this corresponds to increasing the number of chains from 20 to 40. Since the properties deduced from the simulation were independent of system size even for these relatively small systems, we conclude that our results only have minor finite size effects. We will report all quantities in terms of reduced units, which are defined at the end of the paper. We conduct MD simulations using a fifth order Gear algorithm in the microcanonical ensemble with a δt ) 0.001t* ensuring energy conservation to within 0.5% over the whole constant energy simulation run. The starting structures for our runs were constructed by gradually squeezing semidilute solutions to a final density of 0.85 over ∼10 million MD steps at the reduced temperature, T* ) 1. These structures were further run at constant temperature, for an additional 10 million steps to obtain a starting configuration for the constant energy run. After the structure preparation and equilibration, we followed with the constant energy simulation for a minimum of 50 million MD steps for N ) 20, up to 300 million MD steps for N ) 120. The reduced temperature exhibited small (several percent) fluctuations around unity since we did not couple our simulations to a heat bath. It is important to check our simulation results for finite-size effects. Kremer and Grest suggest that for these chain lengths a minimum number of 20 chains is required to avoid such effects. Our systems contain 2400 monomers, i.e., 20 chains for N ) 120. In Table 1 below we have tabulated various quantities like the end-toend vector and diffusivity, for all the chain lengths. In particular, the diffusivity, D, changed from 0.0012 to 0.0011, when the system size was doubled from 2400 to 4800 monomers. Since D is typically the most sensitive to finite size effects, and since these D values closely track the results of Kremer and Grest,5 we conclude that our system sizes are large enough. To check for equilibration we use the autocorrelation function of the end-to-end vector of the chains,4 〈rend(t)rend(0)〉/〈rend〉, where rend(t) denotes the end-to-end distance vector at time t. A similar autocorrelation function can also be defined for the root-mean-square radius of gyration. Both of these autocorrelation functions showed an exponential decay, from which the † Isermann Department of Chemical and Biological Engineering, Rensselaer Polytechnic Institute. ‡ Department of Materials Science and Engineering, Rensselaer Polytechnic Institute. Table 1. a N ) 20 N ) 40 N ) 80 N ) 120 〈Re〉 (σ2) 29.2 62.6 127.5 195.1 D (σ2τ-1) 0.02 0.0074 0.0023 0.0012 a Re is the average end-to-end vector magnitude. D is the diffusivity. BATCH: ma2a46 USER: rjb69 DIV: @xyv04/data1/CLS_pj/GRP_ma/JOB_i03/DIV_ma035487l DATE: December 28, 2004