How Thick Is the Polymer Interphase in Nanocomposites? - American

Article pubs.acs.org/Macromolecules

How Thick Is the Polymer Interphase in Nanocomposites? Probing It by Local Stress Anisotropy and Gas Solubility Evangelos Voyiatzis,*,† Mohammad Rahimi,†,‡ Florian Müller-Plathe,† and Michael C. Böhm† †

Eduard-Zintl-Institut für Anorganische und Physikalische Chemie and Center of Smart Interfaces, Technische Universität Darmstadt, Alarich-Weiss-Strasse 4, D-64287 Darmstadt, Germany ‡ Institute for Molecular Engineering, University of Chicago, Chicago, Illinois 60637, United States S Supporting Information *

ABSTRACT: The influence of the inclusion of a silica nanoparticle on the spatial distribution of the local stresses and the locally resolved excess Helmholtz free energy of sorption of small penetrants (He, H2, O2, and CO2) in a polystyrene matrix is studied by molecular dynamics simulations. The local deviations of these quantities from their bulk averages are correlated with spatial peculiarities in structural and dynamical properties, for instance in the mass density, the segmental orientation, and the mean-square atomic positional fluctuations in the polymer matrix. Relative to the bulk, stress anisotropies are slightly enhanced in the neighborhood of a nanoparticle. This is demonstrated by the spatial variations of the local shear stress and the local von Mises shear stress. For all penetrants considered, the region near the interphase is found to be a preferential sorption site. The two stress estimators are compared against one of the most frequently adopted descriptors, i.e., the radial mass density distribution. We show that the estimated interphase width depends on the quantity considered. The interphase width based on variations of the local stresses is considerably shorter than the estimate obtained using the mass density distribution. Furthermore, we find that the interphase width derived from the locally resolved Helmholtz free energy strongly depends on the size of the inserted molecule. For the smallest particle, the helium atom, a broader interphase is found than for the larger molecular species. In the case of carbon dioxide insertion, we estimate an extension similar to the one derived by stress profiles. By artificially reducing the interactions between the nanoparticle atoms and those of the polymer matrix, an attempt has been made to identify links between different ways to define the interphase thickness. It is shown that the quantities under consideration lead to interphase widths which are independent of each other. mechanical properties of nanocomposites.17−20 Only few of these studies have attempted to detect correlations between the spatial variations of mechanical, thermodynamic, and structural quantities. Such investigations are the prerequisite to explain the microscopic origin of the macroscopically observed peculiarities in nanocomposites. Despite the large number of articles on nanocomposites,1−7,9,11 their description via local or atomic level stress quantities has not been in the focus of recent research activities. Atomic level stresses are defined in an entirely atomic basis without resorting to any theory of continuum mechanics (see section 3.2). Suitable formalisms have been developed several decades ago.21−23 They have found widespread applications in the description of hard matter systems where only pairwise central forces are sufficient to describe the interactions among the different atoms.25,24 The extension of atomic level stress concepts to bonded systems with noncentral forces26,27 has allowed applications to soft matter.23,28 The atomic level stresses characterize the local structure of a material25 and

1. INTRODUCTION As a result of the enormous surface area between a polymer matrix and filler particles in nanocomposites, mechanical and rheological properties are observed that go beyond those of conventional composites.1 A rational design of nanocomposites is still not always possible, as simple extrapolations of design paradigms valid for conventional composites fail to predict their behavior.2 A variety of interesting phenomena with unknown microscopic origin is observed in these materials. Examples are a dramatic reduction in the viscosity of polymer melts after adding a small amount of nanoparticles,3 modifications of the glass transition temperature, Tg, by up to 30 K,4,5 and a significant improvement of the mechanical properties of polymer nanocomposites.2,6 The properties of polymers reinforced with nanoparticles have been studied extensively in the past years by different molecular modeling techniques.6−10 The influence of the shape and size of nanoparticles as well as of their grafting density on the structure and dynamics of the polymer matrix has been the subject of several of these investigations.11−15 Rheological properties of nanocomposites have been interpreted by extending tools for entanglement analysis.16 A variety of techniques have been developed to derive spatially resolved © XXXX American Chemical Society

Received: March 17, 2014 Revised: October 20, 2014

A

dx.doi.org/10.1021/ma500556q | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

the polymer matrix have been reduced artificially. All other structural details of the system, such as the chain length and the radius of the nanoparticle, are the same as in the initial nanocomposite. By monitoring the changes in the estimated interphase lengths obtained by all quantities of interest, we cast light on their interconnections. We also address the question whether it is possible to predict all the interphase widths on the basis of a single quantity. The present analysis serves as a guide to understand the reported differences in the interphase widths encountered in polymer matrices which are in contact with solid surfaces. Nota bene the property profiles employed here to study the interphase are all expected to converge over a very short-range, in some cases below a monomer diameter. In this, our investigation augments and complements previous theoretical43−52 and experimental studies,53−62 by analyzing properties that normalize over considerably longer length scales, such as the chains’ radii of gyration or the mobilities of entire polymer molecules.

provide information on the response of a sample to small elastic deformations.29 They also cast light on the atomic origin of dynamic properties such as the viscosity and various relaxation times. This stems from the fact that the corresponding Green− Kubo expressions can be decomposed into correlation functions between atomic level stresses.31,30 Stress profiles are furthermore useful when combining particle-based simulations with theories from continuum mechanics.32−34 Spatial variations in the gas solubility, or equivalently in the excess Helmholtz free energy of sorption, ΔAex, in nanocomposites have also not been discussed yet. Such calculations can be performed by so-called test particle insertion methods.35,36 They are efficient for small penetrants in not too dense host materials and in the absence of strong and specific host−guest interactions.37 Alternatively, ΔAex can be calculated by removing an existing particle and correcting the results for the excluded volume removed from the system.37,38 The ΔAex values of an inserted molecule are related to the Henry constant.39 Moreover, ΔAex also probes the free volume of a system. Test particle insertion considers the energetic interactions in the system which are completely ignored by geometrical free volume approaches such as the construction of Voronoi cells.40,41 A link between ΔAex measurements and geometric methods already exists.42 The latter techniques map a molecular species into an equivalent one which consists only of hard spheres of different radii. In such a system, the free energy of inserting a hard-sphere test particle is either equal to zero if there is an overlap with an existing particle sphere (atom) or equal to one in the opposite case. In the case of hard spheres, ΔAex is related to the accessible volume. Similar to a system of hard spheres, where the calculated free volume depends on the radius of the inserted hard sphere, the probed volume depends on the interaction parameters of the test particle. The present work is a comparative study of local stress and locally resolved excess Helmholtz free energies of insertion of small molecules. The effect of nanoparticle inclusion on these two quantities, both on the overall mean values and on their spatial distributions, has not been presented up to now in the literature. New insights into the induced modifications of the polymer matrix in nanocomposites are gained by comparing the findings from the local stress and the local excess Helmholtz free energy of sorption with structural and dynamical features of the matrix as captured by the radial mass density, segment orientation, and distribution of atomic spatial fluctuations. We have chosen a nanocomposite with a single silica nanoparticle in an atactic polystyrene matrix. Our approach goes beyond generic models and considers a material with an interesting chemistry, which is analyzed in full atomistic resolution. Nevertheless, certain idealizations concerning for instance the concentration and the shape of the nanoparticle had to be accepted to make such atomistic simulations possible. The spatial distributions of the studied quantities in both the amorphous glass and melt states are presented as a function of the distance from the nanoparticle surface. The interphase width is estimated as the distance from the nanoparticle surface required by the local stress, the excess Helmholtz free energy, the atomic positional fluctuations, or the mass density distributions to reach their bulk averages. The observed differences in the estimated widths are discussed. In an additional step, we consider a fictitious system in which the interactions between the nanoparticle atoms and the atoms of

2. MODELS AND METHODS We perform molecular dynamics (MD) simulations of pure polystyrene (PS) and a polystyrene−silica nanocomposite (NC). In both systems, the polymer matrix consists of 30 monodisperse atactic PS chains of 200 monomers (molecular weight 20 800 amu) which are terminated by methyl groups. The nanocomposite system contains in addition a single silica nanoparticle (NP) with a radius of 2 nm. The NP volume fraction amounts to approximately 3%. This means that the nanoparticle and its periodic images are isolated from each other, approaching the infinite dilution limit. The surface-tosurface distance between the NP and its own periodic image is approximately 6.2 nm. Simulations of polymer films confined between two flat surfaces,63,64 which influence the interphase structure much more than the convex surface of a spherical particle,11,65 have shown that all local properties become bulklike over a range of 2 nm. As the quantities studied in this article are found to converge over even shorter distances, we trust to perform an analysis of polymer quantities under the sole influence of a single NP. The studied quantities are not influenced by a competition between NP−PS and NP−NP interactions. They might appear for NP distances below 3 nm,15 a topic which is beyond the scope of the present contribution. The systems are studied at atmospheric pressure, i.e., P = 101.3 kPa, and at two temperatures, T = 300 and 450 K, corresponding to glass and melt states. The length of the polystyrene chains and the size of the systems under consideration prohibit direct atomistic MD simulations. To derive atomistic data, a multiscale simulation strategy is adopted. The initial configurations are created at a coarse grained (CG) level followed by long CG MD simulations to allow full equilibration. Each styrene monomer is represented by one coarse-grained bead.66 The SiO2 units of the silica nanoparticle are represented by a second type of CG beads.65 After CG equilibration, the resulting structures are backmapped to fully atomistic resolution by inserting all atoms.67 The atomistic force field for polystyrene is based on the OPLS-AA parametrization. 68 Silica interactions are described by the model of Lopes et al.69 A detailed description of the adopted potentials as well as a complete list of all parameters can be found elsewhere.11 The initial NC configurations are created by placing the nanoparticle at the center of the simulation box. Then CG polymer beads are placed around the nanoparticle by using a B


Macromolecules

Article

observer bias when selecting the data points for the high and low temperature fits. The value of Tg is only slightly affected by NP inclusion. The observed 3 K reduction is within the statistical uncertainty. As a first step to quantify modifications in Tg as a function of the separation from the NP surface, we have computed the specific volume of two spherical shells of 1.5 nm thickness around the NP. Their centers are separated by 0.75 and 2.25 nm from the NP surface. The variation of the specific volumes as a function of temperature is shown in Figure 2. The shells are as thick as

self-avoiding walk algorithm. The CG systems are simulated for 20 ns in the isothermal−isobaric (NPT) ensemble at 101.3 kPa and 590 K. The Berendsen thermostat and barostat70 are employed with coupling times of 0.2 and 5 ps. An isothermal compressibility and a time step equal to 10−6 kPa and 5 fs are chosen. The CG simulations are performed using the IBIsCO code.71 The density of both systems converges within less than 1 ns to its final value (approximately 970 kg/m3 for the NC and 940 kg/m3 for PS). The relaxed CG configurations are backmapped to full atomistic resolution using the method of Ghanbari et al.67 They are subjected to an energy minimization using the Polak− Ribière conjugate-gradient method.72 Further equilibration at the atomistic level is performed for 10 ns in the NPT ensemble at 101.3 kPa and 590 K. The Nosé−Hoover thermostat and barostat70 are employed with coupling times of 0.1 and 1.0 ps. A time step of 1 fs using the velocity-Verlet integration scheme is used. A 1.1 nm atom-based cutoff for the summation of van der Waals interactions is employed. All atomistic simulations are performed with the LAMMPS code.73 The systems are cooled in a double-step procedure below their glass transition temperature, Tg. In each double step, the temperature is initially decreased for a time period of 1 ns with a cooling rate of 10 K ns−1, equivalent to an overall temperature drop of 10 K. The employed cooling rate is in line with that of previous work.74,75 After each temperature reduction, the systems are further relaxed at the final temperature by performing a short MD simulation of 1 ns in the NPT ensemble. In total, we employ this double-step procedure 29 times until the temperature reaches 300 K. The variation of the specific volume of PS and the NC as a function of temperature is shown in Figure 1. The Tg of PS is

Figure 2. Temperature dependence of the specific volume of polystyrene in spherical shells with a width of 1.5 nm separated by 0.75 nm (triangles) and 2.25 nm (rhombs) from the nanoparticle surface to derive locally resolved glass transition temperatures, Tg. The solid and the dashed lines are linear least-squares fits in the high and low temperature regimes.

the silica−polystyrene interphase when estimated by segmental structural properties such as the PS radial mass density distribution (see Figure 3). The polymer atoms lying in the

Figure 1. Temperature dependence of the specific volume of pure polystyrene (squares) and the nanocomposite (circles) to derive the glass transition temperature, Tg. The solid and the dashed lines are linear least-squares fits in the high and low temperature regimes.

Figure 3. Polystyrene mass density as a function of the separation from the nanoparticle surface accumulated in 0.05 nm thick bins for the melt (blue line, T = 450 K) and the glass (red line, T = 300 K). The simulations have been performed at P = 101.3 kPa.

found to be 384 ± 23.9 K, which is in sufficient agreement with the experimental value of 380 K,76,77 obtained by differential thermal analysis of pure polystyrene with a weight-average molecular weight of 374 000 amu. The polydispersity index (weight-average molecular weight divided by the numberaverage molecular weight) of the experimental sample was approximately 2.5.77 The molecular weight is high enough to be close to the limit of an infinite chain length. We interpret the experimental result as an indication that Tg does not depend strongly on the polymer length. The rather large error bar in the Tg evaluation is due to a very conservative estimate of the

second shell are influenced only indirectly by the NP inclusion. This influence originates from the propagation of structural modifications in the interphase area to the bulk region. There are no direct energetic interactions between polymer and NP atoms. We observe that the Tg of the first shell has been shifted to a value of approximately 404 ± 25.2 K. Thus, the volumetrically determined transition from the melt to the glass state occurs at a higher temperature than in pure PS. This value is also higher than the one observed for the whole NC C


Macromolecules

Article

sample with its low NP concentration. The Tg value in the second shell is roughly identical to the one of pure PS. Furthermore, we have computed the variation of the specific volume in much thinner spherical shells with a thickness of 0.15 nm. Despite the higher statistical uncertainties, we can conclude that there is no phase transition in the shells where a maximum in the local density is observed (i.e., at distances from the NP of approximately 0.45 and 0.90 nm; see Figure 3). In these two cases, the specific volume data can be fitted to a single straight line. We interpret this finding by the fact that PS in these regimes is always in the same (glass) phase. For the first two layers where a local density minimum is identified (i.e., at NP separations of 0.75 and 1.20 nm) a volumetric phase transition is observed at a temperature of roughly 415 K, which is higher than Tg for pure polystyrene. Our observations are in qualitative agreement with the model proposed by Milner and Lipson.56,57 Some structural data such as the mean-squared radius of gyration, Rg2, the mean-squared end-to-end distance, Ree2, and the density, ρ, are presented in Table 1. They are in good

total stress tensor of the simulation box, S, is recovered by summing over the Vi-weighted atomic level stress tensors σi, i.e., S = ∑Ni=1Viσi/∑Ni=1Vi, where N is the number of atoms in the system. The atomic level stress tensors and the coordinates of all atoms are recorded every 2 ps. Stress profiles in a polymer matrix in contact with a flat surface, i.e., geometrical conditions differing strongly from the present spherical NP, have been reported already some time ago.25 The decisive influence of the shape and radius of a NP on the interphase properties has been commented on in detail in the literature.6−15 However, there has been no study focusing on stress profiles of a spherical interphase. The spatial variation of the excess Helmholtz free energy for insertion of a gas particle, ΔAex(r), is investigated by Widomtype insertions of the test molecule,83 r being the distance from the NP surface. Three linear molecules of different size, i.e., H2, O2, and CO2, and atomic He are used in the insertion procedure. The difference in the potential energy before and after the insertion of the chosen test particle, ΔE(r), is related to changes in ΔAex(r) according to

Table 1. Density ρ, Squared Polymer Radius of Gyration R2g, and Squared Polymer End-to-End Distance R2ee in Pure Polystyrene and the Nanocomposite at a Temperature of 300 K As Derived in the Present MD Simulationsa system

ρ (kg/m3)

R2g (nm2)

R2ee (nm2)

pure polystyrene nanocomposite

1070 ± 1.09 1113 ± 2.11

10.40 ± 4.05 8.84 ± 3.93

52.01 ± 30.52 46.19 ± 33.25

ΔAex (r ) = −kBT ln⟨exp( −ΔE(r )/kBT )⟩

where ⟨...⟩ denotes an ensemble average and kB the Boltzmann constant.39 The sample-averaged excess Helmholtz free energy, ΔA ex, is calculated by summing over the volume weighted ΔAex(r) values ΔAex =

a

For all quantities we have given the standard deviations (±), obtained by calculating mean values per configuration and then computing standard deviations over time. The pressure is 101.3 kPa.

σi ,KL

∑ (ri ,K j≠i

∑ ΔΑex(r)Vsh(r)/ ∑ Vsh(r) shells

(3)

shells

where Vsh(r) is the volume of the spherical shell at distance r from the NP surface. The summation covers all spherical shells. To the best of our knowledge, the work at hand is the first where free energy calculations of gas sorption in polymer nanocomposites are presented in spatially resolved form. We have, however, used a similar setup as described in a calculation of the sorption free energy map of apolar gas molecules near the flat surfaces of cellulose crystals.84 The interaction parameters for the test particles are shown in Table 2. The values proposed in ref 85 for CO2 and in ref 86

agreement with experimental data from the literature.78,79 In particular, the R2g of an atactic polystyrene sample with a weight-average molecular weight of 20 500 amu at 308 K has been measured to be 15.75 nm2.78 Its polydispersity index was approximately 1.02. The Cartesian components of the atomic level stress tensor, σi, of atom i are defined as25 ⎛ 1⎜ 1 = − ⎜mivi ,K vi ,L + Vi ⎝ 2

(2)

⎞ − rmin j ,K )Fi(min j),L⎟⎟ ⎠

Table 2. Force Field Parameters for the Inserted Gas Particles (He, H2, O2, CO2)a (1)

where Vi, mi, vi, and ri are the volume, mass, velocity vector, and position vector of atom i. Summation takes place over all possible i−j atom pairs. rmin j is the position vector of the minimum image of atom j with respect to reference atom i, and Fi(min j) is the corresponding force vector exerted on atom i by the minimum image of atom j. The subscripts K and L label the Cartesian coordinates x, y, and z. More detailed formulas can be found elsewhere.70 The electrostatic contributions to the atomic level stress tensor are calculated by the method of Sirk et al.80 Standard long-range corrections to the atomic level stress tensor for both van der Waals and Coulombic interactions are applied.70 The atomic volume is assumed to be identical to the volume of the Voronoi cell associated with each atom.26 The Voronoi cell is defined by requesting that all its space points are closer to the given atom than to any other atom. The sum of the Voronoi cells is space-filling. Neither overlaps nor empty regions appear; for more details see ref 81. The space tessellation is carried out by adopting the Laguerre− Voronoi method.82 For a system in mechanical equilibrium, the

He H in H2 molecule O in O2 molecule C in CO2 molecule O in CO2 molecule

ε (kJ/mol)

σ (nm)

partial charge

bond length (pm)

0.0848 0.3180 0.3630 0.2390 0.6870

0.2280 0.2318 0.3090 0.2792 0.3000

0 0 0 0.5888 −0.2944

74.1 101.7 116.3

The first two numbers refer to the Lennard-Jones energy (ε) and diameter (σ). In the next column atomic partial charges are given while the adopted bond lengths can be found in the last column. The parameters have been taken from ref 85 for CO2 and from ref 86 for He, H2, and O2. a

for He, H2, and O2 are employed. The structures of the test molecules are always fixed at their equilibrium values. Standard long-range corrections are applied to the energy calculations.87 ΔAex(r) is determined by performing 50 × 103 trial insertions per spherical shell. An exception is CO2, where 70 × 103 trial insertions per spherical shell are attempted. For the three molecules, 10 random orientations per insertion are tested. The D


Macromolecules

Article

in the matrix so that both the real NC and the idealized system containing the hard sphere occupy the same volume with the same overall density. The thickness of the depletion layer can be derived by

normalized autocorrelation function of the difference between the instantaneous value of ΔAex and its ensemble average value approaches zero after 0.1 ns. Thus, the time difference between configurations in the sampling of ΔAex, is set to 0.1 ns. The temperature dependence of the autocorrelation function in both systems is weak. The production runs at both temperatures are 20 ns long. Many results are reported as a function of the distance from the nanoparticle surface. Therefore, it is important to clarify how its position is defined. We follow our established procedure.11 The silica nanoparticle is cut from an α-quartz crystal structure by placing a sphere of 2 nm radius inside it and removing (i) all atoms outside this sphere, (ii) all unconnected atoms near the sphere surface, and (iii) all dangling −SiO3 units. The remaining surface oxygens are saturated by hydrogens; for details, see ref 11. The position of the surface is defined at its radius of 2 nm from its center. With respect to this value the distances from the surface are defined in the analyses described below. They are taken identical in all directions, thereby approximating the nanoparticle as a sphere. Deriving the structure of the nanoparticle from a quartz crystal is one possible choice since little is known about the internal structure of silica nanoparticles at the atomic scale. What is commonly called “amorphous silica” from precipitation, usually refers to amorphicity and porosity of aggregates of these particles on the micrometer scale. We have therefore taken quartz as a simple model to represent the unknown internal structure of the nanoparticle, which in experiment might be crystalline (including other SiO2 modifications), glass-like, or a combination thereof. As we are cutting through several different quartz crystal planes in our procedure, thereby exposing a variety of surfaces to the polymer, we expect some averaging of surface effects.

(R n + δ)3 − R n3 =

∫0

∞

r 2(1 − g (r )) dr

(4)

where g(r) is the pair distribution function between the nanoparticle and a polymer segment around it. The pair distribution function is defined as g(r) = ρ(r)/ρ0, where ρ0 is the bulk polymer density of the system and ρ(r) the radial density at a given distance r from the NP surface. The thicknesses of the depletion layers for the glass and the melt amount to 0.58 ± 0.03 and 0.57 ± 0.03 nm, respectively. The temperature dependence of δ is weak. Already here we want to mention that the interphase dimension derived with the help of δ is much smaller than the above estimates based on the mass density profiles. The spatial dependence of the orientation of the PS bonds and their alignment around the NP is quantified by the angle θ between a bond vector and the vector between the center of the NP and the geometric center of the bond, i.e., the surface normal. The second Legendre polynomial of cos θ, P2(r) = 1 /2(3⟨cos2 θ⟩ − 1), has been adopted to analyze the spatial orientation of the PS bonds. P2 ranges between −0.5 and 1.0, corresponding to a perfectly parallel or perpendicular orientation of the bonds with respect to the NP surface. A random orientation is characterized by P2 = 0. The variation of P2 for the C−C backbone bonds as a function of the distance between the geometric center of the bonds and the NP surface is shown in Figure 4. The blue (red) line corresponds to the

3. RESULTS 3.1. Local Structure and Dynamics. The radial mass density distribution of polystyrene around the nanoparticle is shown in Figure 3. The density is evaluated in approximately 0.05 nm thick spherical shells. The blue (red) line corresponds to the melt (glass). The existence of two pairs of sharp peaks and valleys characterizes the shape of both distributions. They are globally the same as found in a previous melt study at much higher temperatures.11 With the exception of a slightly higher first peak at lower temperature, pronounced differences between the melt and the glass distributions cannot be seen. An approximate interphase width can be estimated as the distance required to reach the average bulk density. The criterion that the difference between the local value of the radial mass density and the bulk density must be statistically different from zero leads to an interphase width of approximately 1.5 nm, corresponding roughly to two monomer diameters. This estimate is in line with previous studies employing both unentangled11,88 and entangled74,89 chains. The observed chainlength independence of the density profiles confirms the local nature of the density. In the present work we will also use the thickness of the depletion layer, δ,90 to estimate the interphase dimension. In contrast to all other quantities where we measure deviations from bulk parameters, δ is a single valued parameter to define the extension of the interphase. The thickness of δ defines a Gibbs equimolar dividing surface in the NC. If the polymer matrix had a uniform density equal to the bulk value, then a spherical hard sphere of radius Rn + δ would have to be inserted

Figure 4. Average order parameter P2 of the aliphatic C−C backbone bonds as a function of their separation from the nanoparticle surface accumulated in 0.05 nm thick bins for the melt (blue line, T = 450 K) and the glass (red line, T = 300 K). P2 refers to the second Legendre polynomial of the angle between the corresponding bond vector and the surface normal. The distance from the surface is determined relative to the geometric center of the two backbone carbons. The simulations have been performed at P = 101.3 kPa. A positive (negative) order parameter P2 indicates a bond orientation perpendicular (parallel) to the nanoparticle surface.

melt (glass). P2 is accumulated in approximately 0.05 nm thick spherical bins. A bond is assigned to a bin if its geometric center resides in the given spherical shell. We find a correlation between maxima (minima) in the radial density distribution and valleys (peaks) in the P2 radial distribution. The influence of temperature on the bond orientation is weak. The predicted E


Macromolecules

Article

spatial distribution of the second Legendre polynomial is in agreement with findings reported for coarse-grained polystyrene−silica composites,65 polystyrene−tightly cross-linked polystyrene,91 and polyethylene-like−silica systems.92 Beyond about 1.5 nm (glass) or 1.9 nm (melt), there is no preferential orientation of PS bonds. Within the MD error bars, we predict a slightly shorter interphase dimensions when measured via the PS mass density than via the bond orientation. This correlation has been expected as peculiarities in the density such as the density enhancement near the NP are a manifestation of the orientation of the PS bonds. To measure the local mobility of the atoms in the glass, their mean-square fluctuations, ⟨(Δr)2⟩, are calculated. A large value of ⟨(Δr)2⟩ indicates a flexible structure. The mean-square fluctuation of the ith particle, ⟨(Δri)2⟩, is given by ⟨(Δri)2 ⟩ =

1 tN

tN

∑ (ri(t j) −

Figure 5. Variation of the mean-square fluctuations ⟨(Δr)2⟩ of the two different types of carbon atoms as a function of the separation from the nanoparticle surface in the glass (T = 300 K, P = 101.3 kPa). The data were accumulated over 20 ns in 0.1 nm thick bins. The continuous horizontal lines indicate the respective bulk values.

ri ̅ )2

tj= 1

(5)

where tN is the total number of configurations in the considered simulation time of 20 ns, ri(tj) the position vector of the ith particle in the jth configuration at time tj, and ri̅ the average position vector of the ith particle when sampled over all configurations. The mean values and the standard deviations of ⟨(Δr)2⟩ for the aliphatic (backbone) and aromatic (phenyl) carbons are given in Table 3. The differences between the PS

is thus in agreement with the one predicted when employing the radial mass density and the bond orientation distributions. This similarity could be expected as a result of the strong correlation between the quantities considered, i.e., mass density and atomic positional fluctuations, on the one hand, and mass density and segment orientation, on the other. 3.2. Local and Atomic Level Stresses. The purpose of the atomic level stress evaluation is not to describe an atom itself but rather its local environment and its surrounding atomic cage.25 Nonzero atomic level stresses are caused by neighbor atoms that are not in positions corresponding to the minimum of the potential energy. An example for vanishing atomic level stresses is a perfect crystal with one atom per unit cell at zero temperature. Atomic level stresses can deviate from zero even for perfectly ordered systems when more than one nonequivalent atom coexist in the unit cell. Three atomic-level quantities that are accessible by atomic level stress calculations are the hydrostatic pressure, pi, the shear stress, σi,sh, and the von Mises stress, σi,vM. Although sometimes misleading, these three stress-based quantities share their names with continuum mechanics analogues. The present labeling has been employed already in previous atomistic studies.25,23,28 Nevertheless, we would like to emphasize that the stress quantities in the present study have been evaluated by conventional particle-based MD simulations without resorting to any kind of continuum mechanics approximation. The above elements are calculated from the components σi,KL of the atomic level stress tensor σi (see eq 1). The atomic level hydrostatic pressure of atom i, pi, is defined as the negative of the arithmetic average of the three diagonal components of the atomic level stress tensor: pi = −(σi,XX + σi,YY + σi,ZZ)/3.26 The atomic level hydrostatic pressure is a measure for local density deviations.25 A positive pi value is related to an average local density above the mean value. On the other hand, a negative value indicates a local density smaller than the mean value and thus a low coordination number. The atomic level shear stress, σi,sh is defined as the arithmetic average of the atomic level off-diagonal stresses. It is the simplest way to quantify anisotropies at the atomic level, i.e., deviations in the local environment of an atom from a perfectly symmetrical one.25 In the case of a perfect isotropic crystal, the σi,sh distribution is a delta function with a value of one at zero stress. In analogy to continuum mechanics,95 an atomic level von Mises shear stress, σi,vM, can be defined which also reflects

Table 3. Sample Averaged Mean-Square Fluctuations and Standard Deviations (±) of the Aliphatic and Aromatic Carbons in the Glass (T = 300 K and P = 101.3 kPa) of Pure Polystyrene and the Nanocompositea atom type aliphatic carbon (backbone) aromatic carbon (phenyl group) a

pure polystyrene (10−3 nm2)

nanocomposite (10−3 nm2)

2.71 ± 0.17

2.68 ± 0.18

5.78 ± 0.28

5.73 ± 0.29

Averaging was performed over 20 ns.

and NC values are negligible as a result of the low NP volume fraction. The aromatic carbons are characterized by a higher ⟨(Δr)2⟩ value than the aliphatic ones which may be attributed to rotations or librations of the phenyl rings around their C2 axis. The ⟨(Δr)2⟩ profile of the two types of carbon atoms as a function of the distance from the NP surface is shown in Figure 5. The continuous horizontal lines correspond to the respective bulk values. The behavior of the aliphatic carbons is different from that of the aromatic carbons. A smooth monotonic increase of ⟨(Δr)2⟩ is observed for the aliphatic carbons until the final bulk value is reached. For the aromatic carbons we observe a sudden increase at a separation of 0.5 nm from the NP surface, i.e., after the first monomer layer. We assume that they correspond to the onset of librational/rotational modes which seem to be hindered in the high density region near the NP. Quite generally, the present results indicate that the atoms close to the nanoparticle have a smaller spatial flexibility than the bulk atoms. These findings are in line with previous studies on nanocomposites where other descriptors for local flexibility have been employed.93−94 The distance necessary to reach the bulk value of ⟨(Δr)2⟩ can also serve as an estimate of the interphase width. It is found to be slightly above 1.2 nm. The estimated interphase length on the basis of the ⟨(Δr)2⟩ profile F


Macromolecules

Article

addition of the nanoparticle shifts the distribution to slightly higher values. The mean values of σi,vM are 13.77 and 13.51 GPa with standard deviations 12.95 and 12.88 for the NC and PS. In analogy to σi,sh, also the two σi,vM numbers indicate that the NP inclusion results in a very small (if any) enhancement of stress anisotropies. The σi,vM distributions in Figure 6c have been obtained by averaging the atomic level von Mises stress distribution of each atom type (i.e., group of atoms which share the same nonbonded interaction parameters) in the two systems. The individual σi,vM distribution for each atom type of pure polystyrene at 300 K is presented in the Supporting Information (Figure S1). The curves for the five atom types are unimodal. The local stress of the atoms increases in the order Harom < Haliph < Carom < Caliph. Based on the notion of the atomic level stress tensor, it is possible to formulate a local stress tensor for a specific region of a material by suitably averaging the σi values of those atoms that reside in the region of interest. The local stress tensor within a spherical shell around the nanoparticle at distance r from the NP surface, σsph(r), can be computed by summation over the Vj-weighted σj tensors of those Nsph atoms that lay within the predefined spherical shells: σsph(r) = ∑j=1NsphVjσj/∑j=1NsphVj.98 By such an approach we can examine the influence of a nanoparticle dispersion on the local stress distribution of the polymer matrix. The hydrostatic pressure, psph, the shear stress, σsph,sh, and the von Mises shear stress, σsph,vM, can be calculated at local resolution using the previous equations where the atomic level stress tensor σi is replaced by the volume averaged stress tensor σsph. It should be noted that the atomic level stresses induce rather long-ranged stress fields around each atom.25 Therefore, the local stress should be interpreted as the sum of two contributions: the first stems from geometrical misfits of an atom with its surrounding while the second one is induced by the long-range stress fields of the neighboring atoms. The variation of the σsph,vM of the aliphatic hydrogens as well as the total psph and σsph,sh as a function of the separation from the NP surface in the melt is displayed in Figure 7. We have studied the spatial variation of the von Mises stress of one atom

the degree of asymmetry in the stress environment of an atom.26 With the help of Cartesian stress components, it reads 1 σi2,vM = [(σi , XX − σi , YY )2 + (σi , YY − σi , ZZ)2 2 + (σi , ZZ − σi , XX )2 + 6(σi2, XY + σi2, YZ + σi2, XZ)] (6) The first three terms on the right-hand side of the above equation measure differences between the normal stresses in the system. The last right-hand term represents the deviation of the atomic level shear stresses from the expected value of zero in a system without overall shearing forces. In a perfectly isotropic environment the von Mises shear stress would be 0. σi,vM is also often used to characterize failure-prone regions in materials, both in continuum mechanics96 and in molecular simulations.97 The spatial regions with the highest von Mises stresses have the highest failure probability. The normalized probability distributions of pi, σi,sh, and σi,vM for PS and the NC at 300 K are given in Figure 6. The

Figure 6. Distribution of (a) the atomic level hydrostatic pressure, pi, (b) the atomic level shear stress, σi,sh, and (c) the atomic level von Mises stress, σi,vM, in pure polystyrene (red) and the nanocomposite (black). An averaging time of 20 ns has been chosen. The simulations have been performed at T = 300 K and P = 101.3 kPa.

contribution of all particles, including the silica atoms, is shown in the plots. The red line corresponds to PS and the black line to the NC. The atomic level stresses have been averaged over 20 ns. The individual pi values differ by orders of magnitude from the overall pressure of the simulation box, which is equal to 101.3 kPa. These differences indicate cancellation effects when summing up the atomic level stresses to compute P via the atomic virial expression. This observation is in agreement with previous studies.18,25,29,26 A symmetric single peak structure in the pi and σi,sh distributions is found both in PS and in the NC. The pi and σi,sh profiles are almost identical, signifying a negligible enhancement of stress anisotropies due to the NP presence. The mean value of the σi,sh distribution is zero since there is no applied external shear while that of pi is close to 101.3 kPa, which is the pressure of the simulation box. We also predict one peak in the σi,vM distributions with an extended tail in the direction of higher stress values. The

Figure 7. Variation of (a) local hydrostatic pressure, psph, (b) local shear stress, σsph,sh, and (c) local von Mises stress, σsph,sh, in spherical shells with a given distance from the nanoparticle surface for the melt (450 K). The error bars correspond to the standard deviation around the local mean values, calculated by ten independent block averages. G


Macromolecules

Article

An estimate of the interphase length can be also obtained by analyzing the spatial stress profiles in Figures 7 and 8. It is the distance from the NP surface that is required to converge the stresses to their bulk values. All prominent features in the local stress profiles occur at distances below 0.5 nm, which is the position of the first density maximum (Figure 3). This indicates that stress perturbations are restricted to the immediate interface; they decay within the first layer of polymer atoms. Smaller deviations in the local density from the bulk value are not coupled to local stress modifications. Additional simulations are required to analyze whether such a correlation between structural properties and stresses is generally valid. The interphase width measured by local stress quantities is smaller than that measured by the polymer local structure and dynamics. Interestingly, though, it agrees with the thickness of the first polymer layer at the NP surface and coincides with the thickness of the depletion layer, δ. 3.3. Gas Sorption. We have also studied the influence of NP inclusion on the solubility of small gas molecules and its spatial variation. We discuss solubility in terms of the excess Helmholtz free energy ΔAex associated with the sorption process. The average values ΔA ex are presented in Table 4. Irrespective of the temperature of the system, ΔA ex decreases in the order He > H2 > O2 > CO2, indicating a solubility enhancement with the particle size. In fact, we observe an approximate linearity between ΔA ex and the sum of the Lennard-Jones parameters ε of the sorbed molecules (not shown), a dependence which has been found before.85 The errors of ΔA ex have been calculated by averaging over ten independent blocks of 2 ns; they increase with penetrant size. The temperature and thus also the phase of the system have a small but non-negligible influence on ΔA ex. The values of ΔA ex for the glass are always above the melt elements, indicating that gas sorption is easier in the melt. At least under the present simulation conditions this might be traced back to a higher polymer flexibility with increasing temperature. The same behavior as predicted for the sorption of a single gas molecule, i.e., increased solubility at higher temperature, is also often found for supercritical and near-supercritical gases.99 The ΔA ex difference between melt and glass is largest for the smallest gas particle, i.e., the He atom where it exceeds 4 kJ/mol. It is lowest for O2 (≈2 kJ/mol). The standard deviations for a given molecule are slightly lower for the melt than for the glass. The NC systems have almost the same overall sorption free energies as the pure PS samples, which seems to be due to the small volume fraction of silica in the system. A larger difference would be expected for higher degrees of filling. The experimental sorption free energies for O2 and CO2 in an atactic polystyrene glass at 300 K are 4.39 and −2.05 kJ/mol, respectively.100 For 450 K a value of −4.19 kJ/mol is obtained for CO2 by the correlations proposed in ref 101. Our

type since the width of the associated von Mises stress distribution is considerably shorter than the width of the overall distribution. The quantities have been accumulated in approximately 0.025 nm thick spherical shells. Atomic level stresses centered on nanoparticle atoms are not shown. For sufficiently large separations from the NP, psph and σsph,sh have the expected bulk values of roughly 101.3 and 0 kPa, since no shearing forces are applied to the studied systems. The parameter σsph,vM reaches a finite nonzero value due to the amorphous structure of the polymer which lacks any specific ordering. The absolute spatial fluctuations of σsph,sh around the bulk value are less intensive than those encountered in psph and σsph,vM. Note the different ordinate scales adopted in Figure 7. One prominent peak and two valleys can be identified in the hydrostatic pressure profile within a 0.75 nm interval around the NP surface. One maximum and minimum pair is found in the σsph,sh profile in the immediate vicinity of the NP surface. An almost flat σsp,vM spatial distribution is predicted. We have observed the same behavior in the spatial variation of the σsp,vM of the other four atom types that exist in our systems. The σsph,vM of the aliphatic hydrogens as well as the total psph and σsph,sh profiles for the NC in the glass are given in Figure 8.

Figure 8. Variation of (a) local hydrostatic pressure, psp, (b) local shear stress, σsp,sh, and (c) local von Mises stress, σsp,vM, with the separation from the nanoparticle surface for the glass (300 K). The error bars correspond to the standard deviation around the local mean values, calculated by ten independent block averages. Note the difference in the y-axis scaling between Figures 7 and 8

The bulk value of σsph,vM for large separations from the NP surface is higher than in the melt. All three distributions show qualitatively the same behavior as in the melt; peaks and valleys occur at the same distances from the surface. Their intensity is, however, stronger than in the melt.

Table 4. Excess Helmholtz Free Energy of Insertion, ΔA ex, of He, H2, O2, and CO2a pure polystyrene (kJ/mol)

a

gas

melt (450 K)

He H2 O2 CO2

8.84 4.61 3.51 −7.34

± ± ± ±

0.06 0.18 0.75 1.36

nanocomposite (kJ/mol)

glass (300 K) 13.11 6.93 5.53 −4.15

± ± ± ±

0.08 0.22 0.85 1.52

melt (450 K) 8.74 4.58 3.58 −7.42

± ± ± ±

0.07 0.20 0.78 1.42

glass (300 K) 12.87 6.93 5.43 −4.12

± ± ± ±

0.09 0.24 0.89 1.59

The simulations have been performed at a pressure of 101.3 kPa. H


Macromolecules

Article

follows the available accessible volume. Much more surprising is the existence of He solubility fluctuations beyond 1.5 nm. At this distance, the polymer density has converged to a bulk-like behavior. We reconcile these facts by assuming that, although the mass density is constant beyond 1.5 nm, the size and distribution of cavities are not. The monomer orientation (Figure 4) supports this view, as some order persists to a distance of 2 nm. The fact that such a wide interphase is only seen in the sorption free energy profile of He is attributed to its smallness. The cavities existing beyond 1.5 nm seem to be too small to be accessible to the studied molecular species, so any variation of them can only be probed by He. The sorption free energy profiles in the glass (Figure 10) do not differ qualitatively from those in the melt (Figure 9).

predictions are in sufficient agreement with the cited experimental estimates. Both the differences between different gas species at a given temperature and the temperature dependence are reproduced by our simulations. The spatial profile of the excess Helmholtz free energies for the sorption of molecules in the vicinity of the nanoparticle is shown in Figure 9 for the melt (T = 450 K). It is sampled in

Figure 9. Profile (blue curve) of the excess Helmholtz energy for gas sorption ΔAex(r) for (a) He, (b) H2, (c) O2, and (d) CO2 as a function of the distance from the nanoparticle surface for the melt (450 K). The red line corresponds to the sample average. The error bars measure the standard deviation around the local mean values, calculated from ten independent block averages.

spherical shells with a width of 0.1 nm around the NP. Three results are obvious. First, for all gas molecules considered, ΔAex is smaller in the immediate vicinity of the nanoparticle (

How Thick Is the Polymer Interphase in Nanocomposites? - American

Recommend Documents