Acceleration and Homogenization of the Dynamics ... - ACS Publications

Oct 19, 2017 - Paul Sotta,. † and Didier R. Long*,†. †. Laboratoire Polymères et Matériaux Avancés, UMR 5268 CNRS/Solvay, 87 avenue des Frèr...
0 downloads 0 Views 4MB Size
Article Cite This: Macromolecules XXXX, XXX, XXX−XXX

Acceleration and Homogenization of the Dynamics during Plastic Deformation Luca Conca,† Alain Dequidt,‡ Paul Sotta,† and Didier R. Long*,† †

Laboratoire Polymères et Matériaux Avancés, UMR 5268 CNRS/Solvay, 87 avenue des Frères Perret, F-69192 Saint-Fons, France Institut de Chimie, Université de Clermont-Ferrand, F-63000 Clermont-Ferrand, France



ABSTRACT: A coarse-grained model has been proposed recently in order to describe plastic properties of glassy polymers with space resolution the scale ξ of dynamical heterogeneities. This model allows for describing plastic flow by assuming that the elastic energy stored at the scale of dynamical heterogeneities reduces by a similar amount the free energy barriers for αrelaxation. The aim of this article is to consider in more details the evolution of the distribution of relaxation times under plastic deformation. This is achieved by taking explicitly into account the so-called facilitation mechanism during plastic flow introduced by Merabia and Long [Eur. Phys. J. E 2002, 9, 195; J. Chem. Phys. 2006, 125, 234901]. This mechanism, which allows for calculating the scale of dynamical heterogeneities, is key for explaining temporal asymmetry regarding the dynamical evolution of a glassy polymer upon cooling or upon heating, as observed experimentally by Kovacs. Upon heating, fast regions appear first and melt the slowest ones. We propose that this same process is at work upon stretching a polymer material. Some subunits get accelerated due to deformation, and their growing number allows for melting the slowest subunits. We show that this facilitation mechanism is responsible for the narrowing of the relaxation time spectrum observed e.g. during plastic strain by using optical probe dynamics measurements by Ediger and co-workers. The physical model is solved numerically in 3D by a method similar to dissipative particle dynamics, in which the strain and stress fields, as well as the local relaxation times, are calculated with a resolution the size of dynamical heterogeneities. The first model for describing plastic deformation has been the Eyring model,7 according to which the motion at the molecular level is biased by the stress. Within this model, the free energy barriers are not affected by the stress, and plasticization is controlled by the so-called activation volume v, an adjustable parameter without clear interpretation.8,9 The idea that plastic flow results from a stress-induced acceleration of the dynamics at the molecular level is supported by several recent experiments.10−14,16−19 Loo et al.10 have studied the dynamics in the amorphous phase of polyamide under stress by NMR. Kalfus et al.19 and Perez-Aparicio et al.20 have studied dielectric relaxation during plastic flow. They observed that secondary relaxations are not significantly affected during deformation, whereas the α-relaxation is modified: an increase of tan δ was indeed observed in the

I. INTRODUCTION Mechanical and dynamical properties of polymers have been intensively studied due to their fundamental and technological importance. They evolve rapidly as the glass transition temperature Tg is approached from above. As temperature decreases below Tg, the elastic modulus G′ increases from about 105 Pa to a few 109 Pa, while the dissipative modulus G″ increases from a low value (which depends on the considered polymer) up to about 108 Pa.1,2 Under stress, glassy polymers exhibit a yield point at a strain of a few percent, characterized by a maximum in the stress−strain curves.3 Typical yield stress values are a few tens of MPa. Yield is followed by the so-called strain-softening regime: the stress drops by a few tens of MPa before reaching a plateau. At even larger strain, depending on molecular weight and cross-linking, the stress may further increase.4,5 This is the so-called strain-hardening regime. These properties have been studied for many years due to their high practical importance.3,6 © XXXX American Chemical Society

Received: June 29, 2017 Revised: October 19, 2017

A

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

Article

Macromolecules

free energy barriers between elementary jumps. It is proposed that the elastic energy stored at the scale of dynamical heterogeneities effectively reduces the free energy barriers present for internal relaxation. In this article, as compared to ref 48, we discuss in more details the impact of mechanical deformation on the distribution of relaxation times. Experiments have shown that the dynamics is more homogeneous during plastic deformation than at rest. We shall see that this narrowing of the distribution of relaxation times can be interpreted as the consequence of the so-called facilitation mechanism introduced in refs 52 and 45 for calculating the scale of dynamical heterogeneities, the diffusion of small probes, and the violation of the Stokes law close to Tg. This facilitation mechanism also provides an explanation for the temporal asymmetry between aging and rejuvenating upon heating. As a consequence, due to similar mechanisms, we shall see that the evolutions of the dynamics upon stretching or upon heating a glassy polymer are very similar. In section II we introduce the notations and discuss the physical basis of the PFVD model.39,48 We emphasize the facilitation mechanism, which is an essential ingredient in this article. In section III we introduce the mechanical model derived from the PFVD model for calculating mechanical properties in both the linear and nonlinear regimes of deformation. By contrast to these earliers works, we extend the PFVD model by taking into accound the Merabia and Long facilitation mechanism, which is a key ingredient for the evolution of the relaxation time distribution under plastic deformation. The numerical method is described in section IV. The mechanical behavior and the evolutions of the distribution of relaxation times under stress are discussed in sections V and VI.

low frequency side of the observation window. They have proposed an explanation for the connection between dynamical heterogeneities and plastic behavior, indicating an enhanced mobility and a more homogeneous dynamics (in comparison to the unstrained system) during uniaxial extension. Ediger and co-workers11−18 have shown that the reorientation dynamics of small molecular probes is accelerated during stretching experiments and that the so-called stretching exponent β increases. They interpreted this result by the fact that the dynamics becomes faster and more homogeneous during plastic deformation. An important feature of glass transition in polymers is the strongly heterogeneous nature of the dynamics close to Tg.21,22 This has been demonstrated experimentally over the past 20 years by NMR,23−25 fluorescence recovery after photobleaching (FRAP),26−30 dielectric hole burning,31 or solvation dynamics.32 The characteristic size ξ has been estimated by NMR24 to be 3−4 nm at Tg + 20 K. Mechanical properties of polymer glasses have also been studied by molecular dynamics (MD) simulations.12,33−38 The possible link between dynamical heterogeneities and mechanical heterogeneities observed by Riggleman et al.12,35 has been discussed in ref 39. This idea has also been considered by Chen and Schweizer within the nonlinear Langevin equation model (NLE).8,9,40−44 This model reproduces many features of plastic deformation. Even though these studies support the fact that the dynamics is enhanced during plastic deformation, a detailed description of the dynamical behavior from the molecular level up to the scale of a few tens of nanometers during plastic deformation is still missing. Stochastic models have been developed by Merabia and Long for studying the evolution of the relaxation time spectrum during cooling and heating45 as well as under applied strain by Medvedev and Caruthers.46,47 To our knowledge, this stochastic approach has been the first one which allows for calculating the evolution of the distribution of relaxation times in out-of-equilibrium conditions.45 However, the spatial aspect of the problem is not explicit in this mean-field description. It lacks an explicit 3D spatial description which was developed by Dequidt et al. and Masurel et al.39,48,49 The physical model described here, which follows f rom earlier works,39,48 is solved numerically in 3D by a method close to dissipative particle dynamics that allows calculating the strain and stress f ields, as well as the local relaxation times, with a resolution the size of dynamical heterogeneities. To the best of our knowledge, this is the first model describing plastic flow in glassy polymers which takes into account dynamical heterogeneities with a 3D spatial description. In particular, it allows observing the appearance of shear bands on the scale 10−15 nm (which was not possible within the earlier stochastic model) and may in principle describe the cascade of relaxation events which lead to these shear bands. In recent publications,48,50 the authors have extended a mechanical model proposed by Dequidt et al.39 which accounts for dynamical heterogeneities in a polymer glass, up to plastic deformation.45,51,52 This is a mesoscale model in which the basic units are the subunits of dynamical heterogeneities of size typically 3−5 nm. When solved by numerical simulations, it allows reaching macroscopic time scales. In this model, denoted “percolation of free volume distribution” (PFVD) model in ref 42, the mechanical properties of polymers close to and below glass transition are calculated within an hypothesis different from that of Eyring: indeed, it was argued in ref 48 that the applied stress not only biases the motion but also lowers the

II. DESCRIPTION OF THE GLASS TRANSITION MODEL Relaxation Processes and Facilitation Mechanism. In this subsection, we summarize discussions introduced in refs 45, 50, and 52. The model derives from the point of view that glass transition in van der Waals liquids exhibits generic features as a consequence of generic thermodynamic properties and relaxation mechanisms at the molecular scale (see refs 53 and 54). Relaxation mechanisms are best qualitatively interpreted within the framework of the free volume model, which assumes that the slowing down is a consequence of the reduction of free volume when cooling the liquid. Throughout the paper, τα denotes the time which dominates the mechanical behavior of the liquid or melt at long times. The variation of τα is given by the empirical Williams−Landel− Ferry (WLF) law:1 ⎛ τ (T ) ⎞ −C1(T − T0) log⎜ α ⎟= C2 + T − T0 ⎝ τα(T0) ⎠

(1)

where T0 is a reference temperature, not to be confused with the Vogel temperature T∞ = T0 − C2. It has long been proposed that α-relaxation is a collective effect which takes place on a cooperative scale ξ.55 The basic assumption of the PFVD model is that density fluctuations at the scale ξ ≈ 3−5 nm generate a wide distribution of relaxation times. Based on the value of the bulk modulus, density fluctuations on scale ξ are of a few ±1%, as discussed in refs 45 and 50−52. Merabia and Long have proposed that two processes compete for allowing individual molecular B

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

Article

Macromolecules

required if the initial local density was maintained fixed by imaginary walls. Merabia and Long have proposed in ref 52 that the size ξ of dynamical heterogeneities results from the competition of these two distinct relaxation mechanisms. By considering the lifetime of density fluctuations (see Figure 2), an expression for ξ was found, expressed by the number Nc of involved monomers.45,52 The scale ξ = aNc1/3 (a is a monomer length) thereby is the smallest scale at which the lifetime of density fluctuations can be equal to or larger than τα. Density fluctuations on smaller scales are irrelevant for the α-relaxation process, being too short-lived. As a consequence of the facilitation mechanism, the relaxation time τα obeys the relation

jumps,45,50,52 as illustrated in Figures 1 and 2: (1) Internal reorganization within a slow subunit at f ixed f ree volume

τα = Nc 2/3τf

(2)

where τf satisfies

∫0 Figure 1. Schematics of a slow (high density) subunit of dynamical heterogeneity embedded in a faster (less dense) surrounding. The subunit can relax either by internal free volume reorganization or by free volume (or, equivalently, monomer) diffusion. The first mechanism is dominant in fast, low-density subunits, while the second one drives relaxation events in slow subunits.

τf

Q (τ ) dτ = qc

(3)

in which Q(τ) is the distribution of relaxation times and the fraction qc of fast subunits is an adjustable parameter, typically a few tens of percent. This relation means that the relaxation of the slowest subunits corresponds both to internal processes and to the melting by the faster environment. In ref 45, Merabia and Long have shown that the rejuvenation dynamics as measured by Kovacs can be correctly described by using the value qc = 0.3. The generic concept of facilitation, i.e., the idea that the mobility at some point may be accelerated by a faster local environment, is not recent (see the review by Ritort and Sollich56 and references therein). It has also been discussed more recently by Chandler and Garrahan.57,58 Merabia and Long propose an explicit description of this mechanism through free volume diffusion from neigbhoring faster subunits.45,52 This process, denoted “facilitation mechanism” in this context by Chen et al.,59 is thus an essential mechanism determining the dominant length scale ξ, of order 3−5 nm close to Tg.45,52 This concept has then been used by Tito et al.60 for describing mobility diffusion in the vicinity of free interfaces of suspended films. As a major outcome of our approach, we have shown that many different features may be explained by the same mechanisms, using identical or similar values for the adjustable parameters, which are essentially determined by the value of Nc. In ref 52, Nc was found to be of order 1000 monomers or equivalently ξ ≃ 3−5 nm. This value was calculated by a saddlepoint (or mini-max) method as noted by Lipson and Milner61 when considering two different relaxation mechanisms in competition. See also the corresponding discussion in the review article by Chen and Schweizer.42 This value is supported by NMR experiments by Spiess et al. (who found 3 nm at Tg + 20 K),24 by the Stokes law violation studied experimentally by Ediger et al.,52 by the kinetics of rejuvenating upon heating,45 and by values of Tg shifts in confinement and related reinforcement in filled elastomers.39,62 All these points are discussed in a detailed and self-contained way in Chapter 9 of the book Polymer Glasses edited by Connie Roth.50 From our point of view, the scale Nc ≈ 1000 monomers at Tg is also consistent with the experiments regarding the yield behavior of glassy polymers.48 Dynamical Heterogeneities and Mechanical Relaxation. The strong heterogeneity of the dynamics in supercooled

Figure 2. The relaxation time τd ∝ N2/3 associated to free volume diffusion (dash-dotted curve) sets the life-time of density fluctuations. On short length scales, one observes comparatively large density fluctuations, which should lead to large relaxation times if one could apply locally the WLF law as a function of the density (solid curve). On the other hand, on short length scales, density fluctuations are short lived. There exists an intermediate length scale on which one has large density fluctuations which are long lived. That corresponds to the scale of dynamical heterogeneities. Typical values have been found of order 1000 close to Tg in reference 52.

f raction. There is no reallocation of free volume on a scale larger than ξ on the considered time scale. This process is dominant in relatively fast subunits, in which there is a large amount of local free volume. (2) The second process is dominant in relatively slow subunits, when the local free volume fraction is too small for allowing an individual molecular jump on the corresponding time scale, at fixed free volume fraction. The process of dissolution by diffusion in faster neighboring subunits is much faster. Then, an individual molecule can jump only when the local free volume fraction has increased through free volume diffusion. This facilitation process occurs in a time that is shorter than the time which would be C

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

Article

Macromolecules

threshold cannot support stress on longer time scales because the environment has relaxed. As a consequence, the effective distribution of relaxation times is given by the “bare” distribution, with an effective cutoff at long times τα defined by (Figure 4)

liquids and/or polymers, with relaxation time distributions spanning many (up to height in some cases) decades, is supported experimentally by dielectric spectroscopy.22,49,63−65 Slow and fast regions coexist on a scale ξ of a few nanometers. Then, one key question is to determine the dominant relaxation time which is measured in a given experiment. We shall see indeed that the probed time scale depends on the considered experiment. In particular, in a mechanical experiment performed in a molten polymer, the dominant relaxation time may be defined by considering the final relaxation of the stress following an applied strain step. How can this dominant relaxation time be identified from the microscopic dynamic picture schematized in Figure 3 and in the distribution of relaxation times plotted in Figure 4?

∫τ



Q (τ ) dτ = pc

(4)

α

The time scale τα effectively controls the final relaxation of the stress. Shorter time scales contribute to the short time relaxation, while longer time scales are not relevant because they are too rare. The time required for melting slow subunits by fast neighboring ones is denoted τd: τd = Nc 2/3τf

(5)

where the relaxation time τf of fast subunits is given by eq 3. According to the above discussion, τd is equal to τα at equilibrium (see eq 2). In out-of-equilibrium conditions, the width of the bare relaxation time distribution changes, and eq 5 is no longer valid. During aging, the distribution of relaxation times is narrower than at equilibrium (see Figure 5). One has then typically Nc2/3τf > τpc. Then, the dominant relaxation time is τpc. Conversely, in rejuvenating experiments, the distribution of relaxation times is broader than at equilibrium, and one has typically Nc2/3τf < τpc. Then, the dominant relaxation τα is equal to Nc2/3τf = τd. In the general situation, the dominant relaxation time is given by

Figure 3. Coarse-grained model for glass transition in polymers, with subunits of size ξ ∼ 3−5 nm having different relaxation times, from fast (gray) to slow (green). The elastic energy stored in slow subunits does not relax by internal slow processes, but by rotation within a faster environment. Left: above Tg, clusters of slow subunits are relatively small. Right: below Tg, slow clusters percolate and determine the macroscopic viscosity.

τα = min[Nc 2/3τf , τp ] c

(6)

This relation expresses in a concise way the facilitation mechanism. The 2/3 power law (eq 5) was introduced in refs 52 and 45. It is related to the mechanism which sets the scale of dynamical heterogeneities. According to this mechanism, the lifetime of density fluctuations plays a key role. We proposed that this lifetime is set by free volume or mobility diffusion. Our point of view is that the mobility can diffuse from the boundary of a slow aggregate toward its interior. The elementary time scale of this diffusion process is set by the fast relaxation at the boundary. Within a slow (dense) aggregate of size ξ embedded in a fluid of relaxation time τf on the molecular scale a, an internal bead will be able to move after a time τf(ξ/a)2, which sets the lifetime τlife = τf(ξ/a)2 of the aggregate. The aggregate does not need to be large for this relation to be correct. By

In the mechanical experiment mentioned above, very fast relaxation times are not relevant, as their contributions to the stress relax very fast. Though they contribute to the response at short time scales, their contribution to the viscosity is very small, being proportional to their own relaxation times. On the other hand, very slow subunits are very rare and surrounded by an ocean of faster subunits. Final relaxation occurs before their own internal relaxation takes place because the stress that they support is transmitted by surrounding ones, which relax in comparatively shorter time scales. Therefore, the internal relaxation time of very slow subunits is not probed in such mechanical experiments. Accordingly, it has been proposed50,51,66,67 that the longest relaxation time probed in such experiments corresponds to a percolation threshold pc. Subunits slower than those corresponding to the percolation

Figure 4. Left: τ ̅ and τα as a function of temperature for PS. Right: distribution of relaxation times for PS, from 370 K down to 350 K. Characteristic times τ̅ and τα at T = 370 K are indicated by blue dashed lines. We define the glass transition temperature Tg the temperature at which the dominant relaxation time is 1 s. Thus, Tg = 370 K. D

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

Article

Macromolecules

reduces the free energy barriers present for internal relaxation. The local relaxation time is accelerated according to ⎛ ξ 3 σ:σ ⎞ τ({ρ}, σ ) = τ({ρ}) exp⎜ − ⎟ ⎝ 2kBT G0′ ⎠

(7)

where τ({ρ}) is the relaxation time that the dynamical heterogeneity would have in the absence of stress and, in principle, depends on the whole density spectrum {ρ}. The stress σ is the local stress calculated for each subunit. The stress-dependent quantity σ:σ is defined as σ:σ =

1 [(σ1 − σ2)2 + (σ2 − σ3)2 + (σ3 − σ1)2 ] 2

(8)

where σ1, σ2, and σ3 are the eigenvalues of the stress tensor σ. Note that for pure shear, where only one eigenvalue σ1 is nonzero, σ:σ = σ12. Note also that eq 8 is similar to a von Mises rupture criterion.69 A yield stress value σy ≈ 2kBTG0′ /ξ 3 of a few 10 MPa can be estimated from eq 7, without any additional adjustable parameter. We proposed to replace the usual Eyring expression by the expression (eq 7) in ref 48. In the Eyring picture, the effect of the stress on the dynamics is controlled by the so-called activation volume v, which is an adjustable parameter without clear interpretation, as pointed out by Chen and Schweizer.8,9 It is usually used for fitting the stress in the regime of plastic flow or at yield (see e.g. refs 4 and 70−72). When so doing, the linear regime cannot be described with the same expression extrapolated at zero stress, as explained in ref 48. There is indeed a discrepancy of 3−4 orders of magnitude if the Eyring expression fitted in the plastic regime is applied in the linear regime. From our point of view, this failure comes from the Eyring model itself, which describes the diffusion motion of a particle trapped between barriers.7 In the presence of an external potential, the motion is biased by a drift which results from the difference ±σv between the free energy levels of two consecutive sites. The detailed balance condition for the exchange of two neighboring positions then results in the Eyring equation. This picture thus compares the energy states in the bottom of the wells, and it would be fine for describing the linear response.73 It does not take into account the modification of the free energy barriers by the applied stress. Hence, the large discrepancy when using the Eyring model over both linear and nonlinear regimes. Only the drift term is present. It is very often interpreted as a change of the free energy barrier in the deeply nonlinear regime, where it is usually applied. However, according to our picture, this change must be quadratic and comes in addition to the biased motion, as explained in ref 48. Shorter time scales results from the lowering of free energy barriers. As discussed above, the Eyring picture does not describes this effect. More direct comparison between the Eyring model and the prediction of our model will be discussed in future work. Distribution of Relaxation Times. To be specific, we assume that the distribution of relaxation times τs at equilibrium can be represented by a log-Gaussian probability distribution of central value τ ̅ and width Δ:48,49

Figure 5. Top graph: Distribution of relaxation times computed after ageing for 104 s and no applied stress (blue curves), as compared to the reference (theoretical) log-Gaussian distribution (red curve), at 355 K ≈ Tg − 20 K. After ageing for 104 s, τd concerns only slowest subunits, negligibly affecting the distribution of relaxation times. Bottom graph: The reciprocal relaxation time τ−1 as a function of the age t of the subunit (continuous line), as estimated from Equation 13. 3 Pa-rameters are qc = 0:3 and N2/3 c = 10 .

definition, the lifetime of dynamical heterogeneties is long enough for the diffusion equation to hold. This idea is supported by theoretical works on how hydrodynamic behavior emerges in an ensemble of microscopic (molecular) objects. Though, to our knowledge, rigorous (exact) solutions do not exist yet, Lebowitz et al. have shown by using mean-field arguments that the density of hard objects with exclusion relaxes according to a standard diffusion equation (even though individual hard spheres do not).68

III. DESCRIPTION OF THE PHYSICAL MODEL The physical ingredients of the model which allows for calculating the dynamical evolution during deformation are summarized here. For more details, we refer to refs 39, 48, and 50 and references therein. Effect of the Stress on the Relaxation Time. The wide distribution of relaxation times is created by density fluctuations on the scale ξ ∼ 3−5 nm.45,52 Each dynamical heterogeneity at the scale ξ is described by a high-frequency elastic modulus G′0 ∼ 109 Pa and a local relaxation time. Entangled or cross-linked polymers also have a low-frequency modulus of order 105 Pa. The model is schematized in Figure 3, which shows subunits with different relaxation times (a similar picture was shown in refs 67 and 39). Mechanical experiments do not probe the internal relaxation times of the slowest subunits because the stress relaxes according to the (faster) relaxation times of their neighbors. It is experimentally known that an applied stress can induce faster relaxation than at equilibrium.3 This implies that the local relaxation time cannot be a function of the local density only. In ref 48, it has been proposed that the elastic energy stored at the scale of dynamical heterogeneities ξ ∼ 3−5 nm effectively

peq (Log(τs)) = E

⎛ Log 2(τ /τ ) ⎞ 1 s ̅ ⎟ exp⎜ − πΔ Δ2 ⎠ ⎝

(9)

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

Article

Macromolecules

The relaxation time associated with a subunit of age tw and local stress tensor σ is thus finally

with Log = log10. We further assume that decreasing temperature slows down the dynamics only by shifting peq toward longer time scales. The central value τ ̅ is computed through the percolation equation

τ(tw , σ , {τs}) = min[τs(tw , σ ), τd({τs})]



and the relaxation probability during dt is given by

∫Log(τ ) peq (Log(τs)) d Log(τs) = pc

(10)

α

dPrup(tw , σ , {τs}) =

with pc ≈ 0.11 the percolation threshold.66 The dependence of peq on temperature is incorporated in τ,̅ via τα, as shown in Figure 4. The constant value Δ = 4 has been chosen, as in ref 48. At each time step dt, a glassy spring relaxes with a probability Prup(t,σ), a function of the age t and of the local stress tensor σ. When a spring relaxes, its age t is set to a lower cutoff value tmin and the glassy stress is released. Otherwise, the spring ages further by dt. As proposed in refs 39 and 48 in the absence of stress, the probability for a spring of age tw to relax is given by ⎛ d ln p (τ ) ⎞ eq ⎟ dt dPs(tw , 0) = −⎜⎜ ⎟ dτ ⎝ ⎠τ= t

=

where peq is given by eq 9. This equation allows for calculating the relaxation time of a spring which has an age t. As discussed in refs 39, 48, and 74, aging and relaxation processes in glassy polymers in the linear regime can be computed in this way. According to eq 7, applying a stress lowers the free energy barrier for the α-relaxation.48 The relaxation probability per unit time thus becomes ⎛ d ln p(τ )eq ⎞ ⎟⎟ e λσ : σ dt dPs(tw , σ ) = dPs(tw , 0)e λσ : σ = −⎜⎜ dτ ⎝ ⎠τ = t w

(12) −3

−2

where λ = ξ /2G0T ≈ 10 MPa (as for example for PS). The deviatoric component of the stress σ appears in the argument of the exponential, which becomes of order one for a stress of a few tens of MPa and a strain of a few percent. The relaxation time τ from (7) depends on τs and τd. τs(tw,σ) is computed according to (12). To summarize, a probability dPs(tw,σ) is associated with τs(tw,σ) according to

τs(tw , σ ) =

dt dPs(tw , σ )

dt τ(tw , σ , {τs})

dt min[τs(tw , σ ), τd({τs})]

(15)

Equations 11−15 altogether contain the effect of aging, the impact of the stress on the acceleration of the dynamics, and the facilitation mechanism (fusion of slow subunits by fast ones). Equations 11 and 12 only, which describe the acceleration of the local dynamics by the applied stress, have been considered in ref 48. To include the facilitation mechanism, which was not implemented in ref 48, the “bare” relaxation time distribution is first calculated from eqs 11 and 12. Then, the time τd, which corresponds to the long time cutoff in this distribution resulting from the facilitation mechanism, is calculated from eq 6. The aim of this article is now to show the impact of the facilitation mechanism on calculated mechanical properties and to discuss this impact in comparison to experimental results. The distributions of relaxation times at equilibrium and after aging 104 s, calculated at T = 355 K ≈ Tg − 20 K, are shown in Figure 5. The values of τf and τd are indicated for both distributions. τd acts as a cutoff at long times for the effective relaxation times. For the distribution aged during 104 s, this cutoff does not change the effective distribution of relaxations times much, as it concerns only the slowest subunits. The continuous line on the bottom graph gives the reciprocal relaxation time τ−1 as a function of the age t of the subunit, as calculated from eq 11 at zero stress, using the parameters qc = 0.3 and Nc2/3 = 103. These two quantities must be considered as adjustable parameters of the model. The choice of qc has been discussed in ref 45. The value chosen for Nc2/3 will be discussed later in the paper. The distributions of relaxation times of systems aged for 104 s at various temperatures below Tg, from Tg − 25 K up to Tg − 5 K, are plotted in Figure 6. This shows that the aging process is indeed observed in the presence of the facilitation mechanism.

(11)

w

3

(14)

(13)

This equation gives a local relaxation time for each dynamical subunit which depends on its age tw and is accelerated by the local stress. A “bare” distribution of local relaxation times at the considered time t, Qσ(τ,t), is thus obtained. Note that these relaxation times are calculated locally, at the level of subunits and for each glassy springs, which have all their own age, stress level, and relaxation times which evolve during the course of the simulations. The “bare” relaxation frequency can be calculated from eq 11 based on the equilibrium distribution of relaxation times at rest. This frequency is accelerated by the effect of the stress according to eq 12. The “bare” distribution, however, only depends on the tw and local stress. On the other hand, according to the facilitation mechanism, subunits may also relax because they are surrounded by fast subunits. The relaxation time τf, and thus also τd = Nc2/3τf, can be calculated from eqs 3 and 5.

Figure 6. Distributions of relaxation times at rest for different temperatures between 345 and 365 K. The systems have aged 104 s at the considered temperature. The facilitation mechanism is taken into account during the aging process. Despite the presence of the cutoff τd, the systems show regions with very slow dynamics at low temperatures. F

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

Article

Macromolecules

distribution of relaxation times of a system approaching equilibrium is reproduced by letting glassy springs (which represent dynamical heterogeneities) age during a time tw = 104 s. Hard core repulsion is described by the following interaction potential:

Note that the results regarding the distribution of relaxation times in e.g. Figure 6 are somehow bumpy. This is due to the finite size of the samples, while relaxation times have very broad distributions. Note also that the peak at very short time scales in the curve at 365 K in Figure 6 comes from the short time cutoff in our simulations and has no physical significance. Note that we define the Tg of our system as the temperature at which the domnant relaxation time is 1 s. This is of course arbitrary. This temperature is 370 K. Thus, the curve at 365 K in Figure 6, which is only 5 K below Tg, indeed reaches the equilibrium distribution, with a dominant relaxation time for mechanical properties close to 10 s, corresponding to the equilibrium relaxation time at Tg − 5 K. Relevant Time Scales. The various time scales introduced for discussing the physics of the model are summarized in Table 1. τK and τc will be defined later in the text. They are

⎛ d ⎞12 UHS(r ) = u0⎜ ⎟ + Ureg(r ) ⎝r⎠

where u0 scales the energy and d ∼ ξ the interaction range. This potential is cut off at a distance rHS = 1.08ξ and regularized through the harmonic term Ureg, so that UHS ″ (rHS) = UHS ′ (rHS) = UHS(rHS) = 0: ⎤ ⎛ d ⎞12 ⎡ ⎛ r ⎞2 r Ureg(r ) = −u0⎜ ⎟ ⎢78⎜ ⎟ − 168 + 91⎥ ⎥⎦ rHS ⎝ rHS ⎠ ⎢⎣ ⎝ rHS ⎠

Table 1. Summary of the Times Used in the Papera τs τ̅

intrinsic local time, function of density ρ (or age), stress σ, temperature T most representative relaxation time

τf τd τpc

fast time relevant for diffusion maximum lifetime of density fluctuations slow time relevant for mechanics

Log τ ̅ = max peq(Log τs) τ ∫ 0fQ(τs) dτs = qc τd = Nc2/3τf ∫∞ τpcQ(τs) dτs = pc

τ τα

true local relaxation time average relaxation time observed in mechanics

τ = min[τs, τd] τα ∼ min[τd, τpc]

τK

KWW time for probe reorientation

τc

average time for probe reorientation

F(t) = exp(−t/ τK)β τc = ∫ ∞ 0 F(t) dt ∼ τα

(16)

freg (r ) = 12

12 ⎞ u0 ⎛ d ⎞ ⎛ r − 14⎟ ⎜ ⎟ ⎜13 rHS ⎝ rHS ⎠ ⎝ rHS ⎠

(17)

(18)

The low-frequency rubbery response is elastic with a low modulus G∞ ′ ∼ 0.1 MPa. The rubbery network is made of permanent rubbery springs of average connectivity 12. For each bead the minimum connectivity is 7, so that the rubbery network is sufficiently homogeneous. The interaction potential of rubbery springs is given by

UR =

τs and τ are local times characterizing a single dynamical heterogeneity. τ,̅ τf, τd, and τpc are times characterizing the distribution of relaxation times. τα, τK, and τc are times determined from experiments or simulations. Note that τpc ≈ τd ∼ τα only at equilibrium. a

1 k∞(r − l0)2 2

(19)

l0 is taken as the average distance to the 12 nearest neighbors in a face-centered cubic lattice of parameter ξ. The rubbery stiffness k∞ = G′∞ξ is the same for all springs. The high-frequency response is modeled by glassy springs described by a potential

related to a correlation function F(t) (defined later) which describes the overall relaxation in the system, based on the distribution of relaxation times. At equilibrium, i.e., after infinite aging time, τpc, τd, and τα are all equal, while, in the general case, τα = min[τd, τpc]. In the situations studied here, under deformation, τd is shorter than τpc.

Ug =

1 k 0( r ⃗ − r0⃗ )2 2

(20)

where r0⃗ is the reference position of the spring and the glassy stiffness k0 = G′0ξ ≫ k∞ is the same for all springs. Glassy springs are created during the preparation of the system, up to an average connectivity of 7. At each simulation time step dt = 10−4 s, glassy springs may relax by breaking. A breaking event occurs randomly with a probability dPrup(t,σ) computed through eq 15. After relaxation, the age of the spring equals a lower cutoff τmin, set to the extremity of the left wing (3%) of the theoretical relaxation time distribution. In addtion, the equilibrium position r0⃗ = r,⃗ so that the glassy stress instantaneously vanishes. Upon deformation, all beads are first displaced affinely. Beads are then iteratively displaced toward the positions of minimum energy during a time step dt. Integration of the equation of motion is performed with the Euler method, with a secondary time step d2t = dt/10. In the overdumped limit, the equation of motion for a bead i reads

IV. NUMERICAL METHOD In this section, we describe how we solve numerically the dynamics for each subunits represented by beads surrounded by glassy springs and rubbery springs. The method is analogous to dissipative particle dynamics and has been described in details in refs 39 and 75−77 and is the same as that used for solving the dynamics in 3D in ref 48. The elementary units are beads interacting through springs and hard core repulsion. Two types of springs are considered: “rubbery” springs model the elastic response of cross-linked or entangled polymers far above the glass transition temperature Tg, and “glassy” springs are responsible for the high modulus below Tg. The unit length in the simulation is naturally defined as ξ = 5 nm. The simulated box has dimensions (15ξ)3 = (75 nm)3. 3D periodic boundary conditions are applied. Before running any simulation, the systems are prepared by iteratively relaxing the velocities of the beads and the deviatoric component of the stress tensor, as detailed in ref 39. Once relaxation has been achieved, the

d2 ri ⃗ =

1 ⃗ 2 Fi d t ζi

(21)

where F⃗i is the sum of the forces acting on the bead. The viscous friction of bead i is controlled by the slowest glassy spring (heterogeneity) connected to i, so that ζi = k0 maxj τα,j. G

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

Article

Macromolecules We calculate thus in 3D the strain and stress fields with the resolution the size ξ of dynamical heterogeneities, with their own local relaxation times, which fluctuate and evolve in time. We obtain then a spatial distribution of relaxation times, which age and relax randomly during the course of the simulations. The values of the parameters used in this paper are summarized in Table 2. The parameters are the same as those considered in ref 48. The elastic modulus calculated from the model with the WLF parameters of polystyrene is given in the Appendix.

To be specific, the chosen C1, C2, and T0 WLF parameters are those of polystyrene (PS), and the chosen Poisson coefficient of 0.3 is representative of glassy polymers.1 Evolution of the Yield Stress. Stress−strain curves obtained either with (continuous curves) or without (dashed curves) the facilitation mechanism (dashed line) are shown in Figure 7a. Below Tg ∼ 370 K, the slopes at small strains are of the order of 1−2 GPa. The linear regime of deformation is negligibly affected by the facilitation mechanism: at equilibrium, the diffusion time τd concerns a small fraction of subunits, which have large relaxation times since τd is not modified significantly by the applied strain. Conversely, the facilitation mechanism has a strong impact at the onset of the plastic regime. At temperatures below Tg, a maximum of the stress of a few tens of MPa is attained at strain values of a few percent. As shown by the comparison of continuous and dashed curves, both the yield stress and yield strain are smaller when the facilitation mechanism is included. In addition, the strain softening is more pronounced when including the facilitation mechanism, resulting in values of the flow stress much lower than without the facilitation mechanism. Stress−strain curves for systems where the facilitation mechanism is included are shown in Figure 7b for various temperatures between Tg − 12 K and Tg. Note that at T = 350 and 370 K equilibrium distributions may be reached in accessible time scales, both experimentally and in our simulations. At and below Tg, the equilibrium relaxation time increases by about 2 decades when lowering temperature by 10 K. Hence, the huge difference between 350 and 370 K may be observed. Below 350 K, we expect not so marked differences as compared to 350 K because, as it is the case experimentally, the samples cannot reach equilibrium and undergo an aging process which does not depend significantly on temperature (see ref 45 and references therein regarding the work by Struik). Then the dynamical state of our samples is controlled by the aging time rather than the temperature itself. Aging time is limited by the accessible time scale, e.g., 106 s of aging. For instance, we expect that our model may describe experiments performed at room temperature on polycarbonate samples (for which Tg = 420 K) aged typically for a few 106 s, such as those described in ref 20. The distribution of relaxation

Table 2. Main Simulation Parameters and Correspondence with Physical Quantitiesa simulation parameter

physical equivalence

numerical value

thermal energy unit length ξ fast percolating fraction qc surface of a domain Nc2/3 glassy stiffness k0 rubbery stiffness k∞ rubbery length at rest l0 average rubbery connectivity minimum rubbery connectivity excluded volume radius d excluded volume energy u0 excluded volume cutoff rHS time step dt secondary time step d2t Poisson coefficient ν

∼Tg 3 nm 30% ∼1000 monomers G0ξ ∼ 3 GPa·ξ G∞ ′ ξ ∼ 105 Pa·ξ (4.8 nm)

∼0.18 1 0.3 103 3000 0.1 1.6 12 7 1 ∼0.01 1.08 10−4 10−5 0.3 12.0 50.0 375.0

WLF parameters

ξ (∼1/20T) (3.26 nm) 10−4/s 10−5/s 0.3 C1 = 12.0 C2 = 50.0 K T0 = 375.0 K

a

Parentheses are put around numerical value used in simulation, but without a particular physical correspondence. More specifically, u0 is set so that the contribution to the stress tensor of the hard-sphere repulsion is much smaller than the glassy one, and l0 is a convenient average distance for the given connectivity, with negligible influence on the behavior of the system close to Tg.

V. MECHANICAL PROPERTIES The mechanical properties based on the 3D model described here are investigated under compression at constant strain rate.

Figure 7. (a) Stress−strain curves obtained during compression tests at constant logarithmic strain rate 0.1 s−1, for simulations with (continuous curves) or without (dashed curves) the facilitation mechanism, at different temperatures. When comparing curves at a given temperature, the stress is systematically lower when the facilitation mechanism is implemented. This is a consequence of the fact that the facilitation mechanism enhances relaxation processes due to the applied stress. (b) Stress−strain curves obtained during compression tests at constant logarithmic strain rate 0.1 s−1, for simulations with the facilitation mechanism, at various temperatures. The facilitation mechanism tends to enhance the stress overshoot as compared to simulations without this mechanism. H

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

Article

Macromolecules

Figure 8. Stress−strain curves at different strain rates: (left) at T = 350 K; (right) at T = 355 K. We observe that the stress depends on the strain rate after yield as it is the case experimentally.

times would not be different from that of a sample aged for the same duration at 390 K, for instance. Stress−strain curves at different strain rates obtained at 350 and 355 K are shown in Figure 8. We observe the standard result that both the yield stress and the stress past yield are an increasing function of the strain rate as can be seen for instance in ref 78. Note that by using the stochastic model introduced by Merabia and Long45 and Medvedev and Caruthers,46 Caruthers and Medvedev obtain a stress during postyield plastic flow which is a constant independent of the strain rate. The evolution of the yield stress σy as a function of the logarithm of the strain rate γ̇ is shown in Figure 9. The curves

is a more intrinsic quantity of the model than the slope itself. In Figure 9, we find that a ≈ 0.1, which is in good agreement with reported data (see Figures 17, 18, and 20 in ref 71, Figure 1 in ref 72, and Figure 4 in ref 70; see also refs 79 and 80). As regards the temperature dependence of the slope in Figure 9, Δσy/σyΔT ≈ 0.03 K−1 in our model, while Δσy/σyΔT ≈ 0.01 K−1 for the data reported in ref 71, Figure 17 for instance. The orders of magnitude are thus similar. Note, however, that this quantity depends on the WLF law of each polymer and that the considered polymers are different. Evolution of the Lower and Upper Cutoffs. The evolution of τd = Nc2/3τf is a measure of the overall number of relaxation events which have taken place. The evolution of τpc gives information about the fraction of slow subunits still present in the system. τd is plotted in Figure 10 as a function of

Figure 9. Yield stress σy as a function of strain rate, at different temperatures, when considering the facilitation mechanism. σy is nearly proportional to the strain rate over several decades. The σy values are smaller than the corresponding ones given in ref 48, for which no facilitation mechanism was taken into account. The slope is also smaller than and comparable to experimental data (see e.g. ref 70).

Figure 10. Evolution of the diffusion time time τd as a function of strain during uniaxial compression at constant logarithmic strain rate 0.1 s−1 at different temperatures between 358 and 370 K.

the strain amplitude. Above Tg, τd remains nearly constant when compressing the system. The amount of fast regions already present is so high that adding further fast regions does not affect significantly τf, as defined by eq 5, because the strain rate is comparable to or smaller than the spontaneous relaxation frequencies of the system. Below Tg, e.g., at 358 or 362 K, τd strongly varies during compression. Starting from relatively high values before deformation, it decreases rapidly under applied strain. This is the consequence of the yielding of many subunits, which acquire faster relaxation times under stress. For strain amplitudes corresponding approximatively to the yield point, at maximum stress, τd reaches a minimum. The lower is the temperature, the higher is the yield stress and the lower τd at yield. In the postyield regime, τd increases again in the systems which show a pronounced strain softening (358 K, red curve, see the corresponding stress−strain curve in Figure 7) as a consequence of the steep drop of the stress.

are nearly linear, which is a classical feature. When compared to results presented in ref 48, for which the facilitation mechanism was not considered, yield stress values are systematically smaller and also the slope of the curve is smaller. This decrease of the slope of σy vs γ̇ obtained when the facilitation mechanism is implemented, as compared to the results in ref 48, shows that the whole relaxation time distribution is shifted toward shorter times by the stress more efficiently when the facilitation mechanism is implemented. Typically, the yield stress increases by 10−15 MPa every 2 decades of shear rate, whereas the corresponding slope in48 was 2−3 times too large. The slope normalized by the yield stress defined by a=

Δσy σy Δ Log γ ̇

(22) I

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

Article

Macromolecules

Figure 11. Distributions of relaxation times τs obtained during uniaxial compression for different strains from 0 to 10%, at a strain rate 0.1 s−1, at two different temperatures, without the facilitation mechanism. Even at 1% strain, the relaxation time distributions are shifted toward shorter time scales, as a consequence of the acceleration induced by the local stress. However, the distributions remain broad even at relatively large deformations.

Figure 12. Distributions of relaxation times τs obtained during uniaxial compression for different strains from 0 to 10%, at a strain rate 0.1 s−1, at two different temperatures. The facilitation mechanism was implemented in the simulations. The long time part of the distribution is suppressed beyond a few percent of deformation and that diffusion time τd acts as a cutoff for the distribution, at about a few seconds in these examples. Note that the distributions are narrower as compared to those in Figure 11.

The results plotted in Figure 10 suggest an interpretation of the mechanical properties observed in Figure 7. At small strains and below Tg, the acceleration of the dynamics induced by the stress affects mainly short and intermediate time scales. At this stage τd remains large, and mobility diffusion (facilitation mechanism) is not enhanced as compared to the linear regime of deformation. Mechanical properties with and without the facilitation mechanism are thus essentially identical. When approaching the plastic regime, the large amount of fast regions which have been created enhances the facilitation mechanism: a higher fraction of faster subunits melts more efficiently slow subunits. τf decreases fast, and τd is shifted toward shorter time scales. The dynamics of all subunits with relaxation times larger than τd is thus accelerated. Therefore, in this case, the facilitation mechanism accelerates the dynamics more efficiently, as compared to the case where it is not implemented. In addition, the increased amount of relaxation events triggers earlier plasticization, as can be seen in Figure 7. Distributions of Relaxation Times. The distributions of relaxation times, computed at various strain amplitudes without the facilitation mechanism, are plotted in Figure 11 for two different temperatures (355 and 360 K). The relaxation time distributions remain broad up to a deformation amplitude of 10%. The shift toward shorter time scales, even at small deformations, is due to the effect of the exponential acceleration induced by the stress. A relatively small fraction of very slow subunits relax under the applied stress, whereas the

dynamics of a majority fraction of them is simply accelerated by the stress. The distribution of relaxation times computed with the facilitation mechanism are plotted in Figure 12 at the same temperatures and range of strain amplitudes as in Figure 11. As compared to Figure 11, the relaxation time distributions are shifted toward shorter times more homogeneously. All regions with relaxation times τ longer than τd are more likely to relax. This process happens in the same way for all slow subunits; they gradually disappear after relaxation upon deformation, and the distributions become nearly flat at long time scales. A cutoff, corresponding roughly to τd as given in Figure 10, appears in the relaxation time distributions. At 360 K, the cutoff is roughly constant and equal to 1 s. At 355 K, the cutoff becomes smaller than 1 s at 2% strain, and then it gradually increases with increasing strain. This is a consequence of the drop of the stress during strain softening and of the subsequent aging process. On the other hand, at 1% strain amplitude, no clear cutoff can be identified because the facilitation mechanism is not yet enhanced. In our model, the evolution of the distribution of relaxation times results from a competition between stress induced relaxation and ongoing aging process. As embodied in eq 11, aging is driven by temperature. Aging may be stopped by the applied stress, which may even lead to rejuvenation according to eq 12. However, as soon as the local stress has relaxed, the aging process starts again according to eq 11. The facilitation mechanism comes also into play. The evolution of the J

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

Article

Macromolecules

For systems at rest, the relaxation function F(t) can be calculated from eq 23 using the theoretical distribution of relaxation times:

distribution of relaxation times results from the interplay of all these processes.

VI. DYNAMICS UNDER APPLIED STRAIN Probing Segmental Dynamics. In experimental studies by Ediger and co-workers, the time correlation function CF(t) of fluorescent probe orientation is measured.11,13,14,16 This function depends on the relaxation time spectrum. Because of their size, the probes interact at any time with relaxation times slower than the fast-time part of the spectrum, which gives a cutoff at short times. Note also that the behavior at very short times is difficult to access experimentally. We thus assume that there is a cutoff at short time scales comparable to the time scale τf introduced above (see Table 1 and Figure 5). For molecules of size close to ξ, this cutoff corresponds roughly to a fraction qc and is set by τfast.45,52 The relaxation function, representative of the average mobility at the scale of dynamical heterogeneities, is thus given by F(t ) = ⟨e−t / τs⟩{τs> τf } =

∫τ

τp

c

e−t / τsp(τs) dτs ≈

f

1 N



Feq(t ) =

∫τ

f

τp

c

e−t / τspeq (τs) dτs

(26)

with Feq(t) the equilibrium value of F(t) and peq(τs) given by eq 9. Equation 26 represents the asymptotic behavior of systems which have reached equilibrium by aging for very long times. The KWW parameters τK and τc corresponding to Feq(t) are plotted in Figure 14. In Figure 14, we consider systems at

e−t / τs

{τs > τf }

(23)

where {τs > τf} is the set of subunits with relaxation time larger than τf and N their number in the system. Note that the longtime behavior, in which we are mostly interested, does not depend on the cutoff at short times. The cutoff τpc at long times has been defined by eq 4. The function F(t) can be fitted with the Kohlrausch function (see e.g. ref 81):

F(t ) ≈ e−(t / τK)

Figure 14. Values of τK and τc obtained by fitting the function F(t) with the KWW function, at the theoretical equilibrium obtained after very long aging. τK and τc are compared to τα and τ.̅ Note that with the symplifying assumption that the width and shape of the distribution of relaxation time are constant (Δ is kept constant and equal to 4) the parameter β remains constant. τK is close to the central value τ̅ of the distribution of relaxation times, while τc is roughly equal to τα. The large difference between τc and τK is due to the particularly small value of β ≈ 0.17.

β

(24)

where the KWW parameters β and τK measure the width of the relaxation time distribution and the average time scale, respectively. The average correlation time τc is defined as τc =

∫0

equilibrium, which can be obtained only after a very long aging time or may also be calculated directly since we know (within the model) the equilibrium distribution of relaxation times. In these systems at equilibrium, the distribution of relaxation times is much broader than that of a system during aging or after deformation. It follows that the fitted value for β = 0.17 is smaller than the values typically obtained in the simulations, in which aging systems, with distributions of relaxation times narrower than at equilibrium, are considered. According to eq 25, this value β = 0.17 corresponds to τc/τK ≈ Γ(β−1)/β ≈ 586hence the large difference between τc and τK observed in Figure 14. Note that τK is comprised between the long time scale represented by τα and the intermediate one, represented by the central value τ.̅ τc, in turn, is very close to τα, which is consistent with the fact that τc is a measure of the mechanical time scale of the system. Homogenization of the Dynamics during Deformation. As seen above, the distribution of relaxation times broadens under deformation in the absence of facilitation mechanism. The evolution of β as a function of strain is plotted in Figure 15a for systems where no facilitation mechanism is considered. β decreases monotonically over almost the full range of deformation. Relaxation of subunits with intermediate relaxation times generates two distinct populations, with fast and slow motions, and the overall distribution of relaxation times broadens. The corresponding results are very different to those observed experimentally by Ediger and co-workers.11−14,16

−1



F(t ) dt ≈ τK

Γ(β ) β

(25)

The β-dependent factor Γ(β−1)/β is equal to one for β = 1 and increases when lowering β, e.g., up to 120 for β = 0.2. We will see later that 0.2 < β < 1 in the simulations, so that τc and τK may differ by up to 2 orders of magnitude. A representative example of the relaxation function F(t) is shown in Figure 13.

Figure 13. A representative example of the relaxation function F(t) as defined by eq 23 (blue crosses), obtained in a system which has been deformed at 355 K at a strain rate 0.03 s−1 up to 10% strain, fitted by a stretched exponential (red curve). Corresponding parameters are β = 0.45 and τK = 0.09 s, which gives τc = 0.22 s. K

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

Article

Macromolecules

Figure 15. Evolution of the stretching exponent β as a function of strain, for various strain rates as indicated. All systems have aged during 105 s prior to deformation. (a): facilitation mechanism not considered. β decreases monotonically from 0.58 down to 0.15. In (b) (T = 355 K), the facilitation mechanism has been implemented. In (b) β first decreases at small deformation, then, after a narrow plateau, it increases at deformation 3 corresponding to yield: the distribution of relaxation time gets narrower. The parameters used here are qc = 0:3 and N2/3 c = 10 .

increase takes place at 5% for the smallest strain rate value, 10−2 s−1. However, the strain rates of our simulations (10−2 s−1 and above) are at least 3 orders of magnitude larger than in experiments (10−5 or 10−4 s−1). It is observed that the larger is the strain rate, the larger is the deformation amplitude at which the β increase takes place. We thus think that both simulation and experiment results are compatible. More systematic comparisons between the results of our model with simulation performed at lower strain rates will help determining the best value for Nc2/3. Note, however, that eq 5 which relates the lifetime of slow subunits to the internal relaxation process of fast neighboring ones tends to underestimate the scale factor between τf and τd. Indeed, as the mobility of the fast subunits diffuses toward the interior of slow subunits, fast subunits simultaneously slow down. This slowing down effect has not been taken into account for establishing eq 5. The fact that a larger Nc value has to be considered for the facilitation parameter Nc2/3 than in the size ξ is a measure of this slowing down process which may be taken into account in future extensions of the PFVD model. Our point of view is that experiments such as those performed by Ediger et al. provide thus a quantitative insight into this facilitation mechanism and its intrinsic kinetics, which is slower than simply described by eq 5. However, at this early stage of the development of the theory, assuming a value 1000 for the facilitation mechanism is consistent with the limitation of the model. The curves in Figures 15b show a small overshoot of order 0.5−0.6 before reaching a steady state, in agreement with experimental data by Ediger et al.11,12 Then, for strain values beyond 8−13%, β remains nearly constant, with values between 0.4 and 0.6 for the slowest deformation rates, again comparable with the results by Ediger et al. within experimental error. For strain rates above 1 s−1, which are beyond the experimentally studied range, the overshoot disappears and β reaches values close to unity. Note that a slight decrease of β is observed in Figure 15b at very small strain. This decrease is not measured experimentally.11,12 Our interpretation is that this behavior is an artifact of the simulation procedure, not of the physical model. Indeed, before applying a deformation, we prepare our samples by letting them age for e.g. 104 or 105 s. During aging, we calculate the evolution of the distributions of relaxation times shown in Figure 6. To accelerate these time-consuming simulations, a cutoff on short time scales, of order 10−3−10−2 s, is applied. The smaller the cutoff, the longer the simulations. Because of

As discussed in ref 48, subunits with intermediate relaxation times relax first, and the amount of fast subunits increases under stress. This is due to the fact that very slow subunits, which are rare, are less induced to relax as compared to their environment: the free energy barrier for making them relax under stress is higher than that of intermediate subunits which surround them and relax first. When the latter have relaxed, very slow subunits are no longer submitted to stress or are submitted to a lower stress value. This is the reason why only intermediate time scales are modified at small deformation amplitudes. This relaxation mechanism is intimately related to the stress softening regime. As subunits with intermediate relaxation times relax, the stress supported by subunits with long relaxation times relaxes as well. This cascade of relaxation processes lead to stress softening. This is also the reason why a long time tail of the distribution still survive when no facilitation mechanism is taken into account, as observed in Figures 11a and 11b at 5% and 10% deformation. This long time tail is the consequence of thepartialrelaxation of the contribution of very slow subunits to the stress, also associated with the stress softening regime. A different behavior is observed in the presence of the facilitation mechanism, as shown in Figure 15b. Upon deformation, β remains constant (apart from a slight initial decrease that will be discussed below) up to a strain value between 5% and 8%, depending on the strain rate. During this plateau regime, the number of fast subunits has not sufficiently increased for melting very slow subunits. Slightly beyond yield, β increases steeply: very slow subunits are starting to melt as a consequence of the increasing number of fast regions. The distribution of relaxation times thus becomes narrower and the overall dynamics more homogeneous. The value of the adjustable parameter Nc2/3 has been set in order to qualitatively reproduce the experiments by Ediger and co-workers.11−13,16 Depending on the particular choice for the “facilitation” parameter Nc2/3, β in Figure 15b increases either too early if Nc2/3 is too small or too late if it is too large. Assuming that the scale of dynamical heterogeneities is 5 nm, Nc is about 1000 and Nc2/3 about 100. If this value is taken in the simulations, the increase of the parameter β appears at an early stage during deformation. The results of Ediger et al. can be recovered by choosing a “facilitation” parameter between 500 and 1000. In Figure 4 of ref 13, β stays constant up to 2% strain and has increased from 0.31 to 0.56 at 4.5% strain. In our paper, the L

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

Article

Macromolecules

Figure 16. τc as a function of strain, for various strain (compression) rates, at two different temperatures: (a) 350 K; (b) 355 K. The facilitation mechanism is implemented in the simulations. τc reaches a constant value in the postyield regime, after dropping by several orders of magnitude at deformations close to yield. The correlation time τc is obtained by measuring the reorientation time of optical probes in the experiments by Ediger and co-workers.13

can be made. A more direct comparison between simulations and experiments performed in extension may be done in future works. Mobility during Deformation. Let us consider how the mobility changes during the compression simulations described in the previous section. The evolution of τc (see eq 25) as a function of strain is plotted in Figure 16 for various strain rates and at two different temperatures. τc decreases by several orders of magnitude as the strain increases up to 6−10%, i.e., slightly above the yield strain. It then reaches a steady value, as was also the case for β. Note that at small deformations τc drops by almost 1 order of magnitude, before increasing and decreasing again. Lee et al.11 mention a fast τc drop immediately after applying deformation. Figure 17 shows the evolution of τc, measured at 15% strain, which corresponds to the constant value attained in the

this cutoff, the contribution of short time scales in the relaxation time distributions computed after aging, as plotted in Figure 6, is underestimated. Subsequent deformations are applied on shorter time scales, a few hundred seconds at most. During deformation, a cutoff on short time scales of order 10−4 s is applied. As a consequence, during the very early stage of the deformation, the distribution broadens at short time scales as a consequence of the relaxation of subunits due to the applied stress. It follows that β first decreases, before reaching a plateau at about 2% deformation. The β value on this transient plateau should thus be considered to be more representative of our model than the very initial β value. This initial, artifical decrease of the parameter β could be suppressed by choosing a smaller cutoff at short time during aging, making aging simulations very long, however. Note that we performed our simulations in compression (and also in shear deformation), whereas experiments by Ediger et al. that we consider for discussing the evolution of the width of the relaxation time distribution have been performed in extension. First of all, the results that we obtained in shear are very similar to those obtained in compression. At the semiquantitative level of description presented here, the same effects or trends are observed, even though the level of stress in shear may be smaller than in compression. Note also that we consider simulations at constant strain rates because most experiments in the literature regarding plastic flow in polymers have been performed this way. Our aim is first to describe plastic properties of glassy polymers by comparing with such mechanical experiments. On the other hand, experiments by Ediger et al. provide a crucial test for the validity of a microscopic model. These experiments are creep experiments, in which the applied force is constant. Conversely, in experiments performed at constant strain rates, the stress varies even beyond yield because of the stress softening regime. However, the stress stays within values of the same order of magnitude. We thus expect our results to qualitatively remain similar to those obtained at constant force. Note also that the stress during some of experiments by Ediger and co-workers varies also because they control the force. Thus, comparisons can only be made at a qualitative or semiquantitative level. Nevertheless, we believe they are crucial for elaborating a model primarily aimed, at this stage, for describing experiments regarding plastic flow of glassy polymers performed at constant strain rates. Our point of view is that both the experimental results of Ediger et al. and our simulations results are qualitatively robust so that a qualitative comparison at least

Figure 17. τc as a function of strain rate at 15% strain amplitude. The results can be fitted by Log τc = −2.0361 − 0.8567 Log γ̇ (at 350 K) and Log τc = −2.2129 − 0.9218 Log γ̇ (at 355 K).

postyield regime, as a function of γ̇. A nearly linear relationship between Log τc and Log γ̇ is observed. Fitting the data yields the relations Log τc−1 = 0.8567 Log γ̇ + 2.0361 at 350 K and Log τc−1 = 0. 9218 Log γ̇ + 2.2129 at 355 K. At both temperatures the slope is very close to one. Experiments,11 theory,44 and simulations12 have given similar results, with slopes comprised between 0.86 and 1. Thus, the average mobility τc−1 is found to be (nearly) linear as a function of the strain rate. Moreover, extrapolating the data to strain rate values γ̇ = 10−4−10−5 s−1 (which is out of our simulation range) gives τc ≈ 102−103 s, which is of the same order as the values obtained by Bending et al.13 M

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

Article

Macromolecules

Figure 18. τK as a function of strain rate for two different temperatures and various strain rates between 10−2 and 3 s−1. The systems have been previously aged during 105 s.

Figure 19. G′ versus temperature measured in compression (right) or in shear (left), with (blue curves) or without (black curves) the facilitation mechanism.

VII. REORGANIZATIONS UNDER STRAIN AND SHEAR BANDING Nonaffine reorganizations play a key role in yield mechanisms. In ref 48 the authors have shown that in the absence of the facilitation mechanism the samples undergo shear banding during deformation below the glass transition temperature. The facilitation mechanism tends to make the dynamics more homogeneous as compared to the model without facilitation mechanism. In principle, it might thus be possible that shear bands would be less present when taking this mechanism into account. This is not the case, however: shear bands similar to those described in ref 48 are observed here as well.

It must be noted that for both experiments and theory the exponent which relates τc to the strain rate, τc ∝ γ̇−ν, is slightly smaller than 1 in absolute value. From our point of view, this difference results from the competition between the applied strain, which tends to impose its own time scale to the sample dynamics, which would lead to an exponent ν = 1, and the ongoing aging process. The evolution of τK as a function of strain is plotted in Figure 18, at the same temperatures as in Figure 17. Prior to deformation, at both temperatures, τK is smaller than the equilibrium value given in Figure 14. This is a consequence of the finite aging time tw = 105 s. A characteristic mechanical relaxation time τm can be defined from the yield stress by the following relation:

τm =

VIII. CONCLUSIONS We have proposed a 3D model to describe relaxation mechanisms in glassy polymers under strain. In a previous article,48 we had proposed that the elastic energy stored at the scale of dynamical heterogeneities (∼5 nm) lowers the local free energy barriers for α-relaxation. In this paper, we have taken another relaxation process, which was proposed in refs 45 and 52 and discussed in details in ref 50, into account. This socalled facilitation mechanism induces relaxation of slow subunits (dynamical heterogeneities) by mobility diffusion from the faster environment. We have considered in details the way in which the various relevant time scales and the width of the relaxation time distribution depend on temperature, strain, and strain rate. We establish a connection between dynamical heterogeneities that is not present in earlier works apart through the so-called Eyring activation volume v, which has no clear interpretation as stressed by Chen and Schweizer.8,9 We proposed in ref 48 that the coupling between stress and dynamics is the consequence

σy G0γ ̇

(27)

Knowing the logarithmic dependence of the yield stress as a function of the strain rate, it follows ⎛ σ0 ⎞ Δy y Log γ ⎟⎟̇ τm = γ −̇ 1⎜⎜ + G0 ⎝ G0 ⎠

(28)

This equation means that the quantity Log τm−1 is expected to depend slightly sublinearly with respect to the strain rate, as a consequence of the logarithmic correction. Fitting Log τm−1 with the data given in Figure 9 gives Log τm−1 ∼ 0.9591 Log γ̇ + 1.73 at 350 K and Log τm−1 ∼ 0.9281 Log γ̇ + 1.8. Indeed, these values are close to those found for τc. The time scale τc is thus representative of mechanical properties. N

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

Article

Macromolecules of the elastic energy stored on the scale ξ of dynamical heterogeneities. Note that the point of view is very different since v ≪ ξ3. The description we propose here for plastic deformation is consistent with earlier works where the size ξ is essential as well, regarding the linear viscoelastic modulus G′ and G″,49 the Stokes law violation,52 aging and rejuvenating kinetics,45,52 and Tg shifts in thin films,39,62 which are all consistent with experimental results if one assumes a length scale for ξ of order 3−5 nm.39,48 Therefore, we stress that our model explains the yield stress behavior semiquantitatively with parameter values perfectly consistent with other aspects of the glass transition and without any new adjustable parameters, such as a so-called activation volume.8,9 The cascade of relaxation events which occur upon increasing the strain (or stress) is as follows. The stress is initially supported by a network of slow subunits which store increasing elastic energy. At some point, some of these subunits, mostly those with intermediate relaxation times, relax because free energy barriers have been reduced. This onset of random relaxation events triggered by the applied strain leads to the stress-softening regime, already observed in ref 48. At a further step, rare, very slow subunits undergo a reduction of the local stress as faster surrounding subunits have relaxed. As a consequence, their relaxation time is not further reduced by the applied stress: they could survive with long relaxation times. This is especially marked because the dependence of the free energy barrier reduction is a quadratic function of the stress. When the facilitation mechanism is not taken into account, these mechanisms result in a constant decrease of the Kohlrausch exponent β (as shown in Figure 15a), which is related to an overall broadening of the relaxation time distribution, as a consequence of the survival of a fraction of slow subunits (as shown in Figure 11). This constant decrease of β is not observed experimentally. Instead, β tends to increase beyond a typical strain of a few percent.11−14 When the facilitation mechanism is implemented in the model, the relaxation time distribution is more effectively shifted toward shorter times and no slow region survives (as shown in Figure 12). It results that the Kohlrausch exponent β increases more or less suddenly at 5% deformation following a plateau regime, in qualitative agreement with experimental observations. The facilitation mechanism is intimately related to the scale ξ of dynamical heterogeneities, which is set by the competition between two relaxation processes: internal reallocation of free volume within a subunit or mobility diffusion from fast to slow subunits.52 The scale ξ plays a double key role in the model described here: (1) for calculating the decrease of free energy barriers and (2) for calculating the dynamic coupling between neighboring subunits. Thus, the evolution of β measured experimentally allows adjusting the scale parameter Nc2/3 which controls the efficiency of the facilitation mechanism in the simulations. The value Nc2/3 = 1000 chosen here leads to an evolution of β during applied strain which is comparable to experimental results.11−14 For Nc2/3 ∼ 100, β would increase too early as compared to experiments, while for Nc2/3 ∼ 10 000, the facilitation mechanism would be less efficient than observed in experiments. Thus, the scale factor Nc2/3 can be adjusted quite precisely by comparison to experimental observations. The value Nc2/3 = 1000 found to reproduce at best the experimental data of Ediger and co-workers may seem quite large as compared to the size ξ of dynamical heterogeneities.

However, it is important to emphasize that this Nc2/3 value is consistent with the size ξ. Indeed, during mobility diffusion, fast subunits undergo a slowing down of their own dynamics as they melt slow neighboring subunits. This slowing down effect has not been taken into account. As a consequence, the facilitation process may be expected to be less efficient than described by the simple relation between τf and τd given by eq 5.52 Thus, the scale of dynamical heterogeneities in the model is absolutely not arbitrary. It is at the core of the relaxation processes introduced in ref 48 and in this paper. The model proposed here takes into account aging and rejuvenating dynamics of subunits and allows for calculating the evolution of the relaxation time distribution during a complex mechanical history. It extends to 3D the description of evolution mechanisms proposed in ref 45 in the context of a Fokker−Planck equation, without explicit 3D description. Equations 11−13 may be used for calculating the evolution of the local relaxation time of subunits upon quenching or heating, in the same way as for mechanical loading. Checking whether thermal histories of glassy polymers studied experimentally by Kovacs can be quantitatively reproduced will also be a key test of this 3D model. The facilitation mechanism is key when considering the rejuvenation kinetics of glassy polymers after an upward temperature jump.45,50 In particular, it explains that melting can be much faster than the initial αrelaxation time of the sample: after heating, the kinetics of fast subunits is made even faster and in turn accelerates the dynamics of slow subunits by mobility diffusion at the scale ξ. Without this dynamic coupling between neighboring subunits, slow regions with local relaxation times larger than the initial τα would survive on a time scale comparable to the initial τα. Studying complex thermomechanical histories based on the present model will be done in future works. Since the facilitation mechanism tends to homogenize the response of the material, an issue was to know whether shear bands were still present, as it was the case in simulations discussed in ref 48. This is still the case. Shear bands are indeed observed on the scale of a few tens of nanometers. This issue will be considered in a forthcoming paper, in which mechanisms leading to strain hardening will be described.78,80,82 Such mechanisms require further developments, as they are not accounted for by the model as presented here.



APPENDIX. THE ELASTIC MODULUS The elastic moduli as obtained from compression and shear simulations, run with (blue curve) or without (black curve) the facilitation mechanism, are plotted in Figure 19. The observed features in compression and shear are nearly identical, except for smaller values observed for the shear modulus, as expected. Close to the glass transition temperature Tg ≈ 370 K, G′ drops to small values, close to G∞. Some tens of kelvin below Tg, the elastic modulus is nearly constant and equal to a few GPa. Note that the transition here is more abrupt than observed in ref 39 for PMMA. This is due to the fact that τα(T) increases much more rapidly for PS than for PMMA when lowering temperature, according to the WLF law. Consequently, the glass transition in PS takes place on a temperature range much narrower in PS than in PMMA. These curves are qualitatively in agreement with data reported in the literature.80 At higher temperatures, G′(T) reaches a plateau, as a consequence of the permanent rubbery network. O

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

Article

Macromolecules



spectroscopy of a stretched polymer glass: heterogeneous dynamics and plasticity. Macromolecules 2016, 49, 3889−3898. (21) Ediger, M. D. Spatially heterogeneous dynamics in supercooled liquids. Annu. Rev. Phys. Chem. 2000, 51, 99−128. (22) Richert, R. Heterogeneous dynamics in liquids: fluctuations in space and time. J. Phys.: Condens. Matter 2002, 14, R703−R738. (23) Schmidt-Rohr, K.; Spiess, H. W. Nature of Nonexponential Loss of Correlation above the Glass-Transition Investigated by Multidimensional NMR. Phys. Rev. Lett. 1991, 66, 3020−3023. (24) Tracht, U.; Wilhelm, M.; Heuer, A.; Feng, H.; Schmidt-Rohr, K.; Spiess, H. W. Length scale of dynamic heterogeneities at the glass transition determined by multidimensional nuclear magnetic resonance. Phys. Rev. Lett. 1998, 81, 2727−2730. (25) Reinsberg, S. A.; Qiu, X. H.; Wilhelm, M.; Spiess, H. W.; Ediger, M. D. Length scale of dynamic heterogeneity in supercooled glycerol near Tg. J. Chem. Phys. 2001, 114, 7299−7302. (26) Cicerone, M. T.; Blackburn, F. R.; Ediger, M. D. Anomalous Diffusion of Probe Molecules in Polystyrene Evidence for Spatially Heterogeneous Segmental Dynamics. Macromolecules 1995, 28, 8224− 8232. (27) Wang, C.-Y.; Ediger, M. D. Enhanced translational diffusion of 9,10-bis(phenylethynyl)anthracene (BPEA) in polystyrene. Macromolecules 1997, 30, 4770−4771. (28) Cicerone, M. T.; Wagner, P. A.; Ediger, M. D. Translational diffusion on heterogeneous lattices: A model for dynamics in glass forming materials. J. Phys. Chem. B 1997, 101, 8727−8734. (29) Fujara, F.; Geil, B.; Sillescu, H.; Fleischer, G. Translational and Rotational Diffusion in Supercooled Orthoterphenyl Close to the Glass Transition. Z. Phys. B: Condens. Matter 1992, 88, 195−204. (30) Hwang, Y.; Inoue, T.; Wagner, P. A.; Ediger, M. D. Molecular motion during physical aging in polystyrene: Investigation using probe reorientation. J. Polym. Sci., Part B: Polym. Phys. 2000, 38, 68−79. (31) Schiener, B.; Böhmer, R.; Loidl, A.; Chamberlin, R. V. Nonresonant spectral hole burning in the slow dielectric response of supercooled liquids. Science 1996, 274, 752−754. (32) Richert, R. Triplet state solvation dynamics: Basics and applications. J. Chem. Phys. 2000, 113, 8404−8229. (33) Yoshimoto, K.; Jain, T. S.; van Workum, K.; Nealey, P. F.; de Pablo, J. J. Mechanical heterogeneities in model polymer glasses at small length scales. Phys. Rev. Lett. 2004, 93, 175501. (34) Leonforte, F.; Boissière, R.; Tanguy, A.; Wittmer, J. P.; Barrat, J.L. Continuum limit of amorphous elastic bodies. III. Threedimensional systems. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 224206. (35) Riggleman, R. A.; Lee, H.-N.; Ediger, M. D.; de Pablo, J. J. Free volume and finite-size effects in a polymer glass under stress. Phys. Rev. Lett. 2007, 99, 215501. (36) Tsamados, M.; Tanguy, A.; Goldenberg, C.; Barrat, J.-L. Local elasticity map and plasticity in a model Lennard-Jones glass. Phys. Rev. E 2009, 80, 026112. (37) Papakonstantopoulos, G. J.; Riggleman, R. A.; Barrat, J.-L.; de Pablo, J. J. Molecular plasticity of polymeric glasses in the elastic regime. Phys. Rev. E 2008, 77, 041502. (38) Riggleman, R. A.; Schweizer, K. S.; de Pablo, J. J. Nonlinear creep in a polymer glass. Macromolecules 2008, 41, 4969−4977. (39) Dequidt, A.; Long, D. R.; Sotta, P.; Sanseau, O. Mechanical properties of thin confined polymer films close to the glass transition in the linear regime of deformation: Theory and simulations. Eur. Phys. J. E: Soft Matter Biol. Phys. 2012, 35, 61. (40) Schweizer, K. S.; Saltzman, E. J. Entropic barriers, activated hopping, and the glass transition in colloidal suspensions. J. Chem. Phys. 2003, 119, 1181−1196. (41) Schweizer, K. S.; Saltzman, E. J. Theory of dynamic barriers, activated hopping, and the glass transition in polymer melts. J. Chem. Phys. 2004, 121, 1984−2000. (42) Chen, K.; Schweizer, K. S. Suppressed Segmental Relaxation as the Origin of Strain Hardening in Polymer Glasses. Phys. Rev. Lett. 2009, 102, 038301.

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected] (D.R.L.). ORCID

Alain Dequidt: 0000-0003-1206-1911 Paul Sotta: 0000-0002-4378-0858 Didier R. Long: 0000-0002-3013-6852 Notes

The authors declare no competing financial interest.



REFERENCES

(1) Ferry, J. D. Viscoelastic Properties of Polymers; John Wiley and Sons, Inc.: 1980. (2) Ediger, M. D.; Angell, C. A.; Nagel, S. R. Supercooled liquids and glasses. J. Phys. Chem. 1996, 100, 13200−13212. (3) Nielsen, L. E.; Landel, R. F. Mechanical Properties of Polymers and Composites; Marcel Dekker: New York, 1994. (4) Meijer, H. E. H.; Govaert, L. Mechanical performance of polymer systems: The relation between structure and properties. Prog. Polym. Sci. 2005, 30, 915−938. (5) van Melick, H. G. H.; Govaert, L. E.; Meijer, H. E. H. Polymer 2003, 44, 2493−2502. (6) Haward, R. N. The Physics of Glassy Polymers; Applied Science Publishers: London, 1973. (7) Eyring, H. Viscosity, plasticity, and diffusion as examples of absolute reaction rates. J. Chem. Phys. 1936, 4, 283−291. (8) Chen, K.; Schweizer, K. S. Microscopic constitutive equation theory for the nonlinear mechanical response of polymer glasses. Macromolecules 2008, 41, 5908−5918. (9) Chen, K.; Schweizer, K. S. Stress-enhanced mobility and dynamic yielding in polymer glasses. Europhys. Lett. 2007, 79, 26006. (10) Loo, L. S.; Cohen, R. E.; Gleason, K. K. Chain mobility in the amorphous region of nylon 6 observed under active uniaxial deformation. Science 2000, 288, 116−119. (11) Lee, H.-N.; Riggleman, R. A.; de Pablo, J. J.; Ediger, M. D. Deformation-Induced Mobility in Polymer Glasses during Multistep Creep Experiments and Simulations. Macromolecules 2009, 42, 4328− 4336. (12) Riggleman, R. A.; Lee, H.-N.; Ediger, M. D.; de Pablo, J. J. Heterogeneous dynamics during deformation of a polymer glass. Soft Matter 2010, 6, 287−291. (13) Bending, B.; Christison, K.; Ricci, J.; Ediger, M. D. Measurement of Segmental Mobility during Constant Strain Rate Deformation of a Poly(methyl methacrylate) Glass. Macromolecules 2014, 47, 800−806. (14) Lee, H.-N.; Paeng, K.; Swallen, S. F.; Ediger, M. D. Direct Measurement of Molecular Mobility in Actively Deformed Polymer Glasses. Science 2009, 323, 231−234. (15) Lee, H.-N.; Paeng, K.; Swallen, S. F.; Ediger, M. D.; Stamm, R. A.; Medvedev, G. A.; Caruthers, J. M. Molecular Mobility of Poly(methyl methacrylate) Glass During Uniaxial Tensile Creep Deformation. J. Polym. Sci., Part B: Polym. Phys. 2009, 47, 1713−1727. (16) Hebert, K.; Bending, B.; Ricci, J.; Ediger, M. D. Effect of Temperature on Postyield Segmental Dynamics of Poly(methyl methacrylate) Glasses: Thermally Activated Transitions Are Important. Macromolecules 2015, 48, 6736−6744. (17) Hebert, K.; Ediger, M. D. Reversing Strain Deformation Probes Mechanisms for Enhanced Segmental Mobility of Polymer Glasses. Macromolecules 2017, 50, 1016−1026. (18) Bending, B.; Ediger, M. D. Molecular Mobility of Poly(methyl methacrylate) Glass During Uniaxial Tensile Creep Deformation. J. Polym. Sci., Part B: Polym. Phys. 2016, 54, 1957−1967. (19) Kalfus, J.; Detwiler, A.; Lesser, A. J. Probing Segmental Dynamics of Polymer Glasses during Tensile Deformation with Dielectric Spectroscopy. Macromolecules 2012, 45, 4839−4847. (20) Pérez-Aparicio, R.; Cottinet, D.; Crauste-Thibierge, C.; Vanel, L.; Sotta, P.; Delannoy, J.-Y.; Ciliberto, S.; Long, D. R. Dielectric P

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

Article

Macromolecules (43) Chen, K.; Saltzman, E. J.; Schweizer, K. S. Segmental dynamics in polymers: from cold melts to ageing and stressed glasses. J. Phys.: Condens. Matter 2009, 21, 503101. (44) Chen, K.; Schweizer, K. S. Theory of Yielding, Strain Softening, and Steady Plastic Flow in Polymer Glasses under Constant Strain Rate Deformation. Macromolecules 2011, 44, 3988−4000. (45) Merabia, S.; Long, D. Heterogeneous dynamics, aging and rejuvenating in van der Waals liquids. J. Chem. Phys. 2006, 125, 234901. (46) Medvedev, G. A.; Caruthers, J. M. Development of a stochastic constitutive model for prediction of postyield softening in glassy polymers. J. Rheol. 2013, 57, 949−1002. (47) Medvedev, G. A.; Caruthers, J. M. Stochastic model prediction of nonlinear creep in glassy polymers. Polymer 2015, 74, 235−253. (48) Dequidt, A.; Conca, L.; Delannoy, J. Y.; Sotta, P.; Lequeux, F.; Long, D. R. Heterogeneous Dynamics and Polymer Plasticity. Macromolecules 2016, 49, 9148−9162. (49) Masurel, R.; Cantournet, S.; Dequidt, A.; Long, D.; Montes, H.; Lequeux, F. Linear Viscoelasticity of polymer in their glass transition domains: 2d finite elements simulations. Macromolecules 2015, 48, 6690−6702. (50) Polymer Glasses; Roth, C., Ed.; Taylor and Francis: Boca Raton, FL, 2016; Chapter Mechanical properties of polymers and nanocomposites close to the glass transition by Dequidt et al. (51) Long, D.; Lequeux, F. Heterogeneous dynamics at the glass transition in van der Waals liquids, in the bulk and in thin films. Eur. Phys. J. E: Soft Matter Biol. Phys. 2001, 4, 371−387. (52) Merabia, S.; Long, D. Heterogeneous dynamics at the glass transition in van der Waals liquids: determination of the characteristic scale. Eur. Phys. J. E: Soft Matter Biol. Phys. 2002, 9, 195−207. (53) Prigogine, I. The Molecular Theory of Solutions; North-Holland Publishing Co.: Amsterdam, 1957. (54) Masnada, E. M.; Julien, G.; Long, D. R. Miscibility Maps for Polymer Blends: Effects of Temperature, Pressure, and Molecular Weight. J. Polym. Sci., Part B: Polym. Phys. 2014, 52, 419−443. (55) Adam, G.; Gibbs, J. H. On temperature dependence of cooperative relaxation properties in glass-forming liquids. J. Chem. Phys. 1965, 43, 139. (56) Ritort, F.; Sollich, P. Glassy dynamics of kinetically constrained models. Adv. Phys. 2003, 52, 219−342. (57) Jung, Y.-J.; Garrahan, J. P.; Chandler, D. Excitation lines and the breakdown of Stokes-Einstein relations in supercooled liquids. Phys. Rev. E 2004, 69, 061205. (58) Garrahan, J. P.; Chandler, D. Geometrical explanation and scaling of dynamical heterogeneities in glass forming systems. Phys. Rev. Lett. 2002, 89, 035704. (59) Chen, K.; Saltzman, E. J.; Schweizer, K. S. Segmental dynamics in polymers: from cold melts to ageing and stressed glasses. J. Phys.: Condens. Matter 2009, 21, 503101. (60) Tito, N. B.; Lipson, J. E. G.; Milner, S. T. Lattice model of mobility at interfaces: free surfaces, substrates, and bilayers. Soft Matter 2013, 9, 3173−3180. (61) Lipson, J. E. G.; Milner, S. T. Percolation model of interfacial effects in polymeric glasses. Eur. Phys. J. B 2009, 72, 133−137. (62) Masurel, R. J.; Gelineau, P.; Cantournet, S.; Dequidt, A.; Long, D. R.; Lequeux, F.; Montes, H. Role of Dynamical Heterogeneities on the Mechanical Response of Confined Polymer. Phys. Rev. Lett. 2017, 118, 047801. (63) Wagner, H.; Richert, R. Dielectric relaxation of the electric field in poly(vinyl acetate): a time domain study in the range 10−3−106 s. Polymer 1997, 38, 255−261. (64) Kremer, F.; Schonhals, A. Broadband Dielectric Spectroscopy; Springer-Verlag: Berlin, 2003. (65) Havriliak, S.; Havriliak, S. J. Result from an unbiased analysis of nearly 1000 s ets of relaxation data. J. Non-Cryst. Solids 1994, 172-174, 297−310. (66) Sotta, P.; Long, D. The cross-over from 2D to 3D percolation. Theory and numerical simulations. Eur. Phys. J. E: Soft Matter Biol. Phys. 2003, 11, 375−388.

(67) Merabia, S.; Sotta, P.; Long, D. Heterogeneous nature of the dynamics and glass transition in thin polymer films. Eur. Phys. J. E: Soft Matter Biol. Phys. 2004, 15, 189−210. (68) Lebowitz, J. L.; Presutti, E.; Spohn, H. Microscopic Models of Hydrodynamic Behavior. J. Stat. Phys. 1988, 51, 841−862. (69) Kanninen, M. F.; Popelar, C. H. Advanced Fracture Mechanics; Oxford University Press: New York, 1985. (70) Holt, D. L. Modulus and yield stress of glassy poly(methyl methacrylate) at strain rates up to 103 in./inch/second. J. Appl. Polym. Sci. 1968, 12, 1653−1659. (71) Ward, I. M. Review - Yield behaviour of polymers. J. Mater. Sci. 1971, 6, 1397−1417. (72) Bauwens-Crowet, C.; Bauwens, J. C.; Homes, G. Tensile Yield Stress Behaviour of Glassy Polymers. Journal of Polymer Science: Part A2 1969, 7, 735−742. (73) Katz, S.; Lebowitz, J. L.; Spohn, H. Phase transitions in stationary nonequilibrium states of model lattice systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1983, 28, 1655−1658. (74) Struik, L. C. E. Physical Aging in Amorphous Polymers and Other Materials; Elsevier: Amsterdam, 1978. (75) Long, D.; Sotta, P. Numerical Simulation for the Mesoscale Deformation of Disoredered and Reinforced Elastomers; IMA Volume in Mathematics and its Applications: Modeling of Soft Matter; Calderer, M.-C. T., Terentjev, E. M., Eds.; 2005; Vol. 141, pp 205−234. (76) Long, D.; Sotta, P. Non-linear and Plastic Behavior of Soft Thermoplastic and Filled Elastomers Studied by Dissipative Particle Dynamics. Macromolecules 2006, 39, 6282−6297. (77) Merabia, S.; Sotta, P.; Long, D. R. A Microscopic Model for the Reinforcement and Non-linear Behavior of Filled Elastomers and Thermoplastic Elastomers (Payne and the Mullins effects). Macromolecules 2008, 41, 8252−8266. (78) Klompen, E. T. J.; Engels, T. A. P.; Govaert, L. E.; Meijer, H. E. H. Modeling of the postyield response of glassy polymers: Influence of thermomechanical history. Macromolecules 2005, 38, 6997−7008. (79) Bauwens-Crowet, C. Compression yield behavior of polymethyl methacrylate over a wide-range of temperatures and strain-rates. J. Mater. Sci. 1973, 8, 968−979. (80) van Melick, H. G. H.; Govaert, L. E.; Meijer, H. E. H. Localisation phenomena in glassy polymers: influence of thermal and mechanical history. Polymer 2003, 44, 3579−3591. (81) Rault, J. Glass: Kohlrausch exponent, fragility, anharmonicity. Eur. Phys. J. E: Soft Matter Biol. Phys. 2012, 35, 9703. (82) Wendlandt, M.; Tervoort, T. A.; Suter, U. W. Non-linear, ratedependent strain-hardening behavior of polymer glasses. Polymer 2005, 46, 11786−11797.

Q

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