Article pubs.acs.org/Macromolecules
Molecular Mechanisms for Conformational and Rheological Responses of Entangled Polymer Melts to Startup Shear Yuyuan Lu,† Lijia An,*,† Shi-Qing Wang,‡ and Zhen-Gang Wang§,† †
State Key Laboratory of Polymer Physics and Chemistry, Changchun Institute of Applied Chemistry, Chinese Academy of Sciences, Changchun 130022, People’s Republic of China ‡ Department of Polymer Science, University of Akron, Akron, Ohio 44325-3909, United States § Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, United States ABSTRACT: In this work, we have carried out Brownian dynamics simulation to describe the detailed characteristics of conformational and rheological responses to startup shear. In addition to the evaluation of the contour length of the primitive chain, the state of chain entanglement and the first normal stress difference as a function of strain, three methods are applied to determine the time-dependent shear stress. Our results show significant stretching of the primitive chain up to many Rouse times, followed by retraction, as the primary origin of stress overshoot for deformation rates lower than the reciprocal of the Rouse time but higher than the reciprocal of the reptation time. The analysis of such results reveals heterogeneous local chain stretching, demonstrating the coupling between stretching and orientation that extends to times considerably longer than the Rouse time. Explicit comparison between the simulation and the theoretical description from the GLaMM theory shows marked differences. For example, the simulation indicates a slower decline in the number of original entanglements than that at equilibrium up to many Rouse times whereas the GLaMM theory predicts a faster decrease. Moreover, contrary to the simulation that depicts a nearly constant slope in the stress−strain relationship during startup shear, the GLaMM theory shows an immediate and precipitous strain softening.
I. INTRODUCTION Deformation and flow of entangled polymers are a fundamental problem in polymer science.1,2 Unlike small-molecule liquids, fluids made of chainlike polymeric molecules involve complicated intermolecular interactions, which affect their structure and dynamics and determine the structure-processing-property relationship in their applications. In the past four decades, extensive theoretical,3−10 experimental,11−22 and simulation23,24 studies have been performed to understand the nonlinear responses of entangled polymers to fast large deformations. The objectives have been to determine the rheological behaviors in both transient and steady states and to depict the structure of polymers in terms of the states of chain orientation, stretching, and entanglement as a function of strain and strain rate. The most notable and widely accepted molecular theory for entanglement dynamics is the reptation/ tube model, as shown in Figure 1. Originally proposed by de Gennes3 to describe the dynamics of entangled polymers under equilibrium conditions, the model was developed by Doi and Edwards (DE)1,4−7 into a full description of linear and nonlinear rheology of entangled polymer melts and concentrated solutions. Over the years, improvements and refinements have been made to the basic DE theory, including explicit description of chain stretching25−28 and additional relaxation mechanisms such as contour length fluctuations (CLF),29,30 constraint release (CR),31−33 and convective constraint release (CCR).34−38 The modified tube model provides a reasonable © XXXX American Chemical Society
Figure 1. Schematic of the tube model, which reduces (a) an intractable interacting many-chain problem to (b) an effective singlechain in a smooth confining tube.
and appealing description of the equilibrium and near equilibrium dynamics of entangled polymer melts and concentrated solutions.39,40 Some of the predictions based on the tube model are also in apparent agreement with macroscopic rheological measurements for large deformations.10,11,41−44 The existence of a confining tube along a test chain as envisioned by the tube model has been confirmed by computer simulation for systems at equilibrium.45−49 Despite these successes, the physical basis for extending the tube model to large deformation is not obvious. While representing the many-chain interactions (i.e., entanglements) by an effective smooth confining tube on a tagged chain seems reasonable under equilibrium conditions, it is not clear to what Received: November 4, 2014 Revised: May 14, 2015
A
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
shear deformation. This decoupling leads to a number of consequences concerning the chain conformation and rheological properties in this rate regime, for instance, (i) the contour length of the primitive chain and stress for step shear should decline substantially on time scales of τR, (ii) the continuous orientation of polymer chain should result in monotonic increase of the radius of gyration during startup shear, and (iii) stress overshoot in this regime results from chain alignment in the confining tube by the shear deformation. In the simplest approximation, the independent alignment approximation (IAA), the correlation between the different “bond” vectors is neglected, i.e., ⟨ui⃗ ·u⃗j⟩ = 0 for i ≠ j, where ui⃗ and uj⃗ are the unit vectors of the ith and jth tube segmental “bond”. For deformation rate such that γ̇ ≪ τR−1, according to DE model,1 the bond length remains unchanged, while the unit vector u⃗i undergoes a uniform rotation by the strain tensor E to the new unit vector
extent the change in the many-chain polymer structure under large deformation can still be represented by deformation of an effective single-chain in a smooth barrier-less tube. Some recent publications have theoretically discussed the weakness of the tube model50−52 and have experimentally demonstrated some discrepancies with the tube model. For example, the time for the onset of time-strain superposibility in experiments has been shown to be much longer than the Rouse time anticipated based on the tube model.2,20−22 Also, in a study by Hassager and co-workers,53 involving extensional deformations of a highly entangled polyisoprene melt, the transient tensile stress difference is significantly above that predicted by the tube model at stretch rates below the reciprocal Rouse time τR. Using particle-tracking velocimetry, Shi-Qing Wang’s group has reported a number of nonlinear rheological phenomena54−56 that suggested the need to establish a new theory to explain these phenomena. The theoretical implications of the reported phenomena have been controversial.52,57−61 More recently, new Brownian dynamics (BD) computer simulations62−64 have emerged to indicate some discrepancies between the tubemodel and the molecular behavior revealed by the BD simulation. Specifically, for startup shear for deformation rates lower than 1/τR but higher than the reciprocal of the reptation time τd, it was found (a) there is considerable chain stretching as measured by the mean-square radius of gyration and the contour length of the primitive chain up to many Rouse times and (b) stress overshoot is primarily due to chain retraction after considerable stretching rather than chain orientation. These findings are in stark contrast to the key features of the tube theory. In this article we further elucidate the molecular characteristics concerning the conformational and rheological responses to startup shear. Our more comprehensive analyses of the new BD simulation results (based on a system of a longer chain length N = 800) have produced additional insights into the nature of the difference between the tube model and the physical picture demonstrated by the simulation. In section II, we describe the pertinent theoretical background that motivates our work. In section III, we describe the model and simulation method. Results and discussion are provided in section IV, where we examine a host of conformation and rheological properties, such as the contour length of the primitive chain, evolution of entanglements, the shear stress and the first normal stress, tube and monomer segmental orientation, coupled effect of chain stretching and orientation on the stress evolution, inhomogeneity of the tube segment length, and the initial slope of stress−strain curve. In section V, we draw some conclusions and provide our overview.
ui⃗ ′ = (E·ui⃗ )/|E·ui⃗ |
(1)
where E = (1,γ,0; 0,1,0; 0,0,1) for simple shear. This results in an unchanged overall mean-square end-to-end distance ⟨Rend2⟩, the mean-square radius of gyration ⟨Rg2⟩, and the contour length of the primitive chain L from their equilibrium values, ⟨Rend2⟩0, ⟨Rg2⟩0 and L0. The stress, apart from the inconsequential contribution from the hydrostatic pressure, is calculated by using the following expression from the primitive chain1,10 intra σαβ =
3G0 Zl0
2
∫0
Z
∂R α ∂Rβ ∂s ∂s
ds (2)
where G0 = ckBT, with c the number of tube segments per unit volume given by c = ρ/Ne, Z = N/Ne is the number of entanglements per chain, l0 is the Kuhn length for the primitive chain (i.e., the size of a tube segment), α or β = x, y, z, and s is the curvilinear coordinate along the primitive chain. The DE model ignores the small amount of chain stretching for WiR = γ̇τR < 1, so the stress is given by DE σαβ ≈ 3G0Sαβ
(3)
where Sαβ is the same-point correlation function for the unit tangent u⃗(s) and is given by Sαβ = Z −1
∫0
Z
ds ⟨uα(s)uβ(s)⟩
(4)
To account for the finite chain stretching for WiR > 1, Pearson et al.16 introduced a decoupling approximation in which the total stress is taken as the orientational part multiplied by a stretching factor (L/L0)2, i.e.
II. THEORETICAL BACKGROUND In order to delineate the objectives and implications of our simulation, it is useful to briefly review the tube theory, including its physical picture, basic assumptions, and predictions.1,2 The tube model simplifies the complex, manybody effects of chain interpenetration and entanglements into a smooth tube-like confinement on a test chain and assumes that the test chain undergoes Rouse dynamics inside the tube. For sufficiently long chains, this simplified physical picture leads to decoupling of chain stretching, which is assumed to occur on time scales of τR, from chain orientation, which occurs on time scales of τd. Specifically, for external deformational rate γ̇ that is higher than τd−1 but lower than τR−1, the tube model perceives little chain stretching, but only chain orientation during the
dec σαβ = 3G0Sαβ(L /L0)2
(5)
The IAA clearly leads to an underestimate of ⟨Rg ⟩. For a more meaningful comparison, we also compare the simulation data with the most sophisticated version of the tube model, the GLaMM (Graham, Likhtman and Milner, McLeish) model,10 which is intended to describe entangled linear chains in the full range of deformation rates from linear to strongly nonlinear flows. The relaxation mechanisms of reptation, CLF, CCR, and retraction are formulated into a stochastic microscopic evolution equation for the dynamics of the space curve describing the tube contour, R⃗ (s,t), where s is a continuous variable labeling the distance along the tube contour and t is 2 64,65
B
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
set to be unity, thus serving as the respective scales for energy and length. The mass m of a monomer is also set to 1. In the energy unit of ϵ, the temperature is kept at kBT = 1. Strictly speaking, the friction should be defined relative to the average moving velocity at the location of interest. However, the average velocity of monomers is about 5 orders of magnitude smaller than the thermal velocity in this work. Furthermore, we have also checked our results using a DPD thermostat that respects momentum conservation and obtained results in agreement with our BD results to within the error bars (of less than 1%). The initial states are prepared by lattice Monte Carlo69 sampling allowing chain crossing, which considerably reduces the time required for system equilibration. We then turn off the lattice grid and run the BD simulation in continuous space to further equilibrate the system for a time on the order of τd. Conformation properties, such as the mean square internal distances ⟨(ri⃗ − rj⃗ )2⟩/|i − j| are evaluated.70 For large |i − j|, ⟨(ri⃗ − rj⃗ )2⟩/|i − j| reaches a plateau value of 1.637 as shown in Figure 2, corresponding to a persistence length lp = 1.319 and
time. From this equation, a partial differential equation for the tube tangent tensor correlation function f(s, s′,t), averaged over all chains f(s,s′,t) =⟨(∂R⃗ (s,t))/(∂s)(∂R⃗ (s′,t))/(∂s′)⟩, can be derived, allowing the contour length of the primitive chain,65 Z L=∑s=1 (Trf(s,s))1/2, the stress tensor,10 σ = ((12Ge )/ Z (5Z))∫ 0 f(s,s) ds + σRouse to be computed for any deformation. Here Ge in ref 10 is related to G0 by a factor of 4/5, i.e., Ge = 5G0/4, and σRouse is the stress contribution from Rouse modes, which can be neglected under the condition τd−1 < γ̇ < τR−1. The full equations and its derivations are described in detail in ref 10. In this work, we choose the same parameters as in ref 65, i.e., αd = 1.15, cν = 0.1, and 9 s = 2.0 for the contour length fluctuations, constraint release and retraction terms, respectively. For additional comparison, we also list the results for the affine deformation. In the case of affine shear deformation, the contour length of the primitive chain L = L0(1 + (γ2/3))1/2, the shear stress σxy = G0γ and the first normal stress N1 = σxx − σyy = G0γ2.
III. MODEL AND SIMULATION METHOD We perform our simulation of entangled polymer melts using BD on the LAMMPS platform. We use the standard Kremer− Grest model45 for our system consisting of M chains with N monomers in a box with volume V. The total number of monomers NM is 2.115 × 105 for N = 500 and is increased to 4.28 × 105 for N = 800. All the monomer pairs interact via a truncated and shifted Lennard-Jones potential defined by ⎡⎛ d ⎞12 ⎛ d ⎞6 ⎛ d ⎞12 ⎛ d ⎞6 ⎤ ULJ(rij) = 4ϵ⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ + ⎜⎜ ⎟⎟ ⎥ ⎢⎣⎝ rij ⎠ ⎝ rij ⎠ ⎝ rc ⎠ ⎝ rc ⎠ ⎥⎦
(6)
where rij is the distance between monomers i and j, and rc denotes the cutoff distance for the interaction chosen to be 21/6d corresponding to pure repulsion between the monomers. We choose the monomer number density to be ρ = NM/V = 0.85d−3. The consecutive monomers on each chain are connected by a finitely extensible nonlinear elastic (FENE) potential of the form ⎡ ⎛ rij ⎞2 ⎤ 1 UFENE(rij) = − kcr0 2 ln⎢1 − ⎜ ⎟ ⎥ 2 ⎣ ⎝ r0 ⎠ ⎦
Figure 2. Mean square internal distances for entangled polymer chains of length N = 500 and N = 800.
the chains have reached equilibrium before turning on shear. All equilibrium properties including the diffusion coefficients are in agreement with known results from Kremer and Grest.45 In the quiescent state, the number of monomers in an entanglement strand is Ne ≈ 36, as determined from the crossover time (from t1/2 to t1/4 scaling) in the monomer meansquare displacement, in agreement with the known literature value.45 The Rouse time τR is determined from τR = 2⟨Rg2⟩0/ (π2DR),1 where DR is the Rouse self-diffusion constant obtained by extrapolation of simulation data of nonentangled polymers. The reptation or disengagement time τd is obtained by the formula τd =2⟨Rg2⟩0/(π2Ds),71 where Ds is the self-diffusion constant obtained by extrapolation of simulation data of shorter entangled polymers.1 For N = 800, the numerical values of the parameters ⟨Rg2⟩0, Ds, DR, τR, and τd are respectively 210, 2.58 × 10−6, 8.12 × 10−5, 5.24 × 105, and 1.65 × 107. We have also estimated the Rouse time by the location of the crossover (from t1/4 to t1/2 scaling) in the monomer mean-square displacement, with the result τR = 7.0 × 105. To be consistent with our earlier publications, we use the Rouse time determined from the center-of-mass self-diffusion. Other measures of the Rouse time are possible and may lead to quantitative differences in the results; the interpretation of our results is
(7)
2
where kc = 30ϵ/d and r0 = 1.5d. The bonded pairs in the same chain thus interact via the sum of the two potentials, U(rij) = ULJ(rij) + UFENE(rij), which has a deep minimum at r = 0.96d. In our simulations the actual average bond length between two monomers on the chain is 0.97d. The additional potential to prevent chain crossing66 used in ref 62 turned out to be unnecessary and was not included since under the strain rates used in the simulation, the probability of chain crossing induced by shear is negligibly small.67,68 For this model, the BD equation of motion for monomer i is 1 1 ri (⃗ ̇ t ) = − ∇Ui + fi ⃗ (t ) μ μ
(8)
Here Ui is the total potential, μ is the friction coefficient, and fi⃗ (t) is a random force related to μ by the fluctuation dissipation theorem ⟨fi⃗ (t)·fj⃗ (t′)⟩ = δijδ(t − t′) 6kBTμ. The equations are integrated with a time step Δt = 0.001τLJ and the friction coefficient is taken to be μ = 0.5τLJ−1, where τLJ is the unit of time defined as τLJ = (mσ2/ϵ) 1/2. The parameters ϵ and d are C
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
chain L with increasing strain. The contour length of the primitive chain is inherently a tube-model concept; however, the tube model does not prescribe a method to identify the primitive chain in terms of the monomer positions. The different PPA methods in the literature67,74−76 are essentially different definitions of the primitive path, and there are no clear criteria of which one is a more faithful representation of the primitive chain in the tube model. The primitive chain obtained by the quenching PPA methods is not a faithful representation of the actual in situ chain conformation. For example, it cannot be used to compute the stress using the intrachain stress expression in ref 1. We have recently proposed an effective way to define the primitive path through a coarse-graining process. The primitive path as envisioned in the original work of Doi and Edwards7 was in essence a coarse-grained view of the polymer chain with a coarse-graining scale of Ne, i.e., the primitive chain is a polymer made of N/Ne coarse-grained segments. We thus define an average contour length of the primitive chain by L = ⃗ ⃗ ⃗ ∑Z−1 s=1 ⟨|Rs+1 − Rs|⟩, where Rs is the center of mass position of a tube segment, obtained by coarse-graining the chain into Z = N/Ncg tube segments each containing Ncg consecutive monomers; here we choose Ncg = Ne. At equilibrium, L = L0. This definition allows a natural discretization scheme to compute the chain orientation and stress. We note this coarse-graining procedure is similar in spirit to a time averaging using time steps of order τe, as first proposed by Read et al.,42 and further developed by Milner and co-workers.77,78 The results for the normalized contour length L/L0 are shown in Figure 3. The nonmonotonic change of the contour length in
subjected to such numerical uncertainties. However, since the validity of the tube model has been well tested experimentally and by computer simulation, we expect the different definitions of the Rouse time to be consistent with each other, i.e., to yield similar numerical values. For most of the results reported in this article, we use periodic boundary condition in the x and z directions and Lees−Edwards boundary condition in the y direction.72 We have verified that the velocity profile is linear in the gradient direction at all stages of the simulation. The data reported are results of averaging over 16 independent samples for N = 500 and 8 independent samples for N = 800. Two different shear rates γ̇ = (1/6)τR−1 (i.e., Rouse-Weissenberg number WiR ≡ γ̇τR = 1/6) and γ̇ = (1/4)τR−1 (WiR = 1/4) have been used. The dimensions of the simulation box are chosen to be Lx=Ly = 6Rg0, and Lz = 4Rg0, where x, y, and z refer to the flow, gradient and vorticity directions, and Rg0 = ⟨Rg2⟩01/2 is the equilibrium radius of gyration. To ensure that our results are not due to finite size effects, we have also carried simulations using two larger box sizes (9 × 6 × 4) and (6 × 9 × 4), where the dimensions are in unit of Rg0. We find that chain conformation properties vary within 3% as the box size varies. Therefore, the larger-box simulations validate the simulation results from the box of Lx = Ly = 6Rg0, and Lz = 4Rg0 and further demonstrate that there is little confinement effect. We have also explored an alternative means to produce startup shear in the BD simulation. Specifically, we have generated shear by displacing the upper wall with a fixed velocity in the x-direction relative to the stationary bottom wall. Monomers within a distance less than d from the two walls are permanently adhered to the walls. Periodic boundary conditions are only imposed in the x- and z-directions. We calculate all the properties using only the free chains in the box. The data reported are results of averaging over 32 independent samples. To quantify the evolution of entanglements during deformation, we adopt the method of primitive path analysis (PPA) introduced by Sukumaran et al.67 PPA applies a simultaneous contour reduction scheme to all chains of a computational sample, such that chain crossing is avoided and the topology of the system of chains is conserved. Operationally, we freeze the conformation of all the chains in the system, and apply the PPA with energy minimization to examine the number of chain strands that hinders the lateral motion of a test chain at any given time. This number is averaged over all chains and is defined as the number of entanglements per chain. This PPA method67 gives an average number of segments between entanglements Nseg ≈ 72 ≈ 2Ne. The factor of 2 arises from different definitions of Ne, which has been discussed in ref 73. Other than setting the modulus G0, the actual values of Ne are not important (within some reasonable range) and do not affect the qualitative conclusions of our work. A difference from previous PPA is that, in our method, we track the identity of the chains that form the entanglements around a given test chain. This allows us to examine the evolution of the survival probability of the original entanglements P, which is obtained by normalizing the number of surviving original entanglements for a tagged chain by the equilibrium value of the number of entanglements per chain.
Figure 3. Evolution of the normalized contour length of the primitive chain as a function of strain. The inset represents the normalized contour length using three different values of Ncg.
the simulation results indicates significant chain stretching followed by retraction. The pronounced peaks, spanning many Rouse times, offer clear evidence of chain stretching and retraction on time scales much longer than τR. In contrast, the GLaMM theory predicts again nearly imperceptible peaks for WiR = 1/6 and 1/4, which occur respectively at γ = 2.64 and 2.66. Thus, the simulation result of significant stretching followed by retraction on time scales of many Rouse times is not captured by the tube model in this rate regime. We have also included in Figure 3 the curve representing affine deformation. Here, the initial increase of the chain contour
IV. RESULTS AND DISCUSSION A. Contour Length of the Primitive Chain. We first examine the evolution of the contour length of the primitive D
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
entanglement network actually tightens up initially, making it more difficult to relax, compared to the loss of entanglements by diffusion in quiescence. The tightening of the entanglements may be related to the hairpins formed by the polymer strands, which have been shown to be central to the theory of stress relaxation in stars.81 While the detailed molecular mechanism for the tightening is currently being investigated, the initial tightening up is consistent with the nearly affine deformation observed for the contour length of the primitive chain (see Figure 3); the evolution of the survival probability relative to its equilibrium behavior provides useful insight into why and how chain deformation starts being affine and eventually becoming nonaffine as a function of the time or strain. C. Stress Response. The stress response can be evaluated in several ways. Below we present the stress vs strain relationship based on three separate methods. The first two are based on simulations using the Lees−Edwards boundary condition, and the third involves the wall-driven simulations, to be discussed respectively as follows. a. Startup Shear with Lees−Edwards Boundary Condition. The shear stress σxy is calculated first using the well-known microscopic expression:1 1 micro σαβ = − ∑ ⟨Fikαrikβ⟩ V k ,i (9)
length from simulation is seen to follow closely the affine deformation curve. Obviously different choices for Ncg or different protocols for defining the primitive chain will yield different values for L. However, qualitatively and statistically they reveal the same information about chain stretching. In the inset of Figure 3, we calculate the normalized contour length using three different values of Ncg ranging from 25 to 50. The feature of a pronounced peak with a width spanning several strain units (and hence many Rouse times) remains unchanged. In fact the peak becomes more pronounced with increasing Ncg (e.g., if we use the Ne ≈ 72 obtained from the PPA). We have also computed the contour length using two alternative methods. One method simply adds the radius of gyration rge of the Ne ≈ 36 monomers, so that L/L0 = rge/rge0. The other method involves locating and connecting the center of mass of the Ne monomers by sliding along the chain going from i = Ne/2 to i = N − Ne/2 such that a nearly smooth curve is obtained. The results for L/L0 from these alternative methods are very similar to those in Figure 3. B. Evolution of Entanglements. Characterizing chain entanglement and elucidating its dynamic evolution under external deformation for entangled polymers is critical to understanding the molecular mechanism of their structure and rheology.79,80 Chain entanglement involves localized intermolecular interactions that may be perceived as network junctions. To describe its dynamic evolution, we have to keep track of the dynamics of these intermolecular coupling junctions, popularly known as entanglement points. A useful global measure of the dynamics of the entanglements is to evaluate how the average number of entanglements per chain evolves with deformation. For startup shear, this number decreases monotonically with the elapsed strain, indicating that deformation causes some degree of disentanglement. Another useful quantity is the average number of surviving original entanglements Z′, which, when normalized by the number of entanglements at equilibrium (i.e., at t = 0) Z0, gives the survival probability of the original entanglements, P = Z′/Z0. As shown in Figure 4, this probability exhibits a slower decline under shear than at equilibrium up to many Rouse times whereas the GLaMM theory predicts a faster decrease according to the differential equation dZ′/dt = −Z′ν, where ν is the constraint release rate obtained by Table I in ref 10. This suggests that the
where Fkiα is the α-component of the total force acting on the monomer i of the kth chain located at rki⃗ . As shown in Figure 5,
Figure 5. Shear stress vs strain during startup shear for the entangled polymeric chain length of N = 800 under the condition of WiR < 1 and Wi > 1.
the tube theory considerably under-predicts the stress compared with the simulation data. Consistent with the experimental observation,15−19 the stress overshoot and the peak strain from the simulation increase appreciably with shear rate. In contrast, the stress calculated from GLaMM shows modest increase with shear rate, but the location of the stress maximum remains nearly independent of the shear rate. For smaller shear rates, GLaMM even predicts the peak strain to decrease with shear rate.12 Second, we can evaluate the stress by the intrachain stress of the primitive chain.1,10 Using the construct discussed in Section A, the shear stress is calculated as
Figure 4. Evolution of the survival probability of the original entanglements P as a function of time in units of τR. E
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules intra σαβ =
3G0 (Z − 1)l0 2
Z−1
deformation is indeed driven by the displacement of the walls, (2) to have a method to study interfacial phenomena such as wall slip, and (3) to compute the stress directly from the force on the walls, which provides an independent, unambiguous means to determine the stress. The shear stress in this case is given straightforwardly as
∑ ⟨(R⃗s + 1 − R⃗s)α (R⃗s + 1 − R⃗s)β ⟩ s=1
(10)
where R⃗ s is the center of mass of the Ncg = Ne ≈ 36 monomers in the unit. As shown in Figure 5, the stress curves produced by eqs 9 and 10 are in good agreement. This agreement implies that the coarse-graining method to define the primitive chain and to evaluate the corresponding shear stress is both adequate and effective. In particular, it confirms the choice of Ne = 36 as the natural coarse-graining unit for defining the primitive chain. In the tube theories, including GLaMM, the dominant contribution to stress comes from orientation in the WiR < 1 regime;63 there is little contribution from chain stretching. On the other hand, by separately examining the orientational and the remaining contributions to stress from the simulation data, we have recently demonstrated that stretching is a significant contribution to stress and that the stress overshoot is primarily due to chain retraction after considerable stretching rather than chain orientation.63 Our new simulation results with N = 800 confirms and strengthens this conclusion. In addition, we find that the universal quantitative master curve in terms of the normalized stress and normalized strain is insensitive to the value of the parameter WiR (as long as WiR < 1) and the chain length; see Figure 1b in ref 63. The agreement between the master curve obtained in simulation and that from experiment indicates that the physical origin for the overshoot in both cases is the same, i.e., chain stretching and subsequent retraction, contradictory to the picture of chain orientation in the tube model. We have also calculated the first normal stress difference N1 = σxx − σyy from eq 9. Consistent with the experimental observation,17,18 a faint overshoot can be seen. The peak value of N1 and the strain at the peak both increase with shear rate as shown in Figure 6. Also consistent with experiment, the
⟨Fx⟩ (11) A where A is the area of wall, and Fx is the x-component of the total force acting on it. Figure 7 shows favorable comparison σxywall =
Figure 7. Comparison of the shear stress from the wall driven method and Lees−Edwards method for the chain length of N = 500. Note that the inset shows the contour length of the primitive chain.
among the three curves, evaluated respectively according to eqs 9 and 11. This demonstrates the ability of the wall-driven shear method to depict the stress response to startup shear. Using this method, the evolution of the contour length is also found to show similar characteristic to that illustrated with the Lees− Edwards boundary condition, as shown in the inset of Figure 7. D. Tube and Monomer Segmental Orientation. Since the tube model attributes stress overshoot primarily to chain alignment in the shear direction under the condition of WiR < 1 and Wi > 1, we examine the behavior of chain orientation. By discretizing the expression on the right-hand side of eq 4, we can obtain Sxy from simulation using Sxy =
1 Z−1
Z−1
∑ s=1
(R⃗ s + 1 − R⃗ s)x (R⃗ s + 1 − R⃗ s)y (R⃗ s + 1 − R⃗ s)2
(12)
As shown in Figure 8a, there is a nearly imperceptible peak for WiR = 1/6, and a weak peak for WiR = 1/4, the peak strain for both being significantly larger than 2.1 predicted by the DE tube model. Thus, for the shear rates used in our study, the evolution of chain orientation does not follow the affine orientational deformation as assumed by the DE tube model. However, the segmental orientation at the monomer level, S′xy = ⟨uxuy⟩ does show a clear overshoot, with peak strains around 2, as seen from Figure 8b. This behavior is consistent with an overshoot of the contour length of the primitive chain: since the tube is made up of blobs of essentially unstretchable monomer segments (for the weak shear considered in the simulation), lengthening of the tube segments must be due to alignment of the monomer segments.
Figure 6. First normal stress difference vs strain during startup shear for the entangled polymeric chain length of N = 800 under the condition of WiR < 1 and Wi > 1.
overshoot of N1 appears at a significantly higher strain than that corresponding to the overshoot of the shear stress. On the other hand, N1 calculated from GLaMM is much lower and shows almost no peak. b. Wall-Driven Startup Shear. Three considerations motivate our use of wall displacement to produce shear: (1) to more directly mimic the shear experiment in which F
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
Figure 9. Comparison of stress during startup shear from the simulation and decoupling approximation. For comparison the corresponding decoupling curve where the initial slope is forced to 1 is also included as the open blue circles. The schematic shows nonuniform stretching of the “bonds”.
simulation data and results calculated from eq 5, where the orientation correlation function and the contour length are obtained from the simulation. Clearly, the decoupling captures some qualitative features of the stress overshoot, but the quantitative agreement is rather poor. In particular, the initial slope is considerably lower. From symmetry considerations, it can be shown that L/L0 is an even function of strain, so the slope in the decoupling approximation is given by the orientational part, which is 3/5. In order to satisfy the requirement of linear viscoelasticity, an argument could be made that, with the decoupling approximation, G0 should be consistently evaluated from the stress calculated with this approximation; i.e., for sufficiently small γ, σxy should coincide with G0γ, thus forcing the initial slope to be 1 in Figure 9. The result of this normalization is shown as the blue open circles. While the initial behavior has been made to agree better with the simulation data, this approach leads to a significant overestimate for the peak stress and steady state stress. Basically, when the stretching of the bonds is inhomogeneous, Z∑Zi=1li2 >(∑Zi=1 li)2 where li is the length of the ith tube segment. Thus, the primary origin of the discrepancy is the inhomogeneity in stretching along the primitive chain; i.e., with a distribution in the “bond” lengths of the tube segments, as illustrated in the inset of Figure 9, it is not a good approximation to factorize out an overall contour length. We now examine this distribution. F. Distribution of the Bond Length of the Primitive Chain. In Figure 10, we show the statistical distribution of the length of tube segments for WiR = 1/6. The equilibrium distribution is Gaussian (multiplied by 4πl2, given by the black curve for γ = 0), but at the peak of the stress overshoot (γ = 1.75), the distribution shifts to larger l; see the blue curve. Also, the distribution noticeably deviates from Gaussian, as measured by a significant non-Gaussian parameter Γ = 3⟨l4⟩/(5⟨l2⟩2) − 1 = 0.31, reflecting the heterogeneous nature of the stretching of the tube segments. This heterogeneity makes it a poor approximation to decouple the orientation and stretching of the chain, because the stress depends nonlinearly on the deformation of each strand of the chain. However, as the strain increases, the distribution of the length of tube segments gradually returns to the Gaussian distribution due to the retraction of the stretched chain, which makes decoupling the orientation and stretching of the chain a better approximation.
Figure 8. (a) Tube and (b) monomer segmental orientation vs strain during startup shear for the entangled polymeric chain length of N = 800 under the condition of WiR < 1 and Wi > 1.
Masubuchi and Watanabe82 recently questioned the validity of our results and conclusions reported in ref 63 using the argument that monotonic behavior reported there83 in the tube segmental orientation implies monotonic behavior in the monomer segmental orientation. By assuming Sxy ′ to be proportional to the shear birefringence, they asserted that the MD simulation results violate the stress-optical rule (SOR). Figure 8 suggests such argument is not valid, i.e., there is no one-to-one correspondence between the monomer level orientation and the tube level orientation. It is also worth mentioning that the ratio of the monomer segmental orientation to the stress in our work is not a constant. Thus, if one would interpret the monomer segmental orientation ⟨uxuy⟩ as being strictly proportional to the experimentally measured optical birefringence, there would be an apparent violation of SOR. The origin of such an apparent violation is not clear and is beyond the scope of this work; we suspect that the issue lies at whether a coarse-grained model such as the Kremer−Grest model can meaningfully be employed to compute the optical birefringence. We note that similar apparent violations can be found in earlier simulations84−87 upon close examination of the data. E. Evaluation of the Decoupling Approximation. The stress eq 5 introduced by Pearson et al.16 was intended to describe deformation rate WiR>1. It includes both chain orientation and stretching, but assumes that these two contributions factorize. It is of interest to see to what extent this decoupling approximation holds in describing the actual stress−strain behavior. Figure 9 shows comparison between the G
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules
predicting the correct initial modulus, but 1/5 of the stress stored in the tube is relaxed after τR by longitudinal modes,88 and the mechanism of chain stretching and retraction in the GLaMM theory is due to a combined effect of convection, CLF and CR (CCR), all taking place from the very beginning of the deformation, with the chain undergoing Rouse motion inside a smooth tube.
V. CONCLUSIONS AND OVERVIEW In summary, we have performed BD simulation under the conditions WiR < 1 and Wi > 1 aimed at verifying the molecular mechanisms envisioned by the tube model during startup shear. We find that (1) our previous results62−64 remain valid for longer chains, larger boxes and different boundary conditions; (2) the orientation and stretching are coupled on time scale considerably longer than the Rouse time for our chain lengths and shear rates; (3) the distribution of the length of tube segments becomes more non-Gaussian around the stress overshoot, reflecting dynamic heterogeneity, and thus a simple meanfield treatment cannot accurately describe the stress response; (4) the initial modulus is in agreement with experiments and close to affine deformation but is in strong disagreement with GLaMM which predicts a precipitous drop; (5) a slowing down of the evolution of initial entanglements is observed during the early stage of startup shear relative to the equilibrium behavior whereas GLaMM predicts accelerated evolution at all times. These clear differences between our simulation results and predictions from the GLaMM model suggest that the basic assumption of a smooth, barrier-less confining tube and spontaneous chain retraction on Rouse time in the tube theories is problematic for startup shear for the chain lengths and shear rates we have studied. In particular, the slower initial decline in the survival probability of original entanglements under shear than at equilibrium up to several Rouse times is a direct contradiction to the relaxation mechanisms envisioned by the tube theories, all of which would predict a faster decline under shear than at equilibrium because of the barrier-less chain retraction. Our BD simulation studies are clearly just the first step toward an explicit elucidation of molecular pictures behind the various remarkable rheological phenomena. Many further studies remain to be carried out in order to obtain a more complete account of the rheological behavior of entangled polymers under large deformation. Examples include, (a) the molecular origin of the 1/3 power scaling for the value of the peak stress or strain in the regime of WiR > 1 in startup shear;13,19,89 (b) the molecular mechanism for elastic yielding during relaxation from large step extension;90−92 and (c) conformational and corresponding rheological responses to large amplitude oscillatory shear.93−95 As our knowledge further accumulates through such computer simulation studies, it can be expected that a more realistic microscopic theory might emerge to better describe the known and unknown nonlinear rheological properties of entangled polymers. We hope that the present and additional (future) simulation results will help identify the necessary ingredients for an improved theoretical description.
Figure 10. Distribution of the length of tube segments in different stages during startup shear at WiR = 1/6.
In other words, we expect decoupling between stretching and orientation to be a reasonable approximation at steady state; this speculation remains to be validated by further examinations. G. The Initial Modulus. To address the fundamental question of when the chain deformation ceases to be affine, i.e., when the transition from dominantly elastic deformation to irrecoverable deformation occurs, we examine the slope of the stress−strain curve (the modulus) in Figure 11. The horizontal
Figure 11. Shear modulus vs strain. The experimental data are taken from ref 19; the shear rate is 0.3 s−1, and the Rouse time is based on τR = τd/(3Z), giving WiR ≈ 0.18.
black line represents the value of the slope for strictly affine deformation. The modulus from our simulation as well as from the experimental data starts from 1 in agreement with the affine deformation, and then gradually decreases, as the deformation is increased. The behavior of the GLaMM prediction is interesting. It starts with the correct initial slope equal to 1 but rapidly drops to ∼0.8. This difference corresponds to different pictures for what is happening to chains during startup shear. In the simulation, chain stretching is caused by a nearly affine deformation of the entangled network, and this affine deformation persists due to the initial “tightening” of the network shown in Figure 4; retraction takes place only when chains undergo significant disentanglement. The GLaMM theory includes Rouse modes and chain stretching, thus
■
AUTHOR INFORMATION
Corresponding Author
*(L.A.) E-mail:
[email protected]. H
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules Notes
(35) Ianniruberto, G.; Marrucci, G. Proc. Int. Congr. Rheol. 2000, 2, 102. (36) Marrucci, G. J. Non-Newtonian Fluid Mech. 1996, 62, 279. (37) Likhtman, A. E.; Milner, S. T.; Mcleish, T. C. B. Phys. Rev. Lett. 2000, 85, 4550. (38) Milner, S. T.; McLeish, T. C. B.; Likhtman, A. E. J. Rheol. 2001, 45, 539. (39) Watanabe, H. Prog. Polym. Sci. 1999, 24, 1253. (40) Mcleish, T. C. B. Adv. Phys. 2002, 51, 1379. (41) Graham, R. S.; et al. Macromolecules 2006, 39, 2700. (42) Read, D. J.; Jagannathan, K.; Likhtman, A. E. Macromolecules 2008, 41, 6843. (43) Sukumaran, S. K.; Likhtman, A. E. Macromolecules 2009, 42, 4300. (44) Likhtman, A. E. J. Non-Newtonian Fluid Mech. 2009, 157, 158. (45) Kremer, K.; Grest, G. S. J. Chem. Phys. 1990, 92, 5057. (46) Kremer, K.; Grest, G. S.; Carmesin, I. Phys. Rev. Lett. 1988, 61, 566. (47) Smith, S. W.; Hall, C. K.; Freeman, B. D. Phys. Rev. Lett. 1995, 75, 1316. (48) Kröger, M.; Hess, S. Phys. Rev. Lett. 2000, 85, 1128. (49) Qin, J.; Milner, S. T.; Stephanou, P. S.; Mavrantzas, V. G. J. Rheol. 2012, 56, 707. (50) Likhtman, A. E.; Graham, R. S. J. Non-Newtonian Fluid Mech. 2003, 114, 1. (51) Larson, R. G. J. Polym. Sci., Part B: Polym. Phys. 2007, 45, 3240. (52) Wang, S.-Q. J. Polym. Sci., Part B: Polym. Phys. 2008, 46, 2660. (53) Nielsen, J. K.; Hassagera, O.; Rasmussen, H. K.; McKinley, G. H. J. Rheol. 2009, 53, 1327. (54) Tapadia, P.; Wang, S.-Q. Phys. Rev. Lett. 2006, 96, 016001. (55) Wang, Y. Y.; Boukany, P.; Wang, S.-Q.; Wang, X. R. Phys. Rev. Lett. 2007, 99, 237801. (56) Ravindranath, S.; Wang, S.-Q.; Olechnowicz, M.; Quirk, R. Macromolecules 2008, 41, 2663. (57) Wang, S.-Q.; Ravindranath, S.; Wang, Y. Y.; Boukany, P. E. J. Chem. Phys. 2007, 127, 064903. (58) Adams, J. M.; Olmsted, P. D. Phys. Rev. Lett. 2009, 102, 067801. (59) Wang, S.-Q. Phys. Rev. Lett. 2009, 103, 219801. (60) Adams, J. M.; Olmsted, P. D. Phys. Rev. Lett. 2009, 103, 219802. (61) Wang, S.-Q.; Wang, Y. Y.; Cheng, S. W.; Li, X.; Zhu, X. Y.; Sun, H. Macromolecules 2013, 46, 3147. (62) Lu, Y. Y.; An, L. J.; Wang, S.-Q.; Wang, Z.-G. ACS. Macro. Lett. 2013, 2, 561. (63) Lu, Y. Y.; An, L. J.; Wang, S.-Q.; Wang, Z.-G. ACS. Macro. Lett. 2014, 3, 569. (64) Lu, Y. Y.; An, L. J.; Wang, S.-Q.; Wang, Z.-G. Macromolecules 2014, 47, 5432. (65) Graham, R. S.; Henry, E. P.; Olmsted, P. D. Macromolecules 2013, 46, 9849. (66) Kumar, S.; Larson, R. G. J. Chem. Phys. 2001, 114, 6937. (67) Sukumaran, S. K.; Grest, G. S.; Kremer, K.; Everaers, R. J. Polym. Sci., Part B: Polym. Phys. 2005, 4, 917. (68) Kröger, M.; Loose, W.; Hess, S. J. Rheol. 1993, 37, 1057. (69) Shaffer, J. S. J. Chem. Phys. 1994, 101, 4205. (70) Auhl, R.; Everaers, R.; Grest, G. S.; Kremer, K.; Plimpton, S. J. J. Chem. Phys. 2003, 119, 12718. (71) Lodge, T. P. Phys. Rev. Lett. 1999, 83, 3218. (72) Lees, A. W.; Edwards, S. F. J. Phys. C 1972, 5, 1921. (73) Zhou, Q.; Larson, R. G. Macromolecules 2006, 39, 6737. (74) Zhou, Q.; Larson, R. G. Macromolecules 2005, 38, 5761. (75) Kröger, M. Comput. Phys. Commun. 2005, 168, 209. (76) Tzoumanekas, C.; Theodorou, D. N. Macromolecules 2006, 39, 4592. (77) Bisbee, W.; Qin, J.; Milner, S. T. Macromolecules 2011, 44, 8972. (78) Cao, J.; Qin, J.; Milner, S. T. Macromolecules 2015, 48, 99. (79) Green, M. S.; Tobolsky, A. V. J. Chem. Phys. 1946, 14, 80. (80) Astarita, G.; Marrucci, G. Principles of Non-Newtonian Fluid Mechanics; McGraw-Hill Inc.: New York, 1974. (81) Levine, A. J.; Milner, S. T. Macromolecules 1998, 31, 8623.
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work is supported, in part, by the National Natural Science Foundation of China (Nos. 21120102037 and 21334007) and is further subsidized by the Special Funds for National Basic Research Program of China (No. 2012CB821500). We thank the anonymous reviewers, whose comments have helped improve the presentation of our work.
■
REFERENCES
(1) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Clarendon: New York, 1986. (2) Dealy, J. M.; Larson, R. G. Structure and Rheology of Molten polymers From Structure to Flow Behavior and Back Again; Hanser Gardner: Cincinnati, OH, 2006. (3) de Gennes, P.-G. J. Chem. Phys. 1971, 55, 572. (4) Doi, M.; Edwards, S. F. J. Chem. Soc., Faraday Trans. 2 1979, 75, 38. (5) Doi, M.; Edwards, S. F. J. Chem. Soc., Faraday Trans. 2 1978, 74, 1818. (6) Doi, M.; Edwards, S. F. J. Chem. Soc., Faraday Trans. 2 1978, 74, 1802. (7) Doi, M.; Edwards, S. F. J. Chem. Soc., Faraday Trans. 2 1978, 74, 1789. (8) Sussman, D. M.; Schweizer, K. S. Phys. Rev. Lett. 2011, 107, 078102. (9) Sussman, D. M.; Schweizer, K. S. Macromolecules 2012, 45, 3270. (10) Graham, R. S.; Likhtman, A. E.; Mcleish, T. C. B.; Milner, S. T. J. Rheol. 2003, 47, 1171. (11) Collis, M. W.; et al. J. Rheol. 2005, 49, 501. (12) Auhl, D.; Ramirez, J.; Likhtman, A. E.; Chambon, P.; Fernyhough, C. J. Rheol. 2008, 52, 801. (13) Boukany, P. E.; Wang, S.-Q.; Wang, X. R. J. Rheol. 2009, 53, 617. (14) Hernadez, A. R.; Detcheverry, F. A.; Peters, B. L.; Chappa, V. C.; Schweizer, K. S.; Muller, M.; de Pablo, J. J. Macromolecules 2013, 46, 6287. (15) Menezes, E. V.; Graessley, W. W. J. Polym. Sci., Polym. Phys. Ed. 1982, 20, 1817. (16) Pearson, D. S.; Kiss, A. D.; Fetters, L. J.; Doi, M. J. Rheol. 1989, 33, 517. (17) Osaki, K.; Inoue, T.; Isomura, T. J. Polym. Sci., Part B: Polym. Phys. 2000, 38, 1917. (18) Osaki, K.; Inoue, T.; Uematsu, T. J. Polym. Sci., Part B: Polym. Phys. 2000, 38, 3271. (19) Ravindranath, S.; Wang, S.-Q. J. Rheol. 2008, 52, 681. (20) Archer, L. A. J. Rheol. 1999, 43, 1555. (21) Sanchez-Reyes, J.; Archer, L. A. Macromolecules 2002, 35, 5194. (22) Inoue, T.; Uematsu, T.; Yamashita, Y.; Osaki, K. Macromolecules 2002, 35, 4718. (23) Baig, C.; Mavrantzas, V. G.; Kröger, M. Macromolecules 2010, 43, 6886. (24) Cao, J.; Likhtman, A. E. Phys. Rev. Lett. 2012, 108, 028302. (25) Marrucci, G.; Grizzuti, N. Gazz. Chim. Ital. 1988, 118, 179. (26) Pearson, D.; Herbolzheimer, E.; Grizzuti, N.; Marrucci, G. J. Polym. Sci., Part B: Polym. Phys. 1991, 29, 1589. (27) Mead, D. W.; Leal, L. G. Rheol. Acta 1995, 34, 339. (28) Mead, D. W.; Yavich, D.; Leal, L. G. Rheol. Acta 1995, 34, 360. (29) Doi, M. J. Polym. Sci., Polym. Phys. Ed. 1983, 21, 667. (30) Milner, S. T.; McLeish, T. C. B. Phys. Rev. Lett. 1998, 81, 725. (31) Milner, S. T. J. Rheol. 1996, 40, 303. (32) Rubinstein, M.; Colby, R. H. J. Chem. Phys. 1988, 89, 5291. (33) Viovy, J. L.; Rubinstein, M.; Colby, R. H. Macromolecules 1991, 24, 3587. (34) Ianniruberto, G.; Marrucci, G. J. Non-Newtonian Fluid Mech. 1996, 65, 241. I
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
Article
Macromolecules (82) Masubuchi, Y.; Watanabe, H. ACS Macro. Lett. 2014, 3, 1183. (83) For the shorter chain length N = 500 studied in ref 63, the tube orientation shows monotonic dependence on the strain for WiR = 1/6 but the monomer orientation shows a clear peak similar to that in Figure 8b. (84) Cao, J. Ph.D. Thesis; see Figure 3.30: http://www.reading.ac. uk/web/FILES/maths/Thesis_Jing_final.pdf. (85) Kröger, M.; Luap, C.; Muller, R. Macromolecules 1997, 30, 526. (86) Yamamoto, R.; Onuki, A. J. Chem. Phys. 2002, 117, 2359. (87) Yamamoto, R.; Onuki, A. Phys. Rev. E 2004, 70, 041801. (88) Likhtman, A. E.; McLeish, T. C. B. Macromolecules 2002, 35, 6332. (89) Divoux, T.; Barentinb, C.; Manneville, S. Soft Matter 2011, 7, 9335. (90) Wang, S.-Q.; Ravindranath, S.; Boukany, P.; Olechnowicz, M.; Quirk, R. P.; Halasa, A.; Mays, J. Phys. Rev. Lett. 2006, 97, 187801. (91) Boukany, P. E.; Wang, S.-Q.; Wang, X. Macromolecules 2009, 42, 6261. (92) Moorcroft, R. L.; Fielding, S. M. Phys. Rev. Lett. 2013, 110, 086001. (93) Ravindranath, S.; Wang, S.-Q. J. Rheol. 2008, 52, 341. (94) Li, X.; Wang, S.-Q.; Wang, X. J. Rheol. 2009, 53, 1255. (95) Hyun, K.; et al. Prog. Polym. Sci. 2011, 36, 1697.
J
DOI: 10.1021/ma502236m Macromolecules XXXX, XXX, XXX−XXX
