Computer Simulation of Particle–Particle Interaction in a Model

Article pubs.acs.org/Macromolecules

Computer Simulation of Particle−Particle Interaction in a Model Polymer Nanocomposite Jonathan Hager,† Reinhard Hentschke,*,† Nils W. Hojdis,‡ and Hossein Ali Karimi-Varzaneh‡ †

Fachbereich Mathematik und Naturwissenschaften, Bergische Universität, D-42097 Wuppertal, Germany Continental Reifen Deutschland GmbH, D-30419 Hannover, Germany

‡

ABSTRACT: We present a fully atomistic simulation study of a filler-to-filler contact embedded in a polymer matrix. The system consists of two silica spheres embedded in polyisoprene. Tests are presented showing that our polymerization algorithm achieves the proper polymer density and characteristic ratio. In this study, the polymer chains are not chemically attached to the particle surfaces and neither are they chemically cross-linked. We apply an external force in the direction of the axis of the system to one of the silica particles and determine the attendant distance to the second particle. The force−distance curves confirm a previously suggested loss mechanism, allowing us to estimate the dissipated energy when a contact is opened and closed due to cyclic loading. We also study the atomic structure and mobility inside the polymer close to the particle surfaces. We do find structural ordering in terms of three appreciable density peaks. The atom mobility inside these density shells is characterized by a fast and a slow process. The latter shows an apparent dependence on the shell’s distance from the particle surface.

1. INTRODUCTION Material properties of polymer matrices are modified routinely by addition of inorganic nanoparticles, i.e., particles with linear dimensions between 1 and 100 nm.1,2 The modifications encompass dynamical mechanical behavior, friction, abrasion, thermal and electrical conductivity, and others. Industrial nanocomposites usually are complex compounds containing not just two but many ingredients combined in an elaborate processing scheme in order to match a broad spectrum of product requirements. Much of the attendant research and development therefore uses an experience-based trial and error approach. One particular difficulty is to understand chemical and physical processes in the “hidden interfaces” between the nanoparticels or between nanoparticles and the surrounding polymer matrix. These interfaces, however, do to large extend govern the aforementioned material properties. In the following we discuss a molecular modeling approach, which is targeting the relation between the above interfaces and dynamical mechanical properties of filled rubbers.2−7 An effect of major importance in this context was extensively studied in the 1960’s A. R. Payne (an early review is ref 8; more recent reviews are refs 6, 9, and 10). The first observation of this effect, however, dates back to the early 1940s.11 The Payne effect, as it is called now, describes the pronounced decrease of the storage modulus, μ′, with increasing strain amplitude, u, in filled rubbers under cyclic loading. Because this does not occur in unfilled rubbers, the cause of the Payne effect must be related to either the rubber−filler or the filler−filler interface(s). Most authors attribute the effect to a breakdown of the filler network in the rubber matrix when the material is strained. On the basis of this concept, we have suggested a dissipative loss mechanism on the molecular scale.12 The idea is that particle aggregates inside a filled elastomer interact with © XXXX American Chemical Society

nearest-neighbor aggregates via short-range van der Waals forces (the long-range Hamaker interaction is not expected to be a major factor9) as well as via (elastic) restoring forces in response to the deformation of the polymer matrix. Assuming local equilibrium of these force contributions allows to spontaneously open and close contacts between aggregates in the branches of a filler network due to external stress. Thereby the contacts dissipate energy into the surrounding polymer matrix. In a previous paper,13 we have studied a simplified variant of this contact model. In this simpler model the van der Waals forces are calculated based on distance-dependent interaction energies of flat silica surfaces covered with alkylsilanes at variable concentration modeled via molecular dynamics simulation. The restoring force due to the surrounding rubber matrix is replaced by a linear spring, whose spring constant is based on the experimental rubber shear modulus. Figure 9 in ref 13 shows the explicit (and quite good) comparison between our model and dynamic mechanical measurements of the storage and the loss modulus at selected strain amplitudes vs silane surface concentration for different alkylsilanes. We like to emphasize that other authors have studied the response, i.e., dynamic moduli as well as stress-induced structural changes, to an entire collection of model filler particles embedded into a polymer matrix, which is subject to dynamic stress (e.g., refs 14−20). Another coarse-grained model for the viscosity of filled entangled polymer melts is ref 21. However, our goal here is different. We focus on a distinct Received: August 24, 2015 Revised: October 30, 2015

A

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

low strain storage modulus, μ′, on filler volume fraction, ϕ. Near the filler percolation threshold this dependence has the form of a power law, Δμ′ ∼ ϕy, where y can be as high as 3.5 (or even higher as shown in ref 24). Δ means that μ′ is calculated relative to its value in the limit of large strain. Our model focuses on a single “breakable” aggregate-toaggregate contact within the filler network as defined in Figure

molecular mechanism, which, in our opinion, should be responsible for much of the Payne effect. The simplifying assumptions made in ref 13 are eliminated in the present work, where we present a fully atomistic simulation study of a filler-to-filler contact embedded in a polymer matrix. The filler in this work is silica, and the matrix polymer is polyisoprene. We start with two spherical silica particles, at close contact, inside liquid isoprene. Using a fast polymerization algorithm, we transform the liquid into about 100 monodisperse chains consisting of 170 monomers each (these numbers are comparable to previous studies of isolated silica particles embedded in polymer; cf. ref 22). Tests are presented showing that the algorithm achieves the proper polymer density and characteristic ratio. In addition, exploratory simulations are presented for particles covered with alkylsilanes, used in practice as simple compatibilizers between the hydrophilic silica and the hydrophobic polymer. We are aware that adding chemical bonding between the polymer and fillers (grafting density) changes the morphology of the filler and therefore influences the moduli (cf. ref 13). For simplicity, however, in the current work the polymer chains are not chemically bonded to the filler surface and neither are they chemically cross-linked. We apply an external force in the direction of the axis of the system to one of silica particles and determine the attendant distance to the second particle. The resulting force−distance curves confirm the above loss mechanism. In addition, they allow to estimate the dissipated energy when a particle−particle contact opens and subsequently closes due to a cyclic load applied to the contact. We also study the atomic structure and mobility of the polymer near the particle surfaces. We do find structural ordering in terms of three appreciable density peaks near the particle surface. The mobility of the atoms in these density shells is characterized by a fast and a slow process. The latter shows an apparent dependence on the shell’s distance from the particle surface. The paper is structured as follows: In section 2 we summarize our model as presented in the aforementioned refs 12 and 13. Section 3 starts with an overview of the molecular model for the contact between silica nanoparticles in polyisoprene on which this paper is based. We then discuss the simulation methods and parametrization including in particular a polymerization algorithm for the polymer matrix. Finally, we describe the measurement of the force between the nanoparticles in our model. Section 4 discusses the results, i.e., the forces as a function of particle-to-particle distance as well as the polymer structure and dynamics in the immediate vicinity of the particles. Section 5 is the concluding discussion.

2. NANOMECHANICAL MODEL FOR DISSIPATIVE LOSS IN FILLER NETWORKS For reasons of clarity we summarize the nanomechanical model for dissipative loss in filler networks developed in refs 12 and 13. This model, which we extend in this work, yields a possible molecular basis for quantitative calculation of the Payne effect. Transmission electron micrographs (TEM) (e.g., ref 3, Figure 7.3 in ref 6, or Figure 3.8 in ref 23) of polymer nanocomposites, under conditions relevant for technical rubber products, do exhibit filler network morphologies consisting of fractal “flocs” of agglomerates. These in turn consist of aggregates composed of primary filler particles. The aggregates are considered “unbreakable”, whereas their superstructures are reversibly stabilized by comparatively weak molecular forces. The picture is supported by the pronounced dependence of the

Figure 1. Top: simulated TEM picture of a filler network inside a rubber matrix (not shown here) generated by a coarse-grained computer model of filled rubber.25,26 Each gray circle is a primary filler particle. Middle: sketch of a inter aggregate contact corresponding to the content of the red box in the previous panel. Here the shaded area is the “rubber collar” referred to in the text. Bottom: cartoon of filler network strands build from contacts depicted in the middle panel. Here each shaded circle is an entire aggregate.

1. The total inter aggregate force in this contact, as shown in Figure 2, is modeled via f (d) = fpp (d) + fmatrix (d)

(1)

Here f pp is the direct force between two adjacent primary particles, each being part of a different aggregate (in principle, B

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 2. Interaggregate force, f, vs interaggregate separation, d (solid line). The dashed lines are the two contributions to f, i.e., f pp and f matrix. The meaning of the path indicated by the red arrows as well as fs and Wloss is explained in the text.

Figure 3. Cutaway view of the central section of the simulated system. Two spherical SiO2 particles embedded in a polyisoprene matrix consisting of at total of 99 monodisperse chains with 170 monomers each. Notice also the blue alkyl moieties of the compatibilizer molecules on the particle surfaces containing eight carbon atoms (in the prequel to this work, i.e.,13 we additionally consider alkylsilanes with alkyl moieties containing 3 and 16 carbon atoms, respectively). The total simulation volume is 13.29 nm × 13.29 nm × 12.65 nm.

there may be more than one direct contact between two aggregates). f matrix is the force in response to straining the “rubber collar” joining the two aggregates. Both forces do depend on the aggregate-to-aggregate separation d. An external force acting on a material as depicted in Figure 1 will give rise to the additional force fs, which does not depend on d. Assuming local equilibrium, the equilibrium aggregate-toaggregate separation, d = d0, is given by the solution of f(d) + fs = 0. If the external macroscopic force is cyclic, as in standard mechanical analysis, fs follows the imposed periodicity. The resulting path in terms of d0 is shown by the red arrows in Figure 2. Notice in particular the horizontal sections of the path, which correspond to spontaneous “jumps” of d0 due to the loop in f(d). It is important to note that the attendant aggregate displacements are quite small, i.e., on the order of several angstroms. If plotted in the fs−d0 plane, the red-arrow path produces a hysteresis corresponding to the energy dissipated during the opening and closing of the contact (jump-out−jump-in), Wloss. With certain additional assumptions, relating fs to the externally applied macroscopic force and the morphology of the filler network, we have shown how storage and loss moduli of rubber nanocomposites can be estimated based on the molecular interactions in the contact (cf. ref 13). In this reference we obtain f pp(d) from molecular force field simulations, whereas f matrix(d) is approximated by simply applying Hooke’s law to the rubber collar. Here we want to show that the proposed mechanism can be obtained from a fully atomistic model of the contact without resorting to linear continuum elasticity.

conjunction with the leapfrog integrator (time step: 0.5 fs). The spherical SiO2 particles, with a diameter of 4.3 nm each, are cut out from β-cristobalite. Truncated valences are saturated with OH groups. Even though industrial silica consists of amorphous particles, there is some justification for the use of β-cristobalite instead. The number density of silanol groups on the surface of our model particle is about 6−8 per nm2, which is in accord with experimental values quoted in the literature.28,29 In addition, there are different types of silanol groups, i.e., geminal, vicinal, and isolated (some authors even distinguish additional types like internal or clustered OH groups and others). In our case the attendant numbers of the three main types are as follows (vicinal: 2.5 nm−2; geminal: 1.4 nm−2; triple: 0.6 nm−2). Experimental information on the relative number of different silanol group types can be obtained from the structure of the adsorption site energy distribution obtained via inverse gas−solid chromatography (e.g., ref 30) or from IR and NMR studies.31,32 However, in practice, these number do strongly depend on manufacturing conditions. Therefore, we have not attempted to study the influence of a variable distribution of the different silanol groups. In our simulations the OH groups may rotate with respect to the Si−O bond. Otherwise, the silica spheres are considered rigid. A number of preliminary simulations are carried out with alkylsilanes distributed on the surface of the spheres (for details of the binding see ref 13). In Figure 3 the alkyl moieties are shown in blue. As explained in detail below, one of the silica spheres can be displaced along the line through the centers of the two particles of Figure 3, keeping the other sphere stationary. During the displacement the force between the two particles as a function of their separation, d, is monitored. Molecular Dynamics Force Field. All MD simulations were carried out using version 4.5.6 of the Gromacs package.33,34 We employ a simple united atom force field (here called IsoSil2015). IsoSil2015 is based on the universal force field35−37 with the following exceptions. The LennardJones (LJ) parameters of (poly-)isoprene are taken from Harmandaris and Doxastakis38 and adjusted for our constant pressure simulations. The adjustments, based on the cohesion energy and the density, consist in the scaling of all σ parameters by a factor 0.96 and all ε parameters by 1.96. Some torsion parameters values are refined using Møller−Plesset quantum calculations applied to selected small molecules. The IsoSil2015 parameter values are compiled in Tables 1−3. The attendant

3. SIMULATION METHODOLOGY Model System Overview. Figure 3 shows a cutaway view of the central section of the simulated system. Two spherical SiO2 particles are embedded in a polyisoprene matrix consisting of 99 monodisperse chains of 170 monomers each. The size of the entire simulation box is 13.29 nm × 13.29 nm × 12.65 nm. Notice that the concept of a spatially well-defined “rubber collar”, introduced as an ad hoc approximation in the simpler original version of the contact model,12,13 based on the assumption that the load-bearing branches of the filler network are composed of contacts involving both the filler aggregates and only the polymer matrix between neighboring aggregates, is now abandoned. In addition, care is taken to include sufficient polymer in order to avoid discernible finite size effects on the results presented here. During our molecular dynamics (MD) simulations the bulk density of polyisoprene is close to 920 kg/m3, the temperature is 300 K, and the pressure is 1 bar. We employ the Berendsen thermostat (τT = 0.1 ps) and barostat (τP = 5 ps)27 in C

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

potential energy group contributions encompass LJ, Coulomb, and simple diagonal valence force field potentials. LJ interactions are calculated via

Table 1. Force Field Parameters for Silica (Si, O, H), Spring (SPR), and Polyisoprene (CA,CD,CE,CH3; cf. Figure 5) σ [nm]

Lennard-Jones Si O H SPR CA CD CE CH3 bonds O−H Si−O Si−SPR CA−CD CD−CE CE−CA CH3−CD CDCD valence angles

b0,ij

0.383 0.312 0.175 0.000 0.385 0.325 0.385 0.385 [nm]

0.110 0.150 0.300 0.1515 0.1500 0.1550 0.1510 0.1338 θ0,ijk [rad]

Si−O−H O−Si−O Si−O−Si Si−Si−SPR O−Si−SPR CA−CDCD CDCD−CE CD−CE−CA CH3−CDCD CE−CA−CD CA−CD−CH3 improper dihedrals

110.000 109.500 180.000 1 1 125.900 125.900 111.650 125.900 111.650 109.000 ζ0,ijkl [rad]

CA−CDCD−CE CH3−CDCD−CE

0 180

ε [kJ/mol]

⎛⎛ ⎞12 ⎛ ⎞6 ⎞ σij σij uLJ(rij) = 4εij⎜⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎟ ⎜⎝ r ⎠ ⎝ rij ⎠ ⎟⎠ ⎝ ij

1.682 0.251 0.770 0.000 0.770 0.823 0.770 1.864 kij [kJ/(mol nm2)]

Cross parameters are obtained using the combination rules σij = (σii + σjj)/2 and εij = (εiiεjj)1/2. The LJ cutoff is 1.2 nm. No LJ long-range corrections are included. Coulomb interactions, calculated using the particle mesh Ewald summation method, are based on atomic partial charges. For the silica particles these charges are calculated using the charge equilibration method and the attendant parameters as described in ref 39. Details of the resulting charge values can be found in Figure 4. Alkylsilane partial charges are shown in Figure 5. We note that (united) atoms within the same molecule do not interact via nonbonded interactions if they are separated by two bonds or fewer. In addition, 1−4 nonbonded interactions are scaled by a factor 1/2. We note also that in order to reduce the computational effort our silica particles are hollow. The thickness of the remaining shell is 0.8 nm in accord with the radial charge distribution in Figure 4. Notice that the integral charge from the neglected atoms is zero. In addition, the dipole moment of the neglected spherical region is small and the overall contribution to the interparticle interaction can be neglected. There is however a dipole layer close to the surface with a width of about 0.5 nm, which is taken into account due to the above thickness of the shell. This was confirmed by separate work.40 We note that this work also shows that in terms of physisorption of small molecules spherical particles cut from crystal structures (saturating cut bonds by hydroxyl groups), the method for particle generation used here, indeed represent the in reality amorphous filler particles quite well. Notice also that the LJ cutoff is slightly larger than the shell thickness. However, we have compared selected results obtained using full particles with the present hollow ones, without finding a discernible difference. The valence force field encompasses bond potentials

5002400 1000000 5000000 278742 278742 297470 297470 297470 kijk [kJ/(mol rad2)]

kijkl

30000 30000 30000 1 1 374.000 374.000 481.160 374.000 481.160 250.000 [kJ/(mol rad2)] 160 160

Table 2. Force Field Parameters for Alkylsilanes (CH1, CH2, CH3silane, Sisilane, Osilane) Lennard-Jones CH1 CH2 CH3silane Sisilane Osilane bonds Sisilane−Osilane Sisilane−CH2 CH2−CH2 CH2−CH1 CH2−CH3silane Osilane−CH2 valence angles Osilane−Sisilane−CH2 Osilane−CH2−CH3silane Sisilane−CH2−CH2 CH2−CH2−CH2 CH2−CH2−CH3silane CH2−CH2−CH1 CH2−CH1−CH2

σ [nm] 0.400 0.400 0.400 0.383 0.312 b0ij [nm] 0.1657 0.1830 0.1506 0.1539 0.1524 0.1421 θ0ijk [rad] 112.900 110.000 115.600 111.700 111.700 111.700 111.700

(2)

ε [kJ/mol]

1 bond kij (rij − b0, ij)2 2 valence angle potentials ubond(rij) =

0.330 0.330 0.330 1.682 0.251 kij [kJ/(mol nm2)]

(3)

1 angle kijk (θijk − θ0, ijk)2 2 improper dihedral potentials 1 improper u improper(ζijkl) = kijkl (ζijkl − ζ0, ijkl)2 2 uangle(rijk) =

326401 201431 297470 278742 287054 425967 kijk [kJ/(mol rad2)]

(4)

(5)

and proper dihedral potentials 5

u proper(Θijkl) =

∑ Cn(cos(Θijkl − π ))n n=0

315.843 587.050 324.023 415.337 415.337 415.337 415.337

(6)

(trans state: Θijkl = 0). Polymerization Algorithm. Our polymerization algorithm is based on the radical like polymerization in ref 41. The authors apply their method to coarse-grained models, whereas here we use it to generate atomistic polymer melts. The polymerization algorithm consists of the steps illustrated in Figure 6. Initially, isoprene monomers are placed on a cubic grid (cf. Figure 6a). Catcher monomers are chosen randomly. D

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Table 3. Proper Dihedral Parameters for Silica (Si, O, H), Spring (SPR), Polyisoprene (CA, CD, CE, CH3; Cf. Figure 5), and Alkylsilanes (CH1, CH2, CH3silane, Sisilane, Osilane) proper dihedrals

C0

C1

C2

C3

C4

C5

Si−O−Si−O Si−O−Si−SPR H−O−Si−SPR CDCD−CE−CA CE−CA−CDCD CD−CE−CA−CD CH2−CH2−CH2−CH2 CH2−CH2−CH2−CH3silane CH2−CH2−CH2−CH1 CH2−CH2−CH1−CH2 Sisilane−Osilane−CH2−CH3silane CH2−Sisilane−Osilane−CH2 Osilane−Sisilane−CH2−CH2 Sisilane−CH2−CH2−CH2

0.100 0.100 0.100 4.603 4.603 −14.036 −0.869 −0.869 −0.869 −0.869 0 0 0 0

0.100 0.100 0.100 1.501 1.501 −20.675 33.868 33.868 33.868 33.868 −52.970 16.990 12.278 27.96

0.100 0.100 0.100 2.835 2.835 1.409 26.52 26.52 26.52 26.52 36.030 12.633 −3.211 21.096

0.100 0.100 0.100 6.216 6.216 23.380 −62.584 −62.584 −62.584 −62.584 63.752 −28.406 −13.992 −64.220

0.100 0.100 0.100 −2.665 −2.665 1.189 3.160 3.160 3.160 3.160 24.201 −2.941 3.543 9.330

0.100 0.100 0.100 −12.742 −12.742 8.334 −0.679 −0.679 −0.679 −0.679 −87.720 −2.363 −3.602 2.223

Figure 4. Radial partial charge distribution in a spherical SiO2 particle as obtained by the charge equilibration method [(a) red: oxygen; (b) black: hydrogen; (c) yellow: silicon]. The horizontal lines indicate the actual respective (lower) charge values used in our simulation. Also shown is the actual SiO2 particle using corresponding colors for the different atom types.

Figure 5. Alkylsilane partial charges (left: AS203; middle: AS208). Right: isoprene monomer notation and geometry.

Their number determines the number of polymer chains in the system. This system is then equilibrated, as shown in Figure 7 (time segment A in the upper panel), via NPT-MD; i.e., the particle number, pressure, and temperature are constant. During the subsequent polymerization each catcher monomer can bind to an unbound monomer, when the monomer is the nearest within a radius of 1.5 nm (cf. panels b and c in Figure 6). The newly acquired monomer becomes the new catcher monomer of the forming chain after a relaxation period of 1 ps. In addition, an extended NPT-MD relaxation of 10 ps is carried out following every 20th catching step. Care is taken to only generate cis-1,4-polyisoprene. This procedure is repeated until the desired polymer distribution is achieved (cf. panel d in Figure 6). Figure 7 shows the development of the system’s density during polymerization (stage B) in a system without filler particles. When filler particles are present, the only difference is that we adjust the density not via the barostat, but via a compression rate adapted to the polymerization rate and

the reduction of occupied volume per newly formed chemical bond. Notice that the density in Figure 7 rises because the bond formation reduces the volume of a bound monomer as compared to a unbound monomer. We note that the efficient formation of chains is favored by the rather large “catching radius” of 1.5 nm in conjunction with the relatively short relaxation time during this phase. This results in the anomalous increase of the density in stage B, which subsequently is relaxed in the last stage of the algorithm. In the last preparation step, the remaining unbound monomers are deleted and the remaining system is equilibrated (cf. stage C in Figure 7). We remark that an extension of the equilibration time to 100 ns did not change the final result. Notice that the densities of the monomer fluid and the final melt of chains are quite close to their experimental values. In addition, we can calculate the characteristic ratio, cn, averaged over the finished polymer chains. Again the result is close to the experiment. It is worth noting that the algorithm also produces E

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Measuring the Force between the Two Embedded Silica Particles. In order to measure the total force between the two silica particles, including enthalpic as well as entropic contributions of the surrounding chains, we employ the concept of a force gauge. The force gauge is realized via a spring between a virtual atom and the right silica particle in Figure 8. Virtual means that the atom is not interacting in any

Figure 8. Principle of the force gauge. The left silica particle is fixed, while the right silica particle moves freely along the line through the centers of the two particles as defined by their initial positions. A harmonic force (spring) is introduced between the right silica particle and a virtual atom shown here as a gray circle. The position of the virtual atom along the aforementioned line and thus the harmonic force due to the spring exerted on the right particle are adjustable.

Figure 6. (a) Isoprene molecules distributed on a cubic grid. Random “catcher” monomers (here shown in red) are selected. (b) A bond is formed between the catcher monomer and another monomer within a radius of 1.5 nm. (c) The newly bonded monomer becomes the catcher monomer for this chain. (d) Repetition of steps b and c until the chains have the required length.

other way with the system. The position of the virtual atom on the line through the centers of the two silica particles in Figure 8 may be changed by an amount Δs giving rise to an attendant force, Δf = kΔs, exerted on the right silica particle due to the spring. This in turn causes the equilibrium distance between the silica particles to change by the amount Δd. By altering s, the position of the virtual atom, we are able to map out the total interaction force between the two particle as a function of their distance. Routinely we use a spring constant, k ≈ 5 × 109 J mol−1nm−2. The virtual atom displacement is roughly Δs = 10−3 nm per 20 ps. These numbers are found by balancing the gauge’s sensitivity and the need for computational speed. The force gauge was checked by comparison of its result to the interparticle force obtained without polymer, which in this case can simply be calculated via the negative gradient of the interaction potential.

4. RESULTS First we conduct a simple test of the effect of polymerization, i.e., the elastic restoring force, on the equilibrium interparticle separation (cf. Figure 9). The two silica particles are initially 0.24 nm apart. Keeping this distance fixed, the polymerization is carried out as described above. Then the interparticle separation is increased with a velocity of 5 cm/s = 0.05 nm/ ns to a certain value. Subsequently, the particles are free to adjust to their equilibrium separation. The black curves in Figure 9 show the time development of the interparticle separation if the spheres are started from two quite different separations (≈ 0.37 nm vs ≈0.55 nm). In both cases the final distance between the spheres is close to 0.47 nm, which, as the next figure shows, correlates well with the zero of the interparticle force curve. In contrast, the red curve shows the result for the larger starting separation if the surrounding medium is nonpolymerized isoprene. There is no appreciable relaxation of the interparticle separation in this case as expected. Figure 10 shows the results obtained during three successive force vs distance cycles for two different sets of partial charge values (cf. below). The displacement velocity is 5 cm/s throughout. Initially the separation is decreased, which results

Figure 7. Top: mass density vs simulation time during the stages A, B, and C of the polymerization. The dashed horizontal lines indicate the experimental densities of the monomers, ρIsoprene, and the polymer, ρPolyisoprene, respectively. Bottom: characteristic ratio, cn, vs the number of monomers, n, between which the mean square end-to-end distance is calculated, of the polyisoprene chains at the end of stage C. The horizontal dashed line is the experimental result from ref 42.

the Flory exponent ν = 1/2 for polymer melts. If for instance ν would be larger than 1/2, then the characteristic ratio cn would continue to rise with increasing n. F

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

scatter, even though a hysteresis remains, distinguishing the force curve during distance reduction from the somewhat lower lying force curve during distance increase. The type of conformation changes, causing the hysteresis are currently unknown, but this will be a subject of future research. In this context we also note that the right particle is moved directly. However, inside a real filler network the force is transmitted by neighboring particles and their direct contacts including associated rubber collars and/or bridges. In our case the direct application of force to the particle causes compression/tension of the rubber matrix on the opposite side of the contact. It is not quite clear how large this contribution exactly is. However, the much larger particle-to-particle separation due to the periodic boundary conditions reduces the strain of the attendant polymer considerably; i.e., the main polymer contribution to the total force should indeed be due to polymer inside the contact. The bottom panel in Figure 10 shows the same result for increased partial charges. Notice that these increased charge values for Si and O are close to the partial charges used in bulk simulations of amorphous silica, i.e., qSi = 2.4 and qO = −1.2.43 Here, despite the hysteresis, we do obtain evidence for a force vs distance curve as assumed previously (cf. Figure 2). Averaging over all paths, with the exception of the initial distance decrease and subsequent increase to full maximum, we obtain the curves shown in Figure 11. This figure again confirms the loss mechanism suggested in previous work in refs 12 and 13 but this time using a fully atomistic model of the contact between two filler particles. The

Figure 9. Interparticle separation equilibration as a function of simulation time. (a) Red curve: the two particles are embedded in nonpolymerized isoprene; black curves: the two particles are embedded in polymerized isoprene. Here (b) and (c) indicate two different starting distances.

Figure 10. Top: total particle−particle force vs interparticle distance for three successive cycles. Black circle: starting distance. Blue curve: decrease (0), increase (1), and again decrease (1) of the interparticle separation. Green (2) and red (3) curves: subsequent cycles. Black line: result in the absence of polymer. Bottom: same as above but with all partial charges increased by a factor of 2.

in an increasingly repulsive interparticle force. Subsequently the displacement direction is reversed until a maximum distance between 0.6 and 0.75 nm is attained. Along this path the force has become attractive during a monotonous decrease. A new reversal of the direction yields a different force curve. However, following this initial reduction of the particle separation from its maximum value all subsequent curves coincide within their

Figure 11. Top: averages over all particle−particle force vs interparticle distance paths in the previous figure. The initial decrease and subsequent increase to full maximum are omitted. (a) Black: normal partial charges; (b) red: increased partial charges. Bottom: hysteresis cycles calculated on the basis of the two curves (a) and (b) in the upper panel. The area enclosed by the hysteresis curves is Wloss (cf. Figure 2). G

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules extend of the energy dissipated certainly depends on the interaction, i.e., the partial charges in the present case. Currently simulation work is under way,40 which is aimed at the refinement of, in particular, the Coulomb interactions on the basis of experiments measuring the site energy distribution of a number of small alkenes on silica using inverse gas chromatography.30 Our current choice of a factor of 2 between the two sets of charges used here is motivated by the present considerable uncertainty regarding these important parameters. If, despite this uncertainty, we calculate Wloss, as explained in the section following the introduction (also cf. Figure 2), we obtain Wloss < 36 kJ/mol and, in the case of increased partial charges, Wloss ≈ 3900 kJ/mol. This is the energy dissipated during one strain cycle due to the opening and closing of one interparticle contact. Wloss will depend on the hysteresis exhibited by the force curves in Figure 10. The hysteresis should be rate dependent. But our current computational means do not permit a significant extension of the overall simulation time to explore this point. Independent of the particle−particle interaction at the small separations, we can analyze the force vs distance curves in Figure 11 on the opposite side, i.e., at large separations, to deduce information on the shear modulus of the matrix surrounding the particles. Hunter44 has calculated that the force, f, needed to displace a macroscopic sphere of radius R embedded in a rubber matrix with shear modulus μ by a small distance Δ is given by f = 6πμRΔ. We may apply this formula, albeit roughly, assuming that the (negative) slope of our force vs distance curves in the limit of large separations (here: >0.5 nm) yields f/Δ. Using our sphere radius of 2.15 nm, we find μ ≈ 500 MPa. We must bear in mind that the molecular simulations take place on the nanoseconds time scale. In terms of frequency this corresponds to the gigahertz range. Experimental information in this range is obtained using the temperature−time superposition principle. A corresponding curve for the storage modulus of natural rubber is shown in Figure 2.72 of ref 5. Comparison with this measurement in the gigahertz range shows that the above simulation value for μ is realistic. The short length of the simulated chains in comparison to the polymers in natural rubber is compensated by the high frequency, which tends to reduce the effects of modes governed by long-range correlations along the chains. We note that the temperature−time superposition principle may be a possible means to map simulation results, like those presented here, into more relevant frequency ranges. We note also that Hudzinskyy et al.18,19 have studied united-atom atactic polystyrene films of chains, whose lengths are comparable to ours, sheared between Lennard-Jones surfaces. Their shear rates are also in the nanoseconds range, and the resulting shear modulus of the polymer is between 500 and 800 MPa, in good accord with our result above. Figure 12 summarizes results regarding the matrix structure and segment dynamics in close proximity of the particles. The upper panel shows the radial density of polymer (averaged over both particles excluding the contact region), where r is measured from the center of each sphere. There is a rather pronounced first layer. But already the second and third layers are much less pronounced. This observation is independent from the separation of the two particles because the contribution of the immediate contact area between the spheres is comparatively small. The bottom panel provides dynamical information. The residence function N(t,θ) is the average number of remaining polymer atoms inside a θ-

Figure 12. Top: radial polymer density (black: interparticle distance is 0.158 nm; red: distance is 0.749 nm). The vertical dashed line indicates the silica surface. Bottom: normalized number of remaining polymer atoms inside circular segments of the first two density shells (solid lines: first shell [2.15 nm < r < 2.65 nm]; dashed lines: second shell [2.65 nm < r < 3.15 nm]) as a function of time and cone angle (blue: θ = 55°; black: 35°; red: 15°; Δθ = 10°). The inset explains the meaning of the cone angle θ. Notice that Δθ is the angular width of the circular surface segments on both particles above which N(t,θ) was determined.

dependent circular segment of either the first or second density layer, as defined in the figure caption, at time t. Notice that a polymer atom, leaving its segment at time t for any extend of time, reduces N(t,θ) by one for all future times. Notice also that large θ means the corresponding circular shell segment samples atoms at a greater distance from the contact. All curves can be fitted rather well with the sum of two exponentials, each with its characteristic residence time. The short residence time in all cases is several tens of picoseconds, whereas the long residence time is several thousand picoseconds. The short residence times can be explained as being due to independent atomic motion, whereas the long residence times are due to the center-of-mass motion of the polymer chain. This interpretation is supported by simulations of polystyrene chains of comparable length near silica particles.45,46 The author determines the chain center-of-mass diffusion near the particle surface. He finds that root-meansquare displacements of a chain’s center of mass, comparable to the present layer widths, require several thousand picoseconds. During the time period governed by the short relaxation time, we do not observe easily distinguishable differences between the two shells or between the different segments inside the same shell. The fact that the normalized number of remaining atoms in a first-shell segment becomes less than the corresponding number in the second shell is due to the significantly smaller spatial size of the segments in the first layer in comparison to their direct neighbors in the second layer. It is H

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

similar to what we find in the case of the “naked” particle. At small distance there appears an additional feature, which is due to a silane molecule obstructing the particle−particle contact. AS208, possessing a longer alkyl chain, has a more pronounced effect. Whereas before we observe the development of a stable cycle after the initial increase of the distance to its maximum, here no such regularity is attained. Inspection reveals that a polymer segment obstructs the particle−particle contact in multiple conformations. Mostly this will cause the force to rise, which prevents the force loop to occur. In turn, this will reduce the dissipative loss caused by the mechanism discussed above. This relation between alkyl chain length and loss, in terms of a reduced loss modulus, is consistent with the experiment (e.g., ref 4). However, Figure 13 is a preliminary result. It is clear that the statistical sampling must be improved with respect to factors like the placement of the silanes on the surface or the number of cycles performed.

worth noting that atoms can leave their segment not only vertically to the particle surface but also through the larger segment surfaces oriented perpendicular to the particle surface. Comparison of N(t,θ) for different θ within the same shell is more straightforward. The data in Figure 12 show that N(t,θ)/ N(t,0) is constantly less for the segment closest to the contact in comparison to the segment farthest from the contact. Notice also that the average long residence time governing N(t,θ) for times greater than 0.5 ns is less in the second shell compared to the first shell. This is consistent with observations of other studies, which investigate the polymer mobility as a function of the distance from filler particles45−47 and find reduced mobility for polymers close to a particle surface. It is worth noting that the presence of fillerappears to affect not only the polymer mobility but also the distribution of entanglements; i.e., increase of strain appears to attract entaglements toward the particle surfaces.48 In practice, silica particles are coated with compatibilizing molecules reducing the surface tension of the particle−matrix interface in order to achieve better dispersion (e.g., ref 4). In our previous work13 we have modeled the concentration dependent effect of alkylsilane compatibilizers on the storage and the loss modulus in silica filled rubbers. Here we apply a random low density coating using a silane with short alkyl moiety (AS203) and one with a longer alkyl moiety (AS208) as specified in ref 13. Figure 13 shows the attendant force vs distance curves. The result obtained with AS203 is rather

5. DISCUSSION AND CONCLUSION We have proposed that the basic unit, responsible for the Payne effect, in a filler network is a contact consisting of two aggregates possessing at least one immediate contact between silica surfaces, which under external stress can be opened and closed. In the section following the introduction we do explain the underlying dissipative mechanism. The ability to describe this mechanism on a molecular scale allows to influence its effects by chemical modification, in this case, of the particle surface. Previously we had treated the polymer in our model of an aggregate-to-aggregate contact as a Hookean spring. In this work we consider a fully atomistic model of a particle-toparticle contact instead. We are able to show that a force loop, which is the key feature of the above dissipative mechanism, can arise depending on the interaction, in particular the Coulomb interaction, between the silica particles. Several critical remarks should be made at this point. It is necessary to improve the description of, in particular, Coulomb interactions between the particles. This requires comparison with experimental data. Probably the best direct comparison can be made by measuring the site adsorption energy distribution of small relevant molecules using for instance inverse gas chromatography30 and industrial silica. Relevant means that the molecules closely represent the monomers in the polymers. In the case of polyisoprene an example is 2,3dimethylbutene. Another relevant molecule is water because the interaction between naked silica particles involves formation and breakage of hydrogen bonds similar to those formed between water and silica. A side benefit from this system may be a better understanding of the role of silica in the improvement of the wet grip of silica filled rubbers as compared to carbon black filled rubbers. Simultaneously this would improve the modeling of the morphology of silica nanoparticles. Even though force fields for amorphous silicon dioxide do exist (e.g., refs 43 and 50), these force fields are parametrized based on bulk properties only. It also is important to improve the description of the silanol groups on the particles surfaces, which thus far are mostly ignored in simulation studies focusing on the interactions between silica nanospheres or between silica nanospheres and their surrounding matrix. The loss mechanism discussed in the section following the introduction is based on an aggregate−aggregate contact rather than the mere contact between primary particles. Realistic sizes of the latter, on the low end of the size distribution, certainly overlap with the size of the particles used in this study.

Figure 13. Total particle−particle force vs interparticle distance in analogy to the upper panel in Figure 10. However, here the particle surfaces are “dressed” with alkylsilanes, which are not present in the case of Figure 10. Top: two cycles for particles containing 0.192 AS03 silanes per nm2. The black circle again marks the starting position start position ((a) solid lines: increasing interparticle separation; (b) dashed lines: reverse direction). Bottom: same as above but for AS08 instead of AS03. I

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

(20) Chen, Y.; Liu, L.; Yang, Q.; Wen, S.; Zhang, L.; Zhong, C. Langmuir 2013, 29, 13932. (21) Masnada, E.; Merabia, S.; Couty, M.; Barrat, J.-L. Soft Matter 2013, 9, 10532. (22) Brown, D.; Marcadon, V.; Mélé, P.; Albérola, N. D. Macromolecules 2008, 41, 1499. (23) Böhm, J. Der Payne-Effekt: Interpretation und Anwendung in einem neuen Materialgesetz für Elastomere Ph.D. Thesis, Universität Regensburg, Regensburg, Germany, 2001. (24) Gurovich, D.; Macosko, C. W.; Tirrell, M. Rubber Chem. Technol. 2004, 77, 1. (25) Xi, H.; Hentschke, R. Eur. Polym. J. 2012, 48, 1777. (26) Meyer, J. Ein coarse-grained Modell für die Simulation dynamisch-mechanischer Eigenschaften gefüllter Elastomernetzwerke. Master Thesis, Wuppertal University, Wuppertal, Germany, 2014. (27) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (28) Krause, M. Untersuchung der Wechselwirkung von Polymer/ Silica-Mischungen mit Festkoerper-NMR. Ph.D. Thesis, Freiburg University, Freiburg, Germany, 2002. (29) Peri, J. B.; Hensley, A. L. J. Phys. Chem. 1968, 72, 2926. (30) Wang, M.-J.; Wolff, S.; Donnet, J.-B. Kautsch. Gummi Kunstst. 1992, 45, 11. (31) Guy, L.; Monton, A.; Jost, P.; Goyard, S.; De Cayeux, S. High Performance Silicas Reinforced Elastomers: A better Understanding of the Interface Nature to Improve Balance Between Hysteresis and Reinforcement. Presented at the Fall 184th Technical Meeting of the Rubber Division of the American Chemical Society, Inc. Cleveland, OH, October 8−10, High Performance Silicas Reinforced Elastomers: A better Understanding of the Interface Nature to Improve Balance Between Hysteresis and Reinforcement2013. (32) Daudey, S.; Guy, L. Silica systems to modulate the compromise reinforcement/hysteresis for the silica filled elastomers; RIEG-IOM, London, UK, 2011. (33) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306. (34) van der Spoel, D.; Lindahl, E.; Hess, B. Gromacs User Manual version 4.5.4; url: www.gromacs.org, 2010. (35) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. (36) Casewit, C. J.; Colwell, K. S.; Rappe, A. K. J. Am. Chem. Soc. 1992, 114, 10046. (37) Casewit, C. J.; Colwell, K. S.; Rappe, A. K. J. Am. Chem. Soc. 1992, 114, 10035. (38) Harmandaris, V. A.; Doxastakis, M. J. Chem. Phys. 2002, 116, 436. (39) Rappe, A. K.; Goddard, W. A. J. Phys. Chem. 1991, 95, 3358. (40) Ballnus, C. Adsorption kleiner polarer und unpolarer Moleküle auf Silica-Nanopartikeln. Master Thesis, Wuppertal University, 2015. (41) Perez, M.; Lame, O.; Leonforte, F.; Barrat, J.-L. J. Chem. Phys. 2008, 128, 234904. (42) Mark, J. E. J. Am. Chem. Soc. 1967, 89, 6829. (43) van Beest, B.; Kramer, G.; van Santen, R. Phys. Rev. Lett. 1990, 64, 1955. (44) Hunter, S. C. Proc. Eding. Math. Soc. 1968, 16, 55. (45) Ndoro, T. V. M Simulation of a Polystyrene Silica Nanocomposite. Ph.D. Thesis, Darmstadt University, Darmstadt, 2011. (46) Ndoro, T. V. M; Böhm, M. C.; Müller-Plathe, F. Macromolecules 2012, 45, 171. (47) Guseva, D. V.; Komarov, P. V.; Lyulin, A. V. J. Chem. Phys. 2014, 140, 114903. (48) Riggleman, R. A.; Toepperwein, G.; Papakonstantopoulos, G. J.; Barrat, J.-L.; de Pablo, J. J. J. Chem. Phys. 2009, 130, 244903. (49) Baeza, G. P.; Genix, A.-C.; Degrandcourt, C.; Petitjean, L.; Gummel, J.; Couty, M.; Oberdisse, J. Macromolecules 2013, 46, 317. (50) Horbach, J.; Stühn, T.; Mischler, C.; Kob, W.; Binder, K. In High Performance Computing in Science and Engineering ’03; Springer: Heidelberg, 2003; p 167.

However, aggregates are composed on several ten primary particles.49 Aggregates are “unbreakable” and act as one large complex particle. Thus, the contact of this study still is a small and simplified model of the contacts inside the filler network, which we think are responsible for the Payne effect. One difference to this model is the much larger amount of rubber forming the “collar” linking the aggregates together (cf. the middle panel in Figure 1). This in turn considerably increases the elastic force, f matrix(d). Here this effect is offset by the high frequency of the cyclic particle displacement, which increases matrix modulus. A possible approach in future simulation studies of this type may be the explicit treatment of the contact on the size scale of the contact between two primary particles, as in this work, combined with a continuum description beyond (an old approach originally introduced in simulations of protein interaction (e.g., ref 51)).

■

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected]; Ph +49 202 439 2628; Fax +49 202 439 3122 (R.H.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Jan Meyer for providing them with the top panel in Figure 1.

■

REFERENCES

(1) Tjong, S. C., Mai, Y.-W., Eds.; Physical Properties and Applications of Polymer Nanocomposites; Woodhead Publishing: Oxford, 2010. (2) Thomas, S., Stephen, R., Eds.; Rubber Nanocomposites: Preparation, Properties, and Applications; John Wiley & Sons: Singapore, 2010. (3) Kraus, G., Ed.; Reinforcement of Elastomers; John Wiley & Sons: New York, 1965. (4) Luginsland, H.-D. A Review on the Chemistry and the Reinforcement of the Silica-Silane Filler System for Rubber Applications; Shaker Verlag: Aachen, 2002. (5) Wrana, C. Introduction to Polymer Physics; LANXESS AG: Leverkusen, 2009. (6) Vilgis, T. A.; Heinrich, G.; Klüppel, M. Reinforcement of Polymer Nano-Composites; Cambridge University Press: New York, 2009. (7) Gil-Negrete, N., Alonso, A., Eds.; Constitutive Models for Rubber VIII; CRC Press/Balkema: Leiden, 2013. (8) Payne, A. R. In Reinforcement of Elastomers, Kraus, G., Ed.; John Wiley & Sons: New York, 1965; p 69. (9) Ouyang, G. B. Kautsch. Gummi Kunstst. 2006, 59, 332. (10) Ouyang, G. B. Kautsch. Gummi Kunstst. 2006, 59, 454. (11) Stambaugh, R. B. Ind. Eng. Chem. 1942, 34, 1358. (12) Hentschke, R. In Constitutive Models for Rubber VIII; GilNegrete, N., Alonso, A., Eds.; CRC Press/Balkema: Leiden, 2013; p 299. (13) Hentschke, R.; Hager, J.; Hojdis, N. W. J. Appl. Polym. Sci. 2014, 131, 40806. (14) Gusev, A. A. Macromolecules 2006, 39, 5960. (15) Merabia, S.; Sotta, P.; Long, D. R. Macromolecules 2008, 41, 8252. (16) Merabia, S.; Sotta, P.; Long, D. R. J. Polym. Sci., Part B: Polym. Phys. 2010, 48, 1495. (17) Raos, G.; Casalegno, M. J. Chem. Phys. 2011, 134, 054902. (18) Hudzinskyy, D.; Michels, M. A. J.; Lyulin, A. V. J. Chem. Phys. 2012, 137, 124902. (19) Hudzinskyy, D.; Michels, M. A. J.; Lyulin, A. V. Macromol. Theory Simul. 2013, 22, 71. J

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules (51) Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and Solutions; John Wiley & Sons: New York, 1991.

K

DOI: 10.1021/acs.macromol.5b01864 Macromolecules XXXX, XXX, XXX−XXX