Primitive Chain Network Simulation of Elongational Flows of

Article pubs.acs.org/Macromolecules

Primitive Chain Network Simulation of Elongational Flows of Entangled Linear Chains: Stretch/Orientation-induced Reduction of Monomeric Friction Takatoshi Yaoita,†,‡,⊥ Takeharu Isaki,‡ Yuichi Masubuchi,† Hiroshi Watanabe,†,* Giovanni Ianniruberto,§ and Giuseppe Marrucci§ †

Institute for Chemical Research, Kyoto University, Kyoto 611-0011, Japan Material Science Laboratory, Mitsui Chemicals Inc., Chiba 299-0265, Japan § Dipartimento di Ingegneria Chimica, Università “Federico II”, P. le Tecchio 80, 80125 Napoli, Italy ‡

ABSTRACT: Well-entangled monodisperse linear polystyrene melts exhibit monotonic thinning of the steady state elongational viscosity with increasing the strain rate ε̇ even beyond the Rouse relaxation frequency, τR‑1. This behavior is quite different from the thinning followed by hardening at ε̇ > τR‑1 observed for entangled semidilute solutions. We attempt to elucidate the molecular origin of this difference by focusing on the concept of stretch/orientation-dependent monomeric friction ζ recently proposed by Ianniruberto and co-workers. Specifically, literature data of the stress relaxation after cessation of transient elongational flow, reported for both PS melts and solutions, are analyzed to evaluate the stretch/ orientation-dependent decrease of ζ. In our working hypothesis, ζ is expressed as a function of the factor Fso = λ2̃ S̅, where λ̃ is the normalized stretch ratio of entangled subchains defined with respect to the fully stretched state, and S̅ is an average orientational anisotropy of the components (polymer plus solvent if any) in the system. The factor Fso was estimated from the stress decay data after flow cessation. The resulting functional form of ζ(Fso) was then used in the primitive chain network (PCN) simulation including finite extensible nonlinear elasticity (FENE) to examine the elongational behavior of melts and solutions. For melts the simulation indicates that ζ decreases significantly under fast elongation because the entangled subchains are short and approach the fully stretched (and fully oriented) limit rather easily. Hence, the steady elongational viscosity ηE follows this decrease of ζ to exhibit the monotonic thinning even at ε̇ > τR‑1. In contrast, for solutions, the simulated ηE exhibits thickening at ε̇ > τR‑1 because the average anisotropy S̅ is governed by the solvent and remains small, thus overwhelming the increase of the subchain stretch λ̃. The simulated results proved to be in satisfactory agreement with the experiments.

1. INTRODUCTION The molecular mechanisms responsible for the viscoelastic behavior of entangled polymeric liquids under fast uniaxial extensional flow have not been fully elucidated yet. For example, the strain-rate (ε̇) dependence of the steady state uniaxial elongational viscosity ηE is quite different for melts1−3 and solutions:4−6 ηE of polystyrene (PS) solutions4−6 decreases with increasing ε̇ up to the equilibrium Rouse relaxation frequency τR‑1 and then increases at higher ε̇ when the entangled subchains begin to be stretched. In contrast, for PS melts,1−3 ηE keeps decreasing monotonically even at ε̇ > τR‑1. The behavior of the semidilute solutions, thinning followed by hardening, is reasonably described by the tube model4,5 that includes the mechanisms of chain stretch and convective constraint release7 (CCR). However, this version of the tube model cannot describe the monotonic thinning seen for PS melts. Marrucci and Ianniruberto8 proposed the interchain tube pressure effect (ICP) that for ε̇ < τR‑1 predicted the −1/2 power law observed in the ηE vs ε̇ curve of PS melts. Wagner et al.9−11 incorporated the ICP mechanism into the version of © 2012 American Chemical Society

tube model called molecular stress function (MSF) model showing that the model did not predict any upturn at ε̇ > τR‑1 and could therefore satisfactorily describe the ηE data of PS melts. However, the tube relaxation time τa included in the MSF-ICP model was not deduced from a clear molecular basis and needed to be adjusted rather arbitrarily (to a value exceeding the terminal relaxation time for the case of low molecular weight PS melts.) In addition, combination of ICP with CCR, the latter being a basic mechanism in the tube model, leads to a failure in describing the ηE data.12 Furthermore, it is not apparent why the ICP mechanism should be included in the tube model only for the case of entangled melts and not also for entangled solutions. Considering the above uncertainty in the molecular understanding of the elongational behavior of entangled polymers, we recently conducted 3D multichain slip link simulations based Received: November 17, 2011 Revised: February 20, 2012 Published: March 6, 2012 2773

dx.doi.org/10.1021/ma202525v | Macromolecules 2012, 45, 2773−2782

Macromolecules

Article

ζ remains constant, i.e., that the friction reduction has no effect on the shear viscosity at experimentally accessible conditions.

on the primitive chain network (PCN) model.13 This model takes into account all known relaxation mechanisms, i.e., reptation,14 thermal constraint release,15 CCR, force balance at the entanglements, and finite extensibility (FENE), but no ICP. The simulated ηE exhibited thinning followed by thickening at ε̇ > τR‑1, which is in quantitative agreement with the data for solutions but not for melts. Kushwaha and Shaqfeh16 reported that their multichain slip link simulation also resulted in thickening at high strain rates under planar extensional flow. It so appears that further improvement is necessary for the PCN-FENE simulation to describe the elongational behavior of both solutions and melts in a consistent way. In this study, we focus on the stretch/orientation-induced reduction of the monomeric (actually Kuhn-segment’s) friction ζ proposed by Ianniruberto et al.17 and incorporate this mechanism in the PCN-FENE simulation. The molecular basis of the proposed mechanism is the notion that elongated objects show a reduced (anisotropic) friction when roughly parallel to one another with respect to the randomly oriented state. Indeed, Uneyama et al.18 analyzed the anisotropic reduction of ζ of a dumbbell model under fast shear flow to show that the viscoelastic relaxation accelerates in fast shear due to friction reduction but the end-to-end vector fluctuation time remains constant, in agreement with the dielectric response. Hua and Yang19 also considered a decrease of ζ in fast flows, though induced in their case by a reduction of the entanglement density due to CCR. However, the prediction of the hardening of ηE at ε̇ > τR‑1 remained unaltered. In our approach, we assume as a working hypothesis that ζ is not affected by the entanglement density but depends only on a stretch/orientation parameter, Fso = λ2̃ S̅, where λ̃ (= λ/λmax) is a reduced stretch ratio of the entangled subchain defined with respect to the fully stretched limit of λ approaching λmax, and S̅ is an average orientational anisotropy of the components (polymer plus solvent) in the system. The flow-induced anisotropic decrease of the friction could result in a tensorial form of the friction coefficient,17,18 but in this study we approximate this coefficient by a scalar quantity, ζ. In order to deduce the dependence of ζ on Fso, we analyze the tensile stress data of PS melts20 and solutions5 after cessation of transient uniaxial elongational flow, and evaluate this dependence from the data5,20 of initial stress decay. Such a decay reflects the Rouse-like chain retraction with the (possibly reduced) monomeric friction ζ(Fso), where Fso is the stretch/orientation parameter reached during the elongational start up just before the flow cessation. The function ζ(Fso) is then used in the PCN-FENE simulation to examine the elongational behavior of PS melts and solutions. It will be shown that the simulated ηE of PS melts exhibits monotonic thinning because the entangled subchains in PS melts are short, and approach the fully stretched and fully oriented limit rather easily. For this case, the reduction of ζ overwhelms the increase of the subchain tension, giving rise to a monotonic thinning of ηE. In contrast, for the semidilute solutions, the simulated ηE exhibits thinning followed by thickening at high ε̇ because the average orientational anisotropy is governed by the solvent and remains small. This effect overwhelms the increase in chain stretch and maintains a small value of the stretch/orientation parameter Fso, thus keeping the friction coefficient close to the equilibrium value. Before closing the paper we will also report on PCN-FENE simulation for PS melts under shear flow in order to show that

2. RESULTS AND DISCUSSION 2.1. Working Hypothesis of Stretch/OrientationInduced Reduction of Friction. The frictional interaction in a melt should weaken when the entangled subchains are stretched and oriented, and overlapping between Kuhn segments is effectively reduced. Conversely, in polymer solutions, since the solvent molecules relax almost instantaneously, the solvent−polymer frictional interaction is expected to be weakly affected by the subchain stretch and orientation, as discussed by Ianniruberto et al.17 For the assumed decrease of the monomeric friction ζ, we may naturally choose the fully stretched and fully oriented state of the whole system (polymer plus solvent if any) as the reference state, where the friction attains the minimum value. Thus, as a working hypothesis, we take ζ to be a function of a stretch/orientation parameter Fso defined with respect to this reference state: ζ = ζ(Fso)

with

2 Fso = λ̃ S ̅

and

λ̃ = λ /λ max

(1)

Here, λ is the (average) subchain stretch ratio defined with respect to the equilibrium state, λ̃ is the reduced stretch ratio normalized to the fully stretched state (λ = λmax), and S̅ is a weighted average orientational anisotropy of the components (polymer and solvent) in the system. Since the solvent stays isotropic under all conditions of our interest, S̅ is assumed to be related to the anisotropy of the polymer, SP, as S ̅ = ϕPSP

(2)

where ϕP is the polymer volume fraction, and the polymer order parameter SP representing the anisotropy of the subchain orientation vector u is obtained from the eigenvalues of the orientation tensor ⟨uu⟩ (with ⟨...⟩ indicating the ensemble average). More specifically, for uniaxial elongational flows SP is given by SP = ⟨ux2⟩ − ⟨uy2⟩, where ux, uy are components parallel and perpendicular to the stretching direction, respectively. In shear flows, Sp is obtained from the eigenvalues associated with the two eigenvectors in the shear plane (see later). The choice of λ2̃ S̅ as the stretch/orientation parameter is suggested by the fact that λ̃2S̅ somehow measures the “monomer” orientation in the elongation direction, where by monomers we mean both the Kuhn segments of the polymer and the solvent molecules. Indeed, the subchain stretch λ̃ becomes unity when all Kuhn segments within the subchains are fully aligned to the subchain orientation u, and the order parameter SP of the polymer subchains becomes unity when the subchains themselves are fully oriented to the elongation direction. Of course, in the case of solutions S̅ is never unity. 2.2. Overview of Stress Data after Cessation of Elongational Flow. The possible decrease of ζ due to chain stretch can be most clearly detected from data of the tensile stress σE‑ after cessation of fast uniaxial elongational flow. The initial decay of σE‑ due to Rouse-like chain retraction reflects relaxation of the stretched/oriented conformation of the entangled subchains reached just before the flow cessation, and it can then be used to verify whether or not ζ has been affected by the flow-induced stretch/orientation of the subchains. Furthermore, the stretch/orientation parameter Fso on which ζ has been assumed to depend is reflected in the σE data just before the flow cessation. On the basis of these 2774


Macromolecules

Article

arguments, we now examine literature data of stress before and after cessation of fast elongational flows. Figure 1 shows the σE‑ data of monodisperse linear PS melt (M = 1.45 × 105; PS145K) at 120 °C reported by Nielsen et

Figure 2. Test of time-rate factorability at long t for the stress decay data reported in Figure 1. Numbers indicate the strain rate.

Figure 1. Stress decay data of PS145K melt20 at 120 °C after cessation of elongational flow (up to εflow = 3). Numbers indicate the strain rate.

al.20 These data were measured after cessation of elongational flow in the transient hardening regime at strain rates ε̇ = 0.001, 0.003, 0.01, and 0.03 s−1; the flow time was chosen to achieve the same macroscopic Hencky strain for all rates, i.e., εflow = 3.0 (elongation ratio =20.1). The Rouse relaxation time of the chain is τR ≅ 1.0 × 103 s,20 and the Weissenberg number based on τR is WiR ≅ 1, 3, 10, and 30 for those ε̇ values. In the experiments, the σE‑ value immediately after the flow cessation is not completely free from artifacts such as the shock (vibration) of the stretching device as well as inertia. Thus, in Figure 1, the σE‑ (t) data are normalized by the data at short time (but still long enough to erase the effect of those artifacts), namely 1 s (≪ τR), and plotted against the time t after the flow cessation. These σE‑ (t)/σE‑ (1) data unequivocally allow us to examine the effect on the stress decay behavior of the elongational rate ε̇ that existed before flow cessation. Figure 1 clearly demonstrates that the increase of ε̇ (>1/ τR) results in an acceleration of the decay of the σE‑ (t)/σE‑ (1) ratio, particularly at short times. For the highest ε̇ (= 0.03 s−1), an appreciable decay not seen at low ε̇ starts at a very short time (t ≅ 2 s). These results already lend qualitative support to the mechanism of stretch/orientation-induced ζ-reduction. However, the finite extensibility (that leads to strain-hardening of the subchains) also contributes to this rapid decay at large ε̇. Thus, a quantitative evaluation of the magnitude of the ζreduction requires a detailed analysis of the σE‑ (t)/σE‑ (1) data, as discussed in the following. To begin with, the terminal relaxation attributable to reptation typically observed in linear viscoelasticity should also dominate the nonlinear σE‑ (t)/σE‑ (1) data at long times. This can indeed be observed in Figure 1 where all curves become roughly parallel to one another at long t. Figure 2 confirms that the terminal relaxation process, occurring well after the Rouse-like chain retraction process, is insensitive to the pre-elongational rate ε̇. Indeed, Figure 2 shows that the normalized ratio hE‑1{σE‑ (t)/σE‑ (0)} for various ε̇ can be superposed at long times with an adequate choice of the damping function hE (with the data for ε̇ = 0.001 s−1 being used as the reference). This time-rate separability is analogous to the time-strain separability14 under step strains.

Figure 3. Same data as in Figure 1 but after subtraction of the terminal relaxation. Solid curves are the result of fitting with a sum of exponential relaxation modes (see Table 1). Dashed curve is the Rouse relaxation curve of linear viscoelasticity, eq 3 with tflow = 3 × 103 s.

Data at long times in Figure 2 are well fitted by the single exponential {σE‑ (t)/σE‑ (1)}terminal = 0.245exp(−t/τd) with τd = 6.8 × 103 s. As can be noted from Figure 2, the terminal relaxation process somewhat contributes to the σE‑ (t)/σE‑ (1) data even at short t, in particular for the case of the smallest ε̇ = 0.001 s−1. We need to subtract this contribution from the σE‑ (t)/σE‑ (1) data for accurate examination of the fast relaxation ‑ ‑ (t)/σE,fast (1) process (Rouse-like chain retraction). The σE,fast ratio for the fast relaxation process obtained after this subtraction is shown in Figure 3 together with the Rouse relaxation curve of linear viscoelasticity after flow cessation (dashed curve), the latter being straightforwardly obtained from the constitutive equation of the Rouse model14 as: σ− Rouse(t ) =Q − σRouse(1) g pR = 2775

⎛ p2 t ⎞ ⎟⎟ g pR exp⎜⎜ − ⎝ τR ⎠ p=1 N

∑

⎧ ⎛ p2 t ⎞⎫ ⎪ flow ⎟⎪ ⎨1 − exp⎜⎜ − ⎬ ⎪ τR ⎟⎠⎪ p2 ⎩ ⎝ ⎭

with

1

(3)


Macromolecules

Article

Equation 3 has been used for comparison with the data at the slowest rate, ε̇ = 0.001 s−1, for which the flow duration is tflow = 3 × 103 s.20 In eq 3, we have used τR = 103 s,20 while Q is a ‑ ‑ (t)/σRouse (1) normalization factor ensuring that the ratio σRouse is unity at t = 1. The number of Rouse segments per chain, N, was evaluated from the molar mass of this segment (= 850 for PS21). Figure 3 shows that the data for ε̇ = 0.001 s−1 are very close to the dashed curve, suggesting that no effect on the friction coefficient is present in such a slow elongational flow. It also shows that the approximation of the terminal relaxation to a single exponential (as used in the subtraction) is indeed acceptable even at short times. ‑ ‑ As seen in Figure 3, the initial decay of the σE,fast (t)/σE,fast (1) ratio due to chain retraction is considerably accelerated for all ε̇ values larger than 0.001 s−1, again lending support to the mechanism of stretch/orientation-induced ζ-reduction. Before reporting on a quantitative estimate of such a reduction (see next subsection), we here briefly consider solution data. Stress decay data after cessation of elongational flow in the transient hardening regime have been reported also for two entangled semidilute PS solutions at room temperature.5 The molecular weight M, concentration c, and solvent, for the two solutions respectively, are as follows: M = 10.2 × 10 6 (PS10.2M), c = 6 wt %, dibutyl phthalate, and M = 3.9 × 106 (PS3.9M), c = 10 wt %, diethyl phthalate.5 The PS10.2M solution was elongated at the rate ε̇ = 1.0 s−1 (WiR = 3.0) up to three values of the Hencky strain εflow = 2.0, 2.5, and 4.0, whereas the PS3.9M solution had ε̇ = 11.7 s−1 (WiR = 3.3), and εflow = 3.3, 4.0, and 4.6. We analyze these data in the same way as for the PS145K melt to evaluate the tensile stress due to the ‑ fast relaxation (chain retraction) process. As before, the σE,fast (t) data were normalized by the stress value at a suitably short time t* (but still long enough to erase possible artifacts): t*=0.1s and t*=0.05s for the PS10.2M and PS3.9M solutions, ‑ ‑ respectively. The σE,fast (t)/σE,fast (t*) ratios thus obtained are shown in Figures 4 and 5, respectively, together with the Rouse relaxation curve (after cessation of steady elongational flow, eq 3 with gpR = 1/p2) properly normalized at t = t*. We note that the data for the solutions are insensitive to εflow and virtually

Figure 5. Stress decay during the fast relaxation process of the PS3.9M solution after cessation of elongational flow at ε̇ = 11.7 s−1. The curve is eq 3 with gpR = 1/p2.

superpose to the Rouse curve. It so appears that the fast relaxation process of the solutions does not exhibit the nonlinear acceleration observed for the melt, despite the fact that the data for solutions and melts covered similar values of WiR (∼3) and εflow (∼3). This result implies that the stretch/ orientation-induced ζ-reduction is negligible for solutions under the experimental conditions examined. 2.3. Evaluation of the Stretch/Orientation-Induced Reduction of the Friction Coefficient. Going back to the PS145K melt case, we observe that the fast relaxation data reported in Figure 3 for different ε̇’s also appear to lie parallel to one another at long times. The time-rate separability is then examined in Figure 6a, where those data are reduced by a suitable damping function hE,fast and semilogarithmically plotted against time. All data sets show unambiguously the same “terminal” relaxation time τ1 (reciprocal of the slope in the semilogarithmic plot at long times), and hence satisfy the separability. In addition, τ1 (=1.1 × 103 s) agrees well with the Rouse time of PS145K, τR ≅ 1.0 × 103 s, reported by Nielsen et al.20 This result reveals that at times of the order of τR the friction coefficient is the same as at equilibrium, which indeed is reasonable because any possible effect of chain stretch has relaxed over that time scale. It now remains to analyze the behavior at shorter times where time-rate separability does not hold, as shown in Figure 6a and even more clearly in Figure 3. At first sight, one may be tempted to evaluate the magnitude of the acceleration due to decrease of the friction coefficient from the initial stress decay rate. However, this evaluation method is not adequate because of the multimode nature of the relaxation at short times. The initial stress decay rate changes even for Rouse relaxation in the linear range: Indeed, although the mode relaxation times τR/p2 are fixed quantities in such a case, the mode intensities gpR change with tflow (see eq 3). We therefore need a multimode analysis of the curves in Figures 3 or 6a. We perform such analysis by fitting the data with a sum of exponentially decaying modes,Σpgp exp(−t/τp), where the slowest one has the relaxation time τ1 (=1.1 × 103 s as determined in Figure 6a for all ε̇ values), and the fastest one has the smallest relaxation time τp experimentally resolved for each ε̇ value. In particular, we have applied the procedure X suggested long time ago by Sobue et al.22 The fitting parameters (times τ and magnitudes g) are reported in Table

Figure 4. Stress decay during the fast relaxation process of the PS10.2M solution after cessation of elongational flow at ε̇ = 1.0 s−1. The curve is eq 3 with gpR = 1/p2. 2776


Macromolecules

Article

reduction, which in turn enables us to evaluate the magnitude of this reduction, as explained later in more detail. In Figure 6b, we also note that the τp’s for the lower order modes remain proportional to p−2 and thus the relaxation of those modes is not accelerated even at large ε̇. This result suggests that the relaxation is accelerated only when the entangled subchains are highly stretched. The stretch relaxes significantly through the higher order modes so that the acceleration vanishes in the time scale of lower order modes. This argument allows us to attribute the accelerated relaxation of the higher order modes to the stretch/orientation-induced reduction of the monomeric friction coefficient ζ (as well as to the FENE effect, as discussed below). In particular, since all relaxation times at time t are proportional to the friction coefficient existing at the same time t, the friction coefficient just before flow cessation can be estimated from the highest order mode resolved experimentally in the subsequent relaxation. “Highest order mode” because all other modes become relevant only at longer times when the stretch/orientation parameter Fso (defined by eq 1) has at least partly relaxed toward zero. We therefore link the ζ(Fso)/ζ(0) ratio existing just before flow cessation to the τp(ε̇)/τp(0) ratio for the highest order mode resolved experimentally (p = 3, 4, and 5 for ε̇ = 0.003, 0.01, and 0.03 s−1; cf. Figure 6b). The τp(ε̇)/τp(0) ratio is represented in Figure 6b by the vertical distance between the data point and the Rouse straight line; i.e., it is the ratio between the smallest τp in Table 1 to the corresponding value from eq 4. Recalling that relaxation times represent the ratio of friction coefficient to spring constant of subchains, we note that τp(ε̇) decreases not only due to the ζ reduction but also because of the FENE effect (increase of the spring constant of the subchain). Hence, the ζ(Fso)/ζ(0) ratio can be expressed in terms of the τp(ε̇)/τp(0) ratio and the FENE factor f FENE as

Figure 6. (a) Same data as in Figure 3 vertically shifted to test timerate separability at long times. (b) Characteristic times of the multiexponential fit of the data in Figure 3. The straight line is eq 4.

τp(ε̇) τp(0)

1, and the resulting curves are drawn in Figure 3, showing an excellent fit. The τp thus determined are double−logarithmically plotted against p in Figure 6b. For the smallest ε̇ = 0.001 s−1, τp is well described by the Rouse relationship shown as a solid line in Figure 6b; τp =

1.1 × 10

1 ζ(Fso) fFENE ζ(0)

(5a)

with fFENE =

1 2 1 − λ̃

(5b)

where, from eq 1, λ̃ = λ/λmax, and λmax is the maximum stretch ratio that is related to the number n of Kuhn segments between consecutive entanglements, i.e.:

3

p2

=

(in s) (4)

λ max =

This result is in harmony with the close coincidence between ‑ ‑ the Rouse curve and the σ E,fast (t)/σ E,fast (1) data for this ε̇ shown in Figure 3. However, Figure 6b also demonstrates that τp for the higher order modes (with p ≥ 3) at larger ε̇ deviates downward from the Rouse scaling of eq 4. This deviation gives the acceleration of the higher mode relaxation due to the ζ

n

(5c)

Notice that, because of flow, n is generally different from the equilibrium value n0 since the number of subchains (sliplinks) decreases due to CCR.16,23 On the other hand, the elongational stress σE(ε̇) just before the flow cessation is related to the stretch/orientation parameter Fso (=λ2̃ S̅) that determines the ζ(Fso)/ζ(0) ratio.

Table 1. Parameters of the Multi-Exponential Relaxation ε̇/s−1

τ1/s

τ2/s

τ3/s

0.001 0.003 0.01 0.03

1100 1100 1100 1100

275 275 273 270

122 75 68 57

τ4/s

15 10

τ5/s

g1

g2

g3

g4

g5

2.5

0.485 0.210 0.096 0.055

0.50 0.37 0.23 0.16

0.05 0.42 0.36 0.215

0.37 0.41

0.30

2777


Macromolecules

Article

asymptotic value of {ζ(Fso)/ζ(0)}/f FENE (={β/(1 + β)}γ for large Fso). Although eq 7 does not directly relate ζ(Fso)/ζ(0) to Fso because the f FENE factor cannot be determined from the experimental data, in the PCN-FENE simulations to be considered in the next subsections both f FENE and Fso are easily calculated at each time, thus allowing direct use of eq 7 to calculate (and then use) the friction reduction ratio ζ(Fso)/ ζ(0). A final comment on the empirical eq 7. In Figure 7, the data points are available only for Fso f FENE smaller than 0.75 and exhibit no leveling-off behavior. Thus, the above β value is obviously highly uncertain. However, the Fso f FENE factor never exceeded 0.8 in our simulations, and hence the uncertainty on β hardly affects the simulation results to be discussed next. 2.4. PCN-FENE Simulation of Elongational Flows. This section tests the hypothesis of the stretch/orientation-induced reduction of friction against literature data with the aid of the primitive chain network (PCN)13 simulation. The simulation code is identical to that used in our previous companion paper23 except for the treatment of the local friction. Polymer chains are described by a sequence of subchains having a molar mass comparable to the entanglement molecular weight Me. The polymer chains are bound through sliplinks to form a network. Dynamics of the system is described by: 3D motion of the sliplinks, 1D monomer transport through the sliplinks, destruction of sliplinks occurring when chain ends pass sliplinks, and creation of new sliplinks when a large number of monomers (=1.5n0)23 slides out of a terminal slip link. The kinetic equations for the 3D motion of the sliplinks and the 1D monomer transport through the sliplinks take the form

Specifically, with the decoupling approximation, σE can be expressed as σE = 3νkBTfFENE λ2SP n S̅ 2 = 3ν0 0 kBTfFENE λ̃ λ 2max n ϕP =

3Gen0 f Fso ϕP FENE

(6)

Here, ν is the number density of entangled subchains and ν0 its equilibrium value (linked by the relationship ν = ν0n0/n), kB and T are the Boltzmann constant and absolute temperature, respectively, and Ge = ν0kBT is the plateau modulus in the linear viscoelastic regime. In deriving eq 6, use has been made of eqs 1, 2, and 5c. As the f FENE factor appearing in eqs 5a and 6 cannot be directly evaluated from the elongational σE(ε̇) data, we will consider the plot of the τp(ε̇)/τp(0) data (={ζ(Fso)/ζ(0)}/ f FENE) against reduced stress data, σEϕP/(3Gen0) (=Fso f FENE), by utilizing the stress value σE(ε̇) at flow cessation and the known values of Ge, ϕP, and n0. Concerning the latter quantity, for the melt n0 is obtained from the ratio Me/MKuhn = 1.1 × 104/720 ≈ 15,23 while for solutions we adopt the scaling n0 = 15ϕP−1.3. This plot for the PS melt20 and the two solutions5 is shown in Figure 7. Clearly, the ζ reduction is observed in the

4

(ζs1 + ζs2)(Ṙ − κ • R) =

∑ Fi + Ff

+ FB (8)

i

ζs ρ

(9) 23

Figure 7. Plot of {ζ(Fso)/ζ(0)}/f FENE (=τp(ε̇)/τp(0) for the highest p available in Figure 6b for the melt, and =1 for solutions) against Fso f FENE (=σEϕP/3Gen0). The fitting curve is eq 7.

where, differently from the previous paper, the subchain friction coefficient ζs = nζ accounts for changes of monomer number n in the subchains due to fluctuations in the number of entanglements, especially important in fast flows because of CCR. (In eq 8 ζs1 and ζs2 each refer to one of the two subchains forming the binary entanglement.) In this way, for a constant value of the monomeric friction coefficient ζ, the overall chain friction coefficient stays constant as it should (and different from the previous paper23). Hence, we first tested the effects of this correction and found that for a constant ζ the results of the present code are numerically close to those obtained in the previous companion paper. The important novelty here is the fact that we also allow for the monomeric friction coefficient ζ (and consequently also for the subchain friction ζs) to change as a consequence of the flow-induced chain orientation/stretch, as discussed in the previous sections. The other symbols in eqs 8 and 9 have the usual meaning, i.e., R and Ṙ are the slip link position and its time derivative, respectively, κ is the strain rate tensor, Fi the subchain elastic force (including the FENE effect at high stretch), Ff the osmotic force, and FB the Brownian random force. In eq 9, ṅ is the rate of change of ni due to sliding of monomers from the (i1)-th to the i-th subchain, ρ the local linear monomer density, Fi the subchain tension, f f the osmotic force and f B the random force. The Brownian forces are related to ζs through the

melt case but not for solutions under the experimental conditions examined. In Figure 7, the solid curve shows the following purely empirical equation that fits the data: ζ(Fso) 1 1 = ζ(0) fFENE (1 + β)γ γ ⎧ 1 − tanh α(F ′so − F ′*so ) ⎫ ⎨β + ⎬ ⎩ ⎭ 2 (7a)

where F′so ≡ Fso f FENE can be calculated from experimental quantities as F ′so = σEϕp/(3Geno) = σEϕpMKuhn /(3ρRT )

n ̇ = (Fi − Fi − 1) + ff + fB

(7b)

In eq 7a, F′*so (=0.14) is the threshold value for the onset of the stretch/orientation-induced reduction of ζ, while α (=20) and γ (=0.15) specify how rapidly ζ decreases with increasing Fso. The remaining parameter, β (=5.0 × 10−9), specifies the 2778


Macromolecules

Article

Table 2. Parameters of PS samples sample

ref

M

Z0

n0

T [°C]

PS3.9M 10 wt % sol. PS100K PS200K PS390K PS200K−S PS145K

5 3 1 1 24 20

3.9 × 106 1.0 × 105 2.0 × 105 3.9 × 105 2.0 × 105 1.45 × 105

33 9 18 35 18 13

160 15 15 15 15 15

21 130 130 130 175 120

fluctuation−dissipation theorem, i.e.,⟨FB(t)FB(t′)⟩ = 2(ζs1+ζs2) kBTδ(t − t′)I (I = unit tensor) and ⟨f B(t)f B(t′)⟩ = 2ζskBTδ(t − t′). In the PCN-FENE simulation the above eqs 8 and 9 are made nondimensional by using the equilibrium distance a between consecutive entanglements as the unit length, and the subchain equilibrium relaxation time τ0 = ζs(0)a2/6kBT as the unit time. These equations are integrated numerically using a simple Euler-type integration scheme. The variable friction coefficient ζs = nζ(Fso) (appearing in the nondimensional equations through the ratio nζ(Fso)/n0ζ(0)) is updated at intervals of unit time (i.e., every τ0), and kept constant in between. (Notice that, although eqs 8 and 9 are integrated with a time step value much smaller than τ0 for obvious numerical reasons, results are significant only every subchain characteristic time τ0 in view of the coarse graining of the PCN model, which is at the level of the distance a between consecutive entanglements.) The friction ratio ζ(Fso)/ζ(0) is calculated through eq 7a after evaluating the stretch orientation parameter Fso (=F′SO/f FENE) from eq 7b. In both equations f FENE is calculated from eq 5b, while the stress appearing in eq 7b is obtained from the simulation as ⎛ rr ⎞ σ Z = fFENE ⎜ ⎟ ⎝n⎠ 3Gen0 Z0

G0 [Pa]

τ0 [s]

τR [s]

× × × × × ×

7.0 × 10−3 5.6 5.6 5.6 5.6 × 10−3 140

3.9 × 10−1 2.3 × 10 9.1 × 10 3.5 × 102 9.2 × 10−2 1.2 × 103

2.0 2.9 2.9 2.9 3.3 2.9

103 105 105 105 105 105

Figure 8. Viscosity growth and relaxation on start-up and cessation of uniaxial elongational flow for PS145K melt at 120 °C (left) and PS3.6M 10 wt % solution at 21 °C (right). Strain rates for PS145K are 0.001, 0.003, and 0.0003 s−1 from left to right and the maximum strain at the flow cessation is 3.0. Strain rate for PS3.6M solution is 11.7 s−1 and the maximum strains at the flow cessation are 3.3, 4.0, and 4.6 from left to right. Dashed and solid (red) curves, respectively, indicate the results of PCN-FENE simulation without and with the stretch/ orientation-induced reduction of friction. Experimental data taken from ref 20 (for PS145K melt) and ref 5 (for PS3.6M solution) are shown with symbols.

elongational viscosity at ε̇ > 1/τR. In contrast, the simulation results with the stretch/orientation-induced reduction of ζ (solid curves) are much closer to the data and reproduce the decrease of ηE cessation with increasing ε̇ (for ε̇ > 1/τR), which lends some support to the hypothesis of the ζ reduction. For the PS3.9M 10 wt % solution (right panel of Figure 8), the simulation results come out insensitive to the ζ reduction mechanism because, due to a low polymer concentration, the average orientation S̅ = ϕPSP is small and, at the same time, the maximum stretch ratio λmax is large (so that λ̃ does not become large enough, see eq 1). The simulation results well mimic the data in this case. Figure 9 shows the elongational viscosity growth function ηE+(t,ε̇) (up to the steady state) for PS100K, PS200K, and PS390K melts. The PCN-FENE simulation without the stretch/orientation-induced reduction of ζ (dashed curves) again overestimates the ηE+(t,ε̇) data1,3 (symbols) at high ε̇ (≫1/τR) and for times such that the Hencky strain ε exceeds the critical value εc ∼ 2, as already noted in our previous paper.23 Incorporation of the ζ reduction in the simulation (solid curves) significantly improves agreement with data. This improvement is most clearly noticeable for the PS390K melt (bottom panel). Figure 10 compares the simulated steady state viscosity ηE(ε̇) with the data for PS100K, PS200K and PS390K melts. Since the applicable range of the simulation is limited as ε̇ < 1/τ0, the predictions are not available in high ε̇ region. The viscosity obtained from the PCN simulation without FENE and ζ reduction (thin dotted curves) diverges at WiR (=ε̇τR) > 1, as noted in our previous paper.23 Incorporation of FENE in the simulation (thick dashed curves) suppresses this divergence but

(10)

where Z/Z0 is the ratio of the current (average) subchain number per chain to its equilibrium value, and r is the nondimensional subchain end-to-end vector. Except for the friction reduction extracted from relaxation data in the nonlinear range, and described by eq 7, all other basic parameters of the PS systems used in this work for the simulation and for comparison with experiments were evaluated with the method described in our previous paper,23 i.e., through a comparison between PCN simulation results and linear viscoelastic data. These parameters are summarized in Table 2 together with the values of M and T. The Rouse relaxation time of the chains, τR, also reported in Table 2, is linked to τ0 through: τR = (Z02/2π2)τ0.14 Further details of the simulation can be found in our previous work.23 Figure 8 compares the simulation results with the elongational stress growth and relaxation data reported for the PS145K melt20 and for the PS3.9M 10 wt % solution.5 (The relaxation data were used in previous sections to examine friction reduction effects.) The dashed and solid curves respectively indicate the results of the PCN-FENE simulation without and with the orientation/stretch-induced reduction of ζ. Here we omit the data for ε̇ > 1/τ0. For the PS145K melt (left panel of Figure 8), the simulation with a constant ζ (dashed curves) clearly overestimates the viscosity data. Furthermore, the simulated viscosity at flow cessation, ηEcessation, increases with increasing ε̇, whereas the experiments show the opposite trend. This overestimation has been observed also in our previous work23 for the steady 2779


Macromolecules

Article

Figure 10. Steady elongational viscosity for PS100K (top), PS200K (middle) and PS390K (bottom) at 130 °C plotted against strain rate. The PCN simulation results without FENE and ζ reduction are shown with thin dotted curves. The PCN-FENE results without and with the stretch/orientation-induced reduction of friction are indicated with thick dashed and solid (red) curves, respectively. Experimental data from ref 3 are shown with unfilled circles.

Moreover ζ should be taken of tensorial form, as mentioned earlier. Modification along these lines might improve agreement between simulation and data. Nevertheless, the simple idea of a reduction of the scalar ζ due to subchain stretch/orientation led to a remarkable overall agreement between the PCN-FENE simulation and the data (Figures 8−10). Before concluding on the extensional flows examined so far, it is informative to examine the relationship between the ζ(Fso)/ζ(0) ratio and the stretch/orientation parameter, Fso = λ̃2S̅. This relationship, underlying the plot shown in Figure 7, could not be specified on a purely experimental basis, as explained earlier. However, the PCN-FENE simulations reported in Figures 8−10 were based on evaluation of λ̃ and S̅ at each simulation step and thus allow us to make a plot of the ζ(Fso)/ζ(0) ratio against Fso that corresponds to the plot in Figure 7. This ζ(Fso)/ζ(0) vs Fso plot is shown in Figure 11 for various melt systems. The results collapse on a master curve, which suggests that the ζ(Fso)/ζ(0) ratio is essentially (though perhaps not rigorously) a function of the single parameter Fso. We also note that the ζ(Fso)/ζ(0) ratio sharply decreases on increasing Fso above a threshold value (≅ 0.15). This relationship between ζ(Fso)/ζ(0) and Fso is qualitatively similar to the empirical relationship between {ζ(Fso)/ζ(0)}/f FENE and Fso f FENE shown in Figure 7. It is also informative to quantify the effect of the ζ reduction on various quantities entering the PCN-FENE simulation. For this purpose, Figure 12 summarizes the viscosity reduction ηE/ ηE,0, the subchain reduced square stretch λ2̃ , the average orientation S̅, the FENE factor f FENE, the entanglement density reduction Z/Z0 (appearing in eq 10), and the ζ(Fso)/ζ(0) ratio, all obtained from the simulation of the PS390K melt in the steady state. These quantities are plotted against the τR-based Weissenberg number, WiR. The solid and dashed curves indicate the results with and without the ζ reduction. The ζ(Fso)/ζ(0) ratio starts decreasing with increasing WiR above unity (cf. bottom panel). This ζ reduction affects S̅ and f FENE

Figure 9. Uniaxial viscosity growth for PS100K, PS200K and PS390K melts at 130 °C. The strain rates are 0.003, 0.01, and 0.03 s−1 for PS100K, 0.001, 0.003, 0.01, and 0.03 s−1 for PS200K and 0.0003, 0.001, 0.003, 0.01, and 0.03 s−1 for PS390K. Experimental data (from refs 1 and 3) are shown with symbols. Dashed and solid (red) curves, respectively, indicate the PCN-FENE simulation results without and with the stretch/orientation-induced reduction of friction. Dotted curve is the linear viscosity growth function.

still leads to thickening of ηE(ε̇) at WiR> 1 where the melt data exhibit monotonic thinning. However, the PCN-FENE simulation with the stretch/orientation-induced reduction of ζ (thick solid curves) well mimics the monotonic thinning of the ηE(ε̇) data. Fair agreement of simulation results with data (Figures 8−10) lends support to the molecular mechanism of orientation/stretch-induced reduction of the local friction. The data for PS melts and solutions are consistently described by this idea. However, we also note a nontrivial deviation from the data, for example, for the PS145K melt (left panel of Figure 8) and for the PS100K melt (top panels of Figures 9 and 10). These deviations may be partly attributable to the numerical integration scheme in our simulation that somewhat delays the ζ reduction with respect to the actual orientation/stretch of the subchains, and hence overestimates the transient viscosity. It is also possible that the empirical procedure adopted to determine eq 7a lacks sufficient precision. It should finally be mentioned that, of course, the starting hypothesis whereby ζ is dependent only on the parameter Fso = λ2̃ S̅ (cf. eq 1) is open to possible revision. Indeed, ζ could be dependent on the stretch λ̃ and the orientation S̅ separately, or on a different combination of them. 2780


Macromolecules

Article

just moderately, but considerably suppresses the subchain stretch (compare the solid and dashed curves in the second top panel). The monotonic thinning of ηE results from a delicate combination of all these effects. (Note that f FENE is slightly larger than unity even in the low WiR region because of the small n0 value for PS melts, see Table 2.) 2.5. PCN-FENE Simulation of Shear Flows. It appears important to verify whether or not the friction reduction effect so far examined for elongational flows is also important for shear flows. In the shear case SP is defined as the positive difference between the two eigenvalues of the orientation tensor ⟨uu⟩ that correspond to the eigenvectors in the shear plane. Hence: SP =

{⟨ux 2⟩ − ⟨uy 2⟩}2 + 4⟨uxuy⟩2

(11)

where x and y are the shear and gradient directions, respectively. All other equations to be used in the simulation remain the same as for the elongational flow case. Figure 13 compares the PCN-FENE simulation results with and without ζ reduction (solid curves and plus symbols,

Figure 11. Reduction factor of the monomeric friction coefficient ζ plotted against the stretch/orientation parameter Fso obtained from the PCN-FENE simulations reported in Figures 8−10. Scales are logarithmic in the top panel and linear in the bottom panel. Different symbols indicate different PS systems.

Figure 13. Shear viscosity growth function for PS200K−S. The PCNFENE simulation results with and without the local friction reduction are shown with the solid curves (red) and the plus symbols, respectively. Data from ref 24 are shown by (gray) unfilled circles. Shear rates are 30, 10, 3, 1, and 0.3 s−1 from bottom to top.

respectively) to the PS200K−S melt data for the shear viscosity growth function η+(t,γ̇) .24 Simulations were made with the parameter values shown in Table 1, and cover the whole experimental range of WiR that extends up to 3.4. At such large WiR, ζ decreases in uniaxial elongational flow, and reduces the elongational viscosity (cf. middle panels of Figures 9 and 10). Conversely, under shear flow conditions, the simulated results with and without ζ reduction are indistinguishable, and in close agreement with the data (see Figure 13). In fact, in shear flow the stretch/orientation parameter Fso remains below the critical value (≅ 0.15; cf. Figure 11) and thus the ζ reduction mechanism incorporated in the simulation is essentially ineffective for WiR up to 3.4. This result shows that the magnitude of ζ reduction is not a function of WiR only, but is also dependent on the type of flow. This result also shows that the previous PCN-FENE simulations for shear flow 25 conducted without accounting for the ζ reduction remain valid even after considering this new effect.

Figure 12. Factors representing the viscosity reduction ηE/ηE,0, the reduced chain stretch squared λ̃2, the average orientation S̅, the FENE factor f FENE, and the local friction reduction ζ/ζ(0) obtained from the PCN-FENE simulation for PS390K in the steady state. These factors are plotted against τR-based Weissenberg number WiR. The simulation results with and without the stretch/orientation-induced reduction of friction are indicated with the solid and dashed curves, respectively. 2781


Macromolecules

Article

3. CONCLUSIONS We have attempted to elucidate the molecular origin of the different behavior in elongational flow between PS melts and solutions by focusing on a concept of stretch/orientationinduced decrease of the local friction ζ, as proposed by Ianniruberto and co-workers.17 Specifically, literature data of the stress relaxation after cessation of transient elongational flow, reported for both a PS melt20and two PS solutions,5 were analyzed to evaluate the magnitude of ζ reduction due to flow. In our working hypothesis, ζ was expressed as a function of an orientation/stretch factor Fso = λ2̃ S̅, where λ̃ is the normalized stretch ratio of the entangled subchains (defined with respect to the fully stretched state), and S̅ is an average orientational anisotropy of the components (polymer plus solvent, if any) in the system. From the initial stress decay data on cessation of elongational flow, the ζ(Fso)/ζ(0) ratio was obtained as an implicit function of Fso. This function, showing a significant decrease of ζ(Fso) upon increasing Fso (>0.15), was utilized in the primitive chain network (PCN) simulation, inclusive of the finite extensible nonlinear elasticity (FENE), to examine the elongational behavior of melts and solutions. For melts, the simulation indicated that ζ decreases significantly under fast elongational flow because the entangled subchains are short in PS melts, and therefore approach the fully stretched and fully oriented limit rather easily. As a consequence, the steady elongational viscosity ηE follows this decrease of ζ to exhibit a monotonic thinning even at ε̇ > τR‑1. In contrast, for solutions, the simulated ηE exhibits thinning followed by thickening at ε̇ > τR‑1 because the average anisotropy S̅ is governed by the solvent and remains small, thereby overwhelming the increase of the subchain stretch λ̃ and limiting the increase of Fso, which remains smaller than 0.15. The simulated results are in satisfactory agreement with the experiments, thus suggesting that the different behavior of ηE in PS melts and solutions reflects the difference of the local orientational field (largely damped by the solvent in semidilute solutions). Further tests for the change of local friction are clearly required, also to examine several other aspects, such as the influence of the molecular weight distribution, of branching, of a different chemistry, etc. For instance, uniaxial elongational experiments on entangled solutions of various polymer concentrations are critical for confirming the 0.15 threshold value of the Fso factor marking the onset of the ζ reduction effect. We are continuing our research along these directions and the results will be reported elsewhere.

■

their useful comments on the basis of which the paper has been improved.

■

REFERENCES

(1) Bach, A.; Almdal, K.; Rasmussen, H. K.; Hassager, O. Macromolecules 2003, 36, 5174. (2) Luap, C.; Müller, C.; Schweizer, T.; Venerus, D. C. Rheol. Acta 2005, 45, 83. (3) Nielsen, J. K.; Rasmussen, H. K.; Hassager, O.; McKinley, G. H. J. Rheol. 2006, 50, 453. (4) Bhattacharjee, P. K.; Oberhauser, J. P.; McKinley, G. H.; Leal, L. J.; Sridhar, T. Macromolecules 2002, 35, 10131. (5) Bhattacharjee, P. K.; Nguyen, D. A.; McKinley, G. H.; Sridhar, T. J. Rheol. 2003, 47, 269. (6) Ye, X.; Larson, R. G.; Pattamaprom, C.; Sridhar, T. J. Rheol. 2003, 47, 443. (7) Marrucci, G. J. Non-Newtonian Fluid Mech. 1996, 62, 279. (8) Marrucci, G.; Ianniruberto, G. Macromolecules 2004, 37, 3934. (9) Wagner, M. H.; Schaeffer, J. J. Rheol. 1992, 36, 1. (10) Wagner, M. H.; Kheirandish, S.; Hassager, O. J. Rheol. 2005, 49, 1317. (11) Wagner, M. H.; Rolón-Garrido, V. H. Korea Austr. Rheol. J. 2009, 21, 203. (12) Dhole, S.; Leygue, S.; Bailly, C.; Keunings, R. J. Non-Newtonian Fluid Mech. 2009, 161, 10. (13) Masubuchi, Y.; Takimoto, J.-I.; Koyama, K.; Ianniruberto, G.; Greco, F.; Marrucci, G. J. Chem. Phys. 2001, 115, 4387. (14) Doi, M.; Edwards, S. The Theory of Polymer Dynamics: Oxford University Press: Oxford, U.K., 1986. (15) Graessley, W. W. Adv. Polym. Sci. 1982, 47, 67. (16) Kushwaha, A.; Shaqfeh, E. S. G. J. Rheol. 2011, 55, 463. (17) Ianniruberto, G.; Brasiello, A.; Marrucci, G. Proc. 7th Annu. Eur. Rheol. Conf. 2011, 61. (18) Uneyama, T.; Horio, K.; Watanabe, H. Phys. Rev. E 2011, 83, 001800. (19) Hua, C. C.; Yang, C. Y. J. Polym. Res. 2002, 9, 79. (20) Nielsen, J. K.; Rasmussen, H. K.; Hassager, O. J. Rheol. 2008, 52, 885. (21) Inoue, T.; Osaki, K. Macromolecules 1996, 29, 1595. (22) Sobue, H.; Murakami, K. Kogyo Kagaku Zasshi 1961, 64, 2005. (23) Yaoita, T.; Isaki, T.; Masubuchi, Y.; Watanabe, H.; Ianniruberto, G.; Marrucci, G. Macromolecules 2011, 44, 9675. (24) Schweizer, T.; van Meerveld, J.; Ottinger, H. C. J. Rheol. 2004, 48, 1345. (25) Masubuchi, Y.; Furuichi, F.; Horio, K.; Uneyama, T.; Watanabe, H.; Ianniruberto, G.; Greco, F.; Marrucci, G. J. Chem. Phys. 2009, 131, 114906.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address ⊥

Fabricated Products Development Division, Mitsui Chemicals Inc., Chiba 299−0265, Japan Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful to Dr. H. K. Rasmussen and to Dr. T. Schweizer for providing us with the extensional flow relaxation data of PS145K and the start-up shear data of PS200K−S, respectively. We are also grateful to the anonymous referees for 2782


Primitive Chain Network Simulation of Elongational Flows of

Recommend Documents