Molecular Simulation of Viscous Dissipation due to Cyclic Deformation

Article pubs.acs.org/Macromolecules

Molecular Simulation of Viscous Dissipation due to Cyclic Deformation of a Silica−Silica Contact in Filled Rubber Jan Meyer,*,† Reinhard Hentschke,† Jonathan Hager,† Nils W. Hojdis,‡ and Hossein Ali Karimi-Varzaneh‡ †

School of Mathematics and Natural Sciences Bergische Universität D-42097 Wuppertal, Germany Continental Reifen Deutschland GmbH D-30419 Hannover, Germany

‡

ABSTRACT: Filled elastomers acquire their mechanical strength through fillers forming spanning branched networks throughout the rubber matrix. Here we study one reversibly breakable filler aggregate-to-filler aggregate contact within a network branch and its contribution to dissipative loss using molecular dynamics simulations. Our model system consists of a pair of spherical silica particles embedded in cis-1,4polyisoprene. Variable chemical options include silanes attached to the particle surfaces at varying density and surface distribution and cross-links between polymer chains as well as chemical bonding of the polymer chains to the silica particles via silanes. We validate key properties of the pure polymer, including density, thermal expansion, and characteristic ratio as well as the dependence of the polymer network density on sulfur content. On the basis of force-vs-inter particle separation curves for cyclic loading, we obtain the energy dissipation in the contact depending on the aforementioned chemical parameters and temperature. We track the segmental relaxation peak in our model system along the temperature axis while moving the filler particle relative to one another at different speeds. From this we are able to show that the peaks in the loss-vs-temperature diagram correspond to the experimental segment relaxation in pure polyisoprene at the attendant frequencies. This in turn allows to relate the high frequency results of our atomistic simulations to the technically relevant frequency range. Finally, we investigate the structure and mobility of the polymer in the vicinity of the particles.

■

INTRODUCTION Elastomer matrices are modified routinely by addition of inorganic nanoparticles, i.e., particles with linear dimensions between 1 and 100 nm.1−4 These nanocomposites are complex compounds containing numerous ingredients combined in an elaborate processing scheme in order to match a broad spectrum of product requirements. The attendant materials development commonly employs an experience-based trial and error approach. However, progress often is hampered by a poor understanding of the molecular mechanisms underlying a material’s macroscopic behavior. Theoretical methods of choice, providing direct access to molecular mechanisms, are computer simulations. Because the literature, even on nanocomposites, is already quite extensive we limit ourselves to examples representing the typical levels of approximation. Most studies employ molecular force fields, retaining chemical detail in a very small system, or develop coarse grained approaches applicable into the micrometerrange. Examples of the former type are ref 5, where the authors compute the interaction of rod-like nanoparticles in an explicit solvent, and ref 6, where the authors investigate the effect of isolated silica nanoparticles on the structure and dynamics of an embedding polymer matrix. In principle, a follow up study, but with considerably more detail, is the work by Ndoro et al.7 Again similarly Guseva et al.8 study non-cross-linked cis-1,4© XXXX American Chemical Society

polyisoprene melts between silica surfaces. These authors observe a strong dependence of the segment dynamics depending on the separation between the surfaces. Examples of the latter type, in the context of filled elastomers, are refs 9−14. Riggleman et al.9 find that the nanoparticles serve as entanglement attractors, particularly at large deformations, altering the topological constraint network that arises in the composite material. Jaber and co-workers10 investigate filler network formation in polymer nanocomposites under shear. Kalathi et al.11 use nonequilibrium molecular dynamics simulations to study the filler effect on the shear viscosity of a polymer melt depending on the filler-to-chain interaction. The authors of refs 12 and 13 study filler dispersion and obtain both the storage and the loss modulus as functions of strain amplitude from nonequilibrium shear experiments on a corse grained nanocomposite for different filler concentrations. Also using a coarse grained model, Masnada et al.14 describe filled entangled polymer melts and study the filler reinforcement on the entanglement network. The authors vary filler volume fraction and distribution (cubic lattice, randomly dispersed, and small clusters randomly dispersed) and consider the presence Received: May 8, 2017 Revised: July 18, 2017

A

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

Article

Macromolecules

Figure 1. Hierarchy of scales within a filled elastomer material. The macroscopic mechanical material properties are governed by the filler network embedded in the elastomer matrix. The basic structural and dynamic element studied in this work is a contact between adjacent filler particles consisting of aggregates of smaller, so-called primary filler particles including surrounding rubber. The highest resolution, depicting our specific model contact, shows two spherical silica particles embedded in 85 polyisoprene chains each consisting of 200 monomer units.

or absence of grafted chains on the fillers. Another class of models simplifies the polymer matrix to a mere network of disordered springs.15−18 The network is subjected to periodic shear and the model is solved by computing the force equilibrium, including viscous forces, as a function of strain amplitude at variable frequency. Again the approach yields the dynamic moduli of the model nanocomposite. This is quite similar to the dissipative dynamics approach. An example is ref 19, where Raos et al. present nonequilibrium dissipative particle dynamics simulations of cross-linked elastomers containing solid filler particles. The authors study the effects of the morphology and of the effective particle−particle interactions. A finite element-type approach to the micromechanical mechanism of reinforcement and loss in filled rubbers is discussed by Gusev.20 Here we focus on energy dissipation associated with the mechanically reversible opening and closing of contacts between filler particles dependent on the contact’s chemical structure. Contacts of the present type are likely to be mechanical key elements within filler network strands traversing the polymer matrix (cf. Figure 1). We follow up on our previous feasibility study in ref 21, where we introduced a method for measuring the dissipated energy due to the relative motion of naked silica particles embedded in a polymer melt in molecular dynamics simulations. Energy dissipation in filled elastomer materials is of major technical importance in a great number of applications (e.g., refs 1−4). However, to the best of our knowledge, it has not yet received much attention in attendant simulation studies on the molecular level. In the following, we investigate a particle−particle contact between silica particles embedded in a polymer matrix subject to cyclic strain. The model is prepared by immersion of two spherical silica particles in liquid isoprene. Subsequently the liquid is transformed into polymer chains by a fast polymerization algorithm based on previous work.21 In addition, sulfur cross-links are introduced between the chains, which can also be linked to a particle’s silanol groups via suitable silanes. The sulfur cross-link density as well as the surface coverage and surface distribution of the silanes are variable model parameters. A variable thermodynamic parameter is temperature. On the basis of force-vs-inter particle separation curves for cyclic loading of the contact during a molecular dynamics simulation we obtain its energy dissipation. We track the segmental relaxation peak in our model system along the temperature axis while moving the filler particles relative to one another at different speeds. From this, we are able to show that the peaks in the loss-vs-temperature diagram correspond to the experimental segment relaxation in pure polyisoprene at the attendant frequencies. This in turn allows us to relate the high frequency results of our atomistic

simulations to the technically relevant frequency range. Finally, we investigate the structure and mobility of the polymer in the vicinity of the particles. We expect our results to add to the elucidation of the socalled Payne effect (e.g. refs 22−25), which is the pronounced decrease of the storage modulus of a filled elastomer material during cyclic loading with increasing strain amplitude. This does not occur in unfilled elastomers and therefore the effect is directly related to the filler particle surfaces and most likely to the particle-to-particle contacts. The Payne effect manifests itself not only by the decrease of the storage modulus but affects the loss modulus, which is directly related to the dissipated energy, as well. In particular it governs the dissipative loss on the high temperature side of the tan δ-peak,26 the ratio of loss to storage modulus, which serves as an important laboratory material performance indicator (e.g., ref 27). Notice however that here we do not present a study of an entire filler network. Thus, we do not obtain the overall shape of the dynamic moduli in terms of strain or temperature. On the other hand, the mechanical reversible opening and closing of a single network contact, which is studied here, yields a model for the energy dissipation on the nanoscale. Subsequently this can be combined with a previous model of the Payne effect (ref 25) to yield a more complete description of the latter. The paper is structured as follows. We begin with an outline of our simulation methodology, including the polymerization scheme and tests thereof, particle modeling, silanization and physical cross-linking as well as the force gauge. Next we discuss our results for the dynamic energy dissipation in the model contact when the particles are subject to cyclic relative motion. This includes the influence of different parameter combinations, e.g., silanes chemically linked to the polymer or not, variable silane surface density, and distribution as well as sulfur cross-links between the polymer chains. In a section on the effects of variable temperature, we discuss the applicability and application of time−temperature-superposition. Finally, we study structural and dynamical effects of the embedded filler particles on the surrounding polymer matrix. The last section is the conclusion.

■

SIMULATION METHODOLOGY Setting up the Model System. The simulation setup follows in part ref 21. Spherical silica particles with fixed diameter, D = 4.2 nm, are cut from β-cristobalite. Subsequently, the particle surfaces are saturated with hydroxyl groups resulting in a silanol surface density of σ ≈ 6.4 nm−2. This number is in accord with experimental data on precipitated silica.28 The ratio of geminal to isolated hydroxyl groups is roughly 2:1. Notice that hydroxyl groups are removed in part B


Article

Macromolecules during silanization with bis[3-(triethoxysilyl)propyl] tetrasulfide (TESPT) (specifically, we attach [3-(triethoxysilyl)propyl] disulfide). In the following step we determine the equilibrium distance of the particles using the force gauge discussed below (see also ref 21). During the subsequent addition of isoprene and its polymerization, the particles are constrained to their initial positions. Figure 2a shows a simulation box containing monomers only, i.e., the system prior to the polymerization. One of the

Figure 2. Illustration of the polymerization process.

monomers is shown in red, which means that this is a catcher monomer. If another monomer is located within 0.8 nm of the catcher monomer the two can “react” and form a chemical bond. In the case of multiple possible monomers, the closest one is selected for reaction. The newly added monomer then becomes the next catcher monomer. Repeated application of this procedure leads to the formation of polymer chains (Figure 2b), whose length and number can be adjusted. Figure 3a illustrates the evolution of the system density during the polymerization process (in this specific example no silica particles are present). Following an initial equilibration of the monomer liquid, the polymerization starts at t = t1 and the density increases linearly. When the chains become long, the number of locally available monomers is reduced, which is indicated by a diminished increase of the density. At t = t2, the polymerization goal is achieved. The remaining free monomers are removed from the simulation box and a final equilibration relaxes the voids. As a check of the final result we calculate the density and its temperature dependence as shown in Figure 3b. The characteristic ratio for chains with length up to 800 monomers is within 10% of the experimental value of 4.7.29 We also note that all polymer configurations in this work are monodisperse. Introducing polydispersity adds another variable quantity to an already complex system, which we wanted to avoid in this study. Optional disulfide bridges between the chains are introduced following the polymerization (cf. below). In addition we link all polymer ends to monomers in neighboring chains to avoid the effects of dangling ends. We also note that there are at least 40 monomers separating cross-links of a chain with itself. This avoids inner-chain cross-links before a chain segment has a significant probability to intersect with itself. Overall the crosslinks are uniformly distributed in the system. This is accomplished by dividing the simulation box into identical subsystems (∼3 nm × 3 nm × 3 nm), each containing the same number of cross-links relative to its polymer volume fraction (notice that a subsystem may overlap with a silica particle). Figure 4 shows the effect of sulfur concentration on the polymer density. We note that the simulation results are in very good accord with the experimental values taken from ref 30. In addition to the polymer−polymer cross-links there are optional disulfide bridges coupling the silanes to the polymers. An

Figure 3. Panel a shows the polymer density evolution during the polymerization in a system without silica particles. Typical values for t1 and t2 are 0.1 and 1.0 ns, respectively. For the final systems (with particles included) these values increase up to 40 ns. ρpoly ∼ 910 kg m−3 is the polymer density which is in accordance with the −3 experimental value ρexp poly ∼ 910 kg m . Panel b shows the equilibrium bulk polymer density vs temperature. The experimental data are taken from ref 31. The two sets of simulation results are based on systems with 20 chains containing 100 monomers and 10 chains containing 200 monomers, respectively. Density values are calculated by averaging over 5 ns of simulation time.

Figure 4. Polymer density vs sulfur concentration in a cross-linked pure polymer system. Density values are calculated by averaging over 5 ns of simulation time. The experimental data are taken from ref 30. phr: parts (by weight) per hundred rubber.

illustration of the entire system including all components and options is depicted in Figure 5. C


Article

Macromolecules

Figure 6. Sketch of the two silica particles and the interactions between them. The dashed line defines the geometric contact area where the cone angle 2θ = 2 × 40°.

It is worth noting that the maximum separation (or amplitude) of the two particles forming the contact not only depends on the maximum displacement of the virtual particle. The spring constant, k, and the chemical composition, influencing the frequency dependent stiffness of the system, also affect the maximum interparticle displacement during a cycle. Because we move the virtual particle using a fixed number of constant displacement steps, the resulting amplitudes can vary. It is useful to to distinguish situations when the bare silica surfaces interact directly from those when the silanes on opposite particles intervene through their mutual interaction. Therefore, we introduce the “contact area”. The contact area of a particle is defined by the cone angle 2θ depicted in Figure 6. The cone angle in turn defines a radially symmetric surface area centered on the axis joining the two particles. The cone angle is chosen according to the separation, s(θ), between the two particles from the rim of one contact area to the rim of the other (parallel to the aforementioned axis). If s(θ) is larger than twice the length of the silane, then silanes grafted to a particle’s surface anywhere outside the contact area do not interact strongly, which here means via steric exclusion. In the following, if not stated otherwise, θ = 40°. Figure 7 is an example illustrating our force measurements. The red up-triangles are measured interparticle force values in

Figure 5. Chemical composition of our contact model.

During the system preparation we apply constant pressure (NPT), whereas the force measurements discussed below are carried out at constant volume (NVT). Long range Coulomb interactions are calculated using the particle mesh Ewald summation with an additional short distance modification (we use the Lennard-Jones parameters and (larger) charges in ref 21). The latter employs Coulomb integrals32 instead of point charges in order to avoid divergencies due to charges at close proximity. The polymer force field is taken from ref 33. Because of the inhomogeneity of our systems, introduced when particles are present, we do not include long-range dispersion corrections and constraints for the internal bonds. Thus, a 2% decrease of all σ values in the aforementioned reference is necessary. For the sulfur bridges we use the parameters of Amber94.34 These were checked by comparing results for small test systems, essentially cross-linked dimers, which were simulated using the force field as well as an ab initio approach.35 Measuring the Inter-Particle Force. We want to measure the force between the silica particles moving on a closed path along the axis through their respective centers. This mimics the attendant motion of particles inside a filler network spanning an elastomer matrix subject to dynamic strain. From the integrated force along the path we calculate the energy dissipated in the contact and its immediate vicinity. The force is determined using the following force gauge during a simulation at constant volume. Figure 6 depicts the different contributions to the particle− particle interaction. It is the configurational entropy contribution of the polymer, which makes it necessary to measure the interparticle force directly. Merely taking the derivative of the potential energy of the contact with respect to the interparticle distance would be a bad approximation. Here the force is measured using a harmonic spring connected to the right particle. The other side of the spring is connected to a virtual particle, which does not interact with any other atom in the system. For every position of the virtual particle there will be a corresponding position of the right particle satisfying the force equilibrium Fspring = Fcontact. By moving the virtual particle periodically back and forth with a certain amplitude, uo, and frequency, ω, we are able to obtain the desired force-vsdisplacement curves.

Figure 7. Example of force curves. Red symbols: Force between two particles when polymer is not present. Green and blue symbols: Forces in the presence of polymer chains (not cross-linked) measured over five complete cycles. The orientation of the triangles indicate the direction of motion. Open spheres are the respective average values. The velocity of the virtual particle is v = 0.5 m s−1. In all cases the temperature is T = 300 K and the pressure is P = 1.013 bar. The shaded area is the dissipated energy. D


Article

Macromolecules the absence of polymer. The quantity r is the distance between the centers of mass of the two particles minus the particle diameter D. The molecular roughness of the surface however yields slight attendant variations of the minimum r-value if different cycles are compared. We obtain the expected van der Waals interaction exhibiting a minimum located, in this case, at r ≈ 0.5 nm. Without polymer the direction of motion has no effect on the force-vs-displacement curve. This changes when polymer is present. The blue left-triangles and the green righttriangles are force values during the forward and the reverse motions, respectively. The circles are the attendant averages over all cycles (note: subsequent figures do show only average force curves). The net result is a hysteresis in the force-vsseparation diagram from which we compute the dissipated energy per cycle via

Wloss =

∮ F (r ) d r

Figure 8. Basic dependencies of interparticle force curves on surface chemistry and cross-linking (polymer-to-polymer and polymer-tosilane). The reference system (a) consists of two silica particles (D = 4.2 nm) embedded in 85 polyisoprene chains with 200 monomers per chain. When silanes are present, then their surface density is about σ ≈ 0.9 nm−2. This corresponds to a total number of 90 silanes. The silane distribution is uniform, except for the contact area which contains no silanes at all. The velocity of the virtual particle is v = 0.5 m s−1. In all cases the temperature is T = 300 K and the pressure is P = 1.013 bar. Note: (pp) polymer-to-polymer coupling (sulfur content is 15 phr); (sp) silane-to-polymer coupling. Notice also that here θ = 20°.

(1)

Here, hysteresis is due to the viscous behavior of the polymer chains. In the presence of polymer a viscous drag force opposes the motion of the particles, which changes sign when the direction of motion is reversed. We remark that during the cyclic motion of the silica particles polymer chains may become trapped in the immediate contact between the two particles. For non-cross-linked polymers, this is observed every second to every fifth cycle, causing repulsion to occur already at larger particle-to-particle separations. In systems containing cross-linked polymers this happens rarely, and if, in addition, the polymer is chemically linked to the particles surface, the effect is not observed at all on the present time scale and at the present cross-link densities. Because of the prohibitively large computational effort required for a statistically significant evaluation, these infrequent occurrences are not included in our current analysis. In addition, we generally observe that the first cycle differs from the subsequent ones, which is an effect known from the experiments as well. Therefore, we omit the first cycle from our calculations of the force curves. Finally we want to comment on the relation between r and the macroscopic strain u. The exact relation does depend on the filler distribution, which we do not consider here. For a simple one-dimensional chain of monodisperse particles embedded in an elastomer matrix the relation is u ≈ (r − r0)/D, where r0 is the particle-to-particle separation at vanishing strain. This point has been discussed in more detail in ref 25.

attendant energy dissipation. A major reason is the exclusion of silane from the particle-to-particle contact. We return to this point below when we investigate the influence of silane coverage in more detail. In addition, at the silane coverage in the present example the effect on the near-surface polymer structure and dynamics is rather small and affects only the very first polymer layer in contact with the particle surface (cf. below). The next system, system c, contains cross-linked polymer chains. Here we observe a reduction in the area of the force hysteresis and an attendant decrease of the dissipated energy. Reduced dissipation is also observed experimentally when, for instance, the sulfur concentration is increased, enhancing the chemical cross-linking between polymer chains (e.g., Figure 34 in ref 36). Dissipation increases, relative to reference a, if instead of cross-linking the polymer chains themselves we merely bond the silanes to the surrounding polymer (system d). This can be understood in terms of an increased viscous drag induced in the polymer matrix by the moving particle. If, however, we combine the bonding of silane to the polymer with chemical cross-linking between polymer chains (system e), then we obtain the overall lowest dissipation. Notice that the dissipated energies for (sp) and (pp) are not simply additive, i.e. Wloss(sp + pp) < Wloss(sp) + Wloss(pp) . We can understand this in terms of viscous loss caused by relative motion of polymer segments with respect to neighboring segments or with respect to the particle surface. The option sp, i.e., polymers bound to the particles via covalent bonds with silanes, causes the largest drag, enhancing relative motion and attendant friction of polymer chain segments due to the moving particles. The pp-coupling, on the other hand, reduces the average relative segment motion. Finally we observe that pp + sp leads to the strongest restoring force, i.e., the steepest slope of the force loop, in comparison to all other systems. Relating Temperature and Frequency. The force curves shown thus far are obtained at very high frequencies, i.e., frequencies in the range of GHz (cf. below), which is typical for

■

SIMULATION RESULTS Selected Basic Effects on Energy Dissipation due to Silanization and Chemical Coupling. In this section we discuss selected basic effects due to silanization and chemical coupling on dissipative loss in comparison to a reference system. The latter consists of two silica particles without silanes embedded in a matrix of non cross-linked polymer chains. Figure 8 compares the force curve and the resulting loss of this system (a) to the same system when the particles are silanized (b), the particles are silanized and the polymer is cross-linked (pp) (c), the particles are silanized and silanes are linked to the polymer (sp) (d), and the particles are silanized and cross-links are introduced in the polymer matrix as well as between silane and polymer chains (sp + pp) (e). Notice that there is no silane in the contact area between the particles. Comparing system b to the reference a, we find no clearly discernible difference regarding the force curves or the E


Article

Macromolecules

minimum, is due to the different contact areas devoid of silane in the two cases. Simultaneously the slope of the force curve at larger distances is decreased with increasing temperature, because of the softening of the polymer (Notice that “softening” here is a subtle issue due to competing enthalpic and entropic effects (cf. thermoelastic inversion).). Figure 9b shows the loss per (simulation) volume versus the temperature for different pulling-velocities. All sets of data exhibit a peak whose position depends on the velocity of the virtual particle and thus on frequency. For further investigation we fit the data using a Cole−Cole-model37 commonly used to describe the temperature dependence of dielectric and mechanical loss:

atomistic computer simulations. Frequencies this high do pose an experimental challenge. Thus, we may want to study their effects using simulations just for this reason. On the other hand and perhaps more importantly, technically relevant frequencies in dynamic elastomer applications are often much lower. A remedy, allowing to extend our present results, at least partially, to significantly lower frequencies, is the time−temperaturesuperposition principle applicable in polymer systems (e.g., ref 1). This requires that we extend our calculations of the dissipative loss to a wide range of temperatures. Because the approach is computationally expensive, we concentrate here on an example, i.e. a fully cross-linked system (case e in Figure 8, but with a silane density of σ = 2.7 nm−2 and a cone angle of 2 × 40°). The contact area remains free of silanes. Selected force curves obtained at three different temperatures are shown in Figure 9a. In accord with experimental data, we observe a

fCC (T ) = G0r(α , T0)1/2 sin(ϕ(α , T0)) + G1 r = X12 + X 2 2 ⎛X ⎞ ϕ = arctan⎜ 2 ⎟ ⎝ X1 ⎠ ⎛ απ ⎞ X1 = sin⎜ ⎟(1.0 + exp{(T0 − T )(1.0 − α)}) ⎝ 2 ⎠ ⎛ απ ⎞ X 2 = cos⎜ ⎟ exp{(T0 − T )(1.0 − α)} ⎝ 2 ⎠

(2)

G0, G1, and α are fit parameters meant to describe the shape of the relaxation peak, whereas T0 indicates the position of the peak maximum. In the following we want to make use of the temperature− frequency-superposition in the form ωτ0 ∼ exp(T − T0)

(3)

where T0 is the temperature at the peak position in Figure 9b. This enables us to apply the resulting f CC to our data. We note that the additional constant G1 improves the fits significantly, which is mainly due to the nonvanishing loss at high temperatures. The underlying mechanism is different from polymer viscosity, as already discussed in previous work,25,38,39 and intimately related to the specific form of the particle-toparticle interaction. A detailed discussion of this dissipative mechanism in the context of the present molecular simulations is forthcoming. We fit eq 2 to the simulation results in Figure 9b requiring all variables to be positive and in addition α < 1.0. The resulting fit-parameters are compiled in Table 1. The overall agreement between the simulated curves and eq 2 is quite satisfactory. Notice that the simulation box is small compared to experimental samples. Simulation results therefore exhibit large fluctuations. Extending the temperature range, which is of course possible only in the simulation, improves the fit quality and thus the measurement of the segment relaxation peak maximum vs frequency. It also improves the fit reliability in the entire range of physically realistic temperatures. Next we want to relate the above to experimental relaxation times for pure polyisoprene. Notice that the experimentally derived loss moduli (cf. refs 40−42), μ″, are related to the dissipated work per volume, w, via w = πμ″uo2, where uo is the strain amplitude, assuming linear viscoelastic behavior (cf. above). Because the displacements of the virtual particle at the different temperatures in Figure 9 are the same, we expect μ″ ∝ Wloss/Vsim. One might object to this because filled elastomers exhibit strong nonlinearity as manifested in the Payne effect as

Figure 9. Panel a shows the effect of selected temperatures on the force curves obtained for system e in Figure 8 at a higher silane density of σ = 2.7 nm−2. Panel b shows the resulting dissipated energy plotted versus temperature at different frequencies. In all cases, the contact area between the two particles does not contain silanes. The data points were fitted with the imaginary part of the Cole−Cole equation (eq 2)37 as described in the text.

reduction of the hysteresis loop with increasing temperature (e.g., Figure 7 in ref 26 showing the temperature dependent loss modulus, μ″. At constant strain amplitude we may assume μ″ ∝ Wloss.). Notice that the minimum exhibited by the green curve in Figure 9a, compared to the green curve, i.e., the curve at the same temperature, in Figure 8, which shown no such F


Article

Macromolecules Table 1. Parameter Values in Equation 2 Used To Describe the Simulation Results Shown in Figure 9ba v/m s−1

ω/109 Hz

0.50 2.50 3.25 5.00

0.12 0.60 0.77 1.19

α

G0/MPa 645.3 927.1 915.6 628.4

± ± ± ±

86.9 212.0 150.0 60.2

0.9857 0.9908 0.9911 0.9891

± ± ± ±

T0/K 0.0017 0.0014 0.0009 0.0008

348.4 385.0 394.0 398.1

± ± ± ±

7.3 10.1 6.7 4.1

G1/MPa 0.38 0.34 0.38 0.89

± ± ± ±

0.20 0.35 0.25 0.12

The errors are the estimated covariances. Notice that v is the velocity of displacement of the virtual particle, while ω is the attendant frequency of the cyclic displacement. a

mentioned in the introduction. However, because w is a welldefined quantity in any case, the above equation for w in terms of μ″ can be used as a definition, regardless of whether the response is linear or not, and usually this is what is done. Here we do not really study a pure polymer system, but a system containing two interacting filler particles. It is well-known that the presence of filler affects the magnitude of the loss modulus significantly (e.g., Figure 6 in ref 43). Nevertheless, in ref 44 (Figure 6), it is shown that the relaxation of segmental modes in filled polyisoprene does not change significantly. For other relaxation modes this does not hold true, where in particular the dispersion of silica plays an important role. In the present simulations, due to the small size of the particles, the loss mainly is related to viscoelastic processes in the pure polymer, albeit induced by the moving particle, and thus we momentarily sidestep the effects of fillers. Figure 10 compares our results for ω vs T0 from Table 1 to measured peak positions of the segmental relaxation (glass

means. However, Figure 10 provides the sought after mapping, which enables us to relate our high frequency results to experimental results at much lower frequencies. For example, at a frequency of 1 Hz the attendant T1Hz 0 is ≈218 K. According to . Thus, if we Table 1, this T0 is about 180 K lower than T1GHz 0 want to obtain dynamical information at for instance 60 °C (≈ 333 K) and 1 Hz, we must adjust our simulation temperature to ≈513 K when the simulation frequency is ≈1 GHz. We should add a note of caution however. Our discussion applies to polymer dynamics. If we study phenomena that couple the polymer dynamics to other temperature dependent effects, e.g., hydrogen bonding between filler particles, the above procedure may not be applicable. Variation of Silane Density and Surface Distribution. Here we address the influence of silane density and distribution (naked and fully covered contact area) based on different systems prepared as follows: In system I, the two particles are close to fully covered by silane. This system is replicated yielding system II. Subsequently the silanes in the contact area of system II are deleted. Following the above procedure polymer is added to both I and II. System II, now including polymer but without silane in the contact area, is replicated again to yield systems II.b and II.c in addition to system II. In the replicas (II.b and II.c), we delete at random 50% and 100% of the silanes, respectively, whereas the original system II remains unaltered. The position of the deleted silanes are saved, however so that the identical deletion procedure can be applied to replicas of system I, which we label I.b and I.c. A pictorial overview of this description is shown in Figure 11. Notice that silane deletion is applied only to silanes outside contact areas and thus I.c ≠ II.c! Finally, in all systems, the polymer is vulcanized and the silanes are linked to the polymer as described above. Notice also that the small particle size makes it difficult to study more detailed silane distributions with statistical significance. It is of course interesting to investigate particularly the effects of silane distribution in the contact on the contact’s contribution to dissipative loss. The resulting force curves of the six systems are shown in Figure 12. Panels a and b in Figure 12 do have in common that an increase in the silane density increases the bonding between the particles, which leads to an increase of the restoring force at large separations. Overall this effect is stronger in Figure 12b, owing to the fact that there are additional silanes due to the silane filled contact areas. One must bear in mind that the contact area is an approximate concept. Silanes near the rim but inside the contact area can be bound to polymer segments which in turn are bound to similarly located silanes in the opposite particle’s contact area. Because of the shortness of the load bearing path thus created, this can significantly enhance the restoring force. When the particle-to-particle separation becomes small, the two types of systems produce quite different force curves. The system type I in panel b, with silane in between the particle

Figure 10. Position of the segmental relaxation peak in the temperature−frequency plane. The simulation results (red downtriangles) are taken from Table 1. Experimental results for pure polyisoprene (not chemically cross-linked) were obtained by different methods and for a range of molecular weights: (i) MW = 97000 g mol−1 from ref 40; (ii) MW = 436000 g mol−1 from ref 41; (iii) MW = 504000 g mol−1 from ref 42.

transition) in pure polyisoprene (molecular weights Mw ≥ 13600 g mol−1 ≅ 200 monomers). In this region the segmental modes are independent of chain length.41 Aside from a comparatively small offset, the reason for which currently is not clear (possible factors include cross-linking45,46 as well as the interface between the particles and the bulk matrix44 and, truncation of long-range hydrodynamic interactions at the box boundaries), our data are in accord with the experimental values, and thus the loss peak in Figure 9b may be identified with the glass-transition. We note that the glasstransition can be measured using Molecular Dynamics simulations for similar polymer models (e.g., ref 47) by other G


Article

Macromolecules

Figure 12. Force curves obtained for different silane densities and distributions. In all cases the rubber is vulcanized and the silanes are coupled to the polymer. The indicated surface densities refer to the surface outside the contact area. The silane density inside the surface area is zero for all II systems in panel a and 2.7 nm−2 for all I systems in panel b.

Figure 11. System preparation schema for the silane study. Silanes inside the contact area are shown in red.

surfaces, exhibits an almost linear increase of the repulsive force. System type II in panel a, with no silanes between the particle surfaces, exhibits a pronounced minimum due to interaction of the hydroxyl groups (In this r-range there is little difference between systems II, II.b, and II.c). Apparently this does not have an influence on the viscous loss. However, we expect that, based on a model presented in refs 25, 38, and 39 (cf. above), there is an independent dissipation mechanism involving short-range inter particle attraction as an essential element when the particle size increases. Notice also that the minimum exhibited by the green curve in Figure 12a, compared to the green curve, i.e., the curve at the same temperature, in Figure 8, which shown no such minimum, is due to the different contact areas in the two cases. Finally we observe that regardless of the presence of silane in the contact the loss decreases with increasing silane coverage (in accord with experimental observations (e.g., Figure 34 in ref 36)). In addition, the loss is diminished by silane in the contact in comparison to systems without silane in the contact area. We attribute this to the “naked” contacts in the II-systems, causing an additional “stickiness” between the particles due to the direct silica−silica-interaction. The attendant interaction potential in conjunction with the restoring force of the polymer matrix gives rise to an additional loss-mechanism during cyclic deformation of the contact. The mechanism as such has been discussed previously38,39 (cf. above). This loss and the viscous loss in the polymer yield the total loss. In a forthcoming publication, we will show how to separate the two loss contributions, and we

have studied their system dependent relative importance. In the present work, the viscous loss is the dominant contribution. Near-Surface Polymer Structure and Mobility. In the introduction we had suggested the present filler−filler contact model as the basic structural element within the filler network. Subsequently we have studied the dynamical viscous loss associated with the contact depending on selected chemical variables. In the following we are interested in the apparent spatial extend of the contact from which the dissipation originates. We address this point via two simple scalar quantities, the (radial) density profile and the mobility in close proximity of the particle surfaces. We focus on the particle−polymer interfaces opposite to the contact between the particle in order to exclude the complicated geometry in the immediate contact. In the following the polymer is not vulcanized and the silanes are not coupled to the polymer. The results are based on simulations of length t = 10 ns, during which the particle positions are kept fixed. Figure 13a shows radial polymer density profiles for different silane densities. Without silanes on the surface (σ = 0 nm−2) three density peaks are discernible. At full coverage (σ = 2.7 nm−2) the silanes have no other option than a side by side arrangement while being perpendicular to the surface. Thus, the position of the first density peak is shifted by about 0.5 nm, which corresponds to the size of one extended silane. In comparison to the naked surface the first density peak vanishes H


Article

Macromolecules

N1 + N2 = Nri(t = 0)

(5)

while τ1 (>τ2) and τ2 are freely adjustable. Note that acceptable fits require Nri(t = 0) > 20 monomers. The attendant residence or relaxation times, τi, depending on the respective shell radius, are shown in Figure 13b. The short relaxation time, τ2, does not exhibit a discernible dependence on the shell radius, whereas the long relaxation time, τ1, exhibits a clearly discernible decrease with increasing shell radius. A significant difference to the bulk values can be found in the same range over which the particles influence the polymer structure. Notice that there are no data points for the highest silane density very close to the particle surface. This is because of the negligible monomer density in this range. A second point is the apparent increase of the relaxation time τ1 for intermediate silane density (i.e., σ = 1.4 nm−2). This can be explained in terms of ’trapping’ of monomers within the silane layer if the latter is not too dense. Nevertheless, there is significant scatter of the data and such interpretations require considerable caution. A final point worth noting concerns the contact subject to external forces, which here are not applied. The attendant hydrodynamic effects induced in the polymer matrix decay algebraically (e.g., ref 51) and thus will be truncated by the box boundaries. However, similar truncation occurs in real highly filled elastomers due to the presence of other filler particles in close proximity of a filler−filler contact.

■

CONCLUSION In this work, we discuss an atomistic model of a filler-to-filler contact in a model elastomer. We expect this type of filler-tofiller contact to be important in the interplay between the morphology of filler networks and the attendant dynamic moduli of the filled rubber. Using a specially designed force gauge, we calculate the energy dissipated in a molecular dynamics simulation of the contact as it opens and closes during a cyclic deformation. We investigate the influence of different chemical parameters on the dissipated energy, e.g., sulfur cross-linking, silane-to-polymer linking, silane distribution, and surface density on the silica nanoparticles. In addition, we study our model contact at variable temperature, which provides a means to map high frequency results to much lower frequencies at increased temperature. Overall, the results are in accord with experimental findings for real filled rubbers, e.g. the reduction of heat build up when the sulfur content is increased, the loss reduction when the temperature is increased (on the high temperature side of the glass transition) or the calculated relaxation-modes for varying temperature and frequency. We are aware, however, that the current model is only one part in a complex system with additional dissipative mechanisms. Nevertheless, we believe that our model can provide the developer with valuable molecular information regarding a technically important hidden interface, i.e., the filler-to-filler interface, which has a pronounced effect on product performance.

Figure 13. Polymer behavior in the vicinity of the particles for three different silane densities. Panel a shows the radial density profile, and panel b shows the corresponding residence time profiles. An additional result obtained for the bulk-polymer system (without silica particles) is included as a reference.

and the second peak is shifted by approximately 0.1 nm. The previously third and least pronounced density peak can no longer be distinguished from the bulk polymer density. A reduction of the silane density to σ = 1.4 nm−2 leads to the vanishing of the first pronounced density peak while the second and third peak remain nearly unaltered. The apparent range of influence of the particles is ≈1.5 nm, which is in accord with previous results from refs.48−50 The near-surface mobility of the polymer is studied here analogous to the approach discussed in ref 21. On the basis of the trajectory at t = 0 ps, we sort all monomers into (half) shells of width 0.2 nm centered on the particle’s centers of mass. During the simulation we check every 10 ps how many of the monomers still do remain inside their original shells. Monomers leaving and remaining outside their shell for more than 100 ps are no longer counted as still residing in this shell. Following this procedure we obtain the number of resident monomers per shell at time t, Nri(t), where ri is the radius of shell i. As in ref 21, we fit Nri(t) by a double exponential function, i.e. ⎛ t ⎞ ⎛ t⎞ Nfit = N1 exp⎜ − ⎟ + N2 exp⎜ − ⎟ ⎝ τ2 ⎠ ⎝ τ1 ⎠

■

AUTHOR INFORMATION

Corresponding Author

*(J.M.) E-mail: [email protected]. Telephone: +49 202 439 3635. Fax: +49 202 439 3122. (4)

ORCID

Reinhard Hentschke: 0000-0003-2739-9552

where I


Article

Macromolecules Notes

(18) Xi, H.; Hentschke, R. The influence of structure on mechanical properties of filler networks via coarse-grained modeling. Macromol. Theory Simul. 2014, 23, 373−382. (19) Raos, G.; Casalegno, M. Nonequilibrium simulations of filled polymer networks: Searching for the origins of reinforcement and nonlinearity. J. Chem. Phys. 2011, 134, 054902. (20) Gusev, A. A. Micromechanical Mechanism of Reinforcement and Losses in Filled Rubbers. Macromolecules 2006, 39, 5960−5962. (21) Hager, J.; Hentschke, R.; Hojdis, N. W.; Karimi-Varzaneh, H. A. Computer Simulation of Particle-Particle Interaction in a Model Polymer Nanocomposite. Macromolecules 2015, 48, 9039−9049. (22) Payne, A. R. Effect of dispersion on the dynamic properties of filler-loaded rubbers. J. Appl. Polym. Sci. 1965, 9, 2273−2284. (23) Ouyang, G. B. Network junction model for carbon black reinforcement: Modulus, hysteresis and the Payne effect. Kautschuk Gummi Kunststoffe 2006, 59, 332−343. (24) Ouyang, G. B. Network junction model for carbon black reinforcement: The Payne effect Part II. Kautschuk Gummi Kunststoffe 2006, 59, 454−458. (25) Hentschke, R. The Payne Effect Revisited. eXPRESS Polym. Lett. 2017, 11, 278−292. (26) Stöckelhuber, K. W.; Svistkov, A. S.; Pelevin, A. G.; Heinrich, G. Impact of Filler Surface Modification on Large Scale Mechanics of Styrene Butadiene/Silica Rubber Composites. Macromolecules 2011, 44, 4366−4381. (27) Zhang, P.; Morris, M.; Doshi, D. Materials Development for lowering Rolling Resistance of Tires. Rubber Chem. Technol. 2016, 89, 79−116. (28) Basic Characteristics of AEROSIL Fumed Silica; Evonik Industries: 1967; Vol. Technical Bulletin Fine Particles 11. (29) Mark, J. Interpretation of random-coil configurations of trans-1, 4-polybutadiene and trans-1, 4-polyisoprene. J. Am. Chem. Soc. 1967, 89, 6829−6835. (30) McPherson, A. T.; et al. Density and electrical properties of the system, rubber-sulphur Part I., Density of rubber-sulphur compounds. Bur. Stand. (U. S.), Sci. Pap. 1927, 22, 385−397. (31) Nemoto, N.; Moriwaki, M.; Odani, H.; Kurata, M. Shear Creep Studies of Narrow- Distribution Poly(cis-isoprene). Macromolecules 1971, 4, 215−219. (32) Rappe, A. K.; Goddard, W. A. Charge equilibration for molecular dynamics simulations. J. Phys. Chem. 1991, 95, 3358−3363. (33) Harmandaris, V. A.; Doxastakis, M.; et al. Detailed molecular dynamics simulation of the self-diffusion of n-alkane and cis-1,4 polyisoprene oligomer melts. J. Chem. Phys. 2002, 116, 436−446. (34) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (35) Alznauer, T.; Becker, J. A. Manuscript in preparation. (36) Luginsland, H. A Review on the Chemistry and the Reinforcement of the Silica Silane Filler System for Rubber Applications; Shaker Verlag: Aachen, Germany, 2002. (37) Cole, K. S.; Cole, R. H. Dispersion and Absorption in Dielectrics I. Alternating Current Characteristics. J. Chem. Phys. 1941, 9, 341− 351. (38) Hentschke, R. In Constitutive Models for Rubber VIII; GilNegrete, N., Alonso, A., Eds.; CRC Press/Balkema: Leiden, The Netherlands, 2013; p 299−314. (39) Hentschke, R.; Hager, J.; Hojdis, N. W. Molecular Modeling Approach to the Prediction of Mechanical Properties of SilicaReinforced Rubbers. J. Appl. Polym. Sci. 2014, 131, 40806. (40) Boese, D.; Kremer, F. Molecular dynamics in bulk cispolyisoprene as studied by dielectric spectroscopy. Macromolecules 1990, 23, 829−835. (41) Abou Elfadl, A.; Kahlau, R.; Herrmann, A.; Novikov, V. N.; Rössler, E. A. From rouse to fully established entanglement dynamics: A study of polyisoprene by dielectric spectroscopy. Macromolecules 2010, 43, 3340−3351.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Prof. Dr. J. A. Becker and Dr. T. Alznauer (University of Hannover) for fruitful discussions and provision of data concerning sulfur cross-links. We also thank F. Fleck (Continental Reifen Deutschland GmbH) for educating us with respect to various experimental aspects of this work.

■

REFERENCES

(1) Wrana, C. Introduction to Polymer Physics; Lanxess AG: Leverkusen, Germany, 2009. (2) Vilgis, T. A., Heinrich, G., Klüppel, M. Reinforcement of Polymer Nano-Composites: Theory, Experiments and Applications; Cambridge University Press: Cambridge, U.K., 2009. (3) Tjong, S. C.; Mai, Y.-W., Eds. Physical Properties and Applications of Polymer Nanocomposites; Woodhead Publishing Series in Composites Science and Engineering; Woodhead Publishing: 2010. (4) Thomas, S.; Stephen, R. Molecular Dynamics Simulation of Macromolecular Interactions in Solution: Poly(y-benzyl glutamate) in Dimethylformamide and Tetrahydrofuran. Rubber Nanocomposites: Preparation, Properties and Applications. Wiley: New York, 2009. (5) Helfrich, J.; Hentschke, R. Molecular Dynamics Simulation of Macromolecular Interactions in Solution: Poly(y-benzyl glutamate) in Dimethylformamide and Tetrahydrofuran. Macromolecules 1995, 28, 3831−3841. (6) Brown, D.; Marcadon, V.; Mélé, P.; Albérola, N. D. Effect of Filler Particle Size on the Properties of Model Nanocomposites. Macromolecules 2008, 41, 1499−1511. (7) Ndoro, T. V. M.; Böhm, M. C.; Müller-Plathe, F. Interface and Interphase Dynamics of Polystyrene Chains near Grafted and Ungrafted Silica Nanoparticles. Macromolecules 2012, 45, 171−179. (8) Guseva, D. V.; Komarov, P. V.; Lyulin, A. V. Molecular-dynamics simulations of thin polyisoprene films confined between amorphous silica substrates. J. Chem. Phys. 2014, 140, 114903. (9) Riggleman, R. A.; Toepperwein, G.; Papakonstantopoulos, G. J.; Barrat, J.-L.; de Pablo, J. J. Entanglement network in nanoparticle reinforced polymers. J. Chem. Phys. 2009, 130, 244903. (10) Jaber, E.; Luo, H.; Li, W.; Gersappe, D. Network formation in polymer nanocomposites under shear. Soft Matter 2011, 7, 3852− 3860. (11) Kalathi, J. T.; Grest, G. S.; Kumar, S. K. Universal Viscosity Behavior of Polymer Nanocomposites. Phys. Rev. Lett. 2012, 109, 198301. (12) Chen, Y.; Liu, L.; Yang, Q.; Wen, S.; Zhang, L.; Zhong, C. Computational Study of Nanoparticle Dispersion and Spatial Distribution in Polymer Matrix under Oscillatory Shear Flow. Langmuir 2013, 29, 13932−13942. (13) Chen, Y.; Li, Z.; Wen, S.; Yang, Q.; Zhang, L.; Zhong, C.; Liu, L. Molecular simulation study of role of polymer-particle interactions in the strain- dependent viscoelasticity of elastomers (Payne effect). J. Chem. Phys. 2014, 141, 104901. (14) Masnada, E.; Merabia, S.; Couty, M.; Barrat, J.-L. Entanglementinduced reinforcement in polymer nanocomposites. Soft Matter 2013, 9, 10532−10544. (15) Merabia, S.; Sotta, P.; Long, D. R. A Microscopic Model for the Reinforcement and the Nonlinear Behavior of Filled Elastomers and Thermoplastic Elastomers (Payne and Mullins Effects). Macromolecules 2008, 41, 8252−8266. (16) Merabia, S.; Sotta, P.; Long, D. R. Unique Plastic and Recovery Behavior of Nanofilled Elastomers and Thermoplastic Elastomers (Payne and Mullins Effects). J. Polym. Sci., Part B: Polym. Phys. 2010, 48, 1495−1508. (17) Xi, H.; Hentschke, R. Dynamic moduli of filled elastomers - A coarse grained computer model. Eur. Polym. J. 2012, 48, 1777−1786. J


Article

Macromolecules (42) Santangelo, P.; Roland, C. Temperature Dependence of Mechanical and Dielectric Relaxation in cis-1,4-Polyisoprene. Macromolecules 1998, 31, 3715−3719. (43) Fritzsche, J.; Klüppel, M. Structural dynamics and interfacial properties of filler-reinforced elastomers. J. Phys.: Condens. Matter 2011, 23, 035104. (44) Fragiadakis, D.; Bokobza, L.; Pissis, P. Dynamics near the filler surface in natural rubber-silica nanocomposites. Polymer 2011, 52, 3175−3182. (45) Hernandez, M.; Ezquerra, T. A.; Verdejo, R.; Lopez-Manchado, M. A. Role of vulcanizing additives on the segmental dynamics of natural rubber. Macromolecules 2012, 45, 1070−1075. (46) Valentín, J. L.; Posadas, P.; Fernández-Torres, A.; Malmierca, M. A.; González, L.; Chassé, W.; Saalwächter, K. Inhomogeneities and chain dynamics in diene rubbers vulcanized with different cure systems. Macromolecules 2010, 43, 4210−4222. (47) Sharma, P.; Roy, S.; Karimi-Varzaneh, H. A. Validation of Force Fields of Rubber through Glass-Transition Temperature Calculation by Microsecond Atomic-Scale Molecular Dynamics Simulation. J. Phys. Chem. B 2016, 120, 1367−1379. (48) Eslami, H.; Behrouz, M. Molecular dynamics simulation of a polyamide-66/carbon nanotube nanocomposite. J. Phys. Chem. C 2014, 118, 9841−9851. (49) Eslami, H.; Rahimi, M.; Müller-Plathe, F. Molecular Dynamics Simulation of a Silica Nanoparticle in Oligomeric Poly (methyl methacrylate): A Model System for Studying the Interphase Thickness in a Polymer-Nanocomposite via Different Properties. Macromolecules 2013, 46, 8680−8692. (50) Eslami, H.; Müller-Plathe, F. How Thick is the Interphase in an Ultrathin Polymer Film? Coarse-Grained Molecular Dynamics Simulations of Polyamide-6,6 on Graphene. J. Phys. Chem. C 2013, 117, 5249−5257. (51) Felderhof, B. U. Hydrodynamic of suspensions. In Fundamental Problems in Statistical Mechanics VII; van Beijren, H., Ed.; NorthHolland: Amsterdam, 1990.

K


Molecular Simulation of Viscous Dissipation due to Cyclic Deformation

Recommend Documents