Determination of the Mechanical Properties of a Poly(methyl

Article pubs.acs.org/Macromolecules

Determination of the Mechanical Properties of a Poly(methyl methacrylate) Nanocomposite with Functionalized Graphene Sheets through Detailed Atomistic Simulations Emmanuel N. Skountzos,† Alexandros Anastassiou,† Vlasis G. Mavrantzas,*,†,‡ and Doros N. Theodorou§ †

Department of Chemical Engineering, University of Patras and FORTH-ICE/HT, GR 26504, Patras, Greece Department of Materials, Polymer Physics, ETH Zürich, HCI H 543, CH-8093 Zürich, Switzerland § School of Chemical Engineering, National Technical University of Athens, GR 15780, Athens, Greece ‡

S Supporting Information *

ABSTRACT: Recent experimental studies have demonstrated that the introduction of oxygen-containing functional groups in graphene sheets can greatly enhance the mechanical properties of their nanocomposites with polar polymers even at extremely low loadings. Motivated by these reports, we determine here the elastic constants of syndiotactic poly(methyl methacrylate) (sPMMA) at small wt % loadings of graphene sheets through atomistic modeling. To carry out a comparative study of the effect of graphene functionalization on the degree of mechanical reinforcement, we address both pure (i.e., unfunctionalized) and functionalized graphene sheets bearing epoxy and hydroxyl groups randomly bound on both sides of their surface in the host sPMMA matrix. The calculation of elastic constants (which involves no adjustable parameters) follows the methodology originally proposed by Theodorou and Suter [Macromolecules 1986, 19, 139], and has been based on the use of the Dreiding all-atom force-field. Our predictions for the elastic constants (which for the pure sPMMA matrix are within the error bars of experimentally computed values) suggest a substantial increase in the elastic constants, especially in the case of functionalized graphene sheets. For example, at just 5.67 wt % loading of the host matrix in functionalized graphene sheets, they indicate an improvement in Young’s modulus E by ∼74%, in the bulk modulus B by ∼19%, and in the shear modulus G by 83%. Our results fully corroborate recent experimental measurements about the unique opportunities that functionalized graphene sheets offer for the design of new, very strong multifunctional materials at low nanofiller content.

1. INTRODUCTION Numerous experimental and theoretical (analytical and computational) studies in the last years have indicated that several properties of polymer matrices (mechanical, electrical, optical, thermal, rheological, and barrier) can be significantly enhanced by the addition of nanosized particles, even at ultralow fractions of the nanofiller.1−24 Such a substantial property improvement is due to the enormous surface area per unit volume characterizing nanoparticles, in distinct contrast to traditional composite structures where larger quantities of the filler are typically needed. The performance of nanocomposites can be enhanced even further by optimizing the nanoparticle (NP) surface chemistry to come up with better dispersion of the NPs, stronger attractive polymer−particle interactions and better polymer adhesion on the NP. This is typically accompanied by a non-negligible increase in the glass transition temperature Tg of the nanostructure which can fundamentally change all physical properties of the material, including processability and thermal stability.2 Another way to improve properties (in particular the phase behavior and the viscoelastic response) is through chain grafting onto the NPs. The key role of interface chemistry and nanoscale morphology in the © XXXX American Chemical Society

performance of polymer nanocomposites has been demonstrated by several studies in recent years.2,5,6,9,10,13−15,18,19,21 Nanoparticles with a variety of shapes can be employed as fillers in polymer nanocomposites: the list includes nanoclays (primarily polymer-layered silicates, PLS), nanoscopic silica particles, nanotubes, nanofibers, and many others. Among the carbon-based nanofillers, and given the prohibitively high cost of carbon nanotubes (CNTs), a material that has attracted considerable attention in the last years is graphite, from which one can produce (with several but nontrivial experimental procedures) expanded graphite (EG), graphitic nanoplatelets (GNP), single-layer graphene sheets (GS), and eventually functionalized graphene sheets (FGS).25−28 According to Schniepp et al.,26 FGS can be produced by a process involving thermal exfoliation of graphite oxide. They are single-layer graphene sheets containing randomly bound epoxy and hydroxyl groups to both of their sides. The epoxy and hydroxyl groups themselves lie ∼0.2 nm above the underlying carbon Received: August 27, 2014 Revised: October 17, 2014

A

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

Macromolecules

Article

nanocomposites. Li et al.17 have further reported that graphene sheets tend to aggregate in the polyethylene matrix, which lowers the degree of their interactions with the host matrix. Miscibility of graphene sheets in PMMA (which is a very important issue for practical applications of the corresponding nanocomposites) has been investigated by Ju et al.31 through the dissipative particle dynamics (DPD) simulation method for different NP volume fractions and different values of the DPD repulsive interaction parameter (to effectively account for surface functionalization). Our goal in this work is to contribute to our understanding of the intimate relationship between atomic-scale structural features and macroscopically exhibited mechanical properties of graphene-based polymer nanocomposites and of the role of surface functionalization through a detailed simulation study that leads directly to the determination of the elastic properties. We wish to offer a more fundamental (at molecular level) understanding of interfacial interactions in FGS-PMMA nanocomposites by investigating in detail PMMA structure and conformation in the bulk and next to graphene sheets, and their connection with elastic constants. In the literature, little has been published on the atomistic simulation of interfacial behavior in graphene-polymer composites to date.32−34 Awasthi et al.32 have carried out atomistic MD simulations using the Consistent Valence Force Field (CVFF) to study nanoscale load transfer between polyethylene and a graphene sheet and characterize the force−separation behavior between CNTs and a polymer matrix. Separation studies were conducted for opening as well as sliding modes, and cohesive zone parameters (such as peak traction) and the energy of separation for each mode were evaluated as a first step toward the development of continuum length-scale micromechanical models for tracking the overall material response incorporating interfacial phenomena. Wang et al.34 have investigated the influence of surface functionalization and layer length on the interfacial load transfer in graphene-polyethylene nanocomposites, using the ab initio polymer consistent force field (PCFF). Their simulations showed that oxygen-functionalized graphene sheets lead to larger interfacial shear force than hydrogen-functionalized or pristine ones during pull-out process. Increasing the oxygen coverage and the layer length enhanced interfacial shear force, but further increase of oxygen coverage to about 7% led to a saturation in the interfacial shear force. In the present work, for the determination of the mechanical properties, we follow a methodology already developed in the past for a simpler class of systems (glassy vinyl polymers such as polypropylene and polystyrene) based on small-strain deformation experiments on the computer of microscopically detailed model structures.35,36 The procedure involves several modeling and mathematical steps and leads to the accurate description of the elastic constants (Young’s modulus E, bulk modulus B, shear modulus G, and Poisson’s ratio ν) of a polymeric glass based on the assumption that vibrational contributions of the hard degrees of freedom are not significant, implying that estimates of the elastic constants can be obtained by computing changes only in the total potential energy of static microscopic structures subjected to simple deformations. For glassy atactic polypropylene for which the method was first developed and implemented by Theodorou and Suter,35 elastic constants were predicted within 15% of the experimental values. In this paper, we extend the method to the more complicated case of graphene-based PMMA nanocomposites.

grid. The carbon atom functionalized by an epoxy or a hydroxyl group is transformed from a planar sp2-hybridized to a distorted sp3-hybridized geometry. Expanded graphite is composed of many graphene sheets held together by van der Waals forces in rigid nanoplatelets that are hundreds of nanometers thick. By using sonication, these rigid nanoplatelets can be broken apart into thinner platelets which then can be dispersed within a solution of the polymer using high-speed shearing methods. Ramanathan et al.2 have reported such a technique to disperse these platelets within a solution of PMMA which led up to a 30 °C increase in the Tg of PMMA at 1−5 wt % loading of the GNP nanofiller, and up to an 80% enhancement of Young’s modulus for a 1 wt % FGS-PMMA nanocomposite material.2 A similar observation has been reported by Li and McKenna16 for a graphene oxide/PMMA nanocomposite material. In the literature, significant improvements in the mechanical (but also electrical, thermal and barrier) properties of graphenebased polymer nanocomposites have been reported for many other glassy polymers, as summarized in several reviews.10,11,18 George and Bhowmick,3 for example, observed a 100% increase in the modulus of poly(ethylene-co-vinyl acetate) reinforced with EG at 4% loading. Fang et al.5 have reported an improvement by 70% for the tensile strength and by 57% for Young’s modulus for graphene−polystyrene nanocomposites with polystyrene chains grafted onto graphene sheets by atomic transfer radical polymerization. Eda and Chhowalla7 studied the electrical properties of solution-processed composite semiconducting thin films consisting of FGS as the filler and polystyrene as the host material; increasing the average size of FGS was found to significantly enhance carrier mobility and thus device performance, thereby demonstrating the possibility of using a commodity plastic to develop low-cost, macro-scale thin film electronics. From a molecular point of view, the experimental studies of Ramanathan et al.2 and Li and McKenna16 have shown that behind the extraordinary mechanical properties of FGS-PMMA nanocomposites (compared to, e.g., EG-PMMA ones) are the enhanced interfacial interactions caused by oxygen functionalities across the surface of graphene with PMMA chains. FGS contain pendant hydroxyl groups on their surfaces which may form hydrogen bonds with the carbonyl groups of PMMA. Additional factors typically invoked2 to explain the extraordinary mechanical properties of FGS-PMMA nanocomposites are the nanoscale surface roughness of FGS, the defects caused during thermal exfoliation of the precursor graphite oxide, and their wrinkled topology at the nanoscale due to their extremely small thickness;26−28 these can enhance mechanical interlocking with the polymer chains and, consequently, lead to better adhesion. The above observations and considerations are in agreement with a recent molecular dynamics (MD) simulation study29,30 which showed strong adhesion of PMMA chains on graphene (especially of the side groups) and considerably slower segmental and chain mobility in the interfacial area. The MD simulations suggest that local mass density, segmental dynamics, and chain terminal relaxation all differ from the bulk behavior for distances extending up to several nanometers from the graphene surface. MD simulations have also been employed by Li et al.17 in their study of the effect of nanoparticle shape on the viscoelastic properties of a polyethylene matrix. They found that it is the surface-tovolume ratio of nanoparticles that plays the dominant role in the structural, dynamical and viscous properties of polyethylene B


Macromolecules

Article

The paper is organized as follows: In section 2, we provide details of the selected systems, of their geometric parameters, and of the potential energy functions describing bonded and nonbonded interactions. In section 3 we describe the statisticalmechanical method (including elements of statistical thermodynamics) for the treatment of the mechanical response to deformation and the estimation of the elastic constants at room temperature. Our predictions for the structural and mechanical properties of the simulated model systems (neat sPMMA and its nanocomposites with graphene) are presented and discussed in section 4. The paper concludes with section 5, summarizing the major conclusions of this work and discussing current and future efforts.

2. SYSTEMS STUDIED, MOLECULAR MODEL, AND SIMULATION METHODS The selected polymer system is syndiotactic PMMA (all chains in the simulation box consist of syndiotactic sequences) at 1 atm and 300 K. PMMA is a strong and transparent thermoplastic, often used as a lightweight or shatter-resistant alternative to glass.37,38 Commercial grades of PMMA (melting temperature Tm = 160 °C, entanglement molecular weight Me = 10 kg/mol) behave as a glass at room temperature (the glass transition temperature Tg ranges from 85 to 165 °C) and have a density between 1.17 to 1.20 g/cm3, which is less than half that of glass. The glass transition temperature Tg of infinite molecular weight sPMMA, in particular, is ∼385 K;39 thus, at the conditions of interest (300 K and 1 atm), the material is glassy, with an experimentally measured specific volume of 0.84 cm3/g.37,38 Clearly, sPMMA can crystallize under ambient conditions. In our simulations, however, all sPMMA model structures were annealed at 500 K for several hundreds of nanoseconds to render them completely amorphous prior to quenching to the conditions of interest; our predictions, therefore, of the elastic constants in the next sections of this paper will reflect the mechanical properties of the amorphous phase of the studied model sPMMA structures. The simulations have been executed with rectangular parallelepiped simulation cells with sPMMA chains subject to three-dimensional periodic boundary conditions. In our simulations (see Table 1 and Figure 1), strictly monodisperse

Figure 1. Typical atomistic configurations of: (A) an sPMMA chain, (B) a nonfunctionalized graphene sheet, and (C) a functionalized graphene sheet from the NPT MD simulations at 500 K. The FGS contains five hydroxyl (OH) groups and three oxygen (O) atoms. In all cases, carbon, oxygen, and hydrogen atoms are represented with gray, red, and white colors, respectively.

sides with oxygen-containing functional groups (hydroxyl −OH, epoxy =O) whose number was chosen (see Table 2) so that the resulting elemental concentration of carbon, oxygen and hydrogen atoms corresponded closely to the elemental analysis and surface area data of the functionalized graphene sheets used in the experimental study of Ramanathan et al.2 Unfunctionalized and functionalized graphene sheets had lateral dimensions of (12 × 12) Å2. Taking into account the terminal hydrogen atoms as well, we find that the concentration of the model GS-sPMMA and FGS-sPMMA nanocomposites in graphene sheets was 5.67 wt % and 6.54 wt %, respectively. The simulations were carried out with the very accurate, explicit-atom Dreiding40 force-field for the calculation of atomistic interactions between PMMA and graphene sheets. All the force-field details (such are schematic representations, potential energy functions and force-field parameters) are presented in the Supporting Information. Initial configurations for all systems were created using the MAPS 3.2 software package41 to build the initial packing cell followed by a static structure optimization using a molecular mechanics algorithm to remove overlaps. The optimization stops when the maximum energy gradient with respect to all

Table 1. Details (Number of Chains and Number of Graphene Sheets) of the Model sPMMA Systems and Their Graphene-Based Nanocomposites Studied in This Work system

sPMMA chains (15 MMA/chain)

graphene sheets

functionalized graphene sheets

1 2 3

27 27 27

− 3 −

− − 3

samples were assumed, each sPMMA chain consisted of 15 monomers (methyl methacrylate, MMA) implying a molecular weight of 1503.75 g/mol, and the simulation cell contained 27 chains. Three model systems were studied (considering computer time limitations): (a) the neat sPMMA matrix (no graphene sheets added), (b) its nanocomposite with three unfunctionalized monolayer graphene sheets (denoted as GSsPMMA in the following), and (c) its nanocomposite with three functionalized monolayer graphene sheets (denoted as FGS-sPMMA in the following). GS and FGS were terminated with hydrogen atoms. FGS were patterned on both of their C


Macromolecules

Article

Table 2. Composition of Graphene Sheets in Hydroxyl and Epoxide Functional Groups and the Corresponding Percentage of Carbon, Oxygen and Hydrogen Atoms number of hydroxyl (−OH) groups

number of epoxide (O) groups

percentage (%) of carbon atoms

percentage (%) of hydrogen atoms

percentage (%) of oxygen atoms

5

3

65

27

8

simulations, we did not explicitly account in our calculation of the tail correction for local variations in the density profile. However, test simulations with various values of the LJ cutoff distance (from small to very large, up to 16 Å) showed that the value of 12 Å used here together with a naı̈ve calculation of the tail correction based on the bulk density value did not indicate any significant changes. Thus, although we plan to carry out a compete treatment of the tail correction problem in a subsequent study, not directly accounting for the inhomogeneous density profile in the present simulations should not be considered as a serious drawback.

microscopic degrees of freedom becomes less than 10 kJ/(mol nm). The resulting minimum-energy configurations for the various model systems were submitted next to an MD simulation in the NPT statistical ensemble at T = 500 K and P = 1 atm, using a Nosé−Hoover thermostat-barostat42,43 with relaxation times equal to 10 and 350 fs, respectively. The simulation cell was subjected to standard periodic boundary conditions in all three dimensions. Electrostatic interactions were computed using the PPPM (particle−particle particle− mesh)44 method with a real space cutoff of 12 Å. For the Lennard-Jones (LJ) interactions, the cutoff radius was set equal to 12 Å. All MD simulation results reported here have been obtained with the LAMMPS software.45 The equations of motion were integrated with the velocity Verlet integrator46−48 with a time step equal to 1 fs. A typical equilibrated atomistic configuration from our simulations with the FGS-sPMMA system is reported in Figure 2.

3. METHODOLOGY FOR THE ESTIMATION OF THE MECHANICAL PROPERTIES To compute the mechanical properties of the simulated polymer nanocomposites, we resorted to the method originally developed by Theodorou−Suter35 for the estimation of the mechanical properties of a glassy vinyl polymer. According to the continuum formulation of thermoelasticity, the polymer is regarded as a collection of material points whose positions are denoted by the vector r0 = (r0,1,r0,2,r0,3) in the reference (undeformed) state and by the vector r = (r1,r2,r3) in the deformed state. In our work, the role of these material points is played by the center of the simulation box (unit cell) containing the molecular model of the material, and by the centers of its periodic images. Deformations are applied by assuming a change in the shape of the simulation box. Thus, if a, b, and c denote the vectors defining the dynamic cell, the state of the system is specified by the tensor h = {a,b,c} whose value in the undeformed state is denoted as h0. The vector s with components si = ri − r0,i is the displacement vector for each material point. The material strain tensor ε is defined as εLM =

Figure 2. Typical atomistic snapshot of the FGS-sPMMA model system. In the simulation cell there exist 27 sPMMA chains and three functionalized graphene sheets with five hydroxyl groups and three oxygen atoms on their surface (6.54 wt % loading). Cell dimensions: (40 × 40 × 40) Å3. Carbon (PMMA), carbon (FGS), and oxygen atoms are represented with gray, yellow, and red color, respectively.

+

3 ⎞ 1 ⎛ ∂s ∂r ∂ri ∂s ⎞ 1⎛ ⎜⎜∑ i − δLM ⎟⎟ = ⎜⎜ L + M ⎟⎟ ∂r0, L ⎠ 2 ⎝ i = 1 ∂r0, L ∂r0, M 2 ⎝ ∂r0, M ⎠

3 ∂s ∂si ⎞ 1⎛ ⎜⎜∑ i ⎟ 2 ⎝ i = 1 ∂r0, L ∂r0, M ⎟⎠

(1)

If we limit ourselves to small deformations, second order terms in the derivatives of the displacement vector with respect to position can be neglected; the sum over i appearing on the right-hand side of eq 1 drops out and ε can be expressed simply in terms of changes in the components of h. The corresponding material stress tensor is denoted by τ and the element τLM is taken equal to the force per unit area acting on an element of surface perpendicular to axis L and along the direction M. The tensors ε and τ are symmetric; thus they can be represented in the Voigt notation as vectors of six components. The two notations are typically used interchangeably by making use of the following correspondence:

From a technical point of view, an important issue to mention is the calculation of the long-range contribution to the nonbonded interactions due to the nonuniform density profile of the PMMA polymer system next to graphene sheets. This issue has been addressed in detail by Mansfield and Theodorou49 and by Daoulas et al.50 For semi-infinite systems, failure to properly account for the long-tail correction might lead to an artificial imbalance between cohesive and adhesive interactions, causing the latter to appear too strong.49 Because of the presence of many graphene layers in the present D


Macromolecules

Article

ε1 = ε11 τ1 = τ11 ε2 = ε22 τ2 = τ22 ε3 = ε33 τ3 = τ33

U = A + TS = U0 + V0 ∑ [τLM + ρ0 cεTγLM ]εLM LM

⎡ ∂C 1 + V0 ∑ ∑ ⎢CLMNK − T LMNK 2 LM NK ⎢⎣ ∂T

ε4 = 2ε23 τ4 = τ23 ε5 = 2ε31 τ5 = τ31 ε6 = 2ε12 τ6 = τ12

The fundamental thermodynamic equation for an elastic system in terms of the stress and strain tensors, at constant chemical composition, is dU = T dS + V0 ∑ τLM dεLM = T dS + V0 ∑ τIdεI I

(3)

In the Helmholtz free energy representation (A = U − TS), the corresponding fundamental thermodynamic expression in differential form reads

U = U0 + V0 ∑ ∑ {τij + ρ0 cεTγij}εij i

dA = −S dT + V0 ∑ τLM dεLM = −S dT + V0 ∑ τI dεI LM

I

This allows one to define the fourth order tensor of isothermal elastic coefficients CLMNK determining the relationship between strain and tension tensors:

Cijkl = T , ε(LM , NK )

(5)

(6)

where μ and λ are the Lamé constants. The tensile (Young’s) modulus E, the shear modulus G, the bulk modulus B, and the Poisson ratio v are related to μ and λ by the following equations: E=μ

3λ + 2μ λ+μ

G=μ B=λ+ v=

2 μ 3

λ 2(λ + μ)

1 ∂ 2U V0 ∂εij∂εkl

T , ε(ij , kl)

(10)

According to this, one can calculate the elastic constants of a glassy polymer by first computing its internal energy U around a minimum in configuration space as a function of the imposed strain tensor ε and then fitting the results with a second-order polynomial; the fitting coefficient of the polynomial’s second order term will provide the values of the elements Cijkl, thus also of the two Lamé constants and consequently of the material’s elastic properties. On the basis of the above considerations, the computational procedure for simulating deformation in glassy polymers can be summarized in the following steps: 1 Start with an undeformed structure in detailed mechanical equilibrium and minimize its total potential energy using a molecular mechanics technique. 2 Choose the type of deformation to be imposed and a small degree of deformation ε (up to 1%). 3 For the given deformation, apply the assumption of affine deformation to define (based on the new shape of the simulation cell) the atomic positions in the deformed state. In our work, the affine deformation is applied to all atoms in the system (both polymer and graphene atoms). 4 Reminimize the total potential energy of the deformed model structure using molecular mechanics calculations. 5 Repeat the procedure (steps 1 to 4) with different initial (minimum energy) structures. 6 Average over all deformed model structures to get estimates of those elastic constants of the simulated system that are involved in the given type of deformation, by fitting the numerical results for the potential energy function with a second order polynomial. 7 Repeat the procedure (steps 1 to 6) for different types of deformation.

As a result of the Voigt symmetry relationships, the above equation can be condensed into a symmetric 6 × 6 matrix C. For an isotropic material (such as an amorphous glassy polymer), the matrix C takes the form ⎡ 2μ + λ λ λ 0 0 0⎤ ⎥ ⎢ λ 2μ + λ 0 0 0⎥ ⎢ λ ⎥ ⎢ λ λ 2μ + λ 0 0 0 ⎥ ⎢ C= ⎢ 0 μ 0 0⎥ 0 0 ⎥ ⎢ ⎢ 0 0 0 0 μ 0⎥ ⎥ ⎢ 0 0 0 0 μ⎦ ⎣ 0

(9)

implying that

2

∂τLM 1 ∂A = ∂εNK V0 ∂εLM εNK

j

1 + V0 ∑ ∑ ∑ ∑ Cijklεijεkl 2 i j k l

(4)

CLMNK =

(8)

where the first terms in the brackets on the right-hand-side are due to the strain derivatives of the Helmholtz energy while the second terms in the brackets originate from the corresponding derivatives of entropy. Theodorou and Suter35 showed that for a typical glassy polymer (e.g., PMMA) at room temperature the internal energy contribution to the isothermal elastic coefficients is considerably more significant than the entropic one. Including the entropic contribution is certainly possible51 but computationally more expensive. In this work, in the interest of keeping computations tractable, we will neglect the entropic contribution to the elastic constants and write

(2)

LM

⎤ ⎥εLM εNK ⎥ ε⎦

(7)

Let us consider an elastic solid subjected to an arbitrary isothermal small deformation. Expanding the internal energy U of the system into a Taylor series around the undeformed state up to second order, and using the definitions introduced above, we can write E


Macromolecules

Article

8 Average again the results for the new deformation experiments to get even more accurate predictions for the elastic constants.

U = U0 − V0(1/3)tr(τ + ρ0 cεTγ )0 ε + (1/2)V0(C1111 + C1122 + C1133 + C2211 + C2222 + C2233 + C3311 + C3322 + C3333)(ε 2 /9)

An important element of the above methodology is the generation of well equilibrated structures which will be representative of the true polymer structure at the conditions of interest. Given that the glass transition temperature of sPMMA is Tg ≈ 385 K,39 at the conditions of interest (T = 300 K and P = 1 atm), sPMMA is a glassy material. To generate well relaxed glassy polymer configurations at these conditions, the following methodology was followed in the present work: • An initial amorphous packing cell was created using the MAPS 3.2 software package41 which was then subjected to a static structure optimization using a molecular mechanics algorithm to remove atom−atom overlaps. • This initial cell was subjected to long MD simulations at a high enough temperature (T = 500 K) well above the melting point of PMMA, for which relaxation was considerably easier to achieve through a direct application of the MD method. • From the MD simulations at the temperature of 500 K, a relatively large number (approximately 15) of independent, thoroughly relaxed atomistic configurations for the three systems (sPMMA, GS-sPMMA, FGS-sPMMA) were selected and cooled down to 300 K where an additional MD simulation was performed until the density of the system and its potential energy reached a steady-state. Clearly, the configuration obtained at 300 K is trapped near a local minimum of the potential energy, equivalent to a glassy structure. The final structure from this last MD step was then subjected to an exhaustive minimization procedure using the Hessian−free truncated Newton algorithm (hf tn). According to this, at each iteration, a quadratic model of the potential energy is solved by a conjugate gradient inner iteration. The Hessian (second derivatives) of the energy is not formed directly, but approximated in each conjugate search direction by a finite difference directional derivative. Close to an energy minimum, the algorithm behaves like a Newton method exhibiting a quadratic convergence rate to high accuracy.45

⎛ αpT ⎞ ⎛2 ⎞ = U0 − V0⎜ −P + ⎟ ε + (1/2)V0⎜ μ + λ⎟ε 2 ⎝ ⎠ k T ⎠0 3 ⎝ (12)

where ρ0 is the mass density in the undeformed state, cε is the specific heat at constant strain, αp is the volume expansivity, κT is the isothermal compressibility, and γ is the so-called Grüneisen tensor. For the uniaxial tension experiment, the degree of deformation ε is defined as a tensile strain along x, y and z. For example, for uniaxial tension along x, the strain tensor is ⎡ε 0 0⎤ ⎢ ⎥ ε = ⎢0 0 0⎥ ⎣0 0 0⎦

while the expression for the potential energy of the system assumes the form U = U0 + V0(τ11 + ρ0 cεTγ11)ε + (1/2)V0C1111ε 2 = U0 + V0(τ11 + ρ0 cεTγ11)ε + (1/2)V0(2μ + λ)ε 2 (14)

To obtain estimates of the elastic constants we proceed as follows. We fit the numerical data for the potential energy function with a second order polynomial, we average the results over all structures and we compute the quantity (C11 + C12 + C13 + C21 + C22 + C23 + C31 + C32 + C33) /9 = (2/3)μ + λ

(15)

for the case of the uniform hydrostatic compression, or the quantity (1/3)(C11 + C22 + C33) = 2μ + λ

(16)

for the case of uniaxial tension. By combining eqs 15 and 16, we obtain the Lamé constants λ and μ, thus also the elastic constants through eqs 7. We close this section by emphasizing that, overall, the proposed method for modeling the deformation of a glassy polymer structure entails the following assumptions: • The 15 glassy structures (for the neat polymer and its graphene nanocomposites) are viewed in configuration space as an ensemble of mutually inaccessible states of microscopic liquid disorder. Each undeformed microstructure and all deformed structures obtained from it satisfy the requirements of detailed mechanical equilibrium. • Entropic contributions to the elastic coefficients are neglected; only potential energy effects are considered. • The analysis is limited to very small degrees of deformation (up to 1%), since the emphasis is exclusively on the elastic response to deformation.

On each one of the cooled (glassy) configurations brought to a minimum of the potential energy, two different types of deformation experiments are performed, hydrostatic compression and uniaxial tension, and the relevant elastic properties are estimated by arithmetically averaging over all individual configurations subject to deformation. For the uniform hydrostatic compression experiment, the degree of deformation ε is defined as the fractional decrease in volume; thus the strain tensor for this deformation is ⎡−ε /3 0 0⎤ ⎢ ⎥ ε=⎢ 0 −ε /3 0⎥ ⎢⎣ ⎥ 0 0 −ε /3⎦

(13)

4. RESULTS A. Validation−Simulation Predictions at T = 500 K. Before employing the all-atom Dreiding40 force-field to carry out small-strain deformation experiments and compute the

(11)

while the potential energy is expressed as F


Macromolecules

Article

species (C, O and H) next to a graphene sheet surface shown in Figure 4. Conformation of Individual PMMA Chains in the Bulk. In Table 4, we report simulation predictions for the average square end-to-end distance ⟨R2⟩ and the average square radius-ofgyration ⟨Rg2⟩ of PMMA chains in all simulated systems. The square root of ⟨R2⟩ provides a measure of the average distance that the main chain spans from its first backbone carbon atom to the last while Rg = (⟨Rg2⟩)1/2 provides a measure of the average size of the polymer chain. The MD predictions for ⟨R2⟩ and ⟨Rg2⟩ reported in Table 4 have been computed by averaging over all chains in the simulation cell and over many statistically uncorrelated system snapshots at 500 K. Our simulations indicate that the presence of graphene sheets does not dramatically affect the conformational properties of the matrix chains: in all cases, the average values of ⟨R2⟩ and ⟨Rg2⟩ in the pure sPMMA melt and in the sPMMA nanocomposites with GS and FGS are the same within the statistical error. This should be attributed to the small concentration of the two nanocomposites in graphene sheets and to the small size (only 12 Å × 12 Å) of the graphene sheets considered in this work. From the data of Table 4, we can further calculate the value of the characteristic ratio Cn of sPMMA chains in their pure melt at 500 K. The characteristic ratio is defined as

mechanical properties of the chosen sPMMA model system and its graphene nanocomposites, we thoroughly validated it by examining MD simulation predictions for the thermodynamic, structural and conformational properties of the various model systems (pure polymer and nanocomposites) considered here at T = 500 K and P = 1 atm and how they compare with available experimental data52−55 and other simulation studies.22,56−62 Density. Our MD predictions for the average density of the simulated sPMMA, GS-sPMMA and FGS-sPMMA nanocomposites are shown in Table 3. For the pure polymer Table 3. Density of all Three Simulated Model Systems at Two Different Temperatures (500 and 300 K) ρ (g cm−3) system

temperature (K)

sPMMA

500 300 500 300 500 300

GS-sPMMA FGS-sPMMA

current work

experimental values

± ± ± ± ± ±

1.07137 1.17−1.2038

1.065 1.158 1.082 1.169 1.091 1.173

0.001 0.005 0.001 0.007 0.001 0.007

matrix, the predicted density values is 1.065 ± 0.001 g/cm3. Experimental data for a pure PMMA sample at 500 K37,38 suggest a value equal to 1.071 g/cm3. At the same temperature and pressure conditions, the corresponding polymer nanocomposite melts are characterized by slightly higher density values. For example, for the 5.67 wt % GS-sPMMA nanocomposite the MD prediction for the density is ρ = 1.082 ± 0.001 g/cm3 while for the 6.54 wt % FGS-sPMMA one the corresponding prediction is ρ = 1.091 ± 0.001 g/cm3. In Figure 3 (where we show the variation of the local PMMA density profile as a function of the distance normal to graphene sheets) we understand that this increase in the density of the polymer nanocomposites is due to adsorption of sPMMA chains on the two faces of graphene sheets, implying a denser packing of polymer mass in the interfacial zones right next to their faces. The adsorption is equally strong for GS and FGS as can be seen from the values of the peaks in the characteristic (slightly oscillatory) profiles of the local mass density next to graphene sheets. The same (slightly oscillatory) behavior is evident in the profiles of the local number densities of the three atomic

Cn =

⟨R2⟩ nl ̅

2

where n denotes the number of MMA units along an sPMMA chain and l ̅ the effective monomer length. Our MD simulations predict Cn = 6.9 for the pure sPMMA melt at 500 K with n = 15. The experimentally measured values for infinite molecular weight PMMA samples are C∞ = 7−8 for sPMMA at T = 300 K.63,64 Microscopic Structure. More detailed information about local packing and structure over short and long length scales in the simulated sPMMA nanocomposites was obtained by analyzing the element pair distribution functions gαβ(r) (α,β ∈{C, O, H}) for all six species pairs (CC, OO, HH, CO, CH, OH). The functions gαβ(r) receive contributions from atom pairs that belong to the same chain (intramolecular part, gintra αβ ) or to different chains (intermolecular part, ginter αβ ). Typical results intra for ginter αβ are shown in Figure 5 and for gαβ in Figure 6. Because of chain connectivity which constrains pair distances within certain limits (depending on the stiffness of the bond stretching

Figure 3. Profiles of the local mass density normal to a graphene sheet as obtained from the present NPT MD simulations (T = 500 K, P = 1 atm) of the sPMMA nanocomposites for the cases of unfunctionalized (a) and functionalized (b) graphene sheets. The blue dashed line at the zero value of the horizontal axis denotes the average position of the midplane of the graphene sheet. G


Macromolecules

Article

Figure 4. Same as with Figure 3 but for the atomic density profiles. Carbon atoms are shown in red, oxygen atoms in blue, and hydrogen atoms in black.

corresponding PMMA-FGS nanocomposites; to emphasize conformational changes induced by the presence of graphene sheets, configurational sampling for the nanocomposites was limited only within a distance approximately equal to 6 Å from a graphene sheet (the so-called “adsorption layers”). The overall shape of the six different types of dihedrals looks the same but clearly one can spot certain differences between the two cases, referring mainly to the population of torsion angles in trans and gauche states in the different dihedrals analyzed. This is explained by the tendency of sPMMA chains to lie on the graphene sheet (see Figure 8a) in a way that maximizes the contact of their skeletal carbon atoms with those of the graphene sheets. In the case of functionalized graphene sheets, the situation is more complicated, as chains tend to lie on the graphene surface in such a way that they can optimize not only carbon−carbon contacts but also the number of hydrogen bonds that can be formed between their carbonyl groups and the functional hydroxyl groups on the graphene sheet (see Figure 8b). More information about the hydrogen bonds that form between the polar sPMMA chains and the surrounding functionalized graphene sheets is provided in section B below. Static Structure Factor S(q). The Fourier transform of the total pair distribution function gtotal(r) leads to the calculation of the static structure factor S(q). Despite a number of complicating factors and the different weightings at different q regimes in the experiments, simulated S(q) curves can be directly compared against experimental S(q) data obtained either from X-ray diffraction or neutron scattering measurements, thereby providing a means for validating the predicted intermolecular pair distribution functions ginter αβ (r). To account for the different atomic species (carbons, hydrogens, and oxygens) present in the simulated systems, the static structure factor S (q) should be calculated as the Fourier transform of the function H(r) defined as

Table 4. Simulation Predictions for the Average Square Endto-End Distance ⟨R2⟩, the Mean Square Radius-of-Gyration ⟨Rg2⟩, and the Characteristic Ratio Cn=15 of All Simulated Model Systems at 500 K system

⟨R2⟩ (Å2)

⟨Rg2⟩ (Å2)

Cn=15

sPMMA GS-sPMMA FGS-sPMMA

476 ± 17 551 ± 30 492 ± 23

77 ± 3 83 ± 1 81 ± 2

6.92 8.02 7.16

and bond bending potentials, as well as on the temperature and pressure conditions considered), intramolecular contributions in the gintra αβ plots displayed in Figure 6 appear as very sharp peaks.35 We further clarify that in the case of the polymer nanocomposites, the terms intermolecular and intramolecular have been used in a broader sense, meaning that each graphene or each functionalized graphene sheet in the sample was considered as a separate molecule (or chain). At long distances, all intermolecular pair distribution functions approach unity, indicating that no long-range structure exists, which is typical of purely amorphous systems. The intermolecular functions ginter αβ (r) rise more or less smoothly with increasing distance exhibiting some subtle structural features up to ∼10 Å but beyond that distance they approach the asymptotic value of one. In all intermolecular pair distribution functions, several maxima and minima appear except for the O−O correlations where practically only a single well-defined maximum shows up in the corresponding plot. The correlation hole effect (humps in the corresponding ginter αβ (r) plots) is evident in all intermolecular correlations and can be explained by considering not only the first coordination shell defined by the sum of the van der Waals radii of the atoms involved in the correlation but also the fact that different skeletal atoms belonging to different chains are kept apart by the substituents surrounding them. Thus, in addition to the first pronounced hump in all intermolecular correlations (except for the O−O ones) at distances between 3 and 4 Å (corresponding to the first coordination shell), several broad humps of intermolecular origin at longer distances (around 5 and 6 Å) are discerned in all distributions. According to Figure 3, variations in the local density of sPMMA in the vicinity of a graphene sheet persist over distances up to 6 Å from the graphene plane. To check the influence of PMMA−graphene interactions on the microscopic structure of PMMA chains, in Figure 7 we report a comparison of all torsional angles between the neat sPMMA matrix and the

3

H (r ) =

3

∑∑ α=1 β=1

xαxβfα fβ 3

(∑γ = 1 xγfγ )2

(g total (r ) − 1) (17)

where xα and xβ denote the number fraction of α-type and βtype atoms in the system, and fα and fβ the respective scattering factors. Values of xα and fα for the three different elements (C, O, H) in the sPMMA melt and in its graphene-based nanocomposites are reported in Table 5. The results for the function H(r) for all model systems are shown in the Supporting Information. By taking the Fourier transform S(q) 2 =∫∞ 0 4πr H (r)exp(−ikr) dr of the function H(r), we computed H


Macromolecules

Article

Figure 5. Predicted intermolecular pair distribution functions ginter αβ (r) for the six different species pairs (CC, OO, HH, CO, CH, OH) in the pure sPMMA and the GS-sPMMA and FGS-sPMMA nanocomposite melts (T = 500 K and P = 1 atm).

the X-ray diffraction patterns shown in Figure 9. The first peak (at small q values) in the three graphs of Figure 9 represents intermolecular correlations while the second and third ones (at large q values) reflect intramolecular packing. For the neat sPMMA, the MD predictions for the location and magnitude of the major intermolecular peak near q = 0.8 Å−1 and of the intramolecular peaks near 2.0 Å−1 and 2.9 Å−1 agree with the S(q) pattern reported for a pure sPMMA melt at 500 K and 1 atm in ref 62 with the OPLS force-field. Additional results for the dynamic properties (segmental dynamics, terminal relaxation, glass transition temperature) of the model systems addressed here (sPMMA, GS-sPMMA, FGS-sPMMA) will be presented in a forthcoming publication. We close this section on the melt simulations by mentioning that although graphene sheets were well distributed in the matrix and always found to be “well-wetted” by PMMA, in some simulations they come close to each other to form a stack

structure (implying the formation of graphene-rich and graphene-poor regions for the macroscopic nanocomposite system). We will discuss this interesting phenomenon in a forthcoming publication but we clarify here that all results reported in this study have been obtained with uniformly dispersed GS or FGS in the polymer matrix. Another important issue is that of the instantaneous shape of graphene sheets in the host sPMMA matrix. Typical atomistic snapshots are shown in Figure 10 indicating that both GS and FGS are very flexible with their shapes undergoing large fluctuations in the course of the simulations. Such a flexibility of their 2-d shape further contributes to the enhancement of mechanical properties since it allows for better contacts with PMMA segments and thus for stronger adhesion (stronger interfacial interactions). Given the inherent stiffness of a PMMA chain (the persistence length of an sPMMA chain at 500 K is 12.6 ± 0.5 Å as predicted by our simulations), if I


Macromolecules

Article

Figure 6. Same as with Figure 5 but for the intramolecular part of the corresponding total pair distribution functions.

polymer in configuration space as an ensemble of mutually inaccessible states of microscopic liquid disorder; by arithmetic averaging over the population of these 15 model structures, one can then obtain estimates of the macroscopic properties. The first property that we calculated through such an arithmetic averaging was density. Our results for the three model systems addressed here (sPMMA, GS-sPMMA, FGSsPMMA) are reported in Table 3. Experimental values37,38 for the density are available in the literature for PMMA of unspecified tacticity and probably refer to commercial PMMA typically characterized by a rather random sequence of m and r diads. Comparison with the available experimental data suggests that the density of the pure matrix (1.158 ± 0.005 g/cm3) is somewhat lower than the reported one (ρexp = 1.17− 1.19 g/cm3). Given that experimental samples consist most likely of the atactic form, the agreement between the theoretical prediction and the experimentally measured value is excellent.

graphene sheets were also stiff, it would not be easy for PMMA chains to bend on the graphene surface to maximize interfacial interactions. Of course, for larger graphene sheets, we expect a wrinkled topology at the nanoscale which should enhance the mechanical properties even further as it can interlock polymer segments. B. Density and Structural Properties at T = 300 K. From the MD equilibration runs carried out at 500 K, a good number (approximately 15) of independent and thoroughly relaxed atomistic configurations of all polymer systems were selected and suddenly cooled down to 300 K for several nanoseconds (typically for ∼200 ns) until their densities reached constant values. Through this procedure, 15 different glassy structures were generated for each system; these glasses (in which chains are in a state of frozen-in liquid disorder) constituted the ensemble of initial microscopic structures for simulating structure and predicting mechanical properties. That is, following Theodorou and Suter,35 we envision the glassy J


Macromolecules

Article

Figure 7. Predicted torsional angle distribution functions in the pure sPMMA melt and in the interfacial areas of the GS-sPMMA and FGS-sPMMA nanocomposite melts (T = 500 K, P = 1 atm).

between the receiver oxygen and the donor hydrogen is less than 2.4 Å, and (c) the angle between the O−O direction and the molecular O−H direction of the donor (H is the hydrogen atom that creates the bond) is less than 30°. According to the second criterion, in a plot of the radial O−H pair distribution function, the H-bonds should correspond to O−H pairs found at a distance less than 2.4 Å.66 The data in Table 6 demonstrate the formation of a non-negligible number of hydrogen bonds between the functional groups (hydroxyls and oxygen atoms) across the surface of a FGS and PMMA segments, indicating a higher affinity of PMMA chains for FGS compared to GS despite the fact that the local sPMMA density profile is practically the same next to the two types of surfaces (the peaks in the density patterns in Figures 11 and 12 are almost identical). A third property that we calculated through arithmetic averaging was the static structure factor S(q) of the simulated

The second property that we calculated through arithmetic averaging was the density of the host sPMMA matrix in the interfacial region next to graphene sheets. The corresponding results are shown in Figures 11 and 12 and, similar to Figures 3 and 4 for the melt case, indicate strong adsorption of PMMA chains on the graphene sheets. According to Ramanathan et al.,2 the strong adsorption on FGS, particularly, is enhanced by the pendant hydroxyl groups and oxygen atoms across the surface of these sheets that mediate the formation of hydrogen bonds with the carbonyl groups of PMMA (an example is shown in Figure 13). To check this, we carried out a detailed search for hydrogen bonds in the simulated FGS-PMMA nanocomposites and we found the results (referring to the average number of hydrogen bonds per surface area of graphene sheets) shown in Table 6. According to Gordillo and Marti,65 a hydrogen bond forms when (a) the distance between two oxygen atoms is less than 3.6 Å, (b) the distance K


Macromolecules

Article

Figure 9. Simulation predictions for the static structure factor S(q) for all three model systems at T = 500 K (P = 1 atm). For the pure sPMMA melt, we also show the S(q) curve obtained with the OPLS force-field (ref 62) at the same temperature and pressure conditions.

Figure 8. Typical snapshots from the atomistic simulations at 500 K showing the arrangement of an sPMMA chain on the surface of (A) a nonfunctionalized graphene sheet and (B) a functionalized graphene sheet. Figure 10. Atomistic snapshots of the GS and FGS in the host sPMMA matrix in the course of the MD simulations at T = 500 K (P = 1 atm).

polymer and polymer nanocomposites systems using the Fourier transform of the function H(r) introduced in eq 17. As in the melt case, the calculation of H(r) requires knowledge of all (six) element pair distribution functions gCC, gHH, gOO, gCH, gCO, and gHO. To obtain these, we computed the six different element pair distribution functions for each of the selected 15 glassy model structures for the pure polymer and its nanocomposites considered here, and then we averaged the results. From the averaged element pair distribution functions we computed the corresponding H(r) function using the values for the number fraction xα and scattering factor fα for the three different elements (C, O, H) reported in Table 5. Plots of the function H(r) obtained from the above calculations for the three model systems are displayed in the Supporting Information. The calculated spectrum in Figure 14 exhibits three peaks at 0.86, 2.13, and 2.95 Å−1. These locations are not far from the experimentally observed peaks at 0.91, 1.98, and 3.02 Å−1 (the other two sharp experimental peaks at 2.7 and 3.1 Å−1 are due to the aluminum sample holder according to Farago et al.55 and

should be ignored). Experimentally, the second peak is the major one while the first is only a shoulder. The major difference between prediction and experiment lies in the intensity of the signal at lowest q that is much stronger in the calculated spectra. We conclude that the diffraction results indicate reasonable agreement between simulated and actual structure. Ramanathan et al.2 have reported XRD patterns for neat PMMA and for a PMMA-FGS nanocomposite using CuKα radiation (λ = 0.15 nm). For both samples, the measured XRD patterns showed peaks at around 2θ = 15° which are indicative of the intrinsic ordering of the amorphous polymer. Furthermore, no graphite layer structure peak at 26° or graphite oxide structure peak at 13° was seen, which suggested that the dispersion of FGS in PMMA in their samples was close to the single-sheet level. In the measured XRD pattern for the FGS-PMMA nanocomposite, the first peak was located at a

Table 5. Values of the Number Fraction xα and Scattering Factor fα for the Three Different Elements (C, O, H) in the Simulated Model Systems scattering factor fα (10−14 m)

number fraction xα system

C

O

H

C

O

H

sPMMA GS-sPMMA FGS-sPMMA

0.3304 0.3477 0.3456

0.1322 0.1267 0.1297

0.5374 0.5256 0.5247

0.665

0.5803

0.667

L


Macromolecules

Article

Figure 11. Profiles of the local mass density normal to a graphene sheet as obtained from the present NPT MD simulations at T = 300 K (P = 1 atm) of the sPMMA nanocomposites for the cases of unfunctionalized (a) and functionalized (b) graphene sheets. The blue dashed line at the zero value of the horizontal axis denotes the average position of the midplane of the graphene sheet.

Figure 12. Same as with Figure 11 but for the atomic density profiles. Carbon atoms are shown in red, oxygen atoms in blue, and hydrogen atoms in black.

Table 6. Number of Hydrogen Bonds Formed in the Simulated FGS-sPMMA Nanocomposites per Unit Area of Graphene Sheet system average number of hydrogen bonds formed per nm2 of graphene sheet

FGS 1 FGS 2 FGS 3 0.29

0.32

0.30

Figure 14. Simulation predictions for the static structure factor S(q) for all three model systems at T = 300 K (P = 1 atm). For the pure sPMMA melt, we also show the experimentally measured X-ray S(q) pattern (ref 55).

Figure 13. Example of a situation where three hydrogen bonds are created between a FGS and the surrounding sPMMA chains. We can observe the formation of two hydrogen bonds with the same chain on the one side of the FGS and of one hydrogen bond with a different chain on the other side of the FGS.

different position (at around 2θ = 180°) than the corresponding one for the neat PMMA (at around 2θ = M


Macromolecules

Article

130°). Our calculations indicate only a slight difference in the location of the first peak between the neat sPMMA matrix and the two nanocomposites (GS-sPMMA and FGS-sPMMA). C. Mechanical Properties at T = 300 K. Our ensemble of 15 equilibrium structures was used next to predict the elastic constants of glassy amorphous sPMMA and its graphene-based nanocomposites by simulating mechanical deformation of various types and degrees. For these computations, the relevant quantity is the potential energy, Upot min, of a local minimum in the system’s configuration space. On the basis of the analysis reviewed in section 3, by monitoring changes in this quantity as a function of deformation one can predict the mechanical properties using as initial structure a member of the ensemble of 15 structures described earlier. An important observation worth mentioning here is that if one starts from a deformed model structure and reverses the deformation one reaches the original undeformed model structure; within the strain range studied here (ε of the order of 10−3), the response of all model systems to deformation was completely reversible (or “elastic”). The detailed mechanical equilibrium requirement seems to have a unique solution around the undeformed state within the range of degrees of deformation studied. This uniqueness characteristic may cease to hold for considerably larger values of ε; discontinuous changes and irreversibilities in the minimum potential energy trajectory in configuration space might indicate the onset of plastic deformation. In Figure 15, we show results for the change in the potential energy as a function of the degree of deformation for several initial structures. Separate results are shown for compression and tension. By fitting (see Supporting Information) the curves with second order polynomials, as suggested by eqs 12 and 14, we computed the values of the two linear combinations of the Lamé constants, (2/3 μ + λ) in compression and (2 μ + λ) in tension, from which we extracted the elastic properties. By averaging over all 15 structures subjected to compression and tension modes of deformation, we obtained the predictions (average values and standard deviations) for the diagonal elastic coefficients and bulk modulus that are listed in Table 7. The large error bars are due to (a) the relatively small size of the model systems, which induces considerable variation of the properties in different directions within a given structure, and (b) to significant variations from structure to structure for the same property in a given direction. But these error bars can become as small as possible by averaging over more model structures or by using additional modes of deformation (e.g., shear) along all possible directions. We further mention that, in all cases, excellent fits with quadratic functions were obtained, indirectly verifying that for all systems the response of the corresponding ensemble of model structures conformed to that of an isotropic elastic solid. In Table 7, we report the predicted values of the four elastic constants (E, B, G, ν) for the neat sPMMA matrix and its GS and FGS nanocomposites. In the literature, there exist several reports of experimentally measured data for the mechanical properties of glassy PMMA.67−73 For example, Mott et al.67 have determined the bulk B and shear modulus G of PMMA as a function of temperature from longitudinal and shear wave speeds at 1 MHz (and from these data, they further calculated the Poisson ratio ν). At 25 °C, they found (see Figure 1 in ref 67): B = 6.13 GPa, G = 1.78 GPa, and ν = 0.35. A value for the elastic modulus E of (neat) PMMA is reported in a recent study by Zhi et al.68 on the mechanical and thermal properties of PMMA-boron nitride nanotube composites: E = 2.07 GPa. Afifi

Figure 15. Plots of the total potential energy divided by twice the volume of the undeformed state (2U/V0) versus fractional increase in the volume of the simulation cell for eight different atomistic structures from the hydrostatic compression-tension experiments. Part A of the figure refers to the pure sPMMA matrix, part B to the GS-sPMMA- nanocomposite, and part C to the FGS-sPMMA nanocomposite (T = 300 K, P = 1 atm).

and Hasan69 have also reported values of the shear G and Young’s modulus E from measurements of ultrasonic velocities (taken at 2 MHz ultrasonic frequency using the pulse echo method) of both longitudinal and shear waves in thermoplastic discs of PMMA, as a function of annealing temperature. At 25 °C, they find (see Figure 4 in ref 69): E = 4.80 GPa and G = 1.80 GPa. Fukuhara and Sampei71 have reported data for Young’s modulus, the shear modulus, the bulk modulus, the N


Macromolecules

Article

Table 7. Predicted Values of the Elastic Constants for all Simulated Systems as Computed from the Present Study lamé constants simulation

elastic constants

system

μ

λ

E (GPa)

B (GPa)

G (GPa)

ν

pure sPMMA GS-sPMMA FGS-sPMMA

1.3 1.4 2.3

3.4 3.7 3.5

3.4 3.8 6.0

4.2 4.7 5.0

1.3 1.4 2.3

0.36 0.36 0.31

Young’s modulus E, of the bulk modulus B, and of the shear modulus G for all three simulated systems on a relative scale. Adding 5.67 wt % nonfunctionalized graphene sheets in PMMA results in an improvement of Young’s modulus E by 12%, of the bulk modulus B by 11%, and of the shear modulus G by 11%. The improvement is more striking in the case of the sPMMA matrix containing functionalized graphene sheets: the FGSsPMMA matrix with a 6.54 wt % concentration in FGS shows an enhancement in Young’s modulus of elasticity E by 75%, in the bulk modulus B by ∼19%, and in the shear modulus G by ∼83%, compared to the corresponding neat matrices. Experimentally, Ramanathan et al.2 have reported an 80% enhancement of Young’s modulus E for a 1 wt % FGS-PMMA nanocomposite material while Li and McKenna16 have measured a 33% enhancement of Young’s modulus E for a 0.005% vol. graphene oxide/PMMA nanocomposite material. The higher concentration for which a similar enhancement is observed in our calculations (6.54 wt % concentration in FGS) should be attributed to the relatively small size of GS considered in our simulations (12 Å × 12 Å only). To check this, we are currently extending our FGS-sPMMA nanocomposite simulations to systems containing considerably larger graphene sheets (compared, e.g., to the size of the matrix PMMA chains). Results from these simulations will be reported in a forthcoming publication.

Lamé parameters, and the Poisson ratio of PMMA as a function of temperature using the ultrasonic pulse method. At 25 °C, they find (see Figures 2 and 3 in ref 71): E = 6.09 GPa, B = 6.17 GPa, G = 2.30 GPa, and ν = 0.34). Additional measurements have been performed by Weishaupt et al.70 (E = 6.20 GPa, B = 5.90 GPa, G = 2.30 GPa, ν = 0.34), Yee and Takemori72 (E = 4 GPa, ν = 0.36), and Kono73 (B = 6.04 GPa, G = 2.08 GPa, ν = 0.35). We have collected all these data and compiled them in the form shown in Table 8. A close look at Table 8. Compilation of Literature Data for the Elastic Constants of PMMA at 25 °C Based on Direct Experimental Measurements (References 67−73) elastic constant ref

Young’s modulus (E) (GPa)

bulk modulus (B) (GPa)

shear modulus (G) (GPa)

Poisson’s ratio (ν)

67 68 69 70 71 72 73

− 2.07 4.80 6.20 6.09 4.00 −

6.13 − − 5.90 6.17 − 6.04

1.78 − 1.80 2.30 2.30 − 2.08

0.35 − − 0.34 0.34 0.36 0.35

the data suggests that the value of E is between 2.07 and 6.20 GPa, B varies between 5.90 and 6.17 GPa, G is between 1.78 and 2.30 GPa, and the Poisson ratio ν is approximately 0.35. We conclude that our predicted values of the elastic constants for sPMMA are in very favorable agreement with the experimentally reported values. The addition of GS, and especially of FGS, to the PMMA matrix greatly improves its elastic properties. This can be seen more clearly in Figure 16, where we have plotted the values of

5. SUMMARY Using detailed MD simulations, we have successfully simulated the structural and conformational properties in the molten state (T = 500 K) of a model sPMMA structure and its nanocomposites with nonfunctionalized and functionalized graphene sheets. In a second step, 15 fully equilibrated and totally uncorrelated structures from the MD simulations at 500 K were chosen and cooled down to 300 K to generate an equal population of independent model glassy structures of the PMMA matrix and its nanocomposites with nonfunctionalized and functionalized graphene sheets. These glassy structures were in all cases amorphous, as could be verified by the corresponding plots of the relevant pair distribution functions; this is particularly important for sPMMA which can crystallize. In a third step, and using techniques already developed in the past for simpler systems (e.g., glassy vinyl polymers such as polypropylene and polystyrene), we were able to compute the mechanical properties of the PMMA nanocomposites. This has been achieved by subjecting a selected number of independent glassy microstructures to small strain deformation following the method originally developed by Theodorou−Suter35 for glassy polypropylene. Our analysis allowed us to carry out a comparative study of the mechanical properties of sPMMA in the absence and presence of uniformly dispersed graphene sheets at relatively low loadings. Our calculations confirmed experimental investigations that ultralow loadings of FGS in PMMA can have a great effect in the mechanical properties of the host matrix. As correctly noted by Ramanathan et al.,2 for PMMA this is due to its stronger adsorption on FGS. FGS

Figure 16. Summary of property improvement for the elastic constants of GS-sPMMA and FGS-sPMMA nanocomposites (T = 300 K, P = 1 atm). All values have been normalized to the values for the neat sPMMA matrix at the same conditions: E = 3.45 GPa, B = 4.22 GPa, Gc= 1.27 GPa, and ν = 0.36. O


Macromolecules

■

contains pendant hydroxyl groups across their surfaces which form hydrogen bonds with the carbonyl groups of PMMA. The surface chemistry results in stronger interfacial interactions with the polymer matrix and thus has a substantially larger influence on the properties of the host polymer. As we will show in a forthcoming publication, PMMA adsorption on graphene sheets is accompanied by a non-negligible increase in the corresponding value of the glass transition temperature of the nanocomposite. Our computations have shown that the addition of small fractions of graphene sheets leads to substantial improvement of these properties, especially when functionalized graphene sheets are used. However, we observed that to get similar levels of enhancement in (e.g.) the Young modulus E as those observed experimentally, a larger concentration of FGS should be used in the simulation cells (compared to experimental data). We have tentatively attributed this to the small size of FGS and the relatively short chain length of polymer chains considered in this work. More simulations are currently underway with significantly larger FGs and longer PMMA chains to address the sensitivity of the predicted elastic constants on the size of graphene sheets and the MW of the polymer matrix.

■

REFERENCES

(1) Ramanathan, T.; Stankovich, S.; Dikin, D. A.; Liu, H.; Shen, H.; Nguyen, S. T.; Brinson, L. C. J. Polym. Sci. Polym. Phys. 2007, 45, 2097. (2) Ramanathan, T.; Abdala, A. A.; Stankovich, S.; Dikin, D. A.; Herrera-Alonso, M.; Piner, R. D.; Adamson, D. H.; Schniepp, H. C.; Chen, X.; Ruoff, R. S.; Nguyen, S. T.; Aksay, I. A.; Prud’homme, R. K.; Brinson, L. C. Nat. Nanotechnol. 2008, 3, 327. (3) George, J. J.; Bhowmick, A. K. J. Mater. Sci. 2008, 43, 702. (4) Villar-Rodil, S.; Paredes, J. I.; Martinez-Alonso, A.; Tascon, J. M. D. J. Mater. Chem. 2009, 19, 3591. (5) Fang, M.; Wang, K. G.; Lu, H. B.; Yang, Y. L.; Nutt, S. J. Mater. Chem. 2009, 19, 7098. (6) Lee, Y. R.; Raghu, A. V.; Jeong, H. M.; Kim, B. K. Macromol. Chem. Phys. 2009, 210, 1247. (7) Eda, G.; Chhowalla, M. Nano Lett. 2009, 9, 814. (8) Kim, H.; Abdala, A. A.; Macosko, C. W. Macromolecules 2010, 43, 6515. (9) Yoonessi, M.; Gaier, J. R. ACS Nano 2010, 4, 7211. (10) Cai, D. Y.; Song, M. J. Mater. Chem. 2010, 20, 7906. (11) Kuilla, T.; Bhadra, S.; Yao, D. H.; Kim, N. H.; Bose, S.; Lee, J. H. Prog. Polym. Sci. 2010, 35, 1350. (12) Potts, J. R.; Dreyer, D. R.; Bielawski, C. W.; Ruoff, R. S. Polymer 2011, 52, 5. (13) Gong, L.; Young, R. J.; Kinloch, I. A.; Riaz, I.; Jalil, R.; Novoselov, K. S. ACS Nano 2012, 6, 2086. (14) Compton, O. C.; Cranford, S. W.; Putz, K. W.; An, Z.; Brinson, L. C.; Buehler, M. J.; Nguyen, S. T. ACS Nano 2012, 6, 2008. (15) Zaman, I.; Kuan, H. C.; Meng, Q. S.; Michelmore, A.; Kawashima, N.; Pitt, T.; Zhang, L. Q.; Gouda, S.; Luong, L.; Ma, J. Adv. Funct. Mater. 2012, 22, 2735. (16) Li, X. G.; McKenna, G. B. ACS Macro Lett. 2012, 1, 388. (17) Li, Y.; Kroger, M.; Liu, W. K. Macromolecules 2012, 45, 2099. (18) Ovid’ko, I. A. Rev. Adv. Mater. Sci. 2013, 34, 19. (19) Crock, C. A.; Rogensues, A. R.; Shan, W. Q.; Tarabara, V. V. Water Res. 2013, 47, 3984. (20) Ganesan, V.; Jayaraman, A. Soft Matter 2014, 10, 13. (21) Yuan, B. H.; Bao, C. L.; Song, L.; Hong, N. N.; Liew, K. M.; Hu, Y. Chem. Eng. J. 2014, 237, 411. (22) Metatla, N.; Soldera, A. Macromol. Theor. Simul. 2011, 20, 266. (23) Rahman, R. J. Appl. Phys. 2013, 113, 24503. (24) Stephanou, P. S.; Mavrantzas, V. G.; Georgiou, G. C. Macromolecules 2014, 47, 4493. (25) Duplock, E. J.; Scheffler, M.; Lindan, P. J. D. Phys. Rev. Lett. 2004, 92, 225502. (26) Schniepp, H. C.; Li, J. L.; McAllister, M. J.; Sai, H.; HerreraAlonso, M.; Adamson, D. H.; Prud’homme, R. K.; Car, R.; Saville, D. A.; Aksay, I. A. J. Phys. Chem. B 2006, 110, 8535. (27) McAllister, M. J.; Li, J. L.; Adamson, D. H.; Schniepp, H. C.; Abdala, A. A.; Liu, J.; Herrera-Alonso, M.; Milius, D. L.; Car, R.; Prud’homme, R. K.; Aksay, I. A. Chem. Mater. 2007, 19, 4396. (28) Geim, A. K.; Novoselov, K. S. Nat. Mater. 2007, 6, 183. (29) Rissanou, A. N.; Harmandaris, V. J. Nanopart. Res. 2013, 15, 1589. (30) Rissanou, A. N.; Harmandaris, V. Soft Matter 2014, 10, 2876. (31) Ju, S. P.; Wang, Y. C.; Huang, G. J.; Chang, J. W. RSC Adv. 2013, 3, 8298. (32) Awasthi, A. P.; Lagoudas, D. C.; Hammerand, D. C. Model. Simul. Mater. Sci. Eng. 2009, 17, 015002 . (33) Montazeri, A.; Rafii-Tabar, H. Phys. Lett. A 2011, 375, 4034. (34) Wang, M. C.; Lai, Z. B.; Galpaya, D.; Yan, C.; Hu, N.; Zhou, L. M. J. Appl. Phys. 2014, 115, 123520. (35) Theodorou, D. N.; Suter, U. W. Macromolecules 1986, 19, 139. (36) Rapold, R. F.; Suter, U. W.; Theodorou, D. N. Macromol. Theor Simul 1994, 3, 19. (37) Brandrup, J.; Immergut, E. H.; Grulke, E. A. Polymer handbook, 4th ed.; Wiley: New York and Chichester, U.K., 2004. (38) Mark, J. E. Physical properties of polymers handbook. AIP Press: Woodbury, NY, 1996.

ASSOCIATED CONTENT

S Supporting Information *

Simulation details, tables of Dreiding force-field parameters, Coulombic potential function parameters, and best-fit values of the 2nd order coefficient, and figures showing the definition of the different atom types and plots of the static function H(r). This material is available free of charge via the Internet at http://pubs.acs.org.

■

Article

AUTHOR INFORMATION

Corresponding Author

*Telephone: +30-6944-602580. E-mail: vlasis@chemeng. upatras.gr (V.G.M.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The research leading to these results has received funding from Greek national funds through the Operational Program “Competitiveness and Entrepreneurship, EPAN II” (General Secretariat of Research and Technology, Ministry of Development, Greece) from the project SYNERGASIA_MEKKA, Code Number 09SYN-42-620. Our work was also supported by the LinkSCEEM-2 project, funded by the European Commission under the Seventh Framework Programme through Capacities Research Infrastructure, INFRA-2010-1.2.3 Virtual Research Communities, Combination of Collaborative Project and Coordination and Support Actions (CP-CSA) under Grant Agreement No. RI-261600. The authors acknowledge that the developments outlined in this paper have been achieved with the assistance of high performance computing resources provided by Cy-Tera/LinkSCEEM-2, based in Cyprus. The assistance of Ms. Thekla Loizou from Cy-Tera/LinkSCEEM-2 in achieving the technical requirements is gratefully acknowledged. P


Macromolecules

Article

(39) Shin, H. S.; Jung, Y. M.; Oh, T. Y.; Chang, T. Y.; Kim, S. B.; Lee, D. H.; Noda, I. Langmuir 2002, 18, 5953. (40) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem-Us 1990, 94, 8897. (41) See http://scienomics.com/ for molecular modeling and simulation software. (42) Nosé, S. Prog. Theor Phys. Supp 1991, 130, 1. (43) Hoover, W. G. Phys. Rev. A 1986, 34, 2499. (44) Eastwood, J. W.; Hockney, R. W.; Lawrence, D. N. Comput. Phys. Commun. 1980, 19, 215. (45) The LAMMPS molecular dynamics simulator can be downloaded from http://www.sandia.gov/∼sjplimp/lammps.html. (46) Verlet, L. Phys. Rev. 1967, 159, 98. (47) Verlet, L. Phys. Rev. 1968, 165, 201. (48) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids. Clarendon Press and Oxford University Press: Oxford, England, and New York, 1987. (49) Mansfield, K. F.; Theodorou, D. N. Macromolecules 1991, 24, 4295. (50) Daoulas, K. Ch.; Harmandaris, V. A.; Mavrantzas, V. G. Macromolecules 2005, 38, 5780. (51) Lempesis, N.; Vogiatzis, G. G.; Boulougouris, G. C.; van Breemen, L. C. A.; Hütter, M.; Theodorou, D. N. Mol. Phys. 2013, 111, 3430. (52) Ward, D. J.; Mitchell, G. R. Phys. Scr. 1995, T57, 153. (53) Lovell, R.; Windle, A. H. Polymer 1981, 22, 175. (54) Miller, R. L.; Boyer, R. F.; Heijboer, J. J. Polym. Sci., Polym. Phys. 1984, 22, 2021. (55) Farago, B.; Chen, C. X.; Maranas, J. K.; Kamath, S.; Colby, R. H.; Pasquale, A. J.; Long, T. E. Phys. Rev. E 2005, 72, 031809. (56) Apel, U. M.; Hentschke, R.; Helfrich, J. Macromolecules 1995, 28, 1778. (57) Soldera, A. Macromol. Symp. 1998, 133, 21. (58) Tsige, M.; Taylor, P. L. Phys. Rev. E 2002, 65, 021805. (59) Jawalkar, S. S.; Adoor, S. G.; Sairam, M.; Nadagouda, M. N.; Aminabhavi, T. M. J. Phys. Chem. B 2005, 109, 15611. (60) Genix, A. C.; Arbe, A.; Alvarez, F.; Colmenero, J.; Schweika, W.; Richter, D. Macromolecules 2006, 39, 3947. (61) Genix, A. C.; Arbe, A.; Alvarez, F.; Colmenero, J.; Farago, B.; Wischnewski, A.; Richter, D. Macromolecules 2006, 39, 6260. (62) Chen, C. X.; Maranas, J. K.; Garcia-Sakai, V. Macromolecules 2006, 39, 9630. (63) Sundarar, P.; Flory, P. J. J. Am. Chem. Soc. 1974, 96, 5025. (64) Cowie, J. M. G. Polymers: chemistry and physics of modern materials; Intertext: Aylesbury, U.K., 1973. (65) Gordillo, M. C.; Marti, J. Chem. Phys. Lett. 2000, 329, 341. (66) Anastassiou, A.; Karahaliou, E. K.; Alexiadis, O.; Mavrantzas, V. G. J. Chem. Phys. 2013, 139, 164711. (67) Mott, P. H.; Dorgan, J. R.; Roland, C. M. J. Sound Vibr. 2008, 312, 572. (68) Zhi, C. Y.; Bando, Y.; Wang, W. L. L.; Tang, C. C. C.; Kuwahara, H.; Golberg, D. J. Nanomater 2008, 2008, 642036. (69) Afifi, H.; Hasan, E. Polym.-Plast. Technol. 2003, 42, 543. (70) Weishaupt, K.; Krbecek, H.; Pietralla, M.; Hochheimer, H. D.; Mayr, P. Polymer 1995, 36, 3267. (71) Fukuhara, M.; Sampei, A. J. Polym. Sci., Polym. Phys. 1995, 33, 1847. (72) Yee, A. F.; Takemori, M. T. J. Polym. Sci., Polym. Phys. 1982, 20, 205. (73) Kono, R. J. Phys. Soc. Jpn. 1960, 15, 718.

Q


Determination of the Mechanical Properties of a Poly(methyl

Recommend Documents