Received 11 March 2016; accepted 26 April 2016; published 29 April 2016

ABSTRACT

Modeling and simulation allow methodical variation of material properties beyond the capacity of experimental methods. The polymers are one of the most commonly used matrices of choice for composites and have found applications in numerous fields. The stiff and fragile structure of monolithic polymers leads to the innate cracks to cause fracture and therefore the engineering applications of monolithic polymers, requiring robust damage tolerance and high fracture toughness, are not ubiquitous. In addition, when “many-parts” cling together to form polymers, a labyrinth of molecules results, which does not offer to electrons and phonons a smooth and continuous passageway. Therefore, the monolithic polymers are bad conductors of heat and electricity. However, it is well established that when polymers are embedded with suitable entities especially nano-fillers, such as metallic oxides, clays, carbon nanotubes, and other carbonaceous materials, their performance is propitiously improved. Among various additives, graphene has recently been employed as nano-filler to enhance mechanical, thermal, electrical, and functional properties of polymers. In this review, advances in the modeling and simulation of grapheme based polymer nanocomposites will be discussed in terms of graphene structure, topographical features, interfacial interactions, dispersion state, aspect ratio, weight fraction, and trade-off between variables and overall performance.

Keywords:

Graphene, Polymer, Nanocomposites, Modeling, Simulation

1. Introduction

A comparison of the properties of polymers with other materials is shown in Figure 1. The polymers show relatively lower values of strength, stiffness, and fracture toughness. Therefore, the applications of monolithic polymers are not that common due to these limitations. However, it was realized that by the time polymers are reinforced, they show a propitious flexibility to manipulate their properties. The Polymer Matrix Composites (PMCs) were the first manufactured composite materials during World War II [1] . The applications of PMCs started on larger scale during 1950s while during 1960s, the applications in consumer sporting equipment extensively increased their market. During the energy crises in 1970s, the demand for PMCs further increased due to their relatively lower cost. During this period, the properties of PMCs, the design methodologies, and manufacturing capabilities were further modified. With the improvement in other properties, the maximum service temperature for PMCs was also increased. The PMCs have replaced aerospace structural alloys such as aluminum alloys as the latter are prone to corrosion and fatigue failure [1] . The high-performance PMCs are increasingly finding demand in the infrastructure commodity market.

Some of the most commonly used polymers include epoxy, polycarbonate, polyester, Polymethyl Methacrylate (PMMA), polyethylene, polypropylene, polyvinyl chloride, Polystyrene (PS), nylon, and teflon etc. The epoxies are most commonly used to produce PMCs. The polymers are reinforced with various kinds of reinforcements. The first manufactured PMCs used glass fibers as reinforcement. Besides glass fibers, some frequently used reinforcements include metallic oxides [2] - [4] , clays [5] - [7] , Carbon Nanotubes (CNTs) [8] - [10] , and other carbonaceous materials [11] - [13] . The trends of the use of various reinforcements in the last decade are shown in Figure 2. The use of titania (TiO₂) is on rise because of its applications in paints, varnishes, and paper industry [14] [15] . However, the number of publications in 2014 with other reinforcements is decreased. Even the state-of-the-art CNTs are being relinquished by the researchers in 2014. It might be suspected that either the global recession has hit the worldwide research or the research has reached its limits. However, Figure 3 provides a quantitative assessment on the popularity of graphene among researchers.

There is an exponential rise in the use of graphene in the last decade. Novoselov et al. [16] experimentally produced Single Layer Graphene (SLG) using scotch-tape method and a low cost method to synthesize graphene on large scale was presented by Stankovich et al. [17] allowing the extension of the applications of graphene to industrial scale. After the groundbreaking experiments on the two-dimensional material graphene by Nobel Laureates, Sir Andre Geim and Konstantin Novoselov [16] from the University of Manchester, graphene came into limelight in research field mainly because of its excellent electrical [18] , thermal [19] , and mechanical properties [20] . Graphene found widespread applications in electronics [21] , bio-electric sensors [22] , energy technology [23] , lithium batteries [24] , aerospace [25] , bio-engineering [26] , and various other fields of nanotechnology [27] . In 2014, 2009 research papers were published on graphene in which about 830 papers were on graphene-epoxy nanocomposites produced using solution casting technique. It shows that epoxy and solution casting technique are still preferred choices as polymer matrix and production method, respectively. It can also be deduced that the researchers are at the moment not interested in introducing some new polymer matrix or to invent a new production method to produce Graphene Based Polymer Nanocomposites (GBPNCs). Instead, research has been carried out to explore the role of graphene in influencing the polymer properties both theoretically and experimentally. To the best of the authors’ knowledge, no article is yet published in which SLG was used in polymer matrices to produce nanocomposites. One reason could be that the SLG is not yet at everyone’s disposal. However, the trends indicate that the graphene would surpass other reinforcements no sooner it becomes easily available. Therefore, extensive theoretical research has been carried out in the last decade to justify the use of graphene as reinforcement in polymers. Although various reviews are available on graphene and polymers, however, no article is yet published which delineates the advances in the theoretical studies of GBPNCs. Therefore, this review article discusses various modeling and simulation approaches employed to correlate between graphene structure, topography, weight fraction, dispersion state, interfacial interactions, and overall performance of GBPNCs.

2. Modeling of Graphene Based Polymer Nanocomposites

Various theoretical and computational approaches have been employed to explore the effect of graphene as reinforcement on the performance of polymer nanocomposites including but not limited to, quantum mechanical-based methods [28] , Continuum Mechanics (CM) [29] , Molecular Mechanics (MM) [30] , Molecular Dynamics (MD) [31] , atomistic modeling [32] , Density Functional Theory (DFT) [33] , and multiscale modeling [34] . The mechanical, thermal, and electrical properties of GBPNCs have widely been studied. Cao [35] has reviewed the atomistic studies on the mechanical properties of graphene and Allegra et al. [36] have reviewedthe modeling of polymer nanocomposites reinforced with spherical nanoparticles or statistically isotropic aggregates. The high strength and stiffness of graphene significantly improve the mechanical properties of polymer nanocomposites. Cho et al. studied the mechanical properties of graphene-epoxy nanocomposites with a combination of MM and Mori-Tanaka Method (MTM) [37] . The properties of nanocomposites are also dependent on the polymer chain networking. Yarovsky and Evans studied the structure and properties of crosslinked polymers using computer simulations [38] . Awasthi et al. used MD simulation to study the load-transfer mechanism at nanoscale in graphene-polyethylene nanocomposites [39] . Due to the semblance of graphene, it is considered as frame-like structure. In the frame, covalent C-C bonds are taken as beams joined together with carbon atoms placed at the joints. Uniaxial beam elements, defined by their cross-sectional area, material properties, and moment of inertia represent the covalent bonds. The parameters of the beam elements are determined by establishing equivalence between structural mechanics and computational chemistry. To generalize the models for GBPNCs, graphene can be wrapped within a “model material” which represents various kinds of polymer matrices as shown in Figure 4. The model material allows systematic manipulation of influencing factors, and accordingly, it provides a tunable elastic platform to determine the overall performance of GBPNCs. The model material allows methodical variation of material properties beyond the capacity of experimental methods.

3. Atomic Structure

A schematic diagram of a hexagonal unit cell of graphene structure is shown in Figure 5. As graphene is a two- dimensional structure, each carbon atom can undergo chemical reaction from the sides. It is one of the reasons for high chemical reactivity of graphene. The carbon atoms on the edge of SLG have three incomplete bonds that impart especially high chemical reactivity to edge carbon atoms. In addition, defects within graphene sheet

are high energy sites and preferable localities for chemical reactants to attack. All these factors make graphene a very highly chemically reactive entity. The carbon atoms are connected through strong covalent bonds. There is sp² orbital hybridization between P_x and P_y that forms σ-bond [40] . The orbital P_z forms π-bond with half-filled band that allows free motion of electrons. These interatomic forces play a crucial part in defining the proficiency of graphene as reinforcement in polymers.

As atoms keep oscillating with characteristic oscillating frequencies, the oscillatory period of carbon atoms for some bounded trajectories reaches terahertz frequencies [41] . To define the trajectory of atoms, Newton’s equation of motion in MD can be used as given in Equation (1) and the relation between force f₁ and potential

energy E is given in Equation (2), where m_I, r_I, and f_I are mass of, position of, and force on I^th atom, respectively, and is potential energy of the N-atom system. The potential energy arises from the bonded

and non-bonded interactions among atoms and molecules. Polymer Consistent Force Field (PCFF) can be used to specify such interactions [42] . The PCFF defines the potential energy components as given in Equation (3) [43] . The first four terms of Equation (3) represent bond, angle, torsion, and dihedral energy, respectively. The other six terms represent cross-terms among various valance energy components, mainly bond and angle. The last two terms define energy arising from Columbic and Van der Walls interactions, respectively. The influence of thermal vibrations on potential energy is addressed by the cross-terms. The total energy of N atoms is represented by Hamiltonian function (H) which can be defined in terms of potential energy (Equation (3)) and is given in Equation (4) while the kinetic energy component is given in Equation (5) [43]

(1)

(2)

(3)

(4)

(5)

where P_I and M_I are momentum and mass of I^th atom, respectively. The energetics of the non-bonding interactions between all pairs of atoms can be defined using Lennard Jones potential (LJ-potential) as given in Equation (6) [44] [45] ,

(6)

where r_ij is the separation distance between atoms i and j, and r_cut is the cut-off distance. For bonding interactions, the bond-stretch and bond-bending potentials are given by Equation (7) and Equation (8), respectively, and the dihedral potential is given by Equation (9) where Q_ijkl is the angle formed by i, j, k, and ℓ bonds [45] .

(7)

(8)

(9)

In continuum-nanoscopic approach, also known as quasi-continuum method, every point in the continuum is considered as an extensive region wherein atomistic degrees of freedom are explicitly taken into account [45] . Therefore, when subjected to deformation, that nanoscopic atomistic structure provides the constitutive response of the material. This nanoscopic atomistic structure is also termed as Nanoscopic Representative Volume Element (NRVE). The standard kinematic assumption to relate the atomistic and continuum deformations is the Cauchy-Born rule [46] [47] as expressed in Equation (10) [48] ,

(10)

where and are the distance vectors before and after deformation between the pair of atoms i and j, respectively. The deformation gradient tensor is given by Equation (11) [49] ,

(11)

where x_ij, y_ij, and z_ij are the displacements of the local nodes in an arbitrary tetrahedral element. The bond length r_ij is given by Equation (12) [45] ,

(12)

where F^T is the transpose of F-matrix. The Green-Lagrange strain tensor E is given by Equation (13) [50] ,

(13)

where I is the identity tensor. From Equation (10) and Equation (12), the bond length in terms of Green-La- grange strain tensor is given by Equation (14) and Equation (15) [50] ,

(14)

(15)

where is the unit vector in the initial bond orientation between atoms i and j. The relationships for non- bonding pair distance are given by Equations (16)-(18) [45] ,

(16)

(17)

(18)

where r_L, r_0L, and are the stretched and un-stretched nonbonding length, and unit vector in the initial bond orientation for bond L, respectively. The cosine of the bond-angle can be written as Equation (19) and Equation (20) [45] ,

(19)

(20)

where θ_ijk is the angle between triplet of atoms i, j, and k. The consine of dihedral can also be written as Equation (21) where r_ij, r_jk, and r_kl are the distances between atoms i and j, j and k, and k and ℓ, respectively, and Q_ijkl is the dihedral angle between the quartet of i, j, k, and ℓ atoms. Equation (21) can take the form of Equation (22) where the values of α, β, M, and N are defined in Equations (23)-(26), respectively [45] .

(21)

(22)

(23)

(24)

(25)

(26)

The graphene has many other stable configurations apart from honeycomb hexagonal lattice structure. For example, different chemisorbed configurations of epitaxial graphene coexist on single crystal Ni (111) such as top-fcc, top-hcp, and top-bridge [51] . In addition, there are structural defects [52] [53] in graphene which significantly influence its physical and chemical properties [54] . Some of the atomic scale defects are schematically shown in Figure 6, which include atomic vacancies, heptagon-pentagon topological defects, adatoms, and dopants. These defects can be inherited by graphene during growth [55] and can be introduced advertently by resonance plasma [56] , chemical treatment [57] , and irradiation [58] . Surface defects are also introduced during functionalization, oxidation, and reduction [59] . For example, in reduced graphene oxide, the graphene layers are found to comprise of a few nanometers interspersed defect areas dominated by clustered heptagons and pentagons [59] . The first-principles calculations show that atomic-scale defects result in both intravelley and intervalley scattering of graphene and Fermi velocity is decreased in the vicinity of the defects due to enhanced scattering [60] . The electronic properties of graphene are strongly influenced by structural defects that deform the original honeycomb lattice [61] [62] .

4. Topographical Features

The exploitation of topographically modified geometries in synthetic and bio-inspired materials is a novel area of research [63] - [66] . The topographically modified carbonaceous materials are produced by various methods and have found numerous applications [67] - [70] . It was shown that the superior electronic properties of graphene are sensitive to topography [71] . Through their crack deflection modeling, Faber and Evans showed that maximum improvement in fracture toughness, among all other nano-reinforcements, can be obtained using graphene mainly because of its better capability of deflecting the propagating cracks [72] [73] . The schematics of various topographies of graphene are shown in Figure 7 [74] . The graphene sheets have coiled structure that helps them to store sufficient amount of energy [75] [76] . The individual sheet and chunk of sheets together are subjected to plastic deformation at the application of external load. The applied energy is utilized in undertaking plastic work that enhances the ability of graphene to absorb more energy [77] . Graphene has shown inclination

for stable folding and bending energy at folds is compensated by intersheet adhesion (Van der Waals interactions) [78] [79] . The individual layers of graphene, under external loadings and thermal stresses, undergo out-of- plane wrapping [80] [81] , rippling [82] [83] , folding [84] [85] , scrolling [86] [87] , and crumpling [88] [89] , making graphene suitable to enhance the toughness of polymers. Wang et al. showed that wrinkle’s wavelength and amplitude are directly proportional to volumetric dimensions of graphene as evident from Equation (27) and eq.(28), where λ is wrinkle wavelength, ν is Poisson’s ratio, L is the length of graphene sheet, t is thickness of graphene sheet, ε is edge contraction on a suspended graphene sheet, and A is wrinkle amplitude [90] .

(27)

(28)

The topographical features of graphene can be observed experimentally and explored using computer simulations [91] - [96] . Abedpour et al. used both analytical and MD-based numerical simulations and correlated the surface topography of graphene with thermal fluctuations [97] . The Monte Carlo-based simulations predicted the formation of spontaneous ripples on graphene surface at finite temperatures. The ripple amplitude is comparable with the interatomic distance of graphene C-C covalent bonds (~0.142 nm) and was found independent of the total length of graphene sheet. It was also observed that the spontaneous ripples had a characteristic wavelength of about 8 nm [98] . Therefore, rippling being the intrinsic feature of graphene sheet is expected to strongly influence the mechanical properties of graphene and graphene based nanocomposites. Liu and Zhang carried out MD simulations to study the influence of thermal fluctuations on bending rigidity of graphene [99] . Montazeri and Rafii-Tabar observed an exponential drop in bending rigidity and lower axial Young’s modulus with increasing temperature in graphene-PMMA nanocomposites [100] . It can be attributed to bond breakage. Graphene has sp² hybridization which is stronger than sp³ hybridization found in diamond. Therefore, graphene shows superior mechanical properties. However, at the application of elevated temperatures, C-C bonds deteriorate and protrude out of the graphene plane. These misaligned bonds do not offer their intrinsic mechanical strength in the axial direction. The graphene sheets are held together by weak Van der Waals forces which can be subdued quite easily making sheet sliding inevitable in Multi-Layered Graphene (MLG). It was also observed that Young’s modulus significantly decreased due to layer sliding at different temperatures [100] .

To model dynamics and energetics of graphene C-C covalent bonds, second-generation Brenner many-body interatomic potential can be used [101] [102] . The Brenner potential is given by Equation (29) [103] ,

(29)

where and are the repulsive and attractive pair-wise interactions, respectively, is the

distance between pairs of nearest-neighbor atoms i and j, and denotes the reactive empirical bond order depending on the local bonding environment [100] . The simulation temperature during MD simulation can be maintained through Nosé-Hoover thermostat [104] . Montazeri and Rafii-Tabar used MD simulations to produce ripples on the surface of graphene to study the influence of ripples on graphene surface due to thermal fluctuations on axial Young’s modulus of graphene-PMMA nanocomposites [100] . The macroscopic behavior of the model was studied using finite element method as shown in Figure 8 [100] . The model was subjected to tensile load to produce 5% strain. The Young’s modulus was calculated using Equation (30) [100] , where F is the total force acting at the end of Representative Volume Element (RVE) (nN), H_o is the initial length of RVE (nm), and ΔH is the elongation. The obtained value of axial Young’s modulus of graphene-PMMA nanocomposite was 59.536 GPa [100] . The experimental results are much lower than the experimentally measured values. It can be attributed to the non-random distribution of graphene in polymer matrix considered in the model. The experimentally produced nanocomposites contain random distribution of graphene

(30)

along with agglomerated graphene. In addition, the proposed model considers perfect structure of graphene. On the contrary, the actual graphene contains surface and structural defects which may deteriorate the mechanical properties of graphene and produced nanocomposites. An additional influential factor could be layer sliding which was not considered in the modeling approach but takes place in reality due to weak Van der Waals interactions between the graphene layers.

5. Interfacial Interactions

The schematic diagram illustrating the interfacial interactions between graphene and polymer matrix is shown in Figure 9. When the nanocomposite is subjected to tensile loading, one possibility is that the graphene C-C bonds get broken. However, the strength of graphene C-C bond as measured by Lee et al. is approximately 130 GPa [105] . On the contrary, the strength of polymers is usually far below 100 MPa and that of nanocomposites also hardly exceeds 100 MPa. Therefore, this phenomenon of C-C bond breakage stretches the limits of plausibility too far. In addition, in non-functionalized carbonaceous reinforcements, usually the reinforcement-matrix bonding consists of weak Van der Waals interactions, physical bonding, and mechanical interlocking. Hu et al. proposed that the Van der Waals forces are responsible for the interfacial interactions in CNTs-polystyrene nanocomposites [106] . These forces are much weaker than C-C covalent bonds. Therefore, it is highly likely that the graphene sheet is pulled-out. The energy is required to drag the graphene out of the polymer matrix to circumvent Van der Waals interactions and physical bonding. This energy dissipation will eventually increase the mechanical properties of nanocomposites. However, this increase would not be commensurate with the abilities of graphene as reinforcement. Therefore, extensive research has been carried out to modify the surface of graphene by functionalization to improve interfacial interactions [107] - [111] .

According to ASM Handbook on composites (vol. 21 [1] ), “when interfacial bond strength is very large, the failure mode changes from interfacial to matrix, and the composite behaves like a brittle material, that is, it becomes ‘notch-sensitive’. Thus, excessive fiber-matrix bond strength may have a detrimental effect on the longitudinal tensile strength of the composite.” However, to the best of our knowledge, no article is yet published on GBPNCs and previously on CNTs reinforced polymer nanocomposites where the interfacial interactions were reported strong enough to have detrimental effect on mechanical properties. Instead, research is still ongoing to further improve the interfacial interactions. Das et al. have shown that mechanical properties of Polyvinyl Alcohol (PVA) significantly increased when reinforced with functionalized few layer graphene [112] .

The improvement in mechanical properties of GBPNCs significantly depends on interfacial interactions as load must efficiently be transferred from the matrix to the reinforcement to alleviate the stress concentration in relatively weaker matrix phase. One of the factors influencing interfacial interactions is interfacial area. The interfacial interactions can significantly be improved with increasing interfacial area. Interfacial area can be increased by reducing reinforcement size and by altering the topography as shown in Figure 10. By texturing the topography, the surface area can significantly be increased even when the size of the sheet is unaltered. With the

increase in surface area, interfacial area with polymer matrix also increases which can significantly improve the mechanical properties of nanocomposites. In addition to that, mechanical interlocking between polymer and reinforcement will take place which can cause a further increase in the mechanical properties of nanocomposites as shown in Figure 11. When reinforcement has smooth surface, there will be poor physical bonding between reinforcement and matrix. However, when reinforcement has corrugated surface, polymer chains will be locked within the spongy structure of reinforcement. This mechanical interlocking can significantly increase the interfacial interactions and concomitant rise in mechanical properties. Surfaces can be made porous or rough to enhance the extent of mechanical interlocking [113] . Karger-Kocsis et al. have studied that hierarchical and hairy fillers have high surface area and capillary wetting by the polymers [114] . The textured fillers also exhibit mechanical interlocking with the polymers and cause local reinforcement of the fiber-matrix interphase [114] . Moon and Jang studied the mechanical interlocking and wetting at the interface between argon plasma treated Ultra-High Modulus Polyethylene (UHMPE) fiber reinforced vinylester resin composite [115] . They observed a significant rise in interlaminar shear strength. It has been shown that plasma etching of UHMPE produces micro-pittings on fiber surface and this spongy surface structure helps improving the mechanical interlocking with the polymer matrix and causes a significant increase in interlaminar shear strength [116] - [119] . Therefore, mechanical interlocking can play a significant role in improving the mechanical properties of GBPNCs.

The graphene-polymer interactions can be defined via LJ-potential wherein spring elements with a nonlinear behavior can be used to couple the graphene and polymer matrix [120] [121] . The graphene sheets bear a major fraction of applied load and as a result, graphene sheets undergo buckling. Parashar and Mertiny developed a multiscale model using finite elements to investigate the phenomenon of graphene buckling in GBPNCs [122] . Lindahl et al. studied the snap-through behavior of convex buckled graphene membranes under the applied electrostatic pressure and reported that bilayer graphene shows a bending rigidity of ~35.5 eV [123] . Depending on the type of interfacial interactions, two main types of buckling takes place: (a) cooperative buckling, and (b) discrete buckling. In cooperative buckling, buckling occurs in the member simultaneously, equally, and in the same direction. Cooperative buckling occurs when interfacial interactions are strong enough to withstand the applied force. Consider a graphene based polymer material as shown in Figure 12(a) subjected to compression loading [124] . When load exceeds the compressive strength, member buckles. When interfacial interactions are strong enough, cooperative buckling takes place as shown in Figure 12(b). The total potential energy, Π₁, of the elastic member subjected to cooperative buckling is given by Equation (31), where b is sample width, D₀ is bending rigidity, is the curvature under buckling, and P is applied axial force. This model allows a first-order (e.g. linear) analysis of graphene based nanocomposites. The first term in Equation (31) belongs to elastic bending energy while second term is measure of work done [124] .

In discrete buckling, also known as delaminated buckling, all layers buckle in disparate way. As indicated by

name, discrete buckling occurs when interfacial interactions are not strong enough to withstand the applied force. When compressive load exceeds the compressive strength and interfacial interactions are poor, delamination and discrete buckling take place as shown in Figure 12(c). The total potential energy of elastic member undergoing delaminated buckling, Π₂, is mentioned in Equation (32), where D₁ and D₂ are bending rigidities of materials 1 and 2, and γ(u) is adhesion energy [124] . The first two terms in Equation (32) belong to elastic bending energy, third term is measure of work done, and fourth term is energetic cost of delamination. This discrete buckling may or may not lead to large-scale delamination or wrinkling of composite sandwich structure. The bending mode depends on the change in effective bending stress which is given by Equation (33). If work done by ap-

plied force is independent of buckling case, then change in potential energy () is given by Equa-

tion (34). Delamination occurs when ΔΠ > 0, and does not occur when ΔΠ < 0 [124] .

To determine whether cooperative buckling occurs or discrete buckling, Ritz method can be used [125] [126] . The approximate displacement function, u(x) is given by, , , and,

where C is an undetermined coefficient and satisfies the boundary conditions for both the cases shown in Figure 12 [124] . It defines the zero rotation at the extremities and is different from Euler buckling case in which zero rotation is enforced at the extremities for classical presumed displacement function. The Ritz method provides a conservative approximation as it is unlikely that exact displacement function, u(x), is constant for both cases. Therefore, the Rtiz method imposes the displacement equivalence. A harmonic approximation is supposed for the change in adhesion energy as given by Equation (35), where k is the stiffness of adhesion (eV/Å⁴). From these values, ΔΠ is given by Equation (36) [124] ,

(31)

(32)

(33)

(34)

(35)

(36)

minimizing with respect to C and setting values as in Equation (37) and Equation (38) [124] ,

(37)

(38)

above equations give the critical condition for ΔD as given by Equation (39) [124] .

(39)

If ΔD < kL⁴/120, the change in stiffness from cooperative to discrete is not sufficient to overcome the adhesion. In this case, cooperative buckling is preferred mode when ΔΠ < 0. In case when ΔD > kL⁴/120, the change in stiffness from cooperative to discrete is sufficient to overcome adhesion. In this case, discrete buckling occurs when ΔΠ > 0. The discrete buckling can be avoided by improving the interfacial interactions [124] .

6. Dispersion State

The graphene is unfortunately produced in entangled form and even if it is somehow dispersed, it tends to re- aggregate due to attractive Van der Waals forces [127] . The graphene agglomeration depends on size, shape, synthesis methods and surface chemistry. The utilization of graphene to enhance the properties of the system to which it is going to reinforce is fairly dependent on its dispersion state [128] . The schematic diagrams of agglomerated and uniformly dispersed graphene are shown in Figure 13. Because of extremely large surface area (~2500 m²/g) and the fact that commercially available graphene is heavily entangled restricts uniform dispersion of graphene in polymer matrix. It is a great barricade in its research and commercial applications. Instead of using agglomerated graphene, it is beneficial for the composite properties to use non-agglomerated graphene, as inter-layer sliding and the fact that aggregates act as stress concentrators can deteriorate the composite strength. The issue of graphene dispersion in polymer matrix has yet not been entirely outwitted, particularly for graphene fraction above a few wt% [129] .

There are two main methods for graphene dispersion. In the first method, the graphene can be separated by applying some external force, such as mechanical mixing and sonication, and is then stabilized by enclosure in a polymer or surfactant complex [130] [131] . It prevents the re-aggregation of graphene. This method can be used to form metastable dispersion. In the second method, the graphene is added in a suitable solvent where the graphene peel off from the bundle layer by layer and is dissolved in the solvent forming a solution [132] .

E.g., super-acids can dissolve graphene in large quantities [133] . In principle, the dispersion of graphene in polymer matrix can be controlled during the production route at two different stages. Firstly, the dispersion state in the uncured state can be tailored with the help of chemical functionalization, mechanical mixing, or by using surfactants. Secondly, dispersionstate significantly alters during curing stage [90] . There are chances of graphene agglomeration if curing requires high temperature. The extent of curing impact on the dispersion state of graphene depends on the curing temperature, time, chemical composition of matrix, and the way of heat input [134] .

The enhancement in the properties of GBPNCs by graphene is still not commensurate with the abilities inherited in the graphene structure [135] [136] . MD simulations showed that uniformly dispersed graphene is more effective in increasing the stiffness of polymers than poorly dispersed graphene [43] . Rahman studied the modeling of three different dispersion states of graphene as shown in Figure 14 [43] . The three dispersion states are, (a) SLG, (b) uniformly dispersed graphene, and (c) agglomerated graphene. Each unit cell contains about 3 wt% graphene with P3m space group symmetry of graphene sheet [137] . The calculations were carried out using open source molecular dynamics code LAMMPS developed at Sandia National Laboratory [138] . To get stress- strain response using MD method, the unit cells were subjected to engineering strain of 1 × 10⁻⁴ under NVT ensemble at every 300 steps up to 1500 steps. The temperature was fixed at 0.1 K. The energy, temperature, and pressure trends obtained during MD equilibration are shown in Figure 15 [43] . It was observed that stiffness of uniformly dispersed graphene-epoxy was higher than agglomerated graphene-epoxy nanocomposites [43] .

7. Aspect Ratio

The schematic diagram of various aspect ratios of graphene is shown in Figure 16. Although surface and interfacial areas are independent of the Aspect Ratio (A.R) of graphene, studies show that the aspect ratio of gra-

phene is an influential factor for the performance of GBPNCs. Rahman carried out modeling of graphene-epoxy nanocomposites with graphene of two different aspect ratios: length to width≥: (a) 5, and, (b) 10 [43] . The graphene sheets were randomly distributed in the unit cell. The molecular models were constructed using XENOVIEW with PCFF and unit cells subjected to Periodic Boundary Condition (PBC). PBC was applied by “wrapped around” fashion. For any atom with spatial coordinate x_i: i = 1, 2, 3., if x_i > L_i, then,

(40)

and if x_i < 0, then,

(41)

here, L_i: i = 1, 2, 3, represents the lengths of unit cell along x, y, z directions [43] .

Rahman observed that uniformly dispersed graphene with high aspect ratio can significantly improve in-plane Young’s modulus [43] . It is because the in-plane Young’s modulus of graphene is as high as 1 TPa.

However, graphene is not effective in increasing the out-of-plane Young’s modulus. It is because the out-of- plane stiffness depends on the modulus of polymer and Van der Waals interactions between graphene-graphene and graphene-polymer. It was observed that experimental Young’s modulus was higher for smaller aspect ratio than the larger one [43] . The highest shear modulus (G) was observed in case of smaller aspect ratio and at 1 wt% graphene. The influence of graphene aspect ratio on G_1C is shown in Figure 17 [95] . When the A.R is 1.4, the increase in G_1C is 18%. And at A.R = 2.8, the increase in G_1C becomes 24%. It shows that the G_1C and fracture toughness increase with increasing aspect ratio of graphene [95] .

Cranford studied the mechanical properties of graphene with different longitudinal lengths and observed that maximum deviation was shown by smallest effective length of 50 Å [124] . It can be attributed to inner graphene sheets which undergo local crumpling prior to buckling. The outer layers undergo axial compression under the influence of applied load. These two factors may add up to raise the local adhesion energy by increasing the molecular density. Therefore, the local stiffness of adhesion is enhanced even when energy per atom is constant. In addition, local fluctuations in structure caused by thermal vibrations yield shear/torsional deformations that

increase the difference in predicted and measured values. Depending on effective stiffness and adhesion strength, either cooperative ordiscrete buckling occurs. A trade-off among critical length (L), stiffness of adhesion (k) and effective rigidity (D) is critical in initial designing of molecular composites. Buckling of graphene layers and interfacial delamination are most commonly observed failure modes in GBPNCs. To avoid such failures, minimum length scale and adhesion energy are necessary to avoid delamination under compressive forces.

8. Weight Fraction

The schematic diagram of varying weight fractions of graphene in polymer matrix is shown in Figure 18. Due to very high surface-to-volume ratio and attractive Van der Waals forces causing agglomeration of graphene, a uniform dispersion of graphene in polymer matrix becomes quite difficult beyond certain weight fraction. In most of the reported literature, the mechanical properties increased up to 1 wt% of graphene and decreased afterwards which can be attributed to agglomeration of graphene and generation of cracks around the graphene agglomerates [90] [139] . The MTM showed increasing values of Young’s modulus with increasing weight fraction of graphene [140] . However, in atomistic simulation, the results showed that Young’s modulus increases up to certain weight fraction of graphene and then starts decreasing depending on the parameters [43] . It is because the micromechanics model such as MTM takes a homogeneous unit cell for calculations. Therefore, it shows a continuously increasing trend of Young’s modulus with increasing weight fraction of graphene. However, in atomistic model, various influential factors contribute toward the outcome such as spatial distribution of atoms, change in molecular energy in the system, and atom density etc.

Rahman developed two different types of models: Fixed Graphene Model (FGrM), and Varying graphene model (VGrM), and studied three different weight fractions of graphene (1 wt%, 3 wt%, and 5 wt%) [43] . The range of unit cell density was 0.7 - 1.0 g/cm³. The weight fraction of graphene was varied by changing the number of epoxy molecules in the unit cell. Each unit cell contained one graphene layer. As graphene layer is fixed in each unit cell, this model was termed as “Fixed Graphene Model (FGrM)”. The alternative way of changing the weight fraction is by changing the number of graphene sheets in the unit cell while number of epoxy molecules is kept constant. In this case, there is no significant change in the total number of atoms in the unit cell. As graphene layer is varied in each unit cell, this model was termed as, “Varying Graphene Model (VGrM)”. The conventional molecular dynamics scheme can be used to equilibrate the model to get the final configuration. The schematic of the unit cell is shown in Figure 19 [43] .

Rahman produced epoxy-graphene nanocomposites to determine the validity of developed models [43] . The graphene was dispersed using a combination of tip sonication and magnetic stirring. It was observed that Young’s moduli as calculated by FGrMs were slightly lower than from VGrMs. This variation can be attributed to change in number of electrons in the unit cell. The total number of atoms reduces as the weight fraction of graphene increases. Also, the total number of atoms in FGrMs is lower than in VGrMs. Moreover, there was large variation in atom density of FGrMs. Therefore, large variation is expected in the elastic constants obtained from FGrMs. However, the average stress is calculated based on all the atoms in the unit cell. It reduces the influence of cell size on overall stress-strain response resulting in negligible variation in Young’s modulus. In both modeling and experimentation, Rahman observed that mechanical properties increase up to a certain weight fraction and decrease afterwards [43] . Rafiee also showed that Young’s modulus increases up to a certain weight fraction of graphene and decreases afterwards [44] . The results obtained indicate that atomistic modeling approach is more accurate than micromechanics approach. The decrease in elastic constants with increasing graphene content can be attributed to the inevitable possibility of agglomeration and void content because of inter- graphene layer gap and graphene-polymer gap [141] .

9. Constitutive Response

The values of critical stress, strain, and Young’s modulus of SLG supported on perforated Si substrate were recorded using Atomic Force Microscopy (AFM) nanoindentation test and values obtained were 130 GPa, 25%, and 1 TPa, respectively [105] . In SLG, there are few slip systems available as no alternative plane is available for glide and Debelak and Lafdi have shown that shear and glide takes place along [001] direction in MLG [142] . The stress tensor σ_ij can be calculated using virial expression as given in Equation (42) [43] ,

(42)

where Vol is volume of the simulation box, M_I is mass of I^th atom, is velocity of Ith atom in I^th direction, N is total number of atoms, ij-component of distance between I^th and J^th atoms, and is ij-component of

force between I^th and J^th atoms. When uniaxial deformation is applied using uniaxial tension or compression, the total potential energy changes which can be used to determine elastic constants. The internal stress tensor can be calculated by taking the first derivative of potential energy with respect to strain. The second derivative of potential energy with respect to strain yields stiffness matrix and deformation in unit cell is given by Equation (43) where the corresponding stresses in the system (and) are calculated using Equation (42) [43] . The stiffness components can be written as given in Equation (44) where σ_i is i^th component of internal stress tensor, and and are stress tensors under tension and compression, respectively, and elastic properties of the system can be calculated using Equations (45)-(47) [43] . By using the MD simulation, the average stress in the system can be calculated using Equation (48) where σ₁₁, σ₂₂, and σ₃₃ are the virial stresses in x, y, and z directions, respectively, and is the average stress of the unit cell [43] . The stress-strain response can be evaluated by applying strain to the system repeatedly followed by MD equilibration. The Young’s modulus can be calculated from the slope of initial linear portion of stress-strain curve.

Ebrahimi et al. used the conjugate pair of the Green-Lagrange strain and the second Piola-Kirchhoff stress tensor to study the stress-strain response of graphene-chitosan nanocomposites under large deformations [45] . The strain energy density and inter-atomic potential, u, of the atoms within the NRVE is given by Equation (49)

where is the sum of potential energies of the atoms in NRVE, and

VRVE is the volume of NRVE [45] . The second Piola-Kirchhoff stress tensor, T, is given by Equation (50) [50] , and can be converted to Cauchy stress tensor using Equation (60) [143] [144] , where J is the Jacobian of the deformation. The constitutive response of nanocomposites is dependent on the ratio of element volume (VEL) to NRVE volume (VNRVE) as in Equation (61) [45] . The variation in Cauchy stress in NRVE of chitosan with uniaxial strain at different values of R is shown in Figure 20. The values of yield stress and tensile modulus are close to experimental values. The variation in tensile modulus with R is shown in Figure 21 which shows that tensile modulus becomes less sensitive to R as value of R increases and finally becomes independent of R [45] . The variation in Cauchy stress with uniaxial strain in chitosan when reinforced with graphene and (10, 10) armchair carbon nanotubes is shown in Figure 22. The carbon nanotubes were supposed to be immovable during simulation and graphene sheets were aligned along the principle stress axis. The graphene showed a higher increase in strength and modulus than carbon nanotubes.

(43)

(44)

(45)

(46)

(47)

(48)

(49)

(50)

(51)

(52)

(53)

(54)

(55)

(56)

(57)

(58)

(59)

(60)

(61)

10. Toughness Enhancement

The ability of a material containing crack to resist fracture, known as fracture toughness, is a simple yet trustworthy indicator of the material’s damage tolerance and hindrance against fracture, and is considered as one of the most important mechanical properties of the engineering materials. The fracture toughness depends on how the pre-existing and newly formed cracks propagate within a material. The plane strain fracture toughness (critical stress intensity factor, K_1C) can be calculated using Equation (62), where P_max is maximum load of load-dis- placement curve (N), f(a/w) is constant related to geometry of the sample and is calculated using Equation (63), B is sample thickness (mm), W is sample width (mm), and a is crack length (should be kept between 0.45 W and 0.55 W according to ASTM D5045) [128] . The critical strain energy release rate (G_1C) can be calculated using Equation (64), where E is the Young’s modulus obtained from the tensile tests (MPa), and ν is the Poisson’s ratio of the polymer. The geometric function f(a/W) strongly depends on the a/W ratio as shown in Figure 23 [145] .

Parashar and Mertiny developed three dimensional RVE based multiscale model in finite element environment to study the influence of graphene on the fracture characteristics of graphene-polymer nanocomposites [95] . The polymer matrix was modeled as a continuum while graphene was modeled in atomistic state. The truss elements were used to simulate Van der Waals interactions between graphene and polymer. Virtual crack closure technique was employed to study the fracture characteristics of the nanocomposites. The bonded interactions between C-C bonds were defined by modified Morse potential [146] [147] . The potential energy (E) of isolated graphene sheet is given by Equations (65)-(67) [95] , where E_S is bond-stretching component of potential energy (J), E_B is angle-bending component of potential energy (J), D_e is dissociation energy (Nm), β is constant controlling the width of the potential (m⁻¹), r is C-C bond length (m), r₀ is C-C equilibrium bond length

(m), K_Ө is force constant for bond bending (N∙m/rad²), Ө is current angle of the adjacent bond (rad), Ө₀ is initial angle of the adjacent bond (rad), and K_sextic is constant in bending term of potential (rad⁻⁴). The continuum polymer matrix, atomistic graphene structure, and interface were modeled using SOLID45, BEAM4, and LINK8, respectively, in ANSYS finite element software [95] . The opening mode crack propagation was simulated by the subsequent release of node constraints under plane stress condition. By differentiating Equation (66), the interatomic force acting between the two C-C atoms can be calculated as Equation (68) and material stiffness values can be assigned to each truss element using Equation (69) [95] .

(62)

(63)

(64)

(65)

(66)

(67)

(68)

(69)

The stress-strain behavior of C-C bonds, configuration of multiscale model, and variation in G_1C with graphene volume fraction are shown in Figures 24-26, respectively [95] . The improvement infracture toughness

can be attributed to variation in stress distribution by graphene in the polymer matrix. Because of high stiffness and space frame structure, graphene absorbs most of the applied load. The graphene also shields the crack tip from crack-driving forces and opening loads. The modeling results showed that when graphene is present on both sides of the crack plane, the crack tip shielding effect is more pronounced.

11. Thermal Properties

Because of superior thermal conductivity of graphene, GBPNCs are promising candidates for high-performance thermal interface materials. The dissipation of heat from electronic devices may also be barricaded when thehigh thermal conductivity of graphene is efficiently utilized. The graphene has shown higher efficiency in increasing the thermal conductivity of polymers than CNTs [148] . It has been found experimentally that the effective thermal conductivity (K_eff) of GBPNCs has a non-linear dependence on graphene weight fraction [149] - [151] . Xie et al. proposed an analytical model to determine effective thermal conductivity of graphene based nanocomposites [152] . Their model proposed very high thermal conductivity values as the model did not take into account the interfacial thermal resistance. Lin et al. developed a model based on Maxwell-Garnett effective medium approximation theory to determine effective thermal conductivity of graphene based nanocomposites [153] [154] . They showed that the enhancement in thermal conductivity is strongly influenced by aspect ratio and orientation of graphene. Hu et al. used molecular dynamics approach to show that the agglomeration of graphene is of major concern in increasing the thermal conductivity of the system [155] .

The theoretical calculations predict thermal conductivity values of polymer nanocomposites far greater than the experimental values [156] . It can be explained on the basis of Thermal Boundary Resistance (TBR, also known as Kapitza resistance) between graphene and polymer chains [157] . The heat transfer in composite system is restricted by interfacial thermal resistance arising from weak phonon-phonon coupling which leads to backscattering of phonons at the boundary region. The interfacial thermal resistance increases with increasing interfacial area. Therefore, this resistance becomes quite significant in case of graphene because of its high surface area. The heat transfer takes place as a consequence of Brownian motion of discrete thermal walkers [158] . In Brownian motion, the thermal walkers exhibit random changes in motion in each time step and can be modeled using random algorithm [159] - [162] . In each space direction, these position changes take values from a normal distribution with a zero mean and a standard deviation depending on the thermal diffusivity (D_m) and the time increment (Δt) as given in Equation (70) [156] .

(70)

The TBR (R_bd) is introduced at graphene-polymer interface by a phonon transmission probability (f_m-GA) which can be measured according to acoustic mismatch theory as given in Equation (71) [163] ,

(71)

where, , and are the density of, specific heat of, and sound velocity in polymer matrix, respec-

tively. The improvement in thermal conductivity can be theoretically described by the percolation theory using a power law relationship as given in Equation (72) [164] - [169] ,

(72)

where and are thermal conductivity and volume fraction of graphene, respectively, is percolation

threshold, and t is universal critical exponent. A comparative infrared microscopy technique can be used to measure thermal properties of nanocomposites [156] . The thermal properties can be approximated by Fourier’s law [170] [171] . The variation in thermal conductivity of graphene-PMMA nanocomposites with varying graphene content is shown in Figure 27 [156] . The thermal conductivity increases up to 3.5 times of pure PMMA at 2.5 vol% of graphene. The increase in thermal conductivity may be attributed to large contact area and strong interface between graphene and matrix. As there are strong π-π interactions between graphene-graphene, the TBR between graphene-graphene can be neglected. However, the TBR between graphene-PMMA must be considered and taken as 1.906 × 10⁻⁸ m²K/W [172] . When graphene is oriented along the heat flux, higher K_eff values are obtained due to efficient heat transfer channel provided by graphene. However, when graphene is oriented

perpendicular to the direction of heat flux, lower K_eff values are obtained due to TBR. The modified effective medium theory to estimate K_eff is given in Equation (73) [156] ,

(73)

where is TBR between graphene and PMMA, H is thickness of Functionalized Reduced Graphene Oxide (m-rGO) sheets, is thermal conductivity of graphene, is thermal conductivity of PMMA, and f is

volume fraction of graphene. The results show that the values predicted by EMT are greater than experimental values. One main reason for this is the particle size distribution which EMT does not take into account. In general, a large particle size distribution can increase the filling ratio as the finer particles can fill in the interstitial sites created by larger particles. The higher filling ratio helps in making strong networks of graphene which can significantly increase the thermal conductivity values compared with narrow size distribution.

Chu et al. proposed an analytical model based on Differential-Effective-Medium (DEM) theory to determine the thermal properties of graphene based nanocomposites [173] . The DEM theory can be used to effectively describe the thermal properties of composites [174] - [176] . Chu et al. considered the Graphene Nanoplatelets (GNPs) as oblate spheroid with very large aspect ratio and considered the isotropic two phase composite system. The composite system was considered with temperature field and heat flux which are given by Equation (74) and Equation (75) [173] , respectively, where d(T) is temperature distribution function.

(74)

(75)

The solution for the temperature fields, T_p and T_m, within the dispersed particles and surrounding matrix, respectively, is given by Equation (76) [173] ,

(76)

where H_j is the depolarization factor of the ellipsoidal particles along j-axis, k_m is the thermal conductivity of matrix, and k_j is thermal conductivity of GNP along j-axis [173] . The depolarization factor (H) obeys the following equation, 2H_x + H_z = 1, where,

(77)

Equation (76) can be used to approximate the dependency of effective thermal conductivity of nanocomposites on the reinforcement volume fraction as given in Equation (78) [174] [176] ,

(78)

when the initial condition is kept k^* = k_m at GNP volume fraction f = 0, Equation (78) can be written after integration as Equation (79) [173] ,

(79)

where K_z and K_x are the thermal conductivity values of GNP along longitudinal and transverse axes, respectively. At a very large aspect ratio (), the values of H_x and H_z can be taken as zero and unity, respectively. Therefore, Equation (79) can be simplified as Equation (80) [173] .

(80)

To consider interfacial thermal resistance, Chu et al. considered GNP as core-shell structure as shown in Figure 28 [173] . The interfacial boundary layer has thickness c and thermal conductivity values of k_s and the effective thermal conductivity can be written as Equations (81)-(83) [154] , where Θ is the dipolar factor of shelled GNPs [177] . When interfacial thermal resistance is incorporated, the effective thermal conductivity of graphene based nanocomposites is given by Equation (84) [173] . Liu et al. also proposed the model based on Effective Medium Approximation (EMA) theory to determine effective thermal conductivity of graphene based nanocomposites as given in Equations (85)-(87) [153] . When interfacial thermal resistance is incorporated, the model takes the form as Equation (88) with comparison of the results shown in Figure 29 [173] .

12. Electrical Properties

Tailoring the electrical properties of graphene can unlock many potential electronic applications of graphene [178] [179] . For example, effective gauge fields are introduced when lattice deformation of graphene takes place. Like the effective magnetic field, the produced effective gauge fields influence the Dirac fermions [180] . The Fermi level in undoped graphene lies at the Dirac point where the minimum conductivity values are achieved [181] . By adding free charge carriers, i.e. dopants, the electrical properties of graphene can be improved and conductivity increases linearly with carrier density [182] [183] . For example, boron as dopants can contribute ~0.5

(81)

(82)

(83)

(84)

(85)

(86)

(87)

(88)

carriers per dopant in graphene sheet [184] . Dopants can be introduced during the synthesis of graphene using Chemical Vapor Deposition (CVD) [185] .

Due to graphitic composition, graphene exhibits excellent electrical conductivity. The charge carrier mobility of graphene at room temperature is as high as 15,000 cm²/Vs [186] . When graphene is uniformly dispersed in the polymer matrix, graphene can provide efficient pathways for conduction which can boost electrical conductivity of nanocomposites. At a certain graphene loading, the electrical conductivity of GBPNCs shows a drastic increase and they lose insulating character and start behaving like conductors. This graphene loading is called as percolation threshold and is usually taken at thirteen orders of magnitude increase in the electrical conductivity of produced nanocomposite [187] - [189] . The percolation threshold should be as low as possible in order to decrease the use of graphene. The dispersion state or degree of agglomeration of graphene is known to have a significant influence on the percolation threshold and electrical conductivity of GBPNCs [190] . The key factors influencing the percolation threshold include, dispersion state, physicochemical interactions between the matrix and the reinforcement, shape of the reinforcement, and method to produce nanocomposite. The conductivity models can be classified into three main categories: (1) structure-oriented [191] , (2) thermodynamic [192] , and (3) statistical [187] [193] . The structure-oriented models consider the microstructure of nanocomposites before and after the processing. As properties are strongly dependent on microstructure, therefore the structure-oriented models give a more realistic approach to predict conductivity behavior of nanocomposites. In certain cases, when achieving the actual microstructure is difficult, the structure-oriented models replace the microstructure with fitting parameters. However, these fitting parameters are not a true replacement of actual microstructure especially at nano-scale.

The electrical conductivity of nanocomposites strongly depends on the microstructure. When polymer is filled with a conducting nano-filler, conduction of the produced nanocomposite depends on how the nano-fillers are aligned [194] . The nano-fillers should make a continuous network to provide a swift conductive pathway for the electrons [195] . Various models have been proposed to approximate percolation threshold for nanocomposites with various nano-fillers such as graphite, CNTs, and graphene. Many proposed models use fitting parameters where achieving actual values is difficult. As microstructure is a key factor in defining the properties of nanocomposites, it should be efficiently incorporated in proposed models to get closer to real-time values. Syurik et al. proposed a modification in McCullough’s model with current maps obtained by Conductive Atomic Force Microscopy (CA-AFM) [192] .

McCullough’s model is a structure-oriented model which assumes that the microstructure and orientation of reinforcement inside the matrix are variables and depend on the reinforcement content [191] . The increase in electrical conductivity can be attributed to the formation of chain-like network of reinforcement. The chain length depends on the size, shape, and volume fraction of reinforcement. Below percolation threshold, this chain-like network is sporadic and all chains are not connected together. Due to cessation of chain continuity, the conduction is obstructed and nanocomposites behave as insulator. By the time percolation threshold is achieved, all chains cling together and construct a continuous network which acts as a highway for electrons and nanocomposites start behaving as conductors. The chain fraction and its dimensions are statistical in nature. The distribution of chain length can be described by a probability density function η(ν_f, a), where a is the effective

aspect ratio and ν_f is the volume fraction of reinforcement. The average value of chain parameter can be calculated using Equation (89) [192] ,

(89)

where λ(a) is chain parameter and is dependent on reinforcement shape. When aspect ratio a of reinforcement exceeds 1, the chain parameter can be calculated using Equation (90) and Equation (91) [192] . The isotropic distribution of reinforcement in the matrix is given by Equation (92) and Equation (93) [192] ,

(90)

where,

(91)

(92)

where,

(93)

where σ_m, σ_f, and σ are electrical conductivities of matrix, reinforcement, and composite, respectively, and ν_f and ν_m are volume fractions of reinforcement and matrix, respectively. The experimental determination of η(ν_f, a) and is quite difficult. Therefore, these factors were taken as fitting parameters. As actual microstructure is not considered, therefore it poses a limitation to McCullough’s model exhibiting a semi-empirical character. Syurik et al. used AFM to determine actual value of [192] . For ideal shapes of reinforcement such as spheres, disks and ropes, the chain parameter can be substituted as given in Equation (94) [192] ,

(94)

where <à²> is mean-square value of conductive chains aspect ratio (a). With the help of AFM in CA-AFM model, the value of <à²> can be calculated as in Equation (95) [192] ,

(95)

where n is the number of conductive chains, ℓ_i is length of i^th chain, d_i is diameter of i^th chain, and

are average length and average diameter of chain, respectively. Therefore, Equation (92) can be reduced to Equation (96), and density of composite can be calculated using rule of mixture as given in Equation (97) [192] ,

(96)

(97)

where ρ, ρ_f, ρ_m, are specific densities of composite, reinforcement, and matrix, respectively. During production of nanocomposites, porosity is inevitable due to, (1) air entrapment, (2) evaporation of volatiles, (3) any shrinkage during curing, and (4) the relative movement of reinforcement and polymer chains. The porosity influences the conductivity which can be considered using an apparent density (ρ_ap) connected with specific density through the pore coefficient as given in Equation (98) [192] .

(98)

When the nanocomposites are produced by a reproducible technology, such as latex technology, in which conductivity trends remain the same when samples are reproduced using the same parameters [196] . The K_ρ can be taken as constant and needs to be measured only once. From Equation (98), Equation (97) can be written as Equation (99) where the coefficient K_ρ depends on the method of composite production and on geometry of the reinforcement [192] . The value of ρ_ap can be defined experimentally if the reinforcement fraction (w), sample’s volume (V) and mass (m) are known as given in Equation (100), and the volume fraction of reinforcement can be calculated using Equation (101) where volume fraction can be converted to weight fraction using Equation (102), and G can be calculated by combining Equation (103) and Equation (104) [192] . The Equation (103) and Equation (104) show the conductivity behavior with respect to reinforcement loading. The CA-AFM image of graphene nanoplatelets reinforced polystyrene (GNPs/PS) nanocomposites is shown in Figure 30(a) and corresponding image of the grains is shown in Figure 30(b) where 368 grains are shown with size distribution of the GNP clusters shown in Figure 30(c) [192] . The GNPs/PS nanocomposites were prepared by latex technology which allows the reproducible

(99)

(100)

(101)

(102)

(103)

(104)

percolation threshold [196] . Most of the GNPs had a thickness of 2 - 3 atomic layers and average surface area 1 - 3 µm². The weight fractions of GNPs were 0, 0.6, 0.9, 1.5, and 2.0 wt%. Electrical conductivity measurements in Direct Current (DC) mode were performed in a direction parallel to the sample top surface using a 2-probe configuration and a Keithley 2602 system source meter. Conductivity was calculated from the obtained I/V characteristics using Equation (105) [192] ,

(105)

where V is applied voltage and I is current value through a cross-section (A) and between the distance (b). The values of for all composites are in the range 0 < λ < 1 and lie in the restrictions of McCullough’s model

[191] . Figure 31 shows the variation in conductivity values measured using: (a) CA-AFM (curve 1), (b) 2-point DC method (curve 2), (c) McCullough’s model (curve 3), and (d) with varying GNP loading (curve 4) [192] .

The relationship between and is shown in Figure 32 [192] . The conductivities by CA-AFM and

2-point DC measurements are in accord. The conductivity of polymer at low GNP loading is close to dielectric polymer matrix as no continuous pathways of GNPs are available for conduction. The percolation threshold of 0.9 wt% is suggested by DC measurements and is corroborated by Scanning Electron Microscopy (SEM) and CA-AFM. At the percolation threshold marked by DC measurements, the conductivity increases sharply by five orders of magnitude. This can be explained by the presence of microscopic conductive sub-networks of GNPs. However, these sub-networks do not create continuous network. Therefore, these sub-networks are not initially detected by DC measurements until a continuous network is established.

13. Overall Microstructure

With the help of optical and electron microscopy, the dispersed reinforcement and polymer matrix can be observed in qualitative and quantitative manners, respectively. However, in this binary system integrated at nanoscale, there are myriad of factors that expound the overall performance of produced nanocomposites. Some of those defining features are summed up in Figure 33. One constituent of the overall microstructure is particle size distribution. The mechanical properties are improved when load is efficiently transferred from matrix to reinforcement. One of the controlling factors for load transfer mechanism is networking of reinforcement. The load can only be transferred from the matrix to the reinforcement if a connected network of reinforcing particles is available. Any disjoint in the network will act as weakest link for load transfer and polymer matrix will be prone to external loading. One of the factors influencing network formation is filling ratio (or packing density). When particle size distribution is narrow, the voids between the particles would not be filled and those empty locations will be a preferred route for the cracks to surmount the reinforcement particles. On the other hand, when reinforcement has wide size distribution, the finer particles can occupy the empty spaces in between large particles. It increases the filling ratio and makes an efficiently connected network of reinforcement. Sohn and Moreland have shown that packing density is dependent on particle size distribution and shows a direct relationship, i.e. packing density increases as the particle size distribution is extended [197] . It was also found that packing density is independent of particle size. They further reported that particle shape also influences packing density. It is obvious as perfect cubes will have 100% packing while voids will be certain in spherical particles which will lower the packing density. Therefore, a wide size distribution is helpful in improving the mechanical properties as strong networking of reinforcing particle can take place because of high packing density.

The properties of nanocomposites are also significantly dependent on filler shape and size. The graphene size, shape, and topography can be controlled simultaneously as shown in Figure 34 [198] . Chen et al. synthesized well-aligned millimeter-sized tetragon-shaped and hexagon-shaped graphene on a polycrystalline copper substrate using low pressure CVD [198] . CVD is an efficient approach to produce high quality and large surface area graphene on metallic substrates [199] [200] . Graphene can be grown epitaxially on single crystal copper substrate [201] - [204] . Graphene shape and orientation strongly depend on crystallographic orientation of copper substrate [205] [206] . The dendritic graphene with multiple branches can be obtained by diffusion-limited growth dynamics [207] . The graphene shape can also be controlled using surface condition of substrate. For example, wet-loaded samples produce tetragonal shaped and dry-loaded samples produce hexagonal shaped graphene [198] . The graphene shape can also be controlled by processing temperature [208] [209] .

It is well established that dispersion state is a crucial factor in defining the properties of nanocomposites and is discussed in detail in Section 6. In addition, the orientation of reinforcement and polymer chains and degree of crosslinking significantly influence the mechanical, thermal, and electrical properties of GBPNCs [210] . Additionally, the porosity is another important influential factor and two kinds of porosity is commonly observed: round porosity which can be because of fluids and air entrapment as fluids exert equal pressure in all sides, and irregular porosity which can be because of reinforcement bridging and relative movement of reinforcement and matrix [211] . One of the factors that dictates the interfacial interactions is the interphase and the properties of GBPNCs are significantly influenced by interphase properties [212] . Another factor that is not an integral part of microstructure but can significantly influence the microstructure is thermal fluctuation [213] . Thermal fluctuations can cause phase transformation and influence topographical features of graphene. Similarly correlations among other variables need to be established. If we consider the variables shown in Figure 33, there may be 156 permutations of the correlations considering two variables at a time (₁₃P²). Although various models have been presented to date regarding GBPNCs in which various aspects have been covered, however, it is just the beginning and myriad of factors are still to be incorporated in the modeling and simulation of GBPNCs.

14. Conclusions

Following conclusions can be drawn based on the critical review of the published literature related to modeling and simulation of GBPNCs:

Ÿ Various theoretical and computational approaches have been employed to explore the influence of graphene as reinforcement on the performance of polymer nanocomposites including but not limited to, quantum mechanical-based methods [28] , continuum mechanics (CM) [29] , molecular mechanics (MM) [30] , molecular dynamics (MD) [31] , atomistic modeling [32] , density functional theory (DFT) [33] , and multiscale modeling [34] .

Ÿ The individual layers of graphene, under external loadings and thermal stresses, undergo out-of-plane wrapping [80] [81] , rippling [82] [83] , folding [84] [85] , scrolling [86] [87] , and crumpling [74] [88] , making graphene suitable to enhance the toughness of polymers.

Ÿ Buckling of graphene layers and interfacial delamination are most commonly observed failure modes in GBPNCs [124] . To avoid such failures, minimum length scale and adhesion energy are necessary to avoid delamination under compressive forces [124] .

Ÿ MD simulation showed that thermal fluctuations produce ripples on graphene surface which influence the mechanical properties of graphene [99] [100] . Spring-based finite element models showed that graphene sheets can slide over each other which can degrade their reinforcing character [214] .

Ÿ With the help of optical and electron microscopy, the dispersed reinforcement and polymer matrix can be observed in qualitative and quantitative manners, respectively. However, in this binary system integrated at nanoscale, there are a myriad of factors that expound the overall performance of produced nanocomposites.

Ÿ If we consider the variables shown in Figure 33, there may be 156 permutations of the correlations considering two variables at a time (₁₃P²). Although various models have been presented to date regarding GBPNCs in which various aspects have been covered, however, it is just the beginning and myriad of factors which are still to be incorporated in the modeling and simulation of GBPNCs.

Acknowledgements

The authors would like to thank the Department of Mechanical and Construction Engineering, Northumbria University, UK for the provision of research facilities, without which the analysis of relevant data was not possible.

Cite this paper

Rasheed Atif,Fawad Inam, (2016) Modeling and Simulation of Graphene Based Polymer Nanocomposites: Advances in the Last Decade. Graphene,05,96-142. doi: 10.4236/graphene.2016.52011

References

List of Abbreviations

A.R: Aspect ratio, AFM: Atomic force microscopy, ASM: American society of metals, ASTM: American society for testing and materials, CA-AFM: Conductive atomic force microscopy, CM: Continuum mechanics, CNTs: Carbon nanotubes, CVD: Chemical vapor deposition, DC: Direct current, DEM: Differential effective medium, DFT: Density functional theory, EMA: Effective medium approximation, FGrM: Fixed graphene model, G1C: Critical strain energy release rate, GBPNCs: Graphene based polymer nanocomposites, GNPs: Graphene nanoplatelets, LAMMPS: Large-scale Atomic/Molecular Massively Parallel Simulator, K1C: Mode-I plane strain fracture toughness, LJ-potential: Lennard-Jones potential, MD: Molecular dynamics, MLG: Multi- layered graphene, MM: Molecular mechanics, m-rGO: Functionalized reduced graphene oxide, MTM: Mori- Tanaka Method, NRVE: Nanoscopic representative volume element, PBC: Periodic boundary condition, PCFF: Polymer consistent force field, PMCs: Polymer matrix nanocomposites, PMMA: Polymethyl methacrylate, PS: Polystyrene, PVA: Polyvinyl alcohol, RVE: Representative volume element, SEM: Scanning electron microscopy, SERR: Strain energy release rate, SLG: Single layer graphene, TBR: Thermal boundary resistance, UHMPE: Ultra-high modulus polyethylene, VGrM: Varying graphene model.

NOTES

^*Corresponding author.