Quantitative predictions of the linear viscoelastic properties of entangled polyethylene and polybutadiene melts via modified versions of modern tube models on the basis of atomistic simulation data

Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

Contents lists available at SciVerse ScienceDirect

Journal of Non-Newtonian Fluid Mechanics journal homepage: http://www.elsevier.com/locate/jnnfm

Quantitative predictions of the linear viscoelastic properties of entangled polyethylene and polybutadiene melts via modified versions of modern tube models on the basis of atomistic simulation data Pavlos S. Stephanou a,⇑, Vlasis G. Mavrantzas b,⇑ a b

Department of Mathematics and Statistics, University of Cyprus, PO Box 20537, 1678 Nicosia, Cyprus Department of Chemical Engineering, University of Patras & FORTH-ICE/HT, Patras, GR 26504, Greece

a r t i c l e

i n f o

Article history: Available online 28 April 2013 Keywords: Modified tube model Leygue et al. model Dual constraint model Molecular dynamics Primitive path Multi-scale modeling

a b s t r a c t We present a hierarchical, three-step methodology for predicting the linear viscoelastic properties of entangled polymer melts. First, atomistic trajectories accumulated in the course of long molecular dynamics simulations with moderately entangled polymer melts are self-consistently mapped onto the tube model to compute the segment survival probability function w(s, t) for primitive paths. Extracted directly from the atomistic simulations, the computed w(s, t) accounts for all possible dynamic mechanisms affecting chain motion in entangled polymers such as reptation, contour length fluctuation, and constraint release. In a second step, the simulation predictions for w(s, t) are compared with modern versions of the tube model, such as the dual constraint model of Pattamaprom et al. and the Leygue et al. model; the comparison reveals ways through which the two models can be improved and parameterized on the basis of the direct molecular simulation data. The key parameters turn out to be the entanglement chain length Ne and the entanglement time se, both of which can be reliably extracted from the simulations. In a third step, the modified versions of the two models are invoked to predict the linear viscoelastic properties of the polymer under study over a broad range of molecular weights. The power of the new methodology is illustrated here for the case of linear polyethylene (PE) and cis- and trans-1,4 polybutadiene (PB) melts for which atomistic molecular dynamics data have already been obtained recently. We present results from the new approach for the zero-shear-rate viscosity g0, and the storage G0 and loss G00 moduli of the three polymers as a function of their molecular weight (MW), and a direct comparison with experimentally measured rheological data. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Dynamics and rheology in entangled polymeric liquids is governed by topological constraints known as entanglements which dominate molecular forces and prevent chains from crossing one the other. Despite the complexity of the corresponding interactions, the pioneering works of Edwards, de Gennes and Doi [1–6] demonstrated that one can account for uncrossability constraints with a simple but powerful mean-field theory known as reptation. The reptation theory is built on the concept of the tube model defining a space around the chain within which this is allowed to execute motion. Over the years, the tube model coupled with mechanisms such as reptation, constraint release, contour length fluctuations and arm retraction in the case of branched molecules ⇑ Corresponding authors. Tel.: +357 22 893911; fax: +357 22 892601 (P.S. Stephanou), tel.: +30 6944 602580; fax: +30 2610 965223 (V.G. Mavrantzas). E-mail addresses: [email protected] (P.S. Stephanou), vlasis@chemeng. upatras.gr (V.G. Mavrantzas). URL: http://www.lstm.chemeng.upatras.gr (V.G. Mavrantzas). 0377-0257/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jnnfm.2013.04.003

has developed to a quantitative theory for the description of the key rheological properties of polymers, especially in the linear viscoelastic regime [7]. Entangled melts of monodisperse and polydisperse linear polymer chains, of star and H-shaped polymers, and of long- and short-chain branched chains can all be described rather quantitatively in the framework of such a model. Reptation involves a geometric reduction from the phase space of the detailed chain to that of its primitive path (PP). The PP is the shortest path that connects the two ends of the chain (which are considered as fixed in space) and follows the main chain contour without violating any topological constraints along the chain from one end to the other. The following two quantities are then introduced: the diameter dt of the tube (effectively representing the strength of the topological constraints on chain dynamics) and the contour length L of the PP. From a mathematical point of view, the most fundamental quantity in the theory is the tube segment probability function w(s, t) [1,2,5]. This is defined as the probability that a segment s along the PP contour will remain inside the initial tube after a time t; it can also be considered as the average probability that the corresponding PP segment has fully relaxed

112

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

the stress initially imposed on it after a time t. If this function is known, all linear viscoelastic (LVE) properties of the system can be obtained in a rather straightforward manner [1]. For chain reptation in the space of fixed obstacles, w(s, t) obeys a typical diffusion equation with appropriate initial and boundary conditions. For chain reptation in a polymer melt, one should further account for the fact that neighboring chains are not fixed obstacles but they move too, thereby continuously destroying and creating new entanglements. The chain itself also changes its conformation as it diffuses through the melt, and this causes fluctuations in its contour length. Thus, for the tube model to quantitatively describe the dynamics of real polymer melts, it has to account further for mechanisms such as contour length fluctuation (CLF) [8] due to the ‘‘breathing’’ motion of chain ends and constraint release (CR) related to the destruction of topological constraints due to the relaxation of chains surrounding the reference chain (or the tube itself) [9]. It is thus not surprising that, since the 1970s when the reptation picture and the tube concept were introduced, the tube model has undergone numerous modifications [7], most of them motivated either by the need to bring model predictions into a quantitative comparison with experimental data for model systems or by the need to describe more complicated chain architectures. Two such typical models are the dual constraint model of Pattamaprom et al. [10] and the Leygue et al. model [11]. These models can indeed offer quantitative predictions of directly measured rheological data for a variety of polymer architectures, but often this comes with the need to assume non-realistic values for some of the fundamental model parameters. We mention, for example, the value of the entanglement time se (strictly defined as the time after which a polymer chain starts to feel the tube constraints) typically chosen to fit experimental data [10,12,13]. For polyethelene (PE) at 177 °C and polybutadiene (PB) at 140 °C, se is computed from direct molecular dynamics (MD) simulations to be about 3 and 2–4 ns, respectively [12–16]. By making use of the temperature shift factor aT reported by Raju et al. [17] for PE and Colby et al. [18] for PB, the corresponding values are 3.75 ns for PE at 190 °C and 44–65 ns for PB at 25 °C which differ appreciably from the values (7 ns for PE at 190 °C and 350 ns for PB at 28 °C) needed for modern tube models to fit experimental data [10]. A remarkable development recently has been the computation of the function w(s, t) describing the survival probability of PP segments directly from atomistic molecular dynamics (MD) data [12– 16]. The methodology (which accounts directly for CLF and CR mechanisms) first maps long atomistic MD trajectories onto time trajectories of primitive chains and then documents primitive chain motion in terms of a curvilinear diffusion in a tube-like region around the coarse-grained chain contour. In the calculations, the effective tube diameter is independently evaluated by observing the effect of tube constraints either on atomistic displacements or on the displacement of primitive chain segments orthogonal to the initial PP. The approach accounts both for chain reptation longitudinally inside the constraining tube as well as for local transverse fluctuations driven from constraint release and regeneration mechanisms, the latter causing parts of the chain to venture outside the average tube surface for certain periods of time. The new approach opened the way to bridging the outcome of direct MD simulations with tube models, since both can be compared on the basis of their prediction for the tube segment probability function w(s, t) of the same system [12–16]. Furthermore, it allows one to utilize the results of the direct PP analysis in order to propose modifications to the models [16]. Our effort in this paper is to capitalize on this work and develop a systematic methodology that will enable the direct computation of the LVE properties of high-MW polymer melts by properly trans-

ferring information from lower-level (atomistic MD) simulations to closed-form tube models having the form of a reaction–diffusionlike partial differential equation. Our paper has been structured as follows: In the following section we provide a concise summary of the modifications (see Ref. [16] for more details) made to the tube models (the dual constraint and the Leygue et al. ones) employed here to describe the rheological properties of long polymer melts. The corresponding original versions of these models have been presented in the Appendix of Ref. [12]. Here, we only note that, as far as the dual constraint model is concerned, wherever in that Appendix we had used 1  cs, here we replace it by s, because in the original model s 2  1c ; 1c with c = 2 for a linear chain, whereas here and in Refs. [12–16] s 2 ½0; 1. In Section 3 we provide some details of the numerical methodology used to solve the two models along with the values of the parameters employed in each one of them and of what experimental data we compared their predictions with. The paper proceeds with Section 4 where the major results are presented and concludes with Section 5 discussing the most important conclusions and future directions. 2. Brief account of the modified dual constraint and Leygue et al. models 2.1. The modified dual constraint model In the modified dual constraint or Pattamaprom et al. model, the reaction–diffusion equation in solved in two stages. In the first stage, the equation incorporates only chain reptation and CLF effects but omits CR (i.e., it initially considers chain motion in a fixed tube):

@  1 @2  1 w ðs; tÞ ¼ 2 w ðs; tÞ   w ðs; tÞ; @t p sd @s2 s ðsÞ where



ð1aÞ

b

searly ðsÞ ¼ se þ 3p4c3=22 Z s2 sR ; slate ðsÞ ¼ sc2R exp½U  ðsÞ;

ð1bÞ



U ðsÞ ¼ 32 Zc s2 ; with the value of the exponent b ¼ ð2bÞ1 taken from the scaling /ðtÞ  t b when se < t < sR , and c = 2 for a linear chain. The partial differential equation (PDE) Eq. (1a) is subject to the following initial and boundary conditions:

w ðs; t ¼ 0Þ ¼ 1 w ðs ¼ 0; tÞ ¼ w ðs ¼ 1; tÞ ¼ expðt=se Þ

:

ð1cÞ

The newly introduced boundary condition in the last equation is obtained, self-consistently, directly from (1a) by considering the case of times t  sd for which the reptation term in (1a) represented by the second derivative may be safely neglected [16]. Then, one obtains as a solution the function w ðs; tÞ ¼ exp½t=sðsÞ. Since close to chain ends early fluctuations prevail, sðsÞ ¼ searly ) sð0Þ ¼ sð1Þ ¼ se , thus w ð0; tÞ ¼ w ð1; tÞ ¼ expðt=se Þ. Such a time dependence has indeed been observed in simulations [12–16]. The expressions for the early- and late-time CLF effects have been borrowed from the works of Doi [8] and Milner and McLeish, [19] respectively (see also the original papers of Pattamaprom et al. [10] for more details). s ðsÞ, on the other hand, is given by

8 for s < C 1 ; > early ðsÞ :  slate ðsÞ for s < C 2 ;

ð1dÞ

where C 1 denotes the segment position close to chain ends corresponding to the first crossover of searly to slate and C 2 the segment position for the second crossover of searly to slate deeper inside the

113

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

tube. The overall tube survival probability is obtained by integrating R1 w ðs; tÞ over s, i.e., as / ðtÞ ¼ 0 dsw ðs; tÞ, and yields the average ⁄ probability U (t):

U ðtÞ ¼



/ ðtÞ

if / ðtÞ > /R ðtÞ;

/R ðtÞ

if / ðtÞ < /R ðtÞ:

ð1eÞ

The meaning of Eq. (1e) is that in order to obtain U⁄(t) one needs to compare / ðtÞ with the following approximate Rouse process [10,20]:

 1=2 t /R ðtÞ ¼ / ðt0 Þ ; t0

ð1fÞ

where t0 is the time after which / ðtÞ starts to decrease faster than t 1=2 . In the second stage, the computed U ðtÞ function is employed in the expression for the activation energy UðsÞ for chain retraction in deep fluctuations, slate ðsÞ, to account for the CR effect. The equations corresponding to Eqs. (1a) and (1b) are

[11]. Proving that indeed the first line of Eq. (3c) does solves Eq. (3a) is a simple exercise. The expression for the curvilinear diffusion coefficient Dc ðsÞ employed in the model [16] is a modification of that proposed by Doi [8,11]:

8 1 2 d > 2 t þX 2 s2 > hLi > > > se > 2 þ4Zs4 > > < p2 sd   2 ~ c ðsÞ  p2 sd Dc ðsÞ ¼ 1d D > 2hLit þX 2 ð1sÞ2 > > > > > 2p2sesd þ4Zð1sÞ4 > > : 1

ð2aÞ

and



3p3=2 Z 2 s 4c2

searly ðsÞ ¼ se þ slate ðsÞ ¼ sc2R exp½UðsÞ; UðsÞ ¼

3 Z 2 c



b

sR ; ð2bÞ

2

U ðtÞs ;

respectively. PDE Eq. (2a) is complemented with the same initial and boundary conditions as in the first stage. In this second stage of the calculations of the Pattamaprom et al. model, searly is given by the same expression as in the first stage, Eq. (1b) above, and the same holds true for s(s) which is obtained in a way completely analogous to that in the first stage:

8 for s < C 1 ; > early ðsÞ : slate ðsÞ for s < C 2 :

ð2cÞ

Despite that the parameters C1 and C2 denote still the crossover points, their values are now time-dependent, since in the second stage U ðtÞ is included in slate . Solving Eq. (2a) for wðs; tÞ allows one to compute eventually the R1 overall probability function WðtÞ ¼ 0 dswðs; tÞ. 2.2. The Leygue et al. model The modified Leygue et al. model is represented by the following modified reaction–diffusion equation:

 @wðs; tÞ @ @wðs; tÞ þ cðt; aÞwðs; tÞ; Dc ðsÞ ¼ @t @s @s

ð3aÞ

with

R1

cðt; aÞ ¼ a

0

  @ ds @s Dc ðsÞ @wðs;tÞ @s d ln WðtÞ ¼a ; R1 dt wðs; tÞds 0

ð3bÞ

subject to the initial and boundary conditions given by Eq. (1c) R1 [here WðtÞ ¼ 0 wðs; tÞds]. Eq. (3a) with the form for cðt; aÞ prescribed by Eq. (3b) leads to a solution for w(s, t) of the form:

wðs; tÞ ¼ w0 ðs; tÞ

Z 0

1

dsw0 ðs; tÞ

a

 @w0 ðs; tÞ @ @w0 ðs; tÞ ; ¼ Dc ðsÞ @t @s @s

; ð3cÞ

indicating that the choice a = 1 corresponds to double reptation [21]. Here w0 (s, t) is the solution of the diffusion equation without the last term in Eq. (3a), i.e. without CR (or, equivalently, a = 0)

ð3dÞ for s > sCLF ; otherwise;

with X and a undetermined, numerical (i.e., fitting) parameters. An interesting point to make here is that in the limit Z  1 we obtain the original Leygue et al. model. The value of the constant sCLF is dic~ c ðsCLF Þ ¼ 1, which leads to tated by the condition that D

2

@ 1 @ 1 wðs; tÞ ¼ 2 wðs; tÞ  wðs; tÞ @t p sd @s2 sðsÞ

for s < sCLF ;

sCLF

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 !ffi3 u u 2 u1 4 2 u 2 s 1 d e t 5: X þ tX 4  16Z 2  ¼t 8Z p sd 4 hLi2

ð3eÞ

Our proposition in Ref. [16] is to adopt the solution prescribed by Eq. (3c) but with a time-dependent dynamic dilution exponent, namely:

aðtÞ ¼ a0 ð1  et=sR Þ;

ð3fÞ

with a0 ¼ 1 (double reptation [21]) or a0 ¼ 4=3 [22]. In this case, the expression for cðt; aÞ is modified to

cðt; aÞ ¼ a



d ln WðtÞ d ln aðtÞ ; þ ln WðtÞ dt dt

ð3gÞ

although this is not needed (since the form of the solution, Eq. (3c), is still applicable with a being time dependent). In the following, the original versions of the dual constraint [10] and Leygue et al. [11] models will be referred to simply as the original dual constraint and original Leygue et al. models, respectively. 3. Numerical details, simulation results, and experimental data We have integrated the partial differential equations representing original and modified versions of the two models using a standard Crank–Nicolson finite-difference scheme [23] with a space step Ds = 5 103 and a time step Dt = 107 (made dimensionless with p2sd). For the parameters entering the modified versions of the two models we have used values directly obtained from the simulations (this is the case for the shorter samples that are amenable to direct MD simulation; see Sections 4.1.1 and 4.2.1) or we have made use (this is the case for the longer samples, whose long time dynamics cannot be accessed by direct atomistic MD simulations; see Sections 4.1.2 and 4.2.2) of the following relationships of the reptation theory between the various characteristic times: sd ¼ 3Z 3 ð1  XZ 0:5 Þse and sR ¼ Z 2 se [5,8], where Z = N/Ne. The latter scaling is also supported by simulation data [24,25]. To obtain the predictions of the original versions of the two models, we have simply used sd ¼ 3Z 3 se (the usual case in the literature). For the Leygue et al. model (both original and modified versions) we have further chosen a0 = 4/3, a value that was found to provide a slightly better agreement with the simulation data [16]. Should we had chosen a0 = 1, we would have needed a slightly smaller value for se for the two sets of data to coincide, and this was our sole criterion in choosing a0 = 4/3 compared to the double reptation value. As a result, our work does not provide a definitive answer as to which of the two values is the most appropriate (since their predictions are very comparable). But what is more important is the fact that the dynamic dilution exponent is a time-dependent quantity.

114

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

We have also chosen X = 0.6 for PE and X = 0.5 for cis-1,4-PB and trans-1,4-PB, as dictated from the predictions of the atomistic MD data for the dependence of the reptation time on the degree of chain entanglement (see Fig. 1). Note here that these values do not differ significantly from those obtained from an analysis of contour length fluctuations (see Fig. 2 of Ref. [12]). In the case of PB which exhibits a larger entanglement length, we noticed that the two shortest cis-1,4 systems (characterized by Z = 2.22 and Z = 2.78) and the shortest trans-1,4-PB system (characterized by Z = 2.69) exhibit longer reptation times than those suggested by the scaling sd ¼ 3Z 3 ð1  XZ 0:5 Þse , which is not surprising since these systems are only marginally entangled. We therefore decided not to present any rheological data for these systems from our methodology, which (strictly speaking) is only valid for truly entangled polymers. The values of all relevant quantities for the tube models adopted here are reported in Table 1. A remark worth making is that the values of se in Table 1 (obtained by monitoring the mean-square displacement of the most inner segments in the atomistic simulations [14]) are consistent with the values obtained through se ¼ ðfN e a2s Þ=ð3p2 kB TÞ [1]. The Rouse time of the PE systems was estimated on the basis of the scaling sR  N2 by using sR data extracted from MD simulations with several shorter (unentangled) PE melts [24–26]; a similar procedure had been followed in a previous study [12]. For trans-1,4-PB, the Rouse times were estimated from those for cis-1,4-PB by making use of the theoret2 ical expression sR ¼ ðfN 2 bb Þ=ð3p2 kB TÞ for the Rouse time [1], implying that sR ðcisÞ=sR ðtransÞ ¼ fðcisÞ=fðtransÞ for the same number of monomers N, where f denotes the monomeric friction coefficient. This allowed us to obtain the Rouse times for all studied PB-trans systems corresponding to the same backbone length by rendering simulation data for f reported by Tsolou and Mavrantzas [27] onto simulation data for sR. We further mention that, although a recent theoretical study [28] has shown that e = 2 in the expression as ¼ edt between the PP step length as and the tube diameter dt for Z  1, we chose e = 1 (as is customarily done in practice). Then, L ¼ ðZ þ 1Þas ðZ þ 1Þdt . A problem we faced in actually mapping MD data onto the dual constraint model was the unavailability of data describing the variation of the exponent bð¼ ð2bÞ1 Þ in the power-law dependence, /ðtÞ  tb , of the segmental mean-square displacement (msd) on time t for times se < t < sR with chain length Z. We overcame this by proposing the following phenomenological expression for the function b ¼ bðZÞ:

b¼

2 ; 1 þ Z n

Z P 1;

ð4Þ

on the basis of the following considerations: (a) it satisfies the two limits, namely b = 1 (i.e., b = 1/2) for Z = 1 (unentangled systems) and b = 2 (i.e., b = 1/4) for very long chains (Z  1) [1,5] and (b) it satisfactorily describes the data obtained from the atomistic MD simulations by choosing n = 0.56 for PE [Fig. 2a], n = 0.68 for trans-1,4-PB, and n = 0.29 for cis-1,4-PB [Fig. 2b]. Actually, and consistent with the data presented in Fig. 1b, the reported n values for the two types of PB systems were obtained by disregarding in the fittings the two shortest cis-1,4-PB systems (characterized by Z = 2.22 and 2.78) and the shortest trans-1,4-PB one (characterized by Z = 2.69). Having solved for w(s, t), one can easily compute all LVE properties of the polymer melt under study by making use of the following expressions of the reptation theory [1]:

WðtÞ ¼

Z

1

wðs; tÞds;

ð5Þ

0

g0 ¼

Z

1 0

GðtÞdt ¼ G0N

G0 ðxÞ ¼ x

Z

1

0

G00 ðxÞ ¼ x

Z 0

1

Z

1

WðtÞdt;

ð6Þ

0

GðtÞ sinðxtÞdt ¼ G0N x

Z

GðtÞ cosðxtÞdt ¼ G0N x

1

WðtÞ sinðxtÞdt;

ð7aÞ

0

Z

1

WðtÞ cosðxtÞdt;

ð7bÞ

0

with G0N denoting the plateau modulus. Our methodology for computing the linear viscoelastic data of entangled polymer melts entails therefore three main steps: Step 1: Map atomistic MD data onto the tube model and extract the PP segment probability function w(s, t). Analyze also the short time dynamics of the simulated systems to obtain se. A further topological analysis allows one to also compute Ne. Step 2: Choose a reliable tube model and parameterize it; if possible, propose modifications to the mathematical form of the model or to its boundary conditions so that (at least for moderately entangled melts) it can satisfactorily reproduce the computed w(s, t) curves for times ranging from se up to several sd values. Step 3: Employ the modified (improved) tube model to actually carry out computations for truly long MW samples of the same polymer and calculate its relevant rheological properties (e.g., g0, G0 and G00 ).

Fig. 1. Reptation time (made dimensionless with 3Z 3 se ) as obtained directly from the MD simulations and its dependence on Z = N/Ne. The values of the entanglement length Ne have been calculated from topological reduction algorithms [14]: for PE, Ne = 76 ± 11, while for cis-PB and trans-PB, we employ Ne = 144 ± 19 and Ne = 119 ± 15, respectively. The entanglement time has been set equal to se = 2.5 ns for PE and equal to 3 ns for both PB-cis and PB-trans. Also shown is how sd =3Z 3 se compares with the 1 quantity 1  XZ 2 .

115

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

Fig. 2. (a) Dependence of the quantity bð¼ ð2bÞ1 Þ on Z for PE (Ne = 76). b is the exponent in the power-law /ðtÞ  tb describing the dependence of the segmental mean-square displacement on time t, for times se < t < sR . The results have been obtained from the atomistic MD simulations presented in Refs. [12,14]. Also plotted is Eq. (4) with n = 0.56. (b) Same as with (a) but now for cis-PB and trans-PB. Also plotted is Eq. (4) with n = 0.68 for trans-1,4-PB and n = 0.29 for cis-1,4-PB.

Table 1 List of quantities needed in the modified and original Leygue et al. and dual constraint tube models: the entanglement time se, the Rouse time sR, the disentanglement time sd, the tube diameter dt, the average contour length L, and the mean square displacement (msd) parameter b. Typical values suggested either directly or indirectly from atomistic MD simulations for three different polymers at several conditions (T = 450 K for the PE and T = 413 K for the PB systems) are also shown. Interested readers are directed to Refs. [12–16] for more details. System

se (ns)

sR (ns)

sd (ns)

PE 320 PE 400 PE 500 PB-trans 400 PB-trans 600 PB-cis 600 PB-cis 800

2.1 ± 0.5 2.9 ± 0.4 2.8 ± 0.4 3.1 ± 0.5 3.2 ± 0.5 2.1 ± 0.8 2.3 ± 0.5

51 ± 8 316 ± 32 33 ± 4 80 ± 14 489 ± 25 33 ± 7 125 ± 20 1042 ± 46 32 ± 6 101 ± 30 355 ± 25 26 ± 3 210 ± 65 1204 ± 50 29 ± 3 77 ± 11 500 ± 70 30 ± 2 138 ± 20 1255 ± 50 31 ± 2

dt (Å)

L (Å)

b

162 ± 4 198 ± 6 251 ± 8 127 ± 5 185 ± 5 154 ± 6 202 ± 9

1.355 ± 0.022 1.468 ± 0.055 1.520 ± 0.069 1.514 ± 0.055 1.477 ± 0.058 1.179 ± 0.020 1.231 ± 0.023

In this work, such a methodology has been applied to three different polymers: linear PE, cis-1,4-PB, and trans-1,4-PB, for which detailed and very long (up to several microseconds) atomistic MD simulations have been realized recently, for chain lengths extending from short (or unentangled) to moderately entangled ones [24–26,29]. To appreciate the improvement over the original versions of the two models, we thoroughly compared our predictions for the three melts with experimentally measured data. For PE, we compared against data for the zero-shear viscosity reported by Pearson et al. [30,31] and by Harmandaris et al. [26] all shifted to 177 °C , and against data for the loss moduli obtained by Raju et al. [17] shifted to 175 °C. At this temperature, the corresponding value of the plateau modulus was set equal to G0N ¼ 2:3 MPa [31]. The entanglement time and entanglement length obtained directly from the simulations [12–16] at 450 K were: Ne = 76 ± 11 and se

2.4 ± 0.8 ns (i.e., between 2 and 3 ns). Experimental measurements yield (443 K) Ne = 82 (Me = 1150 g/mol [32] which agrees with the simulations results; here Me = NeM0 where M0 is equal to 14 g/mol). The value of the plateau modulus, in particular, comes in agreement with the value of 2.17 MPa extracted by using G0N ¼ 45 qMRTe ¼ 45 Mq0RTNe [1] with q = 0.77 g/cm3 and Ne = 76 at T = 450 K [14]. For the corresponding comparisons against rheological data for the two PB polyisomers, we have made use of new experimentally obtained zero-shear viscosity and storage and loss moduli data, properly shifted to 140 °C. For the theoretical calculations, data for se and Ne were borrowed directly from the simulations [12–16]: Ne = 144 ± 19 and se 2.2 ± 0.7 ns for cis-PB and Ne = 119 ± 15 and

se 3.2 ± 0.6 ns for trans-PB, both at 140 °C . The corresponding value for G0N was obtained by using G0N ¼ 45 3

qRT Me

¼ 45

qRT M0 Ne

[1]. At

T = 413 K = 140 °C, q = 0.86 g/cm and Ne = 144 from Ref. [14], yielding a plateau modulus equal to 1.215 MPa for cis-PB. For an experimental sample with composition corresponding practically to pure cis-1,4-PB (cis/trans/vinyl:96.5/1.9/1.6), the corresponding literature values at 298 K are [33,34]: G0N ¼ 0:76 MPa, q = 0.9 g/cm3 and Ne = 174 (Me = 2347 g/mol). By shifting to 413 K using:

G0N ðTÞ G0N ðT 0 Þ

¼

qðTÞ T Me ðT 0 Þ ; qðT 0 Þ T 0 Me ðTÞ

ð8Þ

we obtain G0N ¼ 1:216 MPa, which agrees amazingly well with the computed value. For a PB sample with about 50% trans content (50/40/10 trans/cis/vinyl), the corresponding experimental values at 140 °C are [33]: Ne = 134 (Me = 1815 g/mol) and G0N ¼ 1:25 MPa. For later use, we note that a sample primarily comprised of atactic 1,2-polybutadiene (98% vinyl content) is characterized by a value of Me = 3850 g/mol and G0N ¼ 1:216 MPa at 300 K [33,35]. A simple comparison with the sample that is practically pure cis-1,4-PB at 298 K indicates that the vinyl content causes an increase in both the entanglement length and the plateau modulus. The same trend is expected for the sample containing 50% trans (cis/trans/vinyl:40/ 50/10) [35]. The experimental data reported in Fig. 19 and 24–27 but also in the next Figures of this paper have been generously provided by the Vlassopoulos group. They were obtained from linear viscoelastic measurements in the small amplitude oscillatory mode using an ARES strain-controlled rheometer (TA Instruments). The experimental conditions were identical to those used for other polybutadiene samples [36]. The polybutadienes used were obtained from Polymer Source, Canada, and had predominantly a 1,4-microstructure. The largest molar mass (165000 g/mol) was taken from the literature [37]. Further details will be published in the future [38]. 4. Results and discussion 4.1. Polyethylene 4.1.1. Model predictions for moderately entangled melts In Table 1 we report the values of the set of quantities (sR, sd, L, dt, b) needed in the numerical computations with the modified/improved versions of the two models. These were obtained either directly or indirectly from the atomistic MD simulations and in Table 1 are given as a function of the molecular length of the melt studied (PE, cis-1,4-PB, trans-1,4-PB), and have been used as such

116

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

in the calculations with the modified versions of the two tube models. In contrast, in the calculations with the original versions of the two models, we have followed the common practice in the field, namely to use that sR ¼ Z 2 se and sd ¼ 3Z 3 se . In Figs. 3 and 4, original and modified versions of the Leygue et al. model are compared on the basis of their capability to follow the simulation data for the function w(s, t), for Ne = 76 and se = 2.5 ns. The improvement due to the modifications incorporated in the original model is remarkable, especially at early-to-moderate times. In this time scale, the modified version provides a considerably better prediction of the actual simulation data for the survival probability w(s, t) than the original model. We should mention here that, at early times, the w(s, t) curves obtained from the modified version of the Leygue et al. model differ significantly from those obtained from the original version; this is mainly due to the different boundary condition imposed to the reaction–diffusion equation, Eq. (3c). At later times [see, e.g., the curves corresponding to t = 0.35sd in Fig. 3c], original and modified versions seem to exhibit a similar behavior, which should be attributed to the fact that at long enough times the CLF mechanism becomes irrelevant; an additional reason is that the dynamic dilution exponent, a = a(t), has reached its steady-state value (=a0) in the case of the modified version. As a result of these two reasons, the two versions of the model become practically equivalent. That the modified Leygue et al. model provides a better description of the simulation data for the PP dynamics becomes more clear when the comparison with the original version is made at the level of the function W(t) reported in Fig. 4. In fact, Fig. 4 shows that the modified Leygue et al. model captures the time decay of W(t) accurately irrespective of the molecular length of the melt. Another interesting point to note in the curves of Fig. 4 is the significantly longer relaxation time characterizing the decay of W(t)

in the original model; this is due to the larger reptation times incorporated in this version through the scaling sd ¼ 3Z 3 se (as compared to the smaller value for sd obtained directly from the PP analysis and assumed in the computations with the modified version). The corresponding comparison (between modified and original versions on the one hand and the direct simulation data on the other hand) for the dual constraint model is presented in Figs. 5 and 6. The calculations here have been carried out by using Ne = 76 (exactly as before) but by choosing a slightly larger value for se, equal to 3 ns. This little adjustment of the se value improves somewhat the agreement with the simulation data for w(s, t), and can be justified by the fact that the value extracted from the PP analysis comes with some uncertainty (se is computed to be between 2 ns and 3 ns). Similar to the Leygue et al. case, the modifications made to the dual constraint model improve its predictions mainly at early to moderate times [Fig. 5a–b], especially for the part of the function w(s, t) near chain ends. We had reached a similar conclusion when we studied the implications of a similar set of modifications for the PE 500 system in Ref. [16]. Here, we also observe that the predictions of the modified dual constraint model are more consistent with the simulation results even at later times [see Fig. 5c] for which the function w(s, t) approaches zero for all PP segments s in the interval [0, 1]. Also, and similar to what we observed in the case of the modified Leygue et al. model, the better quantitative agreement of the modified dual constraint model carries on to (and actually is more pronounced) the W(t) function shown in Fig. 6. Additional evidence for the improved capability of the modified Leygue et al. and dual constraint models to describe actual rheological data is provided by their predictions for the zero-shear-rate viscosity. The zero-shear rate viscosity is easily computed from

Fig. 3. Predictions of the modified and original versions of the Leygue et al. model for the segment survival probability function w(s, t) (lines) and comparison with the results of the direct topological analysis (symbols) for several monodisperse PE systems, at three different time instances: (a) t = 0.01sd (early), (b) t = 0.1sd (intermediate), and (c) t = 0.35sd (late). For the original version of the model, the usual assumption has been made that sd ¼ 3Z 3 se . For the modified version, the values of the quantities needed (sR, sd, L, dt) have been obtained directly from the MD simulations (see Table 1). For all systems: Ne = 76 and se = 2.5 ns.

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

Fig. 4. Same as with Fig. 3 but for the average tube survival probability function WðtÞ 

R1 0

117

wðs; tÞds.

Fig. 5. Same as with Fig. 3 but for the original and modified versions of the dual constraint model. For the original version we have used sd ¼ 3Z 3 se . For the modified version the values for the quantities needed (sR, sd, b) have been obtained directly from the MD simulations (see Table 1). Also, Ne = 76 and se = 3 ns.

118

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

Fig. 6. Same as with Fig. 4 but for the original and modified versions of the dual constraint model.

w(s, t) by making use of Eq. (6) of the reptation theory; the data obtained (by taking a value for the plateau modulus equal to 2.3 MPa at 175 °C, as supported from experimental observations [31]) are shown in Fig. 7. They indicate that the modifications incorporated in the two models do not dramatically affect their comparison with the experimental or the direct simulation data. But one has always to bear in mind that although the improvement at the level of g0 is minor (and could be easily adjusted by a slight change in the value of the key parameters, namely G0N and se), the improvement recorded in Figs. 3–6 at the level of the dynamic variables w(s, t) or W(t) onto which a tube model is built is remarkable. 4.1.2. Model predictions for high MW PE melts As mentioned in the Introduction, our effort in this work has primarily been to parameterize and modify currently available tube models on the outcome of direct molecular simulations and then to actually apply these models to industrially relevant samples to compute their key rheological properties without the need to carry out additional simulations at lower scales. We attempt such a levelbridging (or bottom-up) approach in this section and in Section 4.2.2 by employing the modified tube models to predict the rheological data of high-MW PE and PB samples without any further reference to atomistic-level simulations. To this, we abandon the procedure followed in the previous Section to use values for the set of the physically relevant quantities (sR, sd, L, dt, and b) obtained from the molecular simulations but instead we render them functions of only the number Z of entanglements per chain for the given polymer melt. It turns out that what one only needs to specify in this case is the entanglement length Ne and the entanglement time se, which are by default independent of the polymer molecular weight. In turn, this implies that the more accurate the calculation or estimation of these two parameters from lower-level

simulations, the more reliable the predictions of the two tube models will be. Given that Ne is a rather known (fixed) quantity for most polymers, the predictive capability of the chosen tube model will hinge upon the reliable estimation of se. Thus, we now regard the set of quantities (sR, sd, L, dt, b) as being functions of the number of entanglements Z as mentioned in Section 3. As a result, we have to choose values for only two parameters: the entanglement length Ne and the entanglement time se, independent of the molecular weight. We choose to keep the same values as the ones given above [Figs. 3–6], namely Ne = 76 and se = 2.5 ns for the Leygue et al. model, and Ne = 76 and se = 3 ns for the dual constraint model. Then, in Figs. 8–11 we depict the same comparison as in Figs. 3–6. We notice that for both the Leygue et al. model (Figs. 8 and 9) and the dual constraint model (Figs. 10 and 11) the predictions are not changed significantly. This illustrates that the general functions for the set (sR, sd, L, dt, b) capture quite well the values obtained from the molecular simulations. In addition, it proves that these general relations can now be used for the longer systems for which we have not performed, and we could hardly do so, any molecular simulations. We can therefore proceed to calculate the zeroshear-rate viscosity and the loss modulus for larger MW systems and compare the predictions against available experimental data. In Fig. 12a we evaluate the predictive capability of the Leygue et al. model (both original and modified versions) by computing the zero-shear viscosity and comparing it against available experimental data [26,30,31] and our own calculations based on the PP analysis methodology [14] for G0N ¼ 2:3 MPa, Ne = 76 and se = 2.5 ns. We notice that although the quantitative improvement is small, the qualitative one is significant: in fact, it is very encouraging that the modifications cause a shift of the computed curves to the left bringing the modified version of the model closer to the experimental data. And the same is true for the comparison of the model

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

119

Fig. 7. Predictions of the original and modified versions of the Leygue et al. (a) and dual constraint (b) models for the zero-shear-rate viscosity and comparison with experimental [26,30,31] and simulation [14] data for linear PE at 450 K. Parameter values: G0N ¼ 2:3 MPa, Ne = 76 and se = 2.5 ns in (a) and G0N ¼ 2:3 MPa, Ne = 76 and se = 3 ns in (b). In both cases, the values of the set (sR, sd, L, dt, b) needed to solve the two models are those obtained from the simulations (Table 1).

Fig. 8. Same as with Fig. 3 but now for the modified version we have used L/dt = (Z + 1). The values of Ne and se were kept the same.

sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.6, and

(using the same values for Ne and se) against the experimental data of Raju et al. [17] for the loss modulus shown in Fig. 12b. The corresponding comparison for the dual constraint model is shown in Fig. 13a. The numerical results here have been obtained by using G0N ¼ 2:3 MPa, Ne = 76, and se = 3 ns, which are basically the same as those used in the Leygue et al. model (except for se whose value has been slightly increased but still within the bounds suggested by the atomistic simulations: se 2.4 ± 0.8 ns [14]). Overall, the predictions of the modified dual constraint model are better than those of the original version, especially in the case of the moderately entangled systems (Z 2–7). Concerning the comparison for the loss modulus [37], the Leygue et al. model provides a better description than the dual constraint model. In the case of the dual constraint model, for the

less entangled SNPA-1 (Z 18) and SNPA-2 (Z 32) melts, the comparison is in favor of the modified version whereas for the more entangled SNPA-3 (Z 55) and SNPA-4 (Z 73) melts, the comparison is in favor of the original version. This observation came as a surprise to us and led us reconsider our strategy about how we should implement the CR Rouse relaxation process in the model. We found that if we allow it to take place also in the second stage, the results (shown in Fig. 14) improve for the two most entangled systems but worsen for the less entangled ones, SNPA-1 (Z 18) and SNPA-2 (Z 32). It seems that we should include the CR Rouse relaxation process in the second stage only for the larger molecular weight polymers but not for the smaller; unfortunately a crossover is not possible to be given. Alas, this is merely a conjecture and needs to be elaborated further in the future.

120

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

Fig. 9. Same as with Fig. 4, but now for the modified version we have used L/dt = (Z + 1). The values of Ne and se were kept the same.

sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.6, and

Fig. 10. Same as with Fig. 5, but now for the modified version we have used sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.6, and b ¼ 2=ð1 þ Z n Þ for the msd parameter with n = 0.56 [see Fig. 2a]. The values of Ne and se were kept the same.

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

121

Fig. 11. Same as with Fig. 6, but now for the modified version we have used sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.6, and b ¼ 2=ð1 þ Z n Þ for the msd parameter, with n = 0.56 [see Fig. 2a]. The values of Ne and se were kept the same.

Fig. 12. Comparison of the Leygue et al. model predictions (both original and modified versions) for the zero-shear rate viscosity g0 at 177 °C (a) and the loss modulus G00 (b) with the experimental data of Raju et al. [17] for linear PE at 175 °C. The molecular weights of the samples are as follows: SNPA-1 = 19.3 K, SNPA-2 = 33.9 K, SNPA-3 = 58.4 K and SNPA-4 = 77.4 K [17]. Also: G0N ¼ 2:3 MPa, Ne = 76 and se = 2.5 ns.

4.2. Polybutadiene 4.2.1. Model predictions for moderately entangled melts The corresponding set of w(s, t) curves showing the improvement in the predictive capability of the modified Leygue et al. model for the PP segment survival probability function with respect to the direct simulation data are shown in Figs. 15 and 16. The calculations have been carried out by using values for the model parameters obtained directly from the molecular simulations (Table 1). As before, in the case of the original Leygue et al. model, we have used the scalings sR ¼ Z 2 se and sd ¼ 3Z 3 se , the typical choice in the literature. For PB-cis we have used Ne = 144, se = 2.5 ns [12– 16], and X = 0.4 which is close to the value(s) suggested from the

analysis of the direct simulation data for the dependence of sd on MW reported in Fig. 1 (see also Fig. 2 in Ref. [12]). The general conclusion drawn from a careful inspection of Fig. 15 is that again the modified Leygue et al. model provides a better description of the simulation data for cis-PB than the original model. The same is also true for the trans-PB systems [Fig. 16] for Ne = 119 and se = 3.5 ns [12–16]; the corresponding value of X (=0.5) was chosen to provide the best comparison against the simulation results. How the predictions of the modified dual constraint model for the PP segment survival probability function with respect to the direct simulation data improve is analyzed in Figs. 17 and 18. As in the Leygue et al. model, the calculations have been carried out by using values for the quantities needed in the model as obtained

122

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

Fig. 13. Same as with Fig. 12 but for the dual constraint model (both original and modified versions) with G0N ¼ 2:3 MPa, Ne = 76 and se = 3 ns.

Fig. 14. Same as with Fig. 13b but now we have assumed the Rouse relaxation process to be present also in the second stage of the modified dual constraint model analysis. The values of the parameters are the same as in Fig. 13.

directly from the molecular simulations (Table 1). In the case of the original version, we have used again the scalings sR ¼ Z 2 se and sd ¼ 3Z 3 se , the typical choice in the literature. For cis-PB we have used Ne = 144 and se = 3.5 ns [12–16] while for trans-PB we have used Ne = 119 and se = 6 ns; for X we used the same values as in the Leygue et al. model (0.4 for cis-PB and 0.5 for trans-PB). Again, the choice se = 6 ns was made so that it provides the best comparison against the simulation results and is not considerably away from the value proposed from the simulations (see Table 1) [12–16]. Once more, we observe that the modified version of the dual constraint model provides a better description of the simulation data for both cis-PB and trans-PB than the original model. Equally satisfactory is the improvement of the modified version of the two models in their description of experimental data for the zero shear rate viscosity of an almost pure cis-PB melt, shown in Fig. 19. The calculations here have been carried out by using G0N ¼ 1:22 MPa for the plateau modulus of cis-PB (see Section 3). 4.2.2. Model predictions for high MW melts The last stage of our work with PB involved employing generalized scaling relationships for the quantities needed in the model calculations (sR, sd, L, dt, and b) similar to those invoked in the work with PE (i.e., only as a function of Z) in order to obtain predictions for g0, G0 , and G00 for high MW samples. These results are depicted in Figs. 20 and 21 for the Leygue et al. model, and in Figs. 22 and 23 for the dual constraint model. As explained in Section 4.1.2, the use of these scaling laws reduces the number of input data needed from the lower-scale simulations to just two: the number of

entanglements Ne and the entanglement time se. Regarding their values in the calculations, these were the same as in Section 4.2.1, that is: Ne = 144 and se = 2.5 ns for cis-PB, and Ne = 119 and se = 3.5 ns for trans-PB, for the Leygue et al. model, whereas for the dual constraint model we use Ne = 144 and se = 3.5 ns for cisPB, and Ne = 119 and se = 6 ns for trans-PB. The graphs show that the use of the general relationships provides predictions for both models that are not that far compared to the ones obtained by using the values suggested by molecular simulations [compare Fig. 20 with Fig. 15, Fig. 21 with Fig. 16, Fig. 22 with Fig. 17, and Fig. 23 with Fig. 18]. This implies that we can use these general scaling laws rather safely in conjunction with the two models in order to calculate the LVE properties of high-MW polymers. A typical example is shown in Fig. 24 which presents the predictions of the two versions of the Leygue et al. model for the MW-dependence of g0 of an almost pure cis-1,4-PB melt and the comparison with the experimental data. In part (a) of the figure, the calculations have been carried out with the set of values for the parameters G0N , Ne, and se obtained directly from the simulations for model pure cis-1,4-PB melts [12–16]. Part (b) of the figure shows how sensitive the comparison is to the exact value of these parameters. We see that the improvement can be greatly facilitated by choosing: a value for Ne closer to 125 (this is the lowest bound suggested from the simulations, since Ne = 144 ± 19), a relatively larger value for se (equal to 4 ns), and a relatively higher value for G0N (equal to 1.5 MPa). The value se = 4 ns is not far from the value we get from atomistic simulation (between 3 and 4 ns) whereas G0N ¼ 1:5 MPa is slightly larger than the value expected for cis-1,4-PB chains on the basis of a direct topological analysis of atomistically detailed samples (see Section 3). The same holds true also for the G0 and G00 data depicted in Fig. 25: using the values obtained directly from the simulations with the moderately entangled PB systems (i.e., G0N ¼ 1:22 MPa, Ne = 144 and se = 2.5 ns) the comparison favors somewhat the original version [see Fig. 25a and b]. By choosing the optimized values for G0N , Ne, and se that bring the predictions of the modified model and experimental data for g0 on top of each other, one obtains a significantly improved prediction with the modified model also for G0 and G00 [see Fig. 25c and d]. In these, we have also included the high frequency Rouse contribution using the expression proposed by Milner and McLeich [19b] (see also Fig. 14 of Ref. [15]):

"   # Z N X G0N 1 X 2p2 t 2p2 t þ : exp  exp  Z 3 p¼1 sR sR p¼Zþ1

ð9Þ

Finally, we compare the predictions of the two versions of the dual constraint model for the MW-dependence of g0 (Fig. 26) and of G0 and G00 (Fig. 27) as a function of frequency with experimental data. As for the Leygue et al. model, in part (a) of Fig. 26 the

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

123

Fig. 15. Predictions of the modified and original versions of the Leygue et al. model (lines) for the segment survival probability function w(s, t) on the basis of the computed curves from the direct PP analysis (symbols) for the two monodisperse PB-cis systems, at three different time instances: (a) t = 0.01sd (early), (b) t = 0.1sd (intermediate), and R1 (c) t = 0.35sd (late). The average tube survival probability function WðtÞ  0 wðs; tÞds is reported in (d). For the original version of the model, the usual assumption has been made that sd ¼ 3Z 3 se whereas for the modified version the values for the quantities needed (sR, sd, L, dt) have been obtained directly from the simulations (Table 1). The values of Ne and se were taken to be equal to 144 and 2.5 ns, respectively, whereas the CLF parameter X was chosen to be equal to 0.4.

Fig. 16. Same as with Fig. 15 but for the trans-PB systems. The values of Ne and se were taken to be equal to 119 and 3.5 ns, respectively, whereas the CLF parameter X was chosen to be equal to 0.5.

124

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

Fig. 17. Same as with Fig. 15 but for the original and modified versions of the dual constraint model. For the original version, we have used sd ¼ 3Z 3 se whereas for the modified version the values for the quantities needed (sR, sd, b) have been obtained directly from the simulations (Table 1). The values of Ne and se were taken to be equal to 144 and 3.5 ns, respectively.

Fig. 18. Same as with Fig. 17 but for the trans-PB systems. The values of Ne and se were taken equal to 119 and 6 ns, respectively.

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

125

Fig. 19. Predictions of: (a) the Leygue et al. model and (b) the dual constraint model (both original and modified versions) for the zero-shear-rate viscosity, and comparison with experimental and simulation [14] data for cis-PB at 413 K [G0N ¼ 1:22 MPa, Ne = 144, and se = 2.5 ns in (a) or se = 3.5 ns in (b)]. The values of the rest of the quantities needed (sR, sd, L, dt, b) for the modified version were those obtained directly from the simulations (Table 1).

Fig. 20. Same as with Fig. 15 but now for the modified Leygue et al. model we have used X = 0.4, and L/dt = (Z + 1). The values of Ne and se were kept the same.

calculations have been carried out with the set of values for G0N , Ne, se obtained directly from the simulations and used to obtain the results shown in Fig. 19b, whereas in part (b) of this figure we have used Ne = 125 and a larger value for se (equal to 6 ns); the plateau modulus was kept the same. These values seem to provide a better description of the experimental data. And the same is true for the G0 and G00 data depicted in Fig. 27. Parts (a) and (b) of the figure have been obtained by using G0N ¼ 1:22 MPa, Ne = 144 and se = 3.5 ns, which were extracted directly from the simulations with the moderately entangled PB systems; in this case, the comparison with the experimentally measured data is not favorable. In contrast, by using the optimized set for G0N , Ne, and se, the modified version provides a very satisfactory description of the G0 and G00 spectra [see Fig. 27c and d]. But we have to keep in mind that

sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with

the experimental samples are not 100% pure PB-cis and this may explain why we have to deviate slightly from the simulation estimates for Ne and se in order for the modified tube model to provide an accurate description of the experimental data.

5. Conclusions We have presented a systematic and very reliable methodology that allows one to obtain accurate predictions for the entire spectrum of the linear viscoelastic (LVE) properties of polymer melts directly from the chemical constitution of their chains. The methodology is a three-step approach and entails: (a) execution of long atomistic MD simulations for the given polymer with chain lengths

126

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

Fig. 21. Same as with Fig. 20 but for the trans-PB systems. The values of Ne and se were kept the same as in Fig. 16.

Fig. 22. Same as with Fig. 17 but for the modified dual constraint model we take now sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.4, and b ¼ 2=ð1 þ Z n Þ for the msd parameter with n = 0.29 [see Fig. 2b]. The values of Ne and se were kept the same.

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

127

Fig. 23. Same as with Fig. 18 but for the modified dual constraint model we take now sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.5, and b ¼ 2=ð1 þ Z n Þ for the msd parameter with n = 0.68 [see Fig. 2b]. The values of Ne and se were kept the same.

Fig. 24. Same as with Fig. 19a but for the modified version we take now sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.4, and L/dt = (Z + 1). In (a) we have used the same values for the parameters (G0N , Ne, se) whereas in (b) a we have used the values that provide the best fits to the data.

corresponding to moderately (up to well) entangled polymers for times on the order of the longest relaxation time of the melt at the temperature of interest; (b) direct topological analysis using robust geometric or stochastic techniques [39–42] and mapping of the accumulated trajectories of PPs onto modern tube models to parameterize these models, also to suggest improvements based on the outcome of the comparison with the direct PP analysis for the mathematical description of the various dynamic mechanisms, and (c) use of these modified tube models in conjunction with well-tested scaling laws (e.g., validated on the basis of detailed atomistic simulation data with shorter chain length systems) for the key model parameters to calculate the rheological data of the polymer melt but for significantly higher chain lengths extending deeply in the entangled regime.

The approach was designed and illustrated here for three polymers: PE, cis-1,4-PB and trans-1,4-PB, with very promising results. We have also proposed two new versions of two well-known tube models in the literature, the dual constraint and the Leygue et al. models. The most important modification is the use of a non-vanishing boundary condition of the form wðs ¼ 0; tÞ ¼ expðt=se Þ [13,16] in the corresponding reaction–diffusion equation; this brings model predictions for the functions w(s, t) and W(t) close to those obtained directly from the atomistic simulations [12–15]. We have further demonstrated that instead of directly computing the Rouse time sR, the reptation time sd, and the msd parameter b needed in these models from atomistic simulations, one can equally well rely on the use of simple scaling relationships (such as sd ¼ 3Z 3 ð1  XZ 0:5 Þse with X ffi 0:4  0:6, and b ¼ 2=ð1 þ Z n Þ with n = 0.29–0.56). Overall, it

128

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

Fig. 25. Predictions of the modified and original versions of the Leygue et al. model for the storage [(a) and (c)] and loss [(b) and (d)] moduli, and comparison with experimental data. For the modified version, we have used sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.4, and L/dt = (Z + 1). In (a) and (b) we have employed the same values for the parameters as in Fig. 24a whereas in (c) and (d) we have used the parameters values of Fig. 24b. Rouse mode contributions have been accounted according to Milner and McLeish [19b] (see text).

Fig. 26. Same as with Fig. 19b but for the modified version we take now sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.4, and b ¼ 2=ð1 þ Z n Þ for the msd parameter with n = 0.29. In (a) we have used the same values for the parameters (G0N , Ne, se) whereas in (b) a we have used the values that best fit the data.

turns out that the number of parameters actually needed in the calculations boil down to just three: the entanglement length Ne, the entanglement time se, and the plateau modulus G0N . Atomistic simulations [12–16,26,29,42] with relatively short (moderately entangled) systems provide rather straightforwardly reliable estimates for the first two of these parameters, Ne and se. The plateau modulus, on the other hand, can be either estimated through G0N ¼ 45 Mq0RTNe or directly borrowed from the experimental measurements [31–34]. Given the values of these three parameters, the model can be used next to provide reliable predictions for the LVE properties of the polymer under study over the entire regime of chain lengths realized in practical or industrial applications.

In this work, we have used the two modified versions of the two tube models [16] to compute the LVE data of the three polymers studied here (PE, cis-1,4-PB, and trans-1,4-PB) as a function of their MW and to compare against experimentally measured data. Our calculations have demonstrated the importance of the exact value of the entanglement time se for a quantitative comparison of the model predictions with experimentally measured data. Our multi-scale approach (coupling MD simulations with a PP analysis and using elements of the reptation theory) leads self-consistently to the prediction of the key material properties of higher-MW polymers starting from a detailed analysis of the dynamic properties of considerably shorter (but entangled) samples, thereby

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

129

Fig. 27. Predictions of the modified and original versions of the dual constraint model for the storage [(a) and (c)] and loss [(b) and (d)] moduli and comparison against the experimental data. For the modified version, we have used sR ¼ Z 2 se for the Rouse time, sd ¼ 3Z 3 ð1  XZ 0:5 Þse for the reptation time with X = 0.4, and b ¼ 2=ð1 þ Z n Þ for the msd parameter with n = 0.29. In (a) and (b) we have employed the same values for the parameters as in Fig. 26a whereas in (c) and (d) we have used the parameter values of Fig. 26b. Rouse mode contributions have been accounted according to Milner and McLeish [19b] (see text).

avoiding the need to assume non-physical or unjustified values (for example for se). Our methodology has been formulated strictly for linear polymer melts. However, it is in principle extendable to other polymer architectures such as H-shaped [29,43], star [19a,44], even treelike ones [45], as long as a systematic mapping of the atomistic MD data onto a closed-form constitutive model describing their rheology can be made. A typical case is the pom-pom model which proposes two constitutive equations: one for the orientation tensor of the main chain backbone and another for the arm degree of stretching. Our approach relies on the use of a theoretically informed model with input from lower-level (atomistic) simulations and can describe the rheological data of several families of polymers without the need to resort to arbitrary fitting processes with no reference to a theoretical model or no input from detailed molecular simulations. A similar approach has been followed very recently by Li et al. in Ref. [46]. Acknowledgements The authors are thankful to M. Kapnistos, A.E. Likhtman and D. Vlassopoulos for making available their linear viscoelastic data of some polybutadienes, for comparison with our predictive approach. Support provided by the European Commission through the MODIFY (FP7-NMP-2008-SMALL-2, Code 228320) and the Marie-Curie grant via the Project VISCOnanoNET (FP7-PEOPLE-2011CIG, Code 293945) projects are greatly acknowledged. References [1] M. Doi, S.F. Edwards, The Theory of Polymer Dynamics, Clarendon Press, Oxford, 1986.

[2] M. Doi, Introduction to Polymer Physics, Clarendon Press, Oxford, 1995. [3] S.F. Edwards, The statistical mechanics of polymerized material, Proc. Phys. Soc. 92 (1967) 9–16. [4] P.G. de Gennes, Reptation of a polymer chain in the presence of fixed obstacles, J. Chem. Phys. 55 (1971) 572–579. [5] M. Doi, S.F. Edwards, Dynamics of concentrated polymer systems Part 1. – Brownian motion in the equilibrium state, J. Chem. Soc. Faraday Trans. 2 74 (1978) 1789–1801; M. Doi, S.F. Edwards, Dynamics of concentrated polymer systems Part 2. – Molecular motion under flow, J. Chem. Soc. Faraday Trans. 2 74 (1978) 1802– 1817; M. Doi, S.F. Edwards, Dynamics of concentrated polymer systems Part 3. – The constitutive equation, J. Chem. Soc. Faraday Trans. 2 74 (1978) 1818–1832. [6] S.F. Edwards, The theory of rubber elasticity, Polymer 9 (1977) 140–143. [7] H. Watanabe, Viscoelasticity and dynamics of entangled polymers, Prog. Polym. Sci. 24 (1999) 1253–1403; T.C.B. McLeish, Tube theory of entangled polymer dynamics, Adv. Phys. 51 (2002) 1379–1527; E. van Ruymbeke, C.-Y. Liu, C. Bailly, Rheol. Rev. (2007) 53–134. [8] M. Doi, Explanation for the 3.4 power law of viscosity for polymeric liquids on the basis of the tube model, J. Polym. Sci.: Polym. Lett. Ed. 19 (1981) 265–273; M. Doi, Explanation for the 3.4 power law of viscosity for polymeric liquids on the basis of the tube model, J. Polym. Sci.: Polym. Phys. Ed. 21 (1983) 667–684. [9] S.F. Edwards, J.W.V. Grant, The effect of entanglements on diffusion in a polymer melt, J. Phys. A: Math. Nucl. Gen. 6 (1973) 1169–1185; P.G. de Gennes, Reptation of stars, J. Phys. 36 (1975) 1199–1203; P.G. de Gennes, Dynamics of entangled polymer solutions. I. The rouse model, Macromolecules 9 (1976) 587–593; J. Klein, The onset of entangled behavior in semidilute and concentrated polymer solutions, Macromolecules 11 (1978) 852–858. [10] C. Pattamaprom, R.G. Larson, T.J. Van Dyke, Quantitative predictions of linear viscoelastic rheological properties of entangled polymers, Rheol. Acta 39 (2000) 517–531; C. Pattamaprom, R.G. Larson, Predicting the linear viscoelastic properties on monodisperse and polydisperse polystyrenes and polyethylenes, Rheol. Acta 40 (2001) 516–532; C. Pattamaprom, R.G. Larson, A. Sirivat, Determining polymer molecular weight distributions from rheological properties using the dual-constraint model, Rheol. Acta 47 (2008) 689–700.

130

P.S. Stephanou, V.G. Mavrantzas / Journal of Non-Newtonian Fluid Mechanics 200 (2013) 111–130

[11] A. Leygue, C. Bailly, R. Keunings, A differential tube-based model for predicting the linear viscoelastic moduli of polydisperse entangled linear polymers, J. Non-Newtonian Fluid Mech. 133 (2006) 28–34. [12] P.S. Stephanou, C. Baig, V.G. Mavrantzas, Projection of atomistic simulation data for the dynamics of entangled polymers onto the tube theory: calculation of the segment survival probability function and comparison with modern tube models, Soft Matter 7 (2011) 380–395. [13] P. S. Stephanou, Development of Scale-bridging Methodologies and Algorithms Founded on the Outcome of Detailed Atomistic Simulations for the Reliable Prediction of the Viscoelastic Properties of Polymer Melts, Ph.D. Thesis, University of Patras (Greece), 2011. [14] P.S. Stephanou, C. Baig, G. Tsolou, V.G. Mavrantzas, M. Kröger, Quantifying chain reptation in entangled polymer melts: topological and dynamical mapping of atomistic simulation results onto the tube model, J. Chem. Phys. 132 (2010) 124904. [15] C. Baig, P.S. Stephanou, G. Tsolou, V.G. Mavrantzas, M. Kröger, Understanding dynamics in binary mixtures of entangled cis-1,4-polybutadiene melts at the level of primitive path segments by mapping atomistic simulation data onto the tube model, Macromolecules 43 (2010) 8239–8250. [16] P.S. Stephanou, C. Baig, V.G. Mavrantzas, Toward an improved description of constraint release and contour length fluctuations in tube models for entangled polymer melts guided by atomistic simulations, Macromol. Theory Simul. 20 (2011) 752–768. [17] V.A. Raju, G.G. Smith, G. Marin, J.R. Knox, W.W. Graessley, Properties of amorphous crystallizable hydrocarbon polymers I. Melt rheology of fractions of linear polyethylene, J. Polym. Sci. Part B Polym. Phys. 17 (1979) 1183–1195. [18] R.H. Colby, L.J. Feeters, W.W. Graessley, Melt viscosity-molecular weight relationship for linear polymers, Macromolecules 20 (1987) 2226–2237. [19] a S.T. Milner, T.C.B. McLeish, Parameter-free theory for stress relaxation in star polymer melts, Macromolecules 30 (1997) 2159–2166; b S.T. Milner, T.C.B. McLeish, Reptation and contour-length fluctuations in melts of linear polymers, Phys. Rev. Lett. 81 (1998) 725–728. [20] J.L. Viovy, M. Rubinstein, R.H. Colby, Constraint release in polymer melts: tube reorganization versus tube dilution, Macromolecules 24 (1991) 3587–3596; S.T. Milner, Relating the shear-thinning curve to the molecular weight distribution in linear polymer, Macromolecules 40 (1996) 303–315. [21] J.L. Viovy, Constraint release in the slip-link model and the viscoelastic properties of polymers, J. Phys. (France) 46 (1985) 847–853; G. Marrucci, Relaxation by reptation and tube enlangement: a model for polydisperce polymers, J. Polym. Sci., Polym. Phys. Ed. 23 (1985) 159–177; C. Tsenoglou, Viscoelasticity of binary homopolymer blends, Polym. Prepr. (Am. Chem. Soc. Div. Polym. Chem.) 28 (1987) 185–186; C. Tsenoglou, Molecular weight polydispersity effects on the viscoelasticity of entangled linear polymers, Macromolecules 24 (1991) 1762–1767; J. des Cloizeaux, Double Reptation vs. Simple Reptation in Polymer Melts, Europhys. Lett. 5 (1988) 437–442; J. des Cloizeaux, Double Reptation vs. Simple Reptation in Polymer Melts, Erratum: Europhys. Lett. 6 (1988) 475. [22] R.H. Colby, M. Rubinstein, Two-parameter scaling for polymers in h solvents, Macromolecules 23 (1990) 2753–2757. [23] W.H. Press, S.A. Tevkolsky, B.P. Flannery, Numerical Recipes in Fortran 77, Cambridge University Press, Cambridge, 1992. [24] G. Tsolou, V.G. Mavrantzas, D.N. Theodorou, Detailed atomistic molecular dynamics simulation of cis-1,4-poly(butadiene), Macromolecules 38 (2005) 1478–1492. [25] V.A. Harmandaris, V.G. Mavrantzas, D.N. Theodorou, Atomistic molecular dynamics simulation of polydisperse linear polyethylene melts, Macromolecules 31 (1998) 7934–7943. [26] V.A. Harmandaris, V.G. Mavrantzas, D.N. Theodorou, M. Kröger, J. Ramírez, H.C. Öttinger, D. Vlassopoulos, Crossover from the rouse to the entangled polymer melt regime: signals from long, detailed atomistic molecular dynamics simulations, supported by rheological experiments, Macromolecules 36 (2003) 1376–1387.

[27] G. Tsolou, V.G. Mavrantzas, Hierarchical modeling of polymeric systems at multiple time and length scales, in: C.S. Adjiman, A. Galindo (Eds.), Application of Molecular Systems Engineering, vol. 6, Wiley-VCH, 2009, pp. 85–134. [28] H.C. Öttinger, Coarse-graining of wormlike polymer chains for substantiating reptation, J. Non-Newtonian Fluid Mech. 120 (2004) 207–213. [29] N.Ch. Karayiannis, V.G. Mavrantzas, Hierarchical modeling of the dynamics of polymers with a nonlinear molecular architecture: calculation of branch point friction and chain reptation time of H-shaped polyethylene melts from long molecular dynamics simulations, Macromolecules 38 (2005) 8583–8596. [30] D.S. Pearson, G.V. Strate, E. von Meerwall, F.C. Schilling, Viscosity and selfdiffusion coefficient of linear polyethylene, Macromolecules 20 (1987) 1133– 1141. [31] D.S. Pearson, L.J. Fetters, W.W. Graessley, G. Strate, V.E. von Meerwall, Viscosity and self-diffusion coefficient of hydrogenated polybutadiene, Macromolecules 27 (1994) 711–719. [32] L.J. Fetters, D.J. Lohse, S.T. Milner, W.W. Graessley, Packing length influence in linear polymer melts on the entanglement, critical, and reptation molecular weights, Macromolecules 32 (1999) 6847–6851. [33] L.J. Fetters, D.J. Lohse, D. Richter, T.A. Witten, A. Zirkel, Connection between polymer molecular weight, density, chain dimensions, and melt viscoelastic properties, Macromolecules 27 (1994) 4639–4647. [34] J.D. Ferry, Viscoelastic Properties of Polymers, third ed., John Wiley and Sons, New York, 1980. [35] L.J. Fetters, D.J. Lohse, S.T. Milner, William W. Graessley, Packing length influence in linear polymer melts on the entanglement, critical, and reptation molecular weights, Macromolecules 32 (1999) 847–6851. [36] M. Kapnistos, D. Vlassopoulos, J. Roovers, L.G. Leal, Linear rheology of architecturally complex macromolecules: comb polymers with linear backbones, Macromolecules 38 (2005) 7852–7862. [37] J. Roovers, Linear viscoelastic properties of polybutadienes: a comparison with molecular theories, Polym. J. 18 (1986) 153–160. [38] M. Kapnistos, A.E. Likhtman, D. Vlassopoulos, unpublished data. [39] R. Everaers, S.K. Sukumaran, G.S. Grest, C. Svaneborg, A. Sivasubramanian, K. Kremer, Rheology and microscopic topology of entangled polymeric liquids, Science 303 (2004) 823–826. [40] M. Kröger, Shortest multiple disconnected path for the analysis of entanglements in two- and three-dimensional polymeric systems, Comput. Phys. Commun. 168 (2005) 209–232. [41] K. Foteinopoulou, N.C. Karayiannis, V.G. Mavrantzas, M. Kröger, Primitive path identification and entanglement statistics in polymer melts: results from direct topological analysis on atomistic polyethylene models, Macromolecules 39 (2006) 4207–4216. [42] C. Tzoumanekas, D.N. Theodorou, Topological analysis of linear polymer melts: a statistical approach, Macromolecules 39 (2006) 4592–4604. [43] T.C.B. McLeish, R.G. Larson, Molecular constitutive equations for a class of branched polymers: the pom-pom polymer, J. Rheol. 42 (1998) 140–143; H.C. Öttinger, Thermodynamic admissibility for the pompon model for branched polymers, Rheol. Acta 40 (2001) 317–321; H.C. Öttinger, Modelling complex fluids with a tensor and a scalar as structural variables, Rev. Mex. de Física S. 48 (2002) 220–229. [44] E. van Ruymbeke, M. Kapnistos, D. Vlassopoulos, T. Huang, D.M. Knauss, Linear melt rheology of pom–pom polystyrenes with unentangled branches, Macromolecules 40 (2007) 1713–1719. [45] H. Watanabe, Y. Matsumiya, E. van Ruymbeke, D. Vlassopoulos, N. Hadjichristidis, Viscoelastic and dielectric relaxation of a cayley-tree-type polyisoprene: test of molecular picture of dynamic tube dilation, Macromolecules 41 (2008) 6110–6124. [46] Y. Li, S. Tang, B.V. Abberton, M. Kröger, C. Burkhart, B. Jiang, G.J. Papakonstantopoulos, M. Poldneff, W.K. Liu, A predictive mutliscale computational framework for viscoelastic properties of linear polymers, Macromolecules 53 (2012) 5935–5952.