Mr Mr Essa

Mr Mehdipour

Dr Wang

Mrs Umale

Mr Lhadi

Mr Kaboli

Mr Delgado-Rangel

Mr Wen

Mr Barth

Mrs Nierenberger

Bernard

Nanostructured materials are gaining an ongoing demand because of their exceptional chemical and physical properties. Due to complexities and costs of experimental studies at nanoscale, computer simulations are getting more attractive as experimental alternatives. In this PhD work, we tried to use combination of atomistic simulations and continuum modeling for the evaluation of thermal conductivity and elastic stiffness of nanostructured materials. We used molecular dynamics simulations to probe and investigate the thermal and mechanical response of materials at nanoscale. The finite element and micromechanics methods that are on the basis of continuum mechanics theories were used to evaluate the bulk properties of materials. The predicted properties are then compared with existing experimental results.

First, we compared Mori-Tanaka (mean field method), finite element (numerical method) and strong contrast (statistical continuum approach) for the evaluation of thermal conductivity and elastic modulus of composite materials. Then, we developed 3D finite element modeling for the investigation of interphase effects on the effective properties of nanocomposite materials. We also proposed a global statistical approach for the investigation of nanocomposites effective properties filled with spherical particles. The accuracy of the proposed global statistical approach was verified by comparing the modeling results with experimental tests on silica/epoxy nanocomposites. We also developed hierarchical multiscale procedures by the use of molecular dynamics and finite element method for the evaluation of effective elastic modulus and thermal conductivity of expanded graphite/PLA nanocomposites.

Graphene, the planar structure of carbon atoms is the most interesting nanostructured material because of its unique combination of exceptionally high mechanical, thermal and electrical properties. In this thesis, we used molecular dynamics simulations in order to investigate the thermal conductivity and tensile response of graphene. The validity of the developed molecular dynamics models were confirmed by observing close agreements between the simulations results and experimental values for the graphene and graphite. Then, we studied the effect of boron and nitrogen atoms substitution and different types of defects (point vacancy, bivacancy and Stone-Wales defects) on the thermal conductivity and tensile response of graphene. Moreover, we particularly studied the effect of curvature on the thermal conductivity of graphene.

We also proposed a simplified molecular dynamics method in order to construct crosslinked epoxy networks using the OPLS united atoms force field and dynamic crosslinking scheme as well. Finally, we developed a multiscale scheme on the basis of molecular dynamics and finite element methods for the modeling of effective thermal conductivity of graphene epoxy nanocomposites. In this method, we used molecular dynamics simulations for the investigation of thermal conduction in graphene epoxy assembly at atomic scale. Based on the information provided by the molecular dynamics simulation, we developed finite element models for the evaluation of effective thermal conductivity of nanocomposite. The validity of the proposed multiscale approach was confirmed by comparing the predicted results with experimental results in the literature.

Résumé

Les matériaux nanostructurés suscitent un intérêt qui va croissant en raison de leurs propriétés chimiques et physiques exceptionnelles. A cause de la complexité et du coût des développements expérimentaux à l'échelle nano, la simulation numérique devient une alternative de plus en plus populaire aux études expérimentales. Dans ce travail de thèse, nous avons essayé de combiner des simulations à l'échelle atomique avec de la modélisation en milieu continu pour évaluer la conductivité thermique et la réponse élastique de matériaux nanostructurés. Nous avons utilisé des simulations de dynamique moléculaire pour calculer la réponse mécanique et thermique des matériaux sur des volumes à l'échelle nano. Des méthodes de micromécanique et la méthode des éléments finis, qui utilise la mécanique des milieux continus, ont permis d'évaluer les propriétés mécaniques des matériaux. Les résultats obtenus par ces simulations numériques ont été ensuite comparés avec ceux issus de l'expérience.

Dans un premier temps, nous avons comparé la méthode de Mori-Tanaka (basée sur la théorie du champ moyen) avec celle des éléments finis (qui est une méthode numérique) et celle du « strong contrast » (qui est une approche statistique) pour l'évaluation de la conductivité thermique et du module élastique de matériaux composites à deux phases. Afin de pouvoir mieux comparer notre étude avec celles qui existent déjà dans la littérature, notre modélisation a été réalisée pour différentes géométries des charges de renfort, qui vont du cylindre de grande longueur en passant par la sphère et jusqu'au disque de faible épaisseur. Ensuite, nous avons comparé les résultats obtenus par la méthode statistique dite du« strong contrast » avec ceux obtenus par éléments finis et Mori-Tanaka pour différentes géométries et différentes concentration de charges. L'effet de la différence de propriétés entre charges et matrice a également été étudié. Les méthodes de modélisation développées ici peuvent être utilisées pour la prédiction des propriétés effectives de composites à matrice polymère et métallique. Nous avons développé une nouvelle procédure de modélisation pour l'évaluation de la conductivité thermique effective du composite par la méthode des éléments finis. Les résultats ont montré que la méthode des éléments finis est la méthode la plus précise et qu'elle peut donner davantage de détails pour la prédiction des propriétés effectives des composites. Par ailleurs, il est montré que la méthode de Mori-Tanaka est une alternative rapide et valide à celle des éléments finis pour une gamme de géométrie de charge limitée. Nos résultats indiquent que la méthode statistique du « strong contrast » basée sur des fonctions de corrélation à deux points ne peut pas décrire précisément l'effet de la géométrie des charges. Pour les nanocomposites dont les charges sont orientées aléatoirement, les résultats obtenus par éléments finis suggèrent qu'en incorporant des charges en forme de plaquette, on peut obtenir une meilleure amélioration de la conductivité thermique et des propriétés élastiques du composite que dans le cas de charges cylindriques. Ce constat suggère que les nanofeuillets de graphène sont supérieurs aux nanotubes de carbone pour ce qui est de l'amélioration de la conductivité thermique et de la rigidité des nanocomposites, ce qui est en accord avec de nombreux résultats expérimentaux récents.

Bien qu'il soit évident que l'interphase entre les charges et la matrice joue un rôle important dans les propriétés effectives des matériaux nanocomposites, il n'existe aucune information satisfaisante sur l'intensité de l'effet de l'interphase sur l'amélioration de la conductivité thermique et de la rigidité des nanocomposites. D'un point de vue expérimental, il est extrêmement difficile d'obtenir une connaissance complète de l'interphase et donc de son effet sur les propriétés effectives finales du nanocomposite. Pour acquérir des informations générales sur l'intensité de l'effet de l'interphase sur les propriétés effectives du nanocomposite, nous avons développé une modélisation par éléments finis en 3D pour étudier l'effet de la présence d'une interphase sur les propriétés effectives de matériaux nanocomposites. La figure 1 donne des exemples d'un volume élémentaire représentatif (VER) cubique avec différentes géométries de charges orientées aléatoirement et recouvertes d'une couche représentant l'interphase.

L'interphase est alors considérée comme une couche homogène et isotrope d'épaisseur fixée qui couvre toute la surface des charges. Nous observons que pour les nanocomposites comportant des charges avec un facteur de forme effectif élevé (long cylindres et disques fins), l'effet de l'interphase est négligeable sur la valeur calculée des propriétés. A l'inverse, pour les charges sphàériques, l'amélioration des propriétés effectives du nanocomposite est très dépendante de la différence de valeur entre les propriétés de la matrice et de l'interphase mais aussi dépend fortement de l'épaisseur de l'interphase. Nous avons aussi proposé une approche statistique pour évaluer les propriétés effectives de nanocomposites renforcés par des particules sphériques. Comme illustré à la figure 2, nous avons sélectionné dans cette approche une série d'images de microscopie électronique à balayage (MEB) enregistrées sur un échantillon d'un nanocomposite. Chaque image est considérée comme un VER, et la série comme représentative de la microstructure du composite. Ensuite, chaque image a été digitalisée pour distinguer les charges de la matrice tout en respectant la fraction volumique des charges dans la matrice. Par la suite, nous avons utilisée une approche statistique par la méthode du « strong contrast » pour évaluer le module élastique effectif de chaque image MEB. Enfin, nous avons fait appel à un modèle 2D par éléments finis pour évaluer le module élastique effectif de l'ensemble des éléments de nanocomposite représenté par les images MEB sélectionnées. La validité de cette approche a été vérifiée en comparant les résultats issus de la simulation avec des résultats expérimentaux obtenus sur un matériau époxyde renforcé par des particules de silice. Nous avons montré que pour le matériau nanocomposite réel renforcé par des particules sphériques, un calcul par éléments finis en 3D avec des particules parfaitement dispersées dans la matrice peut éventuellement sous-estimer les propriétés du composite par rapport aux résultats expérimentaux. Ainsi, notre étude suggère que, pour ces cas, on doit prendre en compte la microstructure réelle du nanocomposite pour estimer ses propriétés. diminution significative, d'environ la moitié, de sa conductivité thermique. D'ailleurs, une concentration d'atomes d'atomes de bore en substitution de seulement 0,75% entraîne une réduction de la conductivité thermique du graphène de plus de 60%. Nos résultats suggèrent également qu'une concentration proche de 1% d'atomes de bore ou d'azote supprime l'effet de chiralité de la conduction thermique du graphène.

Nos simulations de dynamique moléculaire montrent qu'en augmentant la concentration de défauts, le module élastique décroit graduellement. Dans le cas de défauts de type lacunes doubles et lacunes simples, le module élastique décroît d'environ 20 % pour une concentration de défauts de 2% ; cette décroissance est de seulement 10% pour la même concentration de défaut de Stone-Wales. De plus, la résistance à la traction et la déformation à rupture du graphène diminue brusquement d'environ 20% pour une concentration de défaut de seulement 0.25%. Par ailleurs, nous avons calculé que la baisse de la conductivité thermique du graphène suit une forme exponentielle en fonction de la concentration de défauts. Nos résultats suggèrent qu'une concentration de défauts de seulement 0,25% dans le graphène entraîne une baisse significative de sa conductivité thermique d'environ 50%. Nos avons également montré que l'effet des défauts étudiés sur la conductivité thermique est proche quel que soit le type de défaut.

De plus, nous avons étudié en particulier l'effet de la courbure sur la conductivité thermique du graphène et avons observé que le graphène peut conserver sa haute valeur de conductivité thermique même si sa courbure atteint de grands angles.

Les polymères époxydes se caractérisent par leur excellente adhésion, leur grande rigidité et leur bonne résistance mécanique, leur bonne propriété d'isolation électrique et leur grande résistance aux produits chimiques. Les époxydes sont de nos jours largement impliqués dans de nombreuses applications allant de la microélectronique à l'aérospatiale. Les polymères époxydes sont des thermodurcissables et sont obtenus en mélangeant deux composants distincts, qui sont nommés résine et durcisseur (aussi appelé agent de réticulation). Le choix de la composition optimale de résine et de durcisseur pour atteindre une valeur de propriété donnée est souvent réalisée à partir d'une approche purement empirique, qui est à la fois couteuse en temps passé et en argent. Ainsi, une approche par dynamique moléculaire peut être utilisée afin de prédire quelle est la composition optimale pour atteindre une valeur de propriété donnée. A cette fin, nous proposons une méthode simplifiée qui permette de construire in silico une structure moléculaire constituée de DiGlycidyl Ether de Bisphenol A (DGEBA) réticulé avec de la DiEthylèneTriAmine (DETA). Notre but n'est pas de reproduire la cinétique exacte de la réaction de réticulation, mais de proposer un schéma réactionnel efficace pour obtenir une

Introduction

During the last decade, nanotechnology has developed into a technology with a promising future.

Nanotechnology is one of the important branches of science with specialized research disciplines, such as nanooptics, nanobiotechnology, nanomedicine, nanoelectronics, and nanomaterials. In our time, nanomaterials can be found in many products for daily use, from sun cream to toothpaste with a repair effect, to wall paint which repels dirt, and to new materials for medical implants. In comparison with bulk materials, nanomaterials present considerably higher surface to volume ratio which significantly improve their reactivity. In order to design nanostructured materials for commercial applications, the crucial issue is to obtain comprehensive knowledge about the nanomaterials properties. The main challenge is that the experimental characterizations of materials at nanoscales are complex, time consuming and expensive as well. One solution for this case is to use modeling approaches as experimental alternatives.

Recent advances in the fabrication of nanoscale materials with extraordinary high thermal and mechanical properties such as graphene [1,2], graphite nanosheets, carbon nanotubes [3], nanoclays [4] and metal oxide nanoparticles, has motivated an ongoing demand for the reinforcement of electrical, thermal and mechanical properties of polymer based materials which are called as nanocomposites. For designing of a nanocomposite material with desired effective properties, it is important to select appropriate nanoparticles and also decide about the volume fraction of additives.

We should note that since the nanoparticles are expensive, trial and error experimental tests would be so much expensive and time consuming as well. For these cases, the modeling is inevitable in order to improve the efficiency of nanocomposites design.

The objective of this PhD has been to use combination of atomistic and continuum methods for the evaluation of thermal conductivity and elastic stiffness of nanostructured materials. We used molecular dynamics method to probe and investigate the materials properties at atomic scale. In the next step, the bulk properties were evaluated using finite element, Mori-Tanaka and strong-contrast, continuum approaches. In order to evaluate the accuracy of the developed modeling tools, the modeling results were compared with experimental tests. The details of this PhD thesis are presented in five chapters. We should note that the conducted studies have been already published or submitted in peer-reviewed international journals. Here, we present a brief description of the each chapter of this thesis:

Chapter 1, theoretical background of the methods that have been used in our studies is presented.

First, the basic information about the molecular dynamics method is presented. We describe that how the dynamics and motions of atoms are tracked by the molecular dynamics simulations. We also presented a short and also comprehensive viewpoint about the potential functions that are used in molecular dynamics simulations. Then, in the remaining part of this first chapter, the details of continuum approaches that have been used in our studies are presented. In our works, we used finite element, Mori-Tanaka and strong-contrast methods for the evaluation of bulk properties of nanocomposites. It is worth mentioning that for the evaluation of two-phase composite materials properties, the analytical micromechanics theories such as Mori-Tanaka [5,6] and the Halpin-Tsai [7] have been widely used [8][9][10]. However, these methods cannot accurately consider the interactions between adjacent inclusions and they have limitation in the evaluation of the microstresses involved with individual inclusions. During the last decade, because of considerable improvement of computational tools, numerical methods, such as finite element (FE) method, have been also widely used for the modeling of thermal conductivity and elastic properties of composite materials [10][11][12][13][14][15]. Using the finite element approach, it is possible to evaluate the micro-stresses and we could also consider the adjacent inclusions effects on the effective properties. The main disadvantages of finite element in comparison with Mori-Tanaka and Halpin-Tsai methods is its higher computational costs as well as the complexities for the modeling of high volume concentrations and aspect ratios for inclusions. Statistical continuum mechanics techniques also provide other tools for characterization and reconstruction of heterogeneous materials based on statistical correlation functions [16].

Chapter 2, in this chapter we studied the thermal conductivity and elastic response of nanocomposite materials using the continuum approaches. First, we compared three different groups of modeling techniques for the evaluation of the effective thermal conductivity and elastic modulus of random two-phase composite materials. These three groups are categorized as: mean field methods (by Mori-Tanaka), numerical approaches (by finite element) and statistical continuum methods (by strong contrast). Such a comparison provides a global viewpoint about the different modeling techniques and their advantages and limitations for the evaluation of nanocomposites bulk properties.

It is well known that the interphase between the fillers and matrix plays an important role on the nanocomposites effective properties. However, no adequate information exists about the intensity of interphase effect on the final reinforcement in thermal conductivity and elastic stiffness of nanocomposites. From the experimental viewpoint, it is considerably difficult to obtain comprehensive knowledge on the interphase effect on the final composite effective properties. This way, we studied the interphase effects on the effective thermal conductivity and elastic stiffness of nanocomposite materials using 3D finite element approach. In our works, the interphase was introduced as the homogenous covering layer of the outer surface of fillers within a distinct thickness.

We also developed a global statistical approach for the evaluation of elastic properties of silica/epoxy nanocomposites. In this approach, we selected a series of scanning electron microscopy (SEM) images of a nanocomposite sample as a representation of the microstructure (representative volume element, RVE). In the next step, we used statistical continuum theory by strong contrast method for the evaluation of the effective elastic modulus of each SEM image. In the final step, we used 2D finite element modeling for the evaluation of the effective elastic modulus of the nanocomposite sample consisting of the preselected set of SEM images (RVE). Comparison between the experimental measurements for the nanocomposite elastic modulus and modeling results by the global statistical approach reveals a fine agreement. This observation reveals the importance of nanocomposites microstructure on the effective properties for polymer materials filled with spherical nanoparticles.

Finally, we developed theoretically based multiscale approaches for the investigating of the thermal and mechanical reinforcement of poly lactide (PLA) by addition of expanded graphite (EG) platelets. In the multiscale modeling, first, we evaluated the thermal conductivity and stiffness tensor of graphite by means of molecular dynamics simulations. The molecular dynamics results were found to be in close agreements with experimental results available in the literature. Next, we developed 3D finite element models for the evaluation of thermal conductivity and elastic modulus of the EG/PLA nanocomposites.

Chapter 3, we used classical molecular dynamics simulations for the investigation of tensile response and thermal conductivity of graphene. Graphene [1,2], is the two-dimensional structure of honeycomb carbon atoms. Analogous to other members of carbon family, such as carbon nanotubes (CNTs) and graphite, graphene presents exceptionally high thermal conductivity [17],

mechanical [18] and electrical properties [19]. Recent experimental studies [17,20] have shown that the single-layer graphene (SLG) can present extremely high thermal conductivity of around 4100±500 W/m-K which outperforms all other known materials. This observation proposes the graphene for micro/nanoscale electronic device heat transfer blocks in response to thermal management concerns in the electronics industry. Using atomic force microscope nanoindentation experiments, Lee at al. [18] reported the elastic modulus of 1TPa and breaking strength of 130±10 GPa for defect-free graphene sheets. Combination of high thermal conductivity, elastic properties and surface to volume ratio, also propose the graphene as a unique nanofiller for reinforcement of thermal, mechanical and barrier properties of polymeric materials [21]. The optimized Tersoff potential developed by Lindsay and Broido [22] was used for introducing the carbon atoms bonding interactions. It is shown that the thermal conductivity and mechanical properties predicted by the optimized Tersoff potential are in considerable agreements with experimental results. We particularly studied the effect of boron and nitrogen atoms substitution on the thermal conductivity and tensile response of graphene. These chemically modified graphene structures are expected to have application in the next generation of electronic devices. Boron and nitrogen are the excellent elements for doping of graphene because of their close atomic size and strong valence bonds with carbon atoms.

Similar other known materials, different types of defects exist in experimentally fabricated graphene sheets. Defects in the graphene are usually in the form of point vacancy, bivacancy, Stone-Wales and grain boundaries. Accordingly, we also studied the effect of point vacancy, bivacancy and Stone-Wales defects concentrations on the thermal conductivity and tensile response of graphene.

Moreover, recent atomic-scale microscopy and spectroscopy studies have shown significant deviation of graphene from its planar structures in forms of ripples, kinks, and folds [23,24] which may alter its properties [25]. This way, we particularly studied the effect of curvature on the thermal conductivity of graphene.

Chapter 4, in this chapter we used molecular dynamics simulations for the construction of the crosslinked epoxy networks. We should note that epoxy polymers are known for their excellent adhesion, chemical resistance and its considerable mechanical properties. Epoxy is now widely used in various applications ranging from microelectronics to aerospace. Being in liquid phase before the curing process has made the epoxy a favorable material for the fabrication of nanocomposite materials. The OPLS united-atom force field [26,27] was used to define the bond, angle and dihedral parameters with Lenard-Jones non-bonded interactions. The dynamic crosslinking approach was used to construct the epoxy networks. We also calculated the glass transition temperature of epoxy. It was found that the density of epoxy increases by increasing the crosslinking density.

Chapter 5, we developed a novel multiscale scheme on the basis of molecular dynamics and finite element methods for the modeling of effective thermal conductivity of graphene epoxy nanocomposites. The proposed hierarchical multiscale approach includes three different scales.

First, we used molecular dynamics simulations for the investigation of thermal conduction in graphene epoxy assembly at atomic scale. To this aim, we constructed crosslinked epoxy and graphene assembly molecular models. Molecular dynamics results suggest that the thermal conductivity of single layer graphene, in epoxy matrix, decline by around 30%. Using the developed molecular dynamics model, we also calculated thermal boundary conductance (TBC) between crosslinked epoxy and the graphene sheet. Based on the results acquired by the molecular dynamics method, we developed finite element based representative volume elements (RVE) in order to evaluate local thermal conductivity constants of the nanocomposite. In this modeling, TBC values obtained by the molecular dynamics simulations were used for introducing thermal contact interactions between graphene fillers and epoxy matrix in the finite element modeling. Moreover, graphene thermal conductivity was modified on the basis of molecular dynamics results. At the final step, the nanocomposite effective thermal conductivity was calculated using the finite element approach for homogenization of anisotropic thermal conductivity constants of 3 different microscale RVEs. The proposed finite element methodology in this work could be used as an accurate method for the evaluation of composite effective properties filled with high volume fraction and high aspect ratio of fillers. Comparison between our multiscale modeling results with experimental results [28] for a wide range of graphene concentrations confirm the validity of the proposed multiscale scheme in the evaluation of thermal conduction of polymer nanocomposites.

We also studied the effect of formation of covalent bonds between graphene and epoxy atoms on the effective thermal conductivity of nanocomposites.

Theory 29 In Fig. 1, the different methods that are used in different scales in the multiscale modeling of nanostructured materials are illustrated. We can divide the modeling of nanostructured materials into three-different levels. At the deepest level, the electronic structures of materials are studied using the quantum mechanics (QM) theories. These kinds of calculations are carried out at very short times in the order of femtosecond (10 -15 of a second) and for atomic sized samples. The next step is molecular step, where the molecular dynamics (MD) based methods are used to study the materials properties at nanoscale. In the macroscopic level, the continuum approaches like finite element (FE) are used to evaluate the bulk properties of materials. Since in this PhD work, the electronic structures of materials are not concerned, we focused our research on the higher scales (molecular and continuum scales). However, it is worth mentioning that the potential functions that are used in the molecular dynamics simulations are commonly developed on the basis of quantum mechanics calculations or lattice dynamic theories.

Fig.1, Multiscale modeling

In this chapter, the details of used techniques in the multiscale modeling of thermal conductivity and mechanical response of nanocomposite materials are presented. To this aim, firs, the basic information about the molecular dynamics method is presented. Then, in the remaining part of this chapter the details of continuum approaches are presented. In our works, we used finite element, Mori-Tanaka and strong-contrast methods for the evaluation of bulk properties of nanocomposites.

These three methods could be considered as the representation of three different groups of modeling techniques in the evaluation of bulk properties of nanocomposites. Finite element is the most popular and successful numerical method that has been widely used in the modeling of nanocomposites. The Mori-Tanaka method is the most well known technique among the proposed mean field approaches. The strong-contrast is the modified version of weak-contrast formulation and it is categorized as statistical continuum theories.

Molecular dynamics simulation

Molecular dynamics is a computer simulation technique in which the time evolution of molecules and atoms is tracked by integrating their equations of motion. The molecular dynamics is on the basis of statistical mechanic theories, which is a well established fundamental science in the physics and chemistry. The key factor in the molecular dynamics is the motion, which indicates that how positions, velocities and orientations change with time. Molecular dynamics simulation is the modern realization of an old essentially idea in science; that the behavior of a system can be computed if we have for the system's parts, a set of initial conditions and also the interacting forces [1]. The high computational costs of molecular dynamics simulations impose limits on the number of atoms that are used for the modeling. Conventionally, between 1,000 and 1,000,000 atoms are modeled during the very short time scales of 1 ns to 1 μs (between 10 -9 s and 10 -6 s) and it could be stated that it is physically impossible to develop molecular models at millimeter scale. Despite of size limitation of molecular dynamics simulations, by applying the periodic boundary conditions it is possible to evaluate the bulk properties of materials (infinite number of atoms) using limited number of atoms (several thousands of atoms). In this section, we present a brief description of the molecular dynamics simulation technique, the more detailed information could be found in the molecular dynamics text books.

Equation of motion

In the molecular dynamics simulation method, the motion of individual atoms within an assembly of N atoms or molecules is modeled on the basis of either a Newtonian deterministic dynamics or a Langevin-type stochastic dynamics, given their initial position coordinates and velocities [2]. The energies and dynamics of the atoms are obtained from prescribed two-body or many-body interatomic potentials, U, from which the Newtonian forces experienced by these atoms are derived:

Theory 31
where r ij is the distance between two atoms i and j. The assumption is made that each atom interacts with its nearest neighbors, within a distinct cut-off radius.

The 3N position degrees of freedom set up 3N simultaneous coupled differential equations of motion which are integrated numerically by the use of a variety of finite-difference techniques (for more information refer to Haile [1]). For this purpose, the velocity Verlet algorithm is the most popular method. According to this algorithm, the positions, r i , and velocities, v i , of the atoms of mass m i are update at each simulation time increment, d t , by:

In the molecular dynamics systems, the controlling of temperature and pressure are the crucial points in the evaluation of thermal and mechanical properties. For these purposes, various methods have been developed. The Nose-Hoover dynamics is the most well-known and popular method for the controlling the temperature (thermostat) and pressure (barostat) fluctuations in the molecular dynamics simulations. As an example, the Nose-Hoover thermostat (NVT) method is presented here. In this method, the simulation box is placed in a heat bath at temperature of T. If we assume that the boundaries are allowed to only exchange the energy and not matter, then when the simulation box is equilibrated with the bath, they all attain a common temperature, T. Since, there is a thermal contact between the two, the total energy of the system fluctuates. The velocity verlet formulation of Nose-Hoover thermostat is expressed as follow [2]:

here, Q is a parameter which behaves like the 'mass' associated with the motion of the coordinate and k B is the Boltzman constant. The γ is called the friction coefficient of the bath. This coefficient is not a constant and can take both positive and negative values. From Eq. 7, if the kinetic energy is greater than , the friction coefficient γ, will be positive which prompts a friction inside the bath and accordingly the velocities of the atoms are decreased to lower their kinetic energy to that of the bath. On the other hand, if the kinetic energy is lower than , the friction coefficient γ, will be negative and the bath will be heated and accelerate the motions of atoms. A particular parameterization of the Q is given by: where, τ, is the relaxation time of the heat bath, normally of the same order of magnitude as the simulation time step, dt. The relaxation time, τ, controls the speed, with which the bath damps down the fluctuations in the temperature. This way, is a parameter that should be obtained by trial and error in order to have least deviations from the desired temperature, T. In the Nose-Hoover thermostat and barostat (NPT) method, the simulation box changes to control the pressure of the system (for more information refers to Allen and Tildesley [3]).

Periodic boundary conditions

The macroscopic materials are built by a large number of atoms in the order of 10 23 atoms, which is by a long way far from the molecular dynamics models. The main difference between the nanostructures and bulk materials is the higher ratio of surface atoms in the nanostructures. In another word, the ratio of surface atoms to the total numbers of atoms in the nanostructures would be much larger than bulk materials, causing surface effects to be much more important. Therefore, if the surface atoms effect become negligible, it is possible to evaluate the properties of bulk materials by the means of molecular dynamics simulations. A convincing solution to this problem that is routinely utilized is to use periodic boundary conditions (PBC). When using PBC, particles are enclosed in a box and we can imagine that this box is replicated to infinity by rigid translation in all the three Cartesian directions, completely filling the space. The key point is that now each particle i Theory 33 in the box should be thought as interacting not only with other particles j in the box, but also with their images in nearby boxes [4]. The schematic illustration of PBC is shown in Fig. 2. By the use of PBC the surface atoms effect would be eliminated completely and it is possible to evaluate the bulk properties of materials.

Fig. 2, The schematic illustration of periodic boundary condition.

Inter-atomic potentials

As expressed in Eq. 1, the atomic forces are derived from the inter-atomic potential functions.

These atomic forces are then used to obtain the atomic velocities and update the atomic positions as well. This way, the motions of atoms in the molecular dynamics simulations are strongly dependent on the used potential function. This accordingly means that the predicted properties by the molecular dynamics are intrinsically dependent on the used potential functions for introducing the atomic forces. All mechanical forces in a system of particles can be classified as conservative or nonconservative. Conservative forces are those whose work depends only on positions of the particles, without regard to their instantaneous velocities and trajectories of passage between these positions. The dissipative and gyroscopic forces that are not conservative are called nonconservative [5] (such as mechanical friction). For a system characterized by only conservative forces, there exists The first term represents the energy due to an external force field, such as gravity. The second term

shows potential energy of pair-wise interaction of the particles; the third gives the three-body components, and so on. Accordingly, the function U 1 is also called the external potential, U 2 the pair-wise, and Um at m > 2 a multibody potential.

The Lennard-Jones (LJ), Morse and columbic potentials are the most well known two-bondy interactions. The LJ and Morse potential are categorized as the short-range potentials while the columbic interactions have longe-range nature. The LJ potential is expressed as follows:

Where ε is the depth of the potential well and σ is the equilibrium distance at which the inter-particle potential is zero. The Morse potential has the following form:

The LJ and Morse potentials are commonly used to model the non-bonding interactions in the molecular models. The potential energy and force forces for LJ and Morse potential are shown in Fig. 3. Since in the LJ potential the force and potential energy approaches to zero at 2.5σ, the interactions at higher distances are not computed. The higher-order terms of the potential function (m > 2) are used for the modeling of solids with more complex molecular structures to account for chemical bond formation, topology and spatial arrangement.

In this way, in the three-body potentials, the forces between atoms i and j are not only dependent on the positions of atoms i and j and it is also dependent on the position of atom k. The Tersoff potential [6,7] is the most popular three-body potential. This potential was first developed for the modeling of silicon and carbon atoms bonding interactions. The Tersoff potential is expressed [6,7] as follows:

where, A is a constant, and the three-body coordination of atoms is shown in Fig. 4.

Fig. 4, Three-body coordination of particles [5].

The Tersoff potential mixing rules for the modeling of dissimilar atoms are as follow [8]:

Another special form of a multi-body potential is provided by the embedded atom method (EAM) potentials for metallic systems. One appealing aspect of the EAM potential is its physical picture of metallic bonding, where each atom is embedded in a host electron gas created by all neighboring atoms [5]. In the EAM potentials, the total energy of the system is written as:

Where, ρ h,i is the host electron density at atom i due to all other background atoms in the system, L i [ρ] is the energy to embed the atom i into the background electron density ρ, and ϕ ij is a pair-wise central potential between atoms i and j [2]. The host electron density is a linear superposition of the individual contributions, and is given by: where ρ* is the electron density of atom j as a function of inter-atomic separation.

A sample of many-body and longe-range potential for the modeling of metals is Sutton and Chen (SC) potential. The total energy in the SC potential is written as:

We should note that for the modeling of polymers and bio-molecules, the atomic interactions are modeled by the use of molecular mechanics force fields. In the chapter of 4, a sample of molecular mechanics force field namely OPLS united atom force-field is introduced and used for the modeling of crosslinked epoxy network.

Finite element modeling

Now-a-days, finite element method is considered as a versatile tool for the modeling and simulation of a wide range of engineering problems. Computational limitations and modeling complexities of finite element method impose limits on the maximum number of elements, which are used for introducing the geometries in the model. In this way, the simulations of composite materials are limited to modeling of a representative volume element (RVE) of the system. The finite element simulations in our studies have been performed using ABAQUS (Version 6.10) package. In an attempt to reproduce the RVE in a status closer to those in experimentally fabricated random composites, the 3-Dimensional (3D) inclusions were randomly distributed and oriented in the RVE.

It is worth mentioning that while 2-Dimensional (2D) finite element models could be geometrically acceptable for platelets fillers, they could not accurately describe the geometry for cylindrical and spherical inclusions. Furthermore, Hbaieb et al. [9] has shown that in the case of randomly distributed clay particles, 2D finite element models considerably underestimate the predictions by 3D models. Therefore, the modelings were carried out in ways to take into account different geometries for fillers, as cylinder, sphere and plate. In our studies, the inclusion geometry is introduced by the use of aspect ratio. For inclusions with platelet shape, the aspect ratio is defined as the diameter to thickness ratio of the plate. In the case of cylindrical inclusions, the aspect ratio corresponds to the length to diameter ratio of inclusion. Accordingly, the aspect ratio of one is representative of spherical inclusions. In Fig. 5, samples of 3D cubic RVE with different geometries for inclusions are illustrated. The size of the RVE is adjusted based on the geometry, volume fractions and numbers of perfect particles inside the RVE. To accurately create the desired volume concentration of fillers inside the RVE, the RVE is constructed in a way to satisfy the periodicity criterion. This means that if a filling particle passes one boundary side of the RVE, the remaining part of that particle continues from the opposite side. This accordingly means that if we put these RVE cubes together side by side, no discontinuity will be observed in the constructed sample and all of the particles will have equal geometries. This way, the developed RVEs could also be considered as the representation of a bulk composite material. We should note that we have not applied the periodic boundary conditions in our modeling. The random RVEs were constructed in ABAQUS by developing Python scripts as input files. We developed numerical tools (in C++ language) for creating randomly distributed/oriented hard core inclusions inside the RVE considering the periodicity criterion. The modeling parameters are adjusted in the C++ programs. The provided information by the C++ programs are written in a text file that is later used as an input file for Python scripts for the construction of the final finite element RVEs. We found that the combination of Python and C++ is a faster approach than exclusively relying on Python scripting.

As shown in Fig. 5(d-f), the specimens were meshed using 4-node linear tetrahedron shape elements.

It is worth mentioning that the meshing of the RVEs with cylindrical fillers is considerably more difficult in comparison with RVEs filled with spherical and platelet particles. For the cylindrical fillers, each filler should be partitioned to four symmetric parts (this is also recommended for platelet fillers), otherwise, the meshing of the RVE samples would be almost impossible. The meshing becomes trickier as the fillers length and volume fraction increase and especially when the interphase is also included. For the evaluation of mechanical properties, 3-dimentional stress (C3D4) elements are used. The heat transfer (DC3D4) elements are used in the modeling of thermal conductivity.

For the evaluation of elastic modulus, a small uniform displacement was applied on one of the surfaces along the Z-direction and the remaining surfaces where fixed only in their normal directions. The reaction forces on each surface were calculated in order to obtain the corresponding stresses. Using Hooke's law for an isotropic material, the effective elastic modulus and Poisson's ratio of the RVEs were estimated. For the evaluation of thermal conductivity, a constant heating surface heat flux (+q) was exerted on one of the surfaces along the Z-direction while on the opposite surface the same magnitude cooling surface heat flux (-q) was applied. This introduces a heat flow along the specimen by transferring a constant quantity of heat flux. The total energy of the simulation RVE remains constant and steady state heat conduction is observed. In response to the energy redistribution, thermal energy moves from the hot reservoir to the cold reservoir and a temperature gradient is established in the system. The temperature gradient, , was evaluated by calculating the temperature differences between the hot and cold surfaces. The effective thermal conductivity of the nanocomposite was obtained using the one-dimensional form of the Fourier law, .

Mean field homogenization: Mori-Tanaka approach

Mean field micromechanics methods provide simple tools for the evaluation of the effective thermal and mechanical properties of composites. A typical mean field formulation is based on the inclusion concentration tensors that connect inclusion averaged fields in reinforcements or matrix with the Theory 40 corresponding macroscopic fields. One of these approaches, Mori-Tanaka (MT), approximates the behavior of composites that contain reinforcements at non-dilute volume fractions via dilute inhomogeneities that are subjected to effective matrix fields rather than the macroscopic fields [10].

These effective fields account for the perturbations caused by all other reinforcements in a meanfield sense [11]. Benveniste [12] expressed the MT methods for elastic composites by the following dual relations:

where are the averaged strain and stress tensors of the inclusion and matrix phases, respectively. are macroscopic second order strain and stress tensors, stand for the MT matrices and stand for the inclusion dilute, fourth order strain and stress concentration tensors. The analogous form of the Mori-Tanaka formulation for the thermal conduction has the form [11]:

where and q are the thermal gradient and the heat flux vectors, and stands for the thermal gradient and heat flux concentration second order tensors. The Mori-Tanaka stress concentration tensors of the matrix and inclusions can be expressed by:

where I is the fourth order identity tensor and v f stands for the inclusion volume fractions. The dilute inclusion concentration tensor, , can be obtained in analogy to Hill's [36] expressions as:

where C i and C m stands for the stiffness tensor of the inclusion and matrix respectively and S is the Eshelby tensor [14]. An analogous relation for the inclusion gradient concentration tensors can be obtained by considering the I as the second order identity tensor and replacing the inclusion and matrix stiffness tensors by the conductivity tensors K i and K m , respectively and considering the S as the second order Eshelby tensor for diffusion problem [15].

Due to random orientations of particles inside the RVE, the Mori-Tanaka model results must be representative of the integral of all directions. For the second and fourth order tensors, the volume average can be respectively expressed by:

Here the braces denote the average over all possible orientations, are Euler's angles between the local and global coordinate systems (as shown in Fig. 6) and Q is the corresponding transform matrices. Finally, the Mori-Tanaka macroscopic effective stiffness and conductivity tensors of the composite with randomly oriented inclusions can be expressed by [16]:

Fig. 6, Vector representation in spherical coordinate

Statistical homogenization: Strong contrast method

In our studies, we used statistical continuum theory with strong contrast to predict effective properties of two-phase composite structures. In the statistical continuum theory, correlation functions are used as the statistical description of the microstructure. Generally, probability of occurrence of N-points, which are invariant relative to each other in desired phases, defines the N-point correlation functions. The expression of these functions for a given phase i can be written as follows [17]:

i N 1 2 N 1 2 N P (x , x ,..., x ) Pr obability(x Phase(i) x Phase(i) ... x Phase(i))        (36)
where x i is the position of random points (vectors) in the microstructure.

In multiphase materials, the first order correlation functions represent the volume fractions of different phases and do not describe any information about the distribution and morphology of phases. If N-number of random points are inserted within a given microstructure and the number of points in phase-i is counted as N i , the one-point probability function (P i ) is defined as the volume fraction through the following relation:

M (the total number) is increased to infinity and v i is the volume of phase i.

Here we would like to calculate two-point correlation functions from the two-phase SEM images and computer generated paticles. Now assign a vector, , starting at an arbitrary point in a heterogeneous microstructure. Depending on whether the beginning and the end of these vectors fall within phase-1 or phase-2, there will be four different probabilities (P 11 , P 22 , P 12 and P 21 ) defined as [17]: (38) where N ij is the total number of random vectors that are thrown in the microstructures beginning in phase-i (φ i ) and the end in phase-j (φ j ). Eq. 2 defines a joint probability distribution function for the occurrence of events constructed by two points ( and ) as the starting and end points of a vector when it is randomly inserted in a microstructure. The two-point function can be defined based on two other probability functions in the following form:

The first term on the right hand side is a conditional probability function.

More detailed morphological description of the heterogeneous materials is obtained by using higher order correlation functions.

In our works, we developed C++ codes to calculate two-point correlation functions using Monte-Carlo method. In this method, 3D randomly oriented inclusion (spherical, cylindrical and platelet) are generated and used to calculate the statistical two-point correlation functions. Using the Monte-Carlo approach, we also calculated the two-point correlation functions of the images. Higher order correlation functions can provide more precise morphological detail of heterogeneous systems. In this paper, three-point correlation function

  1 3 23
,, i P x x x is estimated using the two-point correlation function by the following relation [40]:

                    2 1 2 2 1 3 2 1 2 11 2 2 3 33 3 1 2 3 1 2 2 3 2 2 1 3 3 12 13 3 , , , , ,, ,,
i i i i i ii ii i P x x P x x P x x P x x P x x x W W P x P x P x x P x x W Px                        (40)
where the weight coefficients are defined as:

1 2 3 333 1 2 3 1 2 3 1 2 3 1 2 3 ,, R R R WWW R R R R R R R R R        (41)
where R 1 , R 2 and R 3 are the position vectors (as illustrated in Fig. 6). In the strong contrast formulation presented in this study, the subscripts m, i and c stands for the reference phase (matrix), inclusion phase and composite effective properties, respectively.

Strong contrast solution for thermal conductivity

Assuming isotropic properties for the matrix and inclusion, the effective thermal conductivity tensor, K C , for 3-dimensional heterogeneous materials is calculated using the strong-contrast formulation [18]:

    -1 2 1 1 1 1 1 3 2 2 1 1 1 1 1 1 2 1 2 1 2 2 1 1 2 2 1 2 3 1 2 2 3 1 2 2 3 1 2 1 2 3 ( , ) ( ) ( ) 1 I . 2 I I 3 ( , )d ( ) ( ) ( ) ( , , ) ( , ) ( , ) 9( ) ( , ) ( , ) ( ) ( ) ( ) ( ) ( ) C m C m m m P x x P x P x K M x x x P x P x P x P x x x P x x P x x M x x M x x P x P x P x P x P x K K K K K                             i i i m i i i im i i i m m im i i i i i 23 d d .... xx   (42)
In Eq. 42, I is the second-order identity tensor and im  is the polarizability [17] scalar which is calculated using the inclusion and matrix thermal conductivity values, K i and K m , respectively as follows:

im im im - β= +2 KK KK (43)
The second order tensor M m in Eq. 42 is defined as follows:

Theory 44 3 12 12 1 I M ( , ) 3 xx K xx      tt m m ( 44
)
where  is the total solid angle contained in a 3-dimensional sphere and

Strong contrast solution for elastic modulus

Strong-contrast method has been used to determine the effective elastic modulus of two-phase composites [19]. In this method, N-point correlation functions show up in the final equations that characterize the microstructure. The general equation for isotropic composite is written as follows [17,19]:

2 ( ) 2 i im im f h s f n n cm cm k v v B k              I (45)
As a common assumption, the calculations have been performed for the first and second terms of

() i n
B and other terms have been neglected [17,19]:

2 ( ) ( ) 23 ii im im f h s f cm cm k v v B B k            I (46)
Here f v is the inclusions volume fraction, h  and s  are the fourth order hydrostatic projection and shear projection tensors [17].

nm k and nm  are introduced as bulk and shear moduli polarizabilities where the subscript n represent c or i which respectively stand for effective and inclusion properties.

The polarizabilities are expressed in terms of the effective-or-phase bulk and shear moduli, F and G, as follows:

4 3 nm nm nm FF k FG    (47)   3 / 2 4 / 3 2 nm nm m m m n mm GG G F G G FG       (48)
In Eq. ( 46) the tensor coefficients () 2 i B and () 3 i B are given by the following integrals:

( ) ( ) ( ) 2 2 2 1 2 2 1 2 ( , ) ( , ) i m i f B dx U x x P x x v     (49) ( ) ( ) ( ) ( ) 3 2 3 1 2 2 3 3 1 2 3 1 ( , ) : ( , ) ( , , ) i m m i f B dx dx U x x U x x x x x v        (50)
in Eq. 5, the determinant tensor of correlation functions () 3 i  for the inclusion phase [17] is defined as follows:

( ) ( ) 2 1 2 1 2 () 3 1 2 3 ( ) ( ) 3 1 2 3 2 2 3 ( , ) ( ) ( , , ) ( , , ) ( , )
i i i i i P x x P x x x x P x x x P x x  (51)
the U (m) tensor in Eq. 49 and Eq. 50 for the reference phase [17] is based on the position-dependent fourth-order tensor H(r) and the related tensor for reference phase:

   ( ) ( ) ( ) () () ( ) ( ) 5 3 4 ( ) 3( 2 ) 3 5 () 3( 2 )
m m m ij m m m m im im mm m m im mm U r L H r G F G k H r FG G Hr FG                       (52)
H(r) is the symmetrized double gradient tensor [17]. The fourth order tensor L (m) , for the reference phase is expressed as:

  () 5 34 3( 2 ) m m m m im h im s mm G L F G k FG           (53)
In our studies, numerical integration by Monte-Carlo method is used for computing of B tensors (Eq. 50). For this case, due to considerable number of variables and six order integrations, the analytical integration methods are excessively complex. We found that by performing the integration for 20 millions random vectors ( ), the obtained values are completely converged. The latter computational point lies in the evaluation of effective elastic modulus which correlate to k cm and μ cm in Eq. [START_REF]LAMMPS user manual[END_REF]. In this case, we incrementally put various effective elastic modulus and Poison's ratio and try to minimize the error in the two sides of Eq. ( 46). The effective elastic modulus is the one which yields the insignificant differences in the two side of Eq. ( 46).

Introduction

In this chapter, we study the modeling of thermal conductivity and elastic response of nanocomposite materials using the continuum approaches. It is quite well known that the efficiency of reinforcement in the thermal and mechanical properties of nanocomposite materials strongly depends on the additives concentrations, geometry and properties as well. The fundamental understanding of the effective parameters on the final reinforcement is the key issue in the future applications of nanocomposites. Due to difficulties in experimental characterizations of nanocomposites, analytical and numerical simulations are getting more attractive as the experimental alternatives. For the evaluation of two-phase composite materials, the analytical micromechanics theories such as Mori-Tanaka [1,2] and the Halpin-Tsai [3] have been widely used [4][5][6]. These micromechanics methods have been modified to include the effects of interphase on the effective composite properties [7][8][9][10]. These methods could predict the effective properties of composite materials with reasonable accuracy. However, they cannot accurately consider the interactions between adjacent inclusions and they have limitation in the evaluation of the micro-stresses involved with individual inclusions [11]. The numerical methods, such as finite element (FE) have been also

widely used for the modeling of thermal conductivity and elastic properties of composite structures [6][7][8][9][11][12][13][14][15]. By the use of finite element method it is possible to evaluate the micro-stresses and also more elaborately take into account the adjacent inclusions effects on the effective properties.

Statistical characterization can be implemented for different types of random heterogeneous materials [17][18][19][20][21][22]. A wide variety of effective properties of heterogeneous materials can be estimated by exploiting different types of statistical functions. Two-point correlation functions are the simplest order of descriptors to describe the morphology and microstructure properties relationship in statistical continuum mechanics [17]. Two-point correlation functions could be evaluated from different techniques such as microscopy (SEM or TEM), X-Ray scattering and Monte Carlo simulations.

To the best of our knowledge, no study has been performed to compare 3D finite element, Mori-Tanaka and strong-contrast methods for the evaluation of effective thermal conductivity and elastic modulus of random two-phase composite materials. Therefore, in the section 2.1, we elaborately compare finite element, Mori-Tanaka and strong-contrast methods. To provide a better as well as more comprehensive viewpoint in comparison with the existing studies in the literature, the modelings were carried out for different geometries of the fillers, ranging from long cylinders to sphere and thin plates. We should note that in all the studied cases, the filler and matrix materials properties are assumed to be isotropic and the perfect bonding conditions are applied between fillers and matrix.

In term of reinforcement in thermal conductivity and elastic stiffness of polymeric materials, from theoretical point of view, the two-phase nanocomposites materials would not have any priority in comparison with two-phase microcomposite. However, in the case of formation of an interphase between the fillers and matrix with superior material properties than matrix, it would be expected that the nanocomposite materials present higher thermal and stiffness response. This is due to the existence of higher contacting surface of fillers with matrix and consequently higher interphase volume in the nanocomposite structures in comparison with microcomposites at the same fillers loading. Although it is obvious that the interphase also plays an important role in the nanocomposite materials, no adequate information exist about the intensity of interphase effect on the final reinforcement in thermal conductivity and elastic stiffness of nanocomposites. From the experimental viewpoint, it is considerably difficult to obtain comprehensive knowledge on the interphase effect on the final composite effective properties. The experimental characterizations for the evaluation of interphase effect are time consuming and the results are dependent on the availability of nanoparticles with controlled shape and size. Moreover, the quality of experimental results depends on the ability and competencies of experimental scientists in making homogeneous nanocomposites with controlled dispersion of fillers. In this frame, the use of theoretical and numerical approaches sounds promising. As it will be discussed in section 2.1 [11], we conclude that despite of modeling complexities and also computational difficulties, the finite element estimations are more reliable and promising in comparison with Mori-Tanaka and strong contrast methods.

Therefore, in section 2.2, we investigate the interphase effects on the effective thermal conductivity and elastic stiffness of nanocomposite materials using 3D finite element approach. The interphase is introduced as the homogenous covering layer of the outer surface of fillers within a distinct thickness. Then, we studied the effects of fillers geometry, volume fraction and properties contrast and more elaborately the effect of interphase thickness and properties contrast on the effective thermal conductivity and elastic response of the nanocomposite structures. We should note that in all the studied cases, the filler, interphase and matrix materials properties are assumed to be homogeneous and isotropic as well. Moreover, the perfect bonding conditions are applied between fillers, interphase and matrix.

Using the developed continuum approaches for the nanocomposite, we carried out two case studies on the silica/epoxy and expanded graphite/PLA nanocomposites. First, in section 2.3, we propose a global statistical approach for characterization of elastic properties of nanocomposite materials for taking into account the effects of real microstructure of the materials in the effective properties. This approach combines local statistical analysis using strong-contrast method with a 2D finite element analysis used to obtain ensemble average. In section 2.4, we use 3D finite element method for the modeling of expanded graphite/PLA nanocomposites. The transversely isotropic properties of expanded graphite are evaluated using the molecular dynamics approach. In the last section of this chapter, the obtained concluding remarks on the continuum modeling of nanocomposite materials are presented.

Comparison of finite element, Mori-Tanaka and strong-contrast methods

Convergence of finite element results

In the finite element modeling, one should verify the dependency of the reported results to the RVE size. The size of RVE naturally correlates to the number and geometry of the considered particles.

For spherical inclusions, the RVE size is only dependent on the number of spheres and volume fraction as well. For cylindrical and platelet inclusions, by increasing the volume fraction and aspect ratio of the inclusions, we must also check that the RVE size must be higher than the maximum dimension of the particles. Otherwise, there might be the probability that a particle solely connect the two opposite faces together. In such cases, the inclusion is directly involved in the load and heat transfer from one end to the other, which is not an acceptable status for random composite materials.

Fig. 1, shows the predicted elastic modulus and thermal conductivity of the composite for different RVE sizes for spherical (Fig. 1a,b) and platelets inclusion with aspect ratio of 50 (Fig. 1c,d) with 1% volume concentration and inclusion to matrix materials properties contrast ( ) of 500.

The herein reported bounds (variation bars) are based on 6 different simulations with various initial configurations of particles inside the RVE. As an overall trend, by increasing the number of particles, the variation band of the finite element results decreases. In the case of spherical inclusions, the results obtained by RVEs with 30 particles are presenting negligible variations. For these considered RVE sizes, the variation of the predicted properties for the platelet inclusions is considerably higher than spherical inclusions with the same number of particles (see Fig. 1c,d). This 52 reveals that the deviation of the inclusion geometries from spheres, we could expect higher dependency of the finite element results on the RVE size. Development of larger RVE sizes, impose meshing difficulties and higher computational costs. The next factor that should be taken into account is the magnitude of contrast in the inclusion to matrix materials properties. If the properties of the inclusion are considered to be the same as those of the matrix, we could state that the composite effective properties would present no dependency to the RVE size. In this way, we could conclude that the size dependency of the finite element results gets more pronounced by increasing the aspect ratios and contrast in inclusion to matrix materials properties. Since the volume fraction is constant for different RVEs sizes, the variance of the finite element results could be an indication of the RVE anisotropic nature.

Modeling of thermal conductivity by finite element

The thermal conductivity of the composite structure is evaluated using the Fourier law. This law is based on the establishment of a uniform temperature gradient along the material as a result of applied external heat flux. In our finite element modeling, we directly exerted the heat fluxes on the two opposite walls of the RVE and the results are shown in Fig. 2a and Fig. 2b. As it can be clearly seen from Fig. 2b, the established temperature gradient does not show uniform configuration. This observation criticizes the validity of the use of Fourier law in our modeling. Furthermore, for the evaluation of thermal conductivity, we should calculate the temperatures at hot and cold surfaces. In this way, we need to calculate the surface that belongs to each node to obtain the accurate average temperature on the corresponding face. This needs the development of new python script with relatively complex algorithms. To address these issues, we devised new boundary conditions in our finite element modeling. In this method, we merged two thin auxiliary parts to the RVE hot and cold surfaces. As shown in Fig. 2c, these auxiliary parts have the same section size as that of RVE and their thermal conductivity are selected to be one million to that of matrix. In this case, the heat fluxes are directly applied on the outer surfaces of these two highly conductive parts. As illustrated in Fig. 2d, as a result of redistribution of applied heat flux through the auxiliary parts, established temperature gradient throughout the same RVE present considerably more uniform distribution in comparison with initial case (Fig. 2b). Furthermore, as revealed from Fig. 2c, the nodes distributions on the highly conductive parts are completely uniform. It is worth mentioning that the uniformity of the established temperature gradient throughout the specimen correlates to the inclusion's geometry and also their relative thermal conductivity to that of the matrix (see Fig. 2).

Comparison between finite element and Mori-Tanaka methods

In two phase composite materials, the stress and heat flux is transferred between matrix and inclusions through their contacting surfaces. In this way, the inclusions surface to volume ratio plays an important role on the final reinforcement. Since spherical inclusions have minimum surface to volume ratio, we expect to have the minimum reinforcement for this case. The obtained results by FE and MT methods for thermal conductivity and elastic modulus of two-phase random composite structures are compared in Fig. 3. We should note that the Poisson's ratio of the fillers and matrix were assumed to be 0.2. As expected, in all cases, the minimum reinforcement occurs for spherical inclusions.

In the case of inclusions with platelet shape, we could observe a fine agreement between finite element and Mori-Tanaka results when the contrast in materials properties is 50 (Fig. 3a and3c).

Such an observation is also valid for higher contrast in material properties of 500 and 1% volume concentration for inclusions (Fig. 3b and3d). For inclusions with contrast in materials properties of 500 and volume fraction of 3%, the fine agreements in the predictions of FE and MT method is existing only for platelets with low aspect ratios. In this case, we could observe that Mori-Tanaka results start to underestimate the finite element results at higher aspect ratios. We should note that in Continuum modeling of nanocomposites 55 the Mori-Tanaka method, the inclusions are assumed to be placed in an infinite matrix and the reference medium is assumed to be the matrix. In this way, the Mori-Tanaka method does not take into consideration the reinforcement of the reference medium. Therefore, we would expect that if the reinforcements in the composite properties are considerable, the Mori-Tanaka results would underestimate the finite elements results, which is in agreement with our results.

Fig. 3, Comparison between finite elements (FE) and Mori-Tanaka (MT) methods in evaluation of elastic modulus and thermal conductivity of two-phase random composite materials for different inclusions geometries, volume fractions of 1 and 3% and different contrast in material properties of 50 and 500. The finite element results have error bars that are not presented here.

On the other hand, for cylindrical inclusions, the close predictions of finite elements and Mori-Tanaka methods at low aspect ratios vanish completely at higher aspect ratios. In this cases, the finite elements estimations increase by increasing the inclusions aspect ratio up to a critical point, in which the results are no longer dependent on cylinders aspect ratio. This observation has been also reported by Jafari et al. [16]. As it can be observed from Fig. 3, for the higher contrast between inclusion-to-effective composite material properties, the estimated results by finite elements are converged at higher critical aspect ratios.

The results obtained by finite elements and Mori-Tanaka methods suggest strong correlation between composite effective thermal conductivity and elastic modulus to inclusions geometry. For spherical inclusions, the contrast in inclusion to matrix material properties present insignificant effects and the reinforcement is strongly dependent on inclusion volume fraction. The influence of contrast in material properties on the effective properties is getting more pronounced as the inclusion shapes deviate more from spherical geometry. If one considers the finite element modeling as a more accurate method, our results suggest the Mori-Tanaka as a fast and reliable alternative for the ranges that the overall reinforcements in the elastic modulus and thermal conductivity of the matrix materials are lower than approximately 50 and 80 percent, respectively, for platelet inclusions.

For cylindrical fillers, the Mori-Tanaka method is suitable for low aspect ratios. The main disadvantages of finite elements are its computational costs, size dependency of the results and limitations in the development of RVEs with high inclusions volume concentrations and aspect ratios. We should also note that for platelets inclusions, the results by Halpin-Tsai method considerably overestimate the Mori-Tanaka results [5].

The obtained results based on the finite element simulations, suggest that by employing the fillers with platelet geometries, one could achieve higher reinforcement in thermal conductivity and elastic properties in comparison with cylindrical inclusions. Interestingly, in the recent experimental study by Sahil and Balandin [23], it is discussed that using the planar structure of carbon nanostructures (graphene) as fillers inside the epoxy, the reinforcement in thermal conductivity could be higher than epoxy with carbon nanotubes. As it can be observed from the results illustrated in Fig. 3, such a finding could not be concluded from the estimations by Mori-Tanaka method. Especially for the cases that the contrasts in material properties are 500, the Mori-Tanaka method suggests the higher reinforcement percent by incorporation of cylindrical fillers rather than the platelet ones.

Convergence of two-point correlation functions

As discussed in section 1.4, the statistical two-point correlation functions are considered as the description of the microstructure in the strong contrast formulation. Therefore, in the strong contrast formulation, two microstructures with similar two-point correlation functions would present the same elastic modulus and thermal conductivity. In Fig. 4, we have plotted the computed two-point correlation functions (TPCF) for different number of particles (inclusions) and for inclusions with various geometries. As it can be observed, the calculated TPCFs are convincingly 57 independent of the RVE size. However, analogous to the constructed RVEs in the finite element, by having the fillers with higher aspect ratios we could expect higher dependency of the calculated TPCFs on the RVE size.

Fig. 4, The calculated two-point correlation functions (TPCF) for different RVE sizes for (a) cylindrical inclusions with aspect ratio of 50, (b) spherical inclusions and (c) platelets inclusions with aspect ratio of 50. The inclusions volume ratio is 3%.

Strong contrast results

In table 1, the predicted results by statistical strong contrast are compared with those of finite elements and Mori-Tanaka. In the case of spherical inclusions, the reinforcement of elastic modulus predicted by statistical method slightly overestimates the predictions by finite elements and Mori-Tanaka methods. In this case, the strong contrasts estimations for thermal conductivity are in close agreement with those of finite elements and Mori-Tanaka. For cylindrical and platelets inclusions, the strong contrast results show only a slight reinforcement in both thermal conductivity and elastic modulus in comparison with spherical inclusions, while for these cases, finite element and Mori-Tanaka methods present considerable reinforcement effects. This shows the strong limitation of the strong contrast method in describing the inclusions geometry effects. This matter could be also concluded from the results in previous studies [19] at low volume concentrations. We should note that Fullwood et al. [24] also reported the errors and limitations of the weak and strong contrast methods in describing the effect of contrast in properties for high volume concentrations of the fillers. Moreover, we found that two-point correlation functions cannot show the differences when the inclusions aspect ratio exceeds from 20 for both, cylindrical and platelet geometries. We should mention that in the strong contrast formulation, if one selects the inclusion as the reference phase, the inclusion shape will present more significant influence [18,21].

Modeling of interphase using 3D finite element approach

As it was already mentioned in section 1.2, the geometry of the fillers in our works is defined by their aspect ratio. In this work the interphase is modeled as the covering layer of the outer surface of filler within a specific thickness. The interphase thickness for cylindrical and spherical fillers is defined with respect to fillers radius. This way, if the fillers radius is one, the interphase thickness of 0.4 is equivalent with an interphase with outer diameter of 2.8. In a consistent way, for the disc fillers the interphase thickness is introduced with respect to the half of the disc thickness. The hardcore methods are used to generate fillers without intersection or contact (interphases can overlap).

In all the studied cases here, the Poisson's ratio of the fillers, interphase and matrix are assumed to be 0.3. Perfect bonding conditions are used to introduce the contact interactions between fillers, interphase and matrix.

The calculated elastic modulus and thermal conductivity of composite structures by the finite element method for randomly oriented fillers with different geometries and different fillers to matrix properties contrast are shown in Fig. 5. In this study, the effect of fillers to matrix properties contrast is studied for low contrast of 20 and high contrast of 1000. In the conventionally fabricated nanocomposites, silica nanoparticles present low elastic modulus of 70 GPa. This way, if one assumes the elastic properties of conventional polymer matrixes to be around 3.5 GPa (such as epoxy), the fillers to matrix properties contrast of 20 is suitable for silica nanocomposites. The graphene and graphite present high in-plane elastic modulus of around 1000 GPa [25]. if we assume the minimum elastic modulus of conventional polymer matrix to be 1 GPa, the contrast of 1000 could be considered as the upper bound. However, for the thermal conductivity of nanocomposites materials higher values of fillers to matrix contrast in properties still exist (epoxy-graphene nanocomposites). We could observe that for the contrast in materials properties of 20 after a critical aspect ratio for the both discs and cylindrical fillers the further increase of aspect ratio does not influence the effective thermal conductivity and elastic modulus of composite.

In Fig. 6, the effect of interphase on the effective elastic modulus and thermal conductivity of nanocomposite structures are presented for a constant interphase thickness of 0.4. For a better vision for the evaluation of interphase effect on the final nanocomposite properties, we introduced the interphase reinforcement ratio (IRR). This ratio is introduced as the reinforcement of composite structure containing the interphase and fillers with respect to the reinforcement only due to the incorporation of fillers (no interphase). The main reason to use this definition is to exclude the effects of fillers geometry and contrast in material properties from the interphase effect. For a composite property (C), the interphase reinforcement ratio is expressed by the following relation:

here, C if is the composite effective properties including the interphase and fillers, C f is the two-phase composite effective properties (no interphase) and M is the matrix properties which is always set at one in this study (normalizing value). The obtained interphase reinforcement ratios are illustrated for different filer to matrix (FMC) and interphase to matrix (IMC) contrast in material properties. As mentioned earlier, we could expect the priority of nanocomposites in comparison with microcomposites just in the case that the interphase presents higher material properties than matrix. However, in some cases the numerical calculation has revealed that the interphase properties could be lower than that of the matrix [26].

Since the chemical structure of interphase is close to that of the polymer matrix in the case of hard fillers (inorganic, metal oxides and carbon based fillers and etc), we could assume that the interphase properties should be in the same order of matrix. In this study, the interphases to matrix contrast in materials are studied in a range between 0.5 and 3. The results shown in Fig. 6 reveal that the interphase has the maximum effect for spherical fillers. As the filler geometry deviate more from the spherical shape for a constant interphase thickness, the volume of the interphase would be much higher. This way one could expect higher effect for the interphase by the increasing of the fillers aspect ratio which is not the case here. The obtained results suggest that the interphase present the minimum effect in the case of disc fillers. Interestingly, the values of interphase reinforcement ratio for elastic modulus and thermal conductivity are considerably close to each other. Taking into account the reported results in Fig. 5 and Fig. 6, suggest that in order to have higher reinforcement for the nanocomposite materials properties, the effective aspect ratio of the fillers plays the crucial role. The importance of the interphase characteristics is pronounced as the fillers geometries get closer to spherical shape. Since the interphase has the maximum effect for the spherical fillers, in Fig. 7, the interphase effects on the elastic modulus and thermal conductivity of nanocomposite structures filled with spherical particles are shown. As the fillers to matrix properties contrast have the minimum as well as insignificant effect for spherical inclusions, the low fillers to matrix properties contrast of 20 is used for this case. The effects of interphase to matrix properties contrast for constant interphase thickness of 0.5 and different volume loading of spherical fillers inside the matrix are shown in Fig. 7a and Fig. 7b. For the interphase to matrix properties contrast of 0.5, we could observe that the reinforcement due to the addition of fillers is almost neutralized by the low interphase properties.

Further increase of interphase to matrix properties contrast leads in higher reinforcement effects. In the Fig. 7c and Fig. 7d, the effect of interphase thickness at constant interphase properties contrast of 2 is presented. By increasing the interphase thickness, the interphase volume increases by power of 3, that is why we could observe considerably high reinforcement effect. The results shown in Fig. 7 clearly show the importance of interphase for the case of nanocomposite structures filled with spherical particles.

62

In Fig. 8, the finite element RVEs for the evaluation of interphase effects on the longitudinal elastic modulus and thermal conductivity of nanocomposite filled with unidirectional fillers is shown. For this case, the longitudinal thermal conductivity is evaluated using the same procedure described for the randomly oriented fillers. For the evaluation of longitudinal elastic modulus, a small uniform strain is applied on one surface along the z-direction while only the opposite surface is fixed along its normal direction. The resulting stresses are calculated and by using the Hooke's law for unidirectional stress condition, the longitudinal elastic modulus is evaluated. In Fig. 9, the effects of unidirectional fillers on the longitudinal effective elastic modulus and longitudinal thermal conductivity of nanocomposite structures are illustrated. As the first finding in comparison with randomly oriented fillers, similar trends could be observed. In the case of unidirectional fillers, we could observe the priority of the cylindrical fillers in comparison with platelets fillers in the reinforcement of longitudinal composite properties, especially for higher fillers to matrix contrast in material properties. This is due to the fact the maximum dimension of cylindrical fillers at a constant volume is higher than that of the disc fillers. Accordingly, the cylindrical fillers are more engaged in load and heat flux transfer along the longitudinal direction.

The interphase reinforcement ratio for unidirectional fillers is shown in Fig. 9c and Fig. 9d.

63 Interestingly, we could observe that the interphase reinforcement ratio for unidirectional fillers is in the same order of that of randomly oriented fillers. In the same trend, the interphase has the least effect in the case of disc fillers. In this study, we also constructed 3-dimensional (3D) RVE for perfectly dispersed spherical fillers inside the epoxy matrix to compare the obtained results with those obtained by the proposed global statistical approach. It is worth mentioning that the bulk properties of silica nanoparticles and epoxy are used for introducing the materials properties in the strong-contrast and finite element modeling.

We should remind that the details of strong-contrast solution for elastic modulus are described in section 1.4 of this thesis.

Experimental

Materials and processing

The epoxy resin is diglycidyl ether of bisphenol A (DGEBA) provided from Dow chemicals (DER331). The hardener is a diethylenetriamine (DETA) manufactured by Aldrich company. Silica nanoparticles were purchased from Aldrich (medium particle size of 10-20 nm and specific surface area of 140-180 m2/g). Epoxy resin silica nanocomposite was prepared by mixing the epoxy resin and silica at 50°C under magnetic stirring. Quantity of silica nanoparticles is determined so that a composite with 2 wt. % in silica is achieved.

Results and discussions

In Fig. 11, the obtained results by global statistical approach are compared with experimental measurements and finite element results for perfect dispersion of particles. As it can be observed, the global statistical approach yields fine agreement with the experimental measurements from room temperature to around 100 ⁰C. Analysis of DMA curves show that glass transition occurs in a temperature range from 110 to 160°C. Furthermore, a shift of glass transition temperature (Tg) is observed; from 136°C to 140° for epoxy and epoxy nanocomposite, respectively. This increase is assumed to come from lower mobility of epoxy chains close to silica particles. Our modeling results, on the basis of continuum mechanics theories could do not take into account such phenomenon.

Moreover, it should be taken into account that the polymeric materials at high temperatures, close to glass transition temperature may no longer present linear elasticity response. This could be an additional reason why modeling results that are on the basis of linear elasticity start to deviate from experimental results at these temperature ranges. On the other hand, predictions by finite element model for perfect distribution of silica nanoparticles considerably underestimate the experimental results.

Fig. 11, Comparison between the estimations by statistical approach and finite element method with perfect dispersion of particles (random finite element) with experimental measurements for 2 wt. % silica nanopartciles inside epoxy.

We remind that in the two-phase composite materials, the stress and load is transferred between matrix and inclusions through their contacting surfaces. That is why the inclusions surface to volume ratio plays a critical role on the mechanical reinforcement of matrix. Since spherical inclusions have minimum surface to volume ratio compared to other geometries, we expect the minimum reinforcement for this case. In this way, in terms of mechanical reinforcement of twophase nanocomposite materials, having perfect dispersion of spherical fillers do not lead to a noticeable reinforcement. As it can be observed from Fig. 10b, the agglomerations of nano-silica particles result in the deviation of the effective geometries of fillers from spherical shape. These deviations from spherical shape mean higher surface to volume ratio for filling particles which accordingly means higher mechanical reinforcements. This effect could also explain the underestimation of the experimental results by finite element method for perfectly dispersed particles inside epoxy. We note that it has also been shown that the thermal conductivity of nanocomposite materials could be significantly enhanced by the aggregation of spherical nanoparticles [27].

Modeling of expanded graphite/PLA nanocomposites, a case study

In this study, we developed experimental as well as theoretically based hierarchical multiscale procedures for the evaluation of effective elastic modulus and thermal conductivity of expanded graphite (EG)/poly-lactide (PLA) nanocomposites. The incorporation of EG fillers into PLA was carried out by a twin-screw micro-extruder. The dispersion/delamination of EG in PLA was studied using Raman spectroscopy, SEM and TEM. In the multiscale modeling, the thermal conductivity constants and stiffness tensor of EG were first acquired by the means of molecular dynamics (MD) simulations. Using the fillers' properties obtained by the MD, we developed finite elements (FE) models to evaluate the effective thermal conductivity and elastic modulus of EG/PLA nanocomposites. The details of the experimenatl processing and characterization techniques could be found in Ref. [28 and 29].

Modeling

Since the expanded graophite is the multi-layer form of graphene, the molecular dynamics modeling of expanded graphite and the obtained results are presented in section 3.5 of this thesis. We should note that the graphite has the planar structure. The two-phase finite element modeling of isotropic platelet fillers is already explained in section 1.2 of this thesis. However, the graphite has the transversley isotropic nature, which is also shown by the molecular dynamics calculations (see section 3.5). For the modeling of a randomly distributed platelets with transversely isotropic properties, there are two available approaches. The first approach is to define a local coordinate system for each platelet and then assign the transversely isotropic properties for each platelet along

Results and discussions

The comparison between experimental and multiscale predictions of elastic modulus of EG/PLA nanocomposites is shown in Fig. 13. The experimental results are reported for two different EG concentrations of 3 wt. % and 6.75 wt. %, which are equivalent with 1.7 vol. % and 4 vol.%, respectively. As it can be observed, in the case of 3 wt.% concentration of EG at low temperatures (below glass transition), finite element results for the elastic modulus with effective aspect ratio of 120 for EG platelets match fairly well the experimental results. For PLA with 6.75 wt.% concentration of EG and at temperatures below the glass transition, the finite element results with effective aspect ratio of 60 fall close to the experimental elastic modulus. We should remind that for the platelet inclusions, the modeling results are strongly dependent on the effective aspect ratio of the inclusions inside the polymer (see Fig. 3 and Fig. 5 in this chapter). In a fabricated random nanocomposite material, the inclusions shape would not be similar. Therefore, considering similar geometries for all of the EG platelets inside the polymer matrix in our finite element modeling could not be considered as a realistic construction of nanocomposite sample. However, the comparison with finite element results could be used to verify the efficiency of EG platelets distributions and delamination inside the PLA. The agreement of finite element results with aspect ratio of 120 with experimental results could be used as a sign of well distributed and delaminated EG platelets inside PLA. On the other hand, by increasing the volume concentration of EG platelets, the finite element results with effective aspect ratio of 60 match the experiments, accordingly, this could imply the decline of the quality of EG distributions and less delamination of EG platelets inside PLA. This is in accord with the experimental observations [29]. It is worth mentioning that in the finite element modeling presented here we did not include the interphase effect as a third phase. This two-phase finite element modeling could be convincing for herein fabricated EG/PLA nanocomposite, since the TEM images and Raman spectroscopy results [28,29] did not reveal any significant sign for the formation of interphase between EG platelets and PLA matrix. Moreover, as already discussed in section 2.2, the interphase has the least effect in the case of platelet fillers and neglecting the interphase effect is not going have significant effect.

The strong mechanical interaction between carbon atoms of EG and polymer atoms could satisfy our perfect bonding assumption in finite element modeling for elastic modulus. As it is illustrated on Fig. 13, the finite element results are close to the experimental measurements from room temperature to 50°C. At this point, the finite element results deviate from experimental ones. We note that the basic assumption that we made in our finite element modeling is that the nanocomposite structure is a perfectly elastic solid. In the case of polymer materials this assumption is no longer valid if the specimen temperature approaches or gets higher than the glass transition temperature. It is important to underline a shift of the glass transition temperature toward lower values (4-5°C) in the theoretical curve compared to the experimental one. This discrepancy could be explained by the fact that in the theoretical model did not take into account the impact of EG on the polymer chain mobility responsible of the increase of the glass transition temperature of PLA.

Fig. 13, Comparison of the elastic modulus predicted by the modeling with experimental measurements for two EG concentrations of 3 wt.% and 6.75 wt.% in PLA nanocomposites.

Despite strong bonding between EG platelets and the polymer material, perfect heat flux transfer could not be physically acceptable, which is due to the existence of thermal boundary resistance between the fillers and the matrix [23]. Therefore, we could expect the overestimation of thermal conductivity of EG/PLA nanocomposite by the finite element method as shown by our predicted results in Fig. 14. We should remark that the experimental procedure for the evaluation of thermal boundary resistance at the interface between inclusion and matrix at nanoscale is considerably difficult and complex. For these cases, the virtual computer simulations using the molecular dynamics method seems promising. Such a molecular dynamics study is carried out for the graphene/epoxy nanocomposite which is presented in section 4.4 of this thesis. Fig. 9: Thermal conductivity of the nanocomposite predicted by the modeling and experimental measurements for two EG concentrations of 3 wt.% and 6.75 wt.% in PLA.

Concluding remarks

In this chapter, we compared three different groups of modeling techniques for the evaluation of thermal conductivity and elastic modulus of random two-phase composite materials. These three methods are: Mori-Tanaka (as a mean field method), finite element (numerical method) and strong contrast (statistical continuum approach). Moreover, we studied the interphase effect on the thermal conductivity and elastic response of nanocomposites. We also developed a new global statistical approach for the evaluation of nanocomposite properties. We finally used finite element approach in combination with molecular dynamics for studying the thermal and mechanical properties of expanded graphite/PLA nanocomposites. The obtained results could be summarized as follows:

 A new modeling procedure is proposed and discussed for the evaluation of effective thermal conductivity of composite materials by finite element method.

 The dependency of the finite element results on the RVE size increases by increasing the inclusions aspect ratio and contrast in materials properties.

 For platelet inclusions and for the case where the reinforcement in composite elastic and conductivity are lower than approximately 50 and 80 percent, respectively, the Mori-Tanaka results are close to those of finite element simulations.

 For cylindrical inclusions with low aspect ratios, the finite element and Mori-Tanaka results are close.

 The finite element method suggests a critical aspect ratio for both cylindrical and platelet fillers after which further increasing of aspect ratio does not increase the effective properties.

This critical aspect ratio increases by increasing the inclusion to matrix contrast in material properties and decreases by increasing the inclusions volume concentration. The critical aspect ratio for platelet fillers is higher than that for the cylindrical fillers  In the randomly oriented filled nanocomposites, the finite element results suggest that by employing the platelet fillers, one could achieve higher reinforcement in thermal conductivity and elastic properties of the composite in comparison with cylindrical fillers.

 In the nanocomposites filled with unidirectional particles, the finite element results suggest that by employing the cylinderical fillers, the reinforcement in thermal conductivity and elastic properties of the composite is much higher in comparison with platelet fillers. This is due to the fact that the cylindrical fillers are more engaged in load and heat flux transfer along the longitudinal direction.

 The finite element and Mori-Tanaka results show that the effect of inclusion to matrix contrast in materials properties is getting more pronounced as the inclusion geometry deviate more from spherical geometry.

 The estimations by strong contrast method are close to finite element and Mori-Tanaka for spherical inclusions. However, this method could not accurately describe the influence of inclusion geometry for the considered dilute systems.

 The two-point correlation functions cannot accurately capture the effect of high aspect ratio of the fillers.

 For randomly oriented and unidirectional filled composite structures, interphase presents the maximum effect for spherical inclusions. We could observe sharp decrease in the interphase reinforcement ratio as the fillers geometry slightly deviate from spherical shape.

 For the nanocomposite structures filled with fillers with high effective aspect ratios (long cylinders and thin discs), neglecting the interphase effect will present insignificant effect on the predicted estimations. The interphase present minimum effect for the disc fillers.

 For the spherical fillers, the effective reinforcement of nanocomposite properties is strongly dependent on the interphase to matrix contrast in properties and also the interphase thickness.

 A global statistical approach is developed for characterization of elastic properties of nanocomposite for taking into account the effects of nanoparticles agglomeration within the matrix.

 In the experimentally fabricated nanocomposites filled with spherical particles, 3D finite element method for perfectly distributed nanoparticles inside matrix may underestimate the experimental results. For these cases, the real microstructure should be used for the estimation of properties.

 Our study proposes the combination of molecular dynamics and 3D finite element modeling as a characterization technique for the evaluation and prediction of thermal conductivity and elastic modulus of nanocomposite materials.

Introduction

Graphene [1], the two-dimensional structure of honeycomb carbon atoms is an attractive material in fundamental physics and high technology industries due to its extremely high crystal and electronic quality [2][3][4]. Analogous to other members of carbon family, such as carbon nanotubes (CNTs) and graphite, graphene presents exceptionally high thermal conductivity [5], mechanical [6] and electrical properties [7]. Recent experimental studies [8,5] have shown that the single-layer graphene (SLG) can present extremely high thermal conductivity of around 4100±500 W/m-K which outperforms all other known materials. This observation proposes the graphene for micro/nanoscale electronic device heat transfer blocks in response to thermal management concerns in the electronics industry.

Using atomic force microscope nanoindentation experiments, Lee at al. [6] reported the elastic modulus of 1TPa and breaking strength of 130±10 GPa for defect-free graphene sheets.

Combination of high thermal conductivity, elastic properties and surface to volume ratio, also propose the graphene as a unique nanofiller for reinforcement of thermal, mechanical and barrier properties of polymeric materials [9].

In this work, we used classical molecular dynamics (MD) simulations as an alternative of experimental studies [10,11] for the investigation of the effective parameters on the mechanical response and thermal conductivity of graphene. It should be noted that the properties obtained by the MD modeling are strongly dependent on the selection of appropriate interatomic potentials, loading conditions and boundary conditions. In the current study we used the original Tersoff [12,13] and optimized Tersoff [14] potentials for introducing the carbon atoms bonding interactions.

The validity of the presented MD modelings are verified by observing that the acquired Young's modulus, tensile strength and thermal conductivity of free of defect graphene structures, predicted by the optimized Tersoff [14] potential are in close agreement with the latest experimental results.

The Young's modulus and thermal conductivity of bulk graphite is also found to be in agreement with the experimental results reported by Blakslee et al. [15] and Balandin et al. [8], respectively.

The experimentally fabricated graphene sheets are not perfect. Defects, which influence the physical and mechanical properties of a material, also exist in graphene. The defects in the graphene are usually in the form of Stone-Wales (SW), vacancies (point and bivacancy) and substitution atoms [16]. Accordingly, we studied the effects of point vacancy, bivacancy and Stone-Wales defects as the three major types of defects that have been experimentally observed in graphene on the tensile and thermal conductivity response of graphene. It is also worth noting that the chemical doping of materials with foreign atoms is currently considered as an effective method to intrinsically modify materials' properties [17,18]. For the chemical doping of carbon-based materials, boron and nitrogen are excellent elements because of their close atomic size and strong valence bonds with carbon atoms [19]. Experimental and theoretical studies on graphene have shown the possibility of making p-type and n-type semiconducting graphene by substituting carbon atoms with boron or nitrogen atoms [20][21][22]. It was shown that boron and nitrogen dopant atoms can modify the electronic band structure of graphene and open up an energy gap between the valence and conduction bands [23]. The doping of graphene by boron and nitrogen atoms could be also used in various applications such as nanosensors [21,24], spin filter devices [24], fuel cells [25] and Lithium ion batteries [26]. In this way, deep-seated understanding of effect of boron and nitrogen doping on graphene properties is a crucial issue in their future applications. To this aim, we also studied the thermal conductivity and tensile response of the boron and nitrogen doped graphene.

Recent atomic-scale microscopy and spectroscopy studies have shown significant deviation of graphene from its planar structures in forms of ripples, kinks, and folds [27,28] which may alter its properties [29,30]. Such observations have been also studied through different molecular dynamics studies [31][32][33]. Using the flexibility nature of graphene, we studied the curvature effect on the thermal conduction properties of graphene.

Molecular Dynamics Modeling

For performing the MD simulations, LAMMPS (Large-Scale Atomic/Molecular Massively Parallel

Simulator) [34] is used. The bonding interactions between carbon-carbon atoms are modeled using original Tersoff [12,13] and optimized Tersoff [14] potentials. The optimized Tersoff potential could predict the phonon dispersion curves of graphite in better agreement with experiments [14] than the original Tersoff potential. The detailed information about the Tersoff potential is already presented in section 1.1.3. The Tersoff potential parameters for nitrogen and boron atoms where extracted from [35], thus, the carbon-nitrogen and carbon-boron bonding interaction parameters are calculated using equation (1e) in [12] (Eq. 20 in chapter 1). In the modeling of few-layer graphene (FLG), the nonbonding van der Waals interaction between individual carbon atoms at different atomic planes are commonly introduced by Lennard-Jones (LJ) potentials, expressed by: where r refers to the distance between carbon atoms at different atomic layers. ε is the depth of the potential well and σ is the equilibrium distance at which the inter-particle potential is zero. In order to satisfy the bulk density of graphite (2.2 gr/cm 3 ), the equilibrium distance (σ) should be equal to 3.4 Å.

Fig. 1, MD model of uniaxial tension test of graphene with chirality angle of ≈10º in respect to armchair direction. The red and blue atoms are representative of carbon and substituted nitrogen or boron atoms respectively. The loading condition was applies by the extending the periodic simulation box in the X direction. The molecular dynamics simulation model for uniaxial tensile test of graphene with the chirality angle of ≈10º is illustrated in Fig. 1. The chirality angles equal to 0º and 30º are called armchair and zigzag, respectively. As shown in Fig. 1, the red and blue atoms belong to carbon and nitrogen atoms, respectively. The concentrations of nitrogen or boron atoms are defined by the ratio of vacancies or dopant atoms to the total number of available atoms in a free of defect single layer graphene (SLG). In order to prevent the formation of two nitrogen/boron atoms bonds (as in the experiments [24]), the random distributions of nitrogen/boron atoms substitutions were not accepted for two adjacent atoms. In our modeling, the FLG is constructed by putting SLG sheets up together by a spacing distance of 3.4 Å. In all the modeling in the evaluation of tensile properties, the periodic boundary conditions (PBC) was applied along the loading direction (X direction), which is the best way to remove the finite length effect and is the key point for satisfaction of uniaxial stress conditions. Applying the PBC in the Y direction is representative of an infinite graphene sheet. However, free boundary condition in the Y direction is suited for modeling of graphene nanoribbons. The properties of bulk graphite were also evaluated by considering 6 layers of graphene sheets while the PBC was also applied in the graphene sheets' normal direction (Z direction). Due to high computational cost of MD simulations, the tensile tests are usually performed at high strain rates, in the order of 10 9 s -1 . This limitation is due to the small interatomic distances and high velocities of atoms at room temperature. We found that the time increment of 1fs (10 -15 of a second) is computationally efficient and adequately small to accurately describe the deformations of graphene at high loading strain rates.

In Fig. 2, different types of defects that are analyzed in the current study are illustrated. We studied effects of point vacancy, bivacancy and Stone-Wales defects on the graphene effective properties.

Point vacancy (see Fig. 2a) is the simplest defect in any material which is the missing of one atom in the lattice. This defect leads in formation of a five-membered and a nine-membered ring, as highlighted in Fig. 2a. The bivacancy defect (Fig. 2b) can be created either by the coalescence of two point vacancies or by removing two bonding atoms. As highlighted in Fig. 2b, in a bivacancy defect, no dangling bond exists and two pentagons and one octagon are formed instead of four hexagons in ideal graphene. We should note that since the cut-off distance of Tersoff potential for carbon atoms is 2.1 Å, by simply removing carbon atoms, point vacancy and bivacancy defects would not form accurately. Therefore, in these cases, two neighboring atoms of a removed atom should be approached together in order to form covalent bonds (their distance should be below 1.8 Å). The Stone-Wales defect does not involve any removed or added atoms. As highlighted in Fig. 2c, in Stone-Wales defect, four hexagons are transformed into two pentagons and two heptagons by rotating one of CC bonds by 90°. We found that in order to have stabilized defective graphene structures, system should be initially energy-minimized using conjugate gradients method. Similar to experimentally observed defects in graphene, in our modeling defects were randomly distributed along the graphene sheet, observing periodicity of graphene structure. The concentrations of defects are introduced by the ratio of removed (in point vacancy and bivacancy defects) or rotated (in Stone-Wales defects) atoms to the total number of available atoms in a defect-free graphene structure. In all cases, before applying the loading condition, the graphene structure was left to relax to zero pressure along the periodic directions at a temperature of 300 K using constant pressuretemperature (i.e. NPT ensemble) by means of Nosé-Hoover barostat and thermostat method. The loading condition was applied by extending the periodic simulation box size in the X direction under constant engineering strain rates. In order to reduce the intensity of high strain rate effects on the obtained properties, the raising of temperature as a result of high loading strain rates must be damped by performing the constant temperature simulations during the loading step. In the case of bulk graphene sheets, in order to satisfy the uniform uniaxial stress conditions, the periodic simulation box in the Y direction was altered to reach zero stress in this direction, using NPT Nosé-Hoover method. Virial stresses in the loading direction (uniaxial stress) were calculated at each strain level to obtain the stress-strain response. The engineering strain (e e ) at each time step was calculated by multiplying the total time of loading simulation by the applied strain rate, thus, the true strain (e) is obtained by ln(1+ e e ). Finally, the Young's modulus, Y, of the structure can be calculated by: where Us and A are the deformation energy and initial surface of the structure respectively, and t is the thickness of graphene which is equal to 3.4 Å. The strain energy at each strain level is obtained by calculating the potential energy of the strained system minus the potential energy of the relaxed system. By fitting a parabola, 0.5αe 2 + βe + c, to the strain energy at the low strain levels, the Young's modulus can be calculated by .

The MD model developed for the evaluation of thermal conductivity of graphene structures is shown in Fig. 3. In Fig. 3 (a), the MD architecture of perfect and/or nitrogen/boron doped graphene is presented. We fix three rows of carbon atoms at the two ends of the graphene to prevent atoms from sublimating. The periodic boundary condition was applied in the Y direction (see Fig. 3a).

Thus, as mentioned before these simulations are representative of an infinite graphene sheet and not the graphene nanoribbons. In Fig. 3 (b,c) the MD model for the qualitative study of directional dependent thermal conductivity of curved graphene is illustrated. The curvature angle of 0 and 360 degree are respectively equivalent to planar graphene and carbon nanotube. The thermal conductivity is studied along the tangential (circumferential) as well as the axial (longitudinal) direction of the curved graphene. Since the thermal conductivity of graphene is dependent on chirality and length as well, the length (along X direction in Fig. 3b and Fig. 3c) and chirality direction are the same in all curved graphene. In both cases, we fixed 3 layers of carbon atoms at two ends along the X direction and the atoms on the two other edges are left to be free. In this study, thermal conductivity of graphene structures was evaluated using non-equilibrium molecular dynamics method. In this method, we divided the simulation box (excluding the fixed atoms at the two ends) along the longitudinal direction into 21 slabs. First, specimen was relaxed at room temperature (300 K) using constant temperature simulations by the Nosé-Hoover thermostat method (NVT). In the next step, a temperature difference was applied between the first and 21st slabs. In this study, the temperatures at the first (hot reservoir) and 21st (cold reservoir) slabs were set to 310 and 290 K, respectively, using the NVT method, while the remaining slabs are under constant energy (NVE) simulations. In order to keep the temperature differences at the two ends, an amount of energy is added to the hot reservoir and at the same time another amount of energy is removed from the cold reservoir by the NVT method. In Fig. 4a, the adding and removing energy to the hot and cold reservoirs, in the case of 2% concentration of Stone-Wales defects in graphene, are illustrated. As it is illustrated, the amount of energy added to the hot reservoir is equal to removed energy from the cold reservoir which accordingly means that the total energy of the system is accurately conserved. Moreover, the slopes of energy curves in both cases are linear which means that the energy is added or removed with a constant rate. These observations imply that the energy of the simulation box remains constant and a constant heat flux is applied throughout the sample. This heat flux is a constant quantity of kinetic energy transferring from the atoms in the hot reservoir to those in the cold reservoir at each simulation time step. The non-equilibrium steady state heat transfer can be achieved after 500 ps of the exchanging process, when a temperature gradient is established along the longitudinal direction. The applied heat flux by the NVT method is calculated based on the slope of energy curves (Fig. 4a). The thermal conductivity, k, of the graphene could be calculated using the one-dimensional form of the Fourier law:

Here, is the established temperature gradient and Ac is the cross sectional area of the graphene structure. In this work, temperature at each slab is computed as follows:

where T i (slab) is temperature of i th slab, N is number of atoms in this slab, k B is the Boltzmann's constant, m j and p j are atomic mass and momentum of atom j, respectively. Simulations were performed for 3 ns and the averaged temperatures at each slab were computed. As it is shown in Fig. 4b, by neglecting the jumps at two ends, one could find a linear temperature gradient in the middle of the specimen. By computing the established , thermal conductivity could be calculated using Eq. 3.

Tensile response of graphene

As the loading strain rates in MD simulations are commonly in millions of magnitude to those of available experiments, the simulations were performed at different strain rates to verify the intensity of the strain rate effects on the obtained mechanical properties. In Fig. 5, strain energy for armchair SLG is plotted as function of strain for different loading strain rates acquired by optimized and original Tersoff potentials. At low strain levels (up to ≈ 0.05), we could fit a parabola function with high level of accuracy and negligible variance (less than 1%). Assuming this part as the linear elasticity part for the stress-strain response (see later figures), the Young's modulus was obtained using eq. 2. Based on Fig. 5 observations, the Young's modulus of the armchair SLG are independent of the strain rate values. We also verified the Young's modulus consistency for other chirality directions. Moreover, the obtained Young's modulus of SLG, predicted by the optimized and original Tersoff potentials in armchair direction are 0.985 and 1.05 TPa, respectively, which both are in good agreement with experimental Young's modulus of 1 TPa for defect-free SLG reported by Lee et al. [6]. While the obtained Young's modulus is independent of the loading strain rate, the predicted ultimate tensile strength (UTS) which corresponds to the maximum strain energy decreases by decreasing the strain rate. In this way, the UTS should be reported at the corresponding loading strain rate. By calculating the Virial stresses, we found that at the strain rates of 1×10 8 s -1 and In Fig. 6, we compared the deformation behavior of SLG, predicted by the optimized Tersoff (Fig. 6 a-d) with chirality angle of ≈10º with that predicted by the original Tersoff (Fig. 6 e-h) along the armchair direction. It was observed than in both cases and up to relatively high strain levels, the specimen extends uniformly along the loading direction and the structure remains undamaged. The uniform extension continues up to the ultimate tensile strength point, at which voids are initiated in the SLG. For both versions of the Tersoff potentials, the nucleation of these voids takes place during a short time. This indicates that the specimen has reached its maximum load bearing. In the case of optimized Tersoff potential, the voids are initiated by debonding between two adjacent carbon atoms. The growth of these voids seems to occur approximately along the armchair direction with triangle shape of carbon atoms at two ends (ends of ellipsoidal-like void shape). As it can be observed in Fig. 6 (b andc), several voids usually develop in parallel. Since the growth of these voids requires a series of debonding of carbon atoms, we could expect a series of stress raises and falls (see Fig. 7). By the development of several voids along the specimen, damaged SLG presents more flexibility than the original structure (undamaged). As a result of more flexibility of 90 SLG, it goes under considerable deflection to reduce stress concentration at the ellipsoidal void tips (see Fig. 6d). At the final steps, the two part of the specimen are connected together only by several chains of carbon atoms (see Fig. 56). This is due to the coalescence of voids which leads to the final rupture when the last connecting chains fail. In this case, considerable amount of disorder and deflection are persistent in the final structure of SLG. We can therefore conclude that the optimized Tersoff potential predicts a ductile failure behavior for SLG.

On the hand, in the case of original Tersoff potential, the nucleation of voids (Fig. 6f) and subsequent failure (Fig. 5h) occur at considerably close strain levels. We also observe that the nucleation, growth and coalescence are highly localized (see Fig. 6f-h). These observations suggest that the deformation behavior of SLG presents a high level of void sensitivity which accordingly means brittle failure behavior. Furthermore, as the rupture takes place, the undamaged parts of the graphene sheet (Fig. 6h) tend to retrieve the original structure and contract rapidly.

It should be noted that the first-principles density-functional theory study, by Topsakal and Ciraci [36], predicted the deformation behavior of graphene nanoribbons closer to that we observed by optimized Tersoff rather than original Tersoff. Observing the high accuracy of optimized Tersoff potential in the prediction of mechanical properties and deformation response of graphene, the remaining results are reported on the basis of optimized Tersoff potential.

In Fig. 7, the acquired stress-strain response of SLG with different loading chirality angles, ranging from 0⁰ (armchair) to 30⁰ (zigzag), at the strain rate of 1×10 8 s -1 is shown. At the first glance, the stress-strain relations are almost coinciding which suggests a limited chirality dependency. However, as an overall trend, both Young's modulus and UTS of SLG increase slightly by increasing the chirality angle from 0 to 30º and the Young's modulus of zigzag SLG is around 10 GPa higher than that of the armchair. In Fig. 8, the effect of graphene nanoribbon width on the stress-strain response is shown along with the stress-strain response of bulk graphene. As it can be observed, the Young's modulus of graphene nanoribbon decreases by decreasing the graphene width. This observation suggests that the edge atoms have a softening effect on the Young's modulus of graphene which turns to be more significant when the nanoribbons widths are lower than 3 nm. In all cases, the UTS of graphene nanoribbons are lower than the bulk SLG. However no particular trend as function of varying width can be derived.

Fig. 8, Effect of nanoribbon width on stress-strain response of zigzag graphene compared with graphene sheet, the strain rate is equal to 1×10 8 s -1 .

For the fabrication of graphene, it is difficult to get only to SLG sheets, and we usually obtain a combination of single and few-layer graphene structures. It should be noted that the latest experimental and theoretical study [8] suggest significant and continues decline of the in-plane thermal conductivity of graphene as one goes from single layer to 4 layers graphene structures. Thus, in our work, we examine the thickness effects on the mechanical properties of graphene by considering few-layer graphene (FLG) sheets. There are wide variety of parameterization available in the literature for the depth of well (ε), in Lennard-Jones potentials for nonbonding interactions for carbon atoms. To best of our knowledge, the most popular ones are: ε=2.4 m-eV [37] followed by ε=12 m-eV [38]. Using these two LJ potentials, we compared the effect of the number of atomic planes on the Young's modulus of FLG in Fig. 9. As the first finding, for both LJ potentials, the obtained in-plane Young's moduli of bulk graphite are of close values and are in a fine agreement with the experimental value of 1.02±0.03 TPa [15]. We also found that the ultimate tensile strength and its corresponding strain of FLGs and bulk graphite vary within a maximum of 5% from those of 93 SLG. However, the FLG's could undergo higher strain levels before they completely lose their load bearing ability (failure). For both LJ potentials, we could say that using the optimized Tersoff potential, the Young's modulus of FLG is acceptably independent of the number of atomic planes.

Such a finding can be interesting for multiscale modeling of graphene nanocomposites. Thus, the elastic properties of few-layer graphene can be accurately approximated by using those of bulk graphite.

Fig. 9, The effect of number of atomic planes on the Young's modulus of few-layer graphene along the armchair direction, using two different parameterization for interlayer Lennard-Jones potentials depth of well, ε=2.4 m-eV [36] and ε=12 m-eV [37].

In Fig. 10, we studied the effects of Nitrogen doping and point vacancies concentrations on the Young's modulus (Y) and ultimate tensile strength of SLG. As it can be observed, by increasing the concentration of nitrogen atoms up to 6%, the Young's modulus of SLG is not changing. However, the tensile strength decreases by more than 35% only by having 2% concentration of nitrogen atoms. On the other hand, we observe that the Young's modulus of SLG is very sensitive to point vacancies. It decreases by increasing vacancy concentration particularly when the concentration is more than 2%. Furthermore, the tensile strength decreases considerably by around 25% percent only by 0.5% point vacancies. Further increasing of concentration of vacancies up to 6% decreases the tensile strength of SLG in an approximate linear pattern to around 60% of that of the defect-free SLG structure. In Fig. 11, we plot the acquired engineering stress-strain response of defect-free graphene and 4% boron doped graphene. We have performed the uniaxial tension tests for six various concentrations of boron atoms in graphene. We observed that by increasing the boron atoms concentrations, the tensile strength and elastic modulus of graphene decrease slightly. We found that in the case of 4% boron doped graphene, the elastic modulus and tensile strength declined by approximately 5% and 8%, respectively. Our results suggest that boron atoms substitutions present negligible effects on the mechanical response of graphene. 95 Fig. 11, Engineering stress-strain response of defect-free and 4% boron doped graphene.

The effect of nitrogen atoms concentration on the tensile behavior of SLG at different strain levels is illustrated in Fig. 12. In a similar trend to that of perfect SLG, the specimen uniformly extends along the loading direction and voids initiation is the sign that the specimen has reached its ultimate tensile strength. Unlike the free of defect SLG, where the initial voids are localized and the remaining parts of the SLG are undamaged (see Fig. 6b), in nitrogen doped SLG, the initial voids seem to be distributed along the entire sheet. These initial voids are the consequences of debonding between two adjacent carbon and nitrogen atoms along the loading direction. In nitrogen doped graphene, the first debondings between carbon-nitrogen atoms promote debonding between the closest pair atoms (see Fig. 12b). This extends and forms a weak region along the specimen width (see Fig. 12c). At this point, the initiation of voids in the remaining parts of the specimen is stopped.

Subsequently, void growth occurs in the formed weak region of the specimen (see Fig. 12c). As shown in Fig. 12d, the coalescence of the voids results in the rupture of the specimen. Similar to ideal graphene sheets, the formation of carbon chains helps the specimen to keep its load bearing ability at higher strain levels after the tensile strength. However, the existence of substituted nitrogen atoms in the carbon chains reduces the stretching limit of these chains. We can therefore conclude that the nitrogen doping of graphene decreases its failure ductility. In Fig. 13, the deformation process of 4% doped graphene at various strain levels is shown. As it is clear from this figure, the specimen extends uniformly and remains undamaged up to the tensile strength. The tensile strength is a point at which the first debonding occurs between two adjacent carbon and boron atoms (Fig. 13b). The initial void is formed in the place where we have the maximum density of the doped boron atoms. Shortly after the initiation of the first void, more debonding happens between the closest pair atoms (Fig. 13c). As shown in Fig. 13d, the coalescence of the voids results in the final rupture of the specimen. As it is clear, the nucleation of initial voids (Fig. 13b) and subsequent failure (Fig. 13d) occur at considerably close strain levels. Accordingly, deformation of the boron doped graphene presents void sensitivity which means brittle failure behavior. In this way, the boron atoms substitutions in graphene result in changing the failure mechanism form ductile to brittle. Fig. 13, Tensile deformation process of single layer graphene with 4% concentration of boron atoms at various stages of loading.

In Fig. 14, effects of point vacancy, bivacancy and Stone-Wales defects concentrations on the graphene effective elastic modulus, tensile strength and strain at failure are illustrated. In all cases, by increasing the defects concentration, elastic modulus decreases gradually. Since in the Stone-Wales defects the number of atoms is the same as that of the defect-free graphene, we could observe the minimum effect in this case. Decreasing of elastic modulus is more significant in the case of bivacancy and point vacancy defects. In these cases, elastic modulus decreases by around 20% for 2% concentration of defects which is almost two times higher than the same concentration of Stone-Wales defects. In the case of bivacancy defects, we could observe that at low concentrations (up to 0.5%) elastic modulus is close to that of the Stone-Wales defects. In this case, by increasing the defects concentration the elastic modulus gets closer to that of the point vacancy defects. On the other hand, we could observe that effect of Stone-Wales and bivacancy defects on the tensile strength and strain at failure of graphene are considerably close to each other. In all studied cases, we could observe a sharp decline by around 20% in the tensile strength and strain at failure of graphene for only 0.25% concentration of defects. In the case of Stone-Wales and bivacancy defects, further increasing of defects concentration does not present considerable decreasing effect on the tensile strength and strain at failure. On the other hand, in the case of point vacancy defects, by increasing the defects concentration an overall reduction trend could be observed for both of tensile strength and strain at failure. In comparison with boron and nitrogen doped graphene, we could observe considerably higher destructive effects on the graphene mechanical response by point vacancy, bivacancy or Stone-Wales defects.

In Fig. 15, deformation process of defect-free graphene is compared with that of the Stone-Wales and bivacancy defects at various strain levels. As it is shown, in the case of defect-free graphene, specimen extends uniformly and remains undamaged up to the tensile strength point (Fig. 15b). The tensile strength is a point at which the first deboning occurs between two adjacent carbon atoms.

Shortly after, we could observe the growth of the initial void and extension of this defective zone along the specimen (see Fig. 15c). At this step (Fig. 15c), initial voids are localized and the remaining part of the specimen remains undamaged. The main characteristic of these initially formed voids is the formation of mono-atomic carbon chains. As it is shown in Fig. 15c, mono-atomic carbon chains are formed along the armchair direction and initial debondings are also occurring along armchair direction. By increasing the strain levels, we could observe the initiations of more voids and also elongation of mono-atomic carbon chains that tend to keep the graphene sheet connected.

The final rupture occurs when the last connecting carbon chain fails (Fig. 15d). As shown in Fig. 15d, in the ruptured graphene sheets, edges are mainly along the zigzag direction. Moreover, in this case the ruptured graphene sheet is considerably damaged (Fig. 15d). In the case of defective graphene sheets, in a similar trend of that of defect-free graphene, we could observe the formation of initial voids, extension of these voids and final rupture of the specimen when last connecting mono-atomic carbon chain fails. In the case of bivacancy defects, at strain levels lower than tensile strength limit, we could observe the formation of ten-membered carbon rings, in places where two bivacancy defects are connected together by only one CC bond (see Fig. 15f highlighted rings).

This way, the connecting bond of two bivacancy defects fails and forms a ten-membered ring connecting two bivacancy octagons. However, similar to defect-free graphene sheet, the tensile strength happens at a point in which debonding of a CC bond occurs along the armchair direction (Fig. 15f). This debonding promotes more debondings of carbon bonds which result in the formation of mono-atomic carbon chains (Fig. 15g) and leading to specimen rupture (Fig. 15h).

Fig. 15, Deformation process of single-layer graphene sheets at various engineering strain levels.

Molecular dynamics study of thermal and mechanical response of graphene 101 Unlike defect-free graphene sheet, in the case of bivacancy defects damage of graphene is more localized and the rest of graphene sheet remains almost undamaged. In the case of Stone-Wales defects, at strain levels below the tensile strength point, we could observe formation of several parallel connected ten-membered rings as a result of debonding of CC bonds between pentagons and heptagons (see Fig. 15l). As it is shown in Fig. 6l, in some cases, one of the ten-membered rings extends and form fourteen-membered ring which is in parallel with previously formed tenmembered ring. These formed defects are the starting point for the extension of more parallel voids at tensile strength point (see Fig. 15m). In a similar trend, the extension of voids and failure of formed mono-atomic carbon chains lead to the specimen rupture (Fig. 15n).

Fig. 16, Instantaneous angles distributions of graphene sheets at different stages of deformation.

In Fig. 16, we compared the instantaneous angles distribution of graphene sheets at various deformation stages. In a honeycomb planar structure, angle between two linked bonds is 120°.

Therefore, in a perfectly planar graphene structure, angles distribution would present only one peak at 120°. This way, broadening of angles distribution peak means waviness of the structure and deviation from ideally planar structure. From the results illustrated in Fig. 16, we could conclude that in defective graphene sheets, the structure deviate more form planar structure which is concluded on the basis of more broadening of angles distribution peaks in defective graphene sheets in comparison with defect-free graphene sheet. Results of Fig. 16 at equilibrium stage suggest that point vacancy defects result in a maximum waviness of graphene structure while the Stone-Wales defects present least effect for this case. By stretching the honeycomb structure, two angles increase and the other angle decreases. Consequently, as it is shown in Fig. 16, under uniaxial tensile strain, we could observe presence of two wider peaks for angles distribution. This way, the probability of angles (angles distribution) for the peak occurring at higher degree is almost two times of that occurring at lower degree. We note that relative degrees in which these peaks are occurring change as strain level changes. In all cases, after the occurrence of rupture in graphene sheets, angles distributions approximately match that at zero strain level. Such an observation suggests that the level of planarity of ruptured graphene sheets is acceptably the same as that of the undamaged graphene structure.

Thermal conductivity of graphene

In Fig. 17, the inverse of thermal conductivity of armchair SLGs as a function of the inverse of the graphene length is illustrated. Due to high computational cost of NEMD simulations, the maximum length of the SLGs in the current study was limited to 1000 nm. As it can be observed, by increasing the graphene length the thermal conductivity of increases. The length dependency of the obtained thermal conductivity of SLG can be eliminated by extrapolation of the obtained results for infinite length. For this, a simple approach is to define the effective phonon mean free path ( ) as

. Since the thermal conductivity is proportional to , the thermal conductivity of the infinite system can be obtained by extrapolating to 1/L=0, [38]. In this way, the size independent thermal conductivity of SLG is obtained to be around 3,000 ± 100 W/m-K. It should be noted that the obtained thermal conductivity of SLG by extrapolation of 6 data sets between 50 to 200 nm also falls within the herein reported range. This estimate is close, but underestimates, the experimental results reported [8]. It is worth mentioning that in a similar study in [39], by using the second generation REBO potentials [40], the thermal conductivity of SLG was predicted to be around 350 W/m-K [39]. Using the original version of the Tersoff potential [12], Wei et al. [41] obtained the thermal conductivities of around 850 W/m-K for SLG. This shows the considerable accuracy of the optimized Tersoff potential in analyzing the thermal conductivity of carbon structures.

103 Fig. 17, Measured and extrapolating curves of inverse of thermal conductivity as function of inverse of single-layer graphene length.

In Fig. 18, the effects of (a) nitrogen atoms and (b) boron atoms concentration on the thermal conductivity of SLG are studied at a finite length. The simulations were performed along the armchair and zigzag directions and the reported values are normalized based on the thermal conductivity of zigzag at this length. As a first finding, the thermal conductivity of zigzag SLG takes a higher value than that in the armchair direction. Such a finding is in agreement with previous MD studies [42][43] at room temperature. This was explained by the fact that in the zigzag direction there are more particles to conduct phonons than along the armchair direction [44]. Fig. 14a results suggest that the substitution of just 1% of nitrogen atoms in graphene results in a decline of the thermal conduction properties of graphene by more than 50% and in a vanishing superiority of zigzag heat conductivity. The situation is even worse in the case of boron doped graphene. We can observe that only 0.75% concentration of boron atoms in graphene leads to more than 60% reduction in thermal conductivity of graphene along both chirality directions. In the both cases of boron and nitrogen doped graphene we could observe an approximate saturation limit in the reduction of thermal conductivity at concentrations higher than 2%. A recent MD study [45] on the nitrogen doped CNTs has demonstrated that the point vacancies have more destructive effect on the CNTs thermal conductivity than the substituted nitrogen atoms. Taking into consideration the 104 similar atomic size and strong bonding interactions between carbon and nitrogen/boron atoms, we can conclude that the thermal conduction properties of graphene are excessively sensitive to defects and doping.

Fig. 18, Effect of (a) nitrogen and (b) boron atoms substitutions on the thermal conductivity of single-layer graphene along the armchair and zigzag chirality directions.

In Fig. 19, normalized thermal conductivity of graphene as a function point vacancy, bivacancy or Stone-Wales defects concentrations are illustrated. In Fig. 19, we also compared obtained thermal conductivity of defective graphene with that of nitrogen doped graphene reported in. As a first finding, we could observe similar exponential trends in the reduction of graphene thermal conductivity due to defects and dopant atoms concentrations. Interestingly, we found that the effects of studied defects on the thermal conductivity of graphene are considerably close to each other and a distinct difference exists in comparison with nitrogen doped graphene. Accordingly, thermal conduction in graphene is much more sensitive to studied defects in comparison with chemical doping. Reported results in Fig. 19 suggest that only 0.25% concentration of defects in graphene could result in significant reduction of the thermal conduction property of graphene by around 50%. In all studied cases here, we could observe the maximum and minimum destructive effects on graphene thermal conductivity for the case of point vacancy and Stone-Wales defects, respectively. In an analogous trend like nitrogen doped graphene, sharp reduction in graphene thermal conductivity occur up to 0.5% concentration of defects followed by a plateau trend for higher defects concentrations. We note that since the energy of the simulation cell does not conserve accurately, we could not report the effect of point vacancy defects at higher concentrations 105 than 1%. However, for this case the reported results for bivacancy defects could be used as close upper bound estimates.

Fig. 19, Normalized thermal conductivity of graphene as a function of defects concentrations.

In order to find the mechanism that results in the reduction of thermal conductivity of nitrogen doped SLG, we plotted a carbon atom longitudinal phonon density of states (DOS) in Fig. 15. The longitudinal DOS can be obtained by:

where ω is the frequency, v x is the velocity of atoms along the established heat flow . In order to report size independent DOS, the fixed atoms were removed and periodic boundary condition was also applied in the longitudinal (X direction in Fig. 2a) direction. The DOS was plotted for a carbon atom with no doping (Fig. 15a) and when it is bonded with a nitrogen atom (Fig. 20b). As it can be observed from Fig. 20, the nitrogen atom does not change the relative frequency of the maximum number of accessible phonons in the DOS. However, we can observe some broadening of the DOS peaks in the frequency range. These broadening of the peaks in DOS can be explained as an indication of phonon scattering. The scattering of phonons by higher concentration of nitrogen or boron atoms can explain the vanishing of chirality effect. The simulations were conducted for different nanoribbons widths. The widths of graphene nanoribbons were set to 20 nm, in which obtained thermal conductivity of graphene nanoribbons falls only within 5 percent difference with that of infinite width graphene. In Fig. 16, the effects of curvature angle on the longitudinal and circumferential heat transfer of a single layer graphene are illustrated. The results are normalized based on the thermal conductivity of planar graphene nanoribbons. We believe that any kind of deviation from the natural structure of carbon atoms, planar (graphene) and cylindrical (CNT) structure makes the phonon transportation phenomena more difficult. This will accordingly results in lower thermal conduction performance. As illustrated in Fig. 21, by increasing the curvature angle up to the 180 degree, the thermal conductivity in the circumferential direction drops slightly while this drop is more considerable in the longitudinal direction. Interestingly, the longitudinal conductivity increases for higher angles and the minimum thermal conductivity of curved graphene in this direction occurs for 180 degree. Since the thermal conductivity of 270 and 360 degrees curved graphene are close to those of 90 and 0 degree respectively, it can be concluded that an approximate symmetric relation exist. For further increase of the curvature angle from 180 degree, the circumferential heat transport in graphene declines considerably. For a 350 degrees curvature structures it drops to around 75% of the original graphene sheet. Our simulations suggest that the curvature effect has less destructive impact on the thermal conductivity of graphene in comparison with chemical doping and defects. Such findings 107 propose graphene as a flexible and reliable block for fast heat dissipation in micro/nanoscale electronic devices.

Fig. 21, Effects of longitudinal and circumferential curvature angle on the thermal conductivity of a graphene sheet at room temperature

Thermal conductivity and elastic constants of graphite

In Fig. 22, the equilibrated MD model for the evaluation of graphite properties is shown. The graphite was constructed by putting the graphene sheets up together with a spacing distance of 3.4

Å. In the MD simulations in this study, we applied the periodic boundary conditions in the planar directions (X and Y directions in Fig. 22). In the evaluation of bulk graphite properties, we also applied the periodic boundary conditions along the thickness direction (Z directions in Fig. 22). It is worth noting that in all of the simulations in this study, prior to the step for properties evaluation, the structure was allowed to relax to zero pressure at 300 K in periodic directions using constant pressure-temperature (i.e. NPT ensemble) simulations by Nosé-Hoover barostat and thermostat method. Then the final equilibrated structure was obtained by performing constant energy (i.e. NVE ensemble) simulations. The time increment of the simulations was set to 0.5 fs. The optimized 108

Tersoff potential [14] was used to introduce the bonding interactions between carbon atoms. The non-bonding interactions were modeled using the LJ potential proposed in Ref. [36].

We developed a six layer graphene structure (as illustrated in Fig. 17) for the evaluation of expandedgraphite elastic constants at room temperature. Calculation of elastic constants at finite temperature is always challenging. One way to do this is to measure the change in average stress tensor at constant temperature NVT simulations when the simulation box undergoes finite deformations [START_REF]LAMMPS user manual[END_REF].

By performing the MD simulations, the elastic constants of expanded graphite compared to different experimental studies are listed in Table 1. As it can be observed, the molecular dynamics results are in close agreement with different experiments. We also verified that the EG presents transversely isotropic behavior. This means that the elastic responses along the X and Y directions are the same. It should be mentioned that due to the uncertainties observed in the evaluation of C44, we preferred not to report any value.

Fig. 22, The constructed six-layers graphene structure for evaluation of thermal conductivity and stiffness tensors of expanded graphite by molecular dynamics. The equilibrium molecular dynamics (EMD) method is used for the evaluation of thermal conductivity of expanded graphite. The EMD method relies on relating the ensemble average of the heat current auto-correlation function (HCACF) to the thermal conductivity K, via the Green-Kubo expression:

where ɑ denotes the three Cartesian coordinates, K B is the Boltzmann's constant, V and T are the volume and temperature of the system respectively. The auto-correlation functions of the heat current (HCACF) can be calculated using the heat current as expressed by [START_REF]LAMMPS user manual[END_REF]:

here, e i and v i are the total energy and velocities of atom I, respectively, f ij and r ij are respectively the interatomic forces and position vector between atoms i and j. The in-plane thermal conductivity of EG is obtained by averaging the thermal conductivity along the X and Y directions. We found that the thermal conductivities along the X and Y directions are considerably close to each other, thus, the thermal conductivity of EG also presents transversely isotropic nature. It is worth noting that in order to minimize the statistical errors, the thermal conductivity is calculated based on the average HCACF, calculated from eight independent MD runs with different initial atomic velocities.

Fig. 23, depicts the in-plane and along the thickness normalized HCACF obtained from MD simulation for bulk graphite. It shows a decay behavior over the entire time range. The thermal conductivity is calculated by performing the integration of HCACF using Eq. 6. By direct integration of HCACF over the time range, the in-plane and out of plane thermal conductivity of graphite are obtained to be 1200±100 W/m-K and 2.2±0.2 W/m-K, respectively. The size independency of the 110 obtained results was verified by performing the simulations for different length, ranging between 3 to 8 nm and different number of layers from 6 to 8. We also performed the same study on 5 layers graphene sheets (no PBC along the Z direction) and we found that the in-plane thermal conductivity falls within the obtained range for bulk graphite. It is worth mentioning that the latest experimental results [8] suggest a continuous decline of in-plane thermal conductivity of graphene from singlelayer (4100±540 W/m-K) to four-layers graphene (1230±160 W/m-K) structures. Further increase of graphene number of atomic planes does not change the reported in-plane thermal conductivity value. The in-plane thermal conductivity of graphite has also been reported to be in the range between 800-2000 W/m-K [8]. As it can be concluded, using the optimized Tersoff potential the predicted thermal conductivity of graphite is in a considerable agreement with experimental values.

We should note that using the original Tersoff potential [12,13], the thermal conductivity of graphite is predicted to be 500 W/m-K [41] which is far from the experimental range.

Concluding remarks

We used molecular dynamics simulations in order to probe the thermal conductivity and tensile response of graphene as the best advanced material with unique combination of extraordinary high thermal and mechanical properties. The optimized Tersoff potential developed by Lindsay and Broido [14] is used for introducing the carbon atoms bonding interactions. It is shown that the thermal and mechanical properties predicted by the optimized Tersoff potential are in considerable agreements with experimental results. We studied the effect of boron and nitrogen atoms substitution on the thermal conductivity and tensile response of graphene. These chemically modified graphene structures are expected to have application in the next generation of electronic devices. Boron and nitrogen are the excellent elements for doping of graphene because of their close atomic size and strong valence bonds with carbon atoms. Moreover, we studied the effects of point vacancy, bivacancy and Stone-Wales defects as the major types of defects that have been experimentally observed on the thermal conductivity and tensile response of graphene. We particularly studied the effect of curvature on the thermal conductivity of graphene. Finally, we studied the thermal conductivity and elastic constants of graphite. The observed results could be summarized as follows:

 The optimized Tersoff potential predicts the thermal conductivity of graphene in a much better agreement with the experiments than the predictions based on the original Tersoff and REBO potentials.

 The graphene presents highly ductile failure behavior  The Young's modulus of graphene is acceptably independent of loading chirality direction.

 The edge atoms in graphene nanoribbons present reduced mechanical properties than those in the center.

 The Young's modulus and tensile strength of few-layer graphene are independent of the number of atomic planes.

 The existence of point vacancies in graphene structure decreases considerably the tensile strength and Young's modulus of graphene.

 The Young's modulus of nitrogen doped graphene is almost independent of nitrogen atoms concentration (up to 6%). On the other hand, substituted nitrogen atoms decrease the tensile strength and ductile failure behavior of graphene as well.

 The boron atoms substitution has limited effects on mechanical response of graphene.

 The boron atoms concentration changes the failure behavior of graphene form ductile to brittle.

 By increasing the defects concentration, elastic modulus decreases gradually. In the case of bivacancy and point vacancy defects elastic modulus decreases by around 20% for 2% concentration of defects which is almost two times higher than the same concentration of Stone-Wales defects.

 Tensile strength and strain at failure of graphene sharply decrease by around 20% for only 0.25% concentration of defects.

 Small quantities (1%) of nitrogen atoms in graphene structure result in a noticeable decline the thermal conductivity of graphene by almost one half of that of ideal sheet.

 The thermal conductivity of graphene is considerably sensitive to boron doping. Only 0.75% concentration of boron atoms reduces the thermal conductivity of graphene by more than 60%.

 Around 1% concentration of nitrogen or boron atoms removes the chirality effect on thermal conduction of graphene.

 We found exponential trends in reduction of graphene thermal conductivity as results of defects concentrations. Our results suggest that only 0.25% concentration of defects in graphene result in significant reduction of thermal conduction property of graphene by around 50%. We also found that the effects of studied defects on the thermal conductivity of graphene are considerably close to each other.

 The high thermal conductivity of graphene is considerably sensitive to any kind of chemical doping and defects.

 The graphene could keep its high thermal conduction properties, even if it goes under high curvature angles.

 The graphite present transversely isotropic elastic stiffness and thermal conductivity response. The predicted properties by optimized Tersoff potential for graphite are in close agreement with experiments.

Introduction

Polymer materials are now widely used in many applications because of their low density and their ability to be relatively easily processed to the desired final shapes. Moreover, the mechanical and thermal properties of polymeric materials could be enhanced and modified effectively by incorporation of micro or nanosized fillers in order to provide composite materials with advanced properties for engineering applications (electronic devices, building, automotive parts, aircraft applications, etc.).

Epoxy polymers are well known for their excellent adhesion, considerable elastic stiffness and mechanical strength, high electrical insulating properties and good chemical resistance. This way, the epoxy polymers are now employed in a wide range of applications for various industries such as:

paints and coatings, industrial tooling, composite materials, electronics and adhesives [1]. Epoxy is a copolymer, which means that it is formed from two different chemicals. These are referred as the resin and the hardener (crosslinking/curing agent). The resin consists of monomers or short chain polymers with an epoxide group at either end and the hardener consists of polyamine monomers.

When the resin and crosslinking agent are mixed together, the amine groups react with the epoxide groups to form a covalent bond. Each NH group can react with an epoxide group to form a rigid and strong crosslinked structure. In this way, the final properties of epoxy could be manipulated by selecting the appropriate resin and hardeners materials and ratios.

For the design of optimized composition of resin and hardener mixture to reach the desired properties, trial-and-error experimental approaches are time consuming and expensive as well. For this purpose, modeling and computer simulation techniques have excited academic and industrial interests as alternative tools of experimental studies [2,3]. For the modeling of physical and mechanical properties of crosslinked epoxy system, classical molecular dynamics (MD) simulations have been used during the last two decades. Due to computational concerns, the accuracies of the molecular dynamics studies have been enhanced by continuing the advances in computers processing abilities. This way, first molecular dynamics studies on epoxy systems were carried out using simple bead-spring [4,5] models. The simplified bead-spring models could not accurately consider the details of the molecular structures and their influence on the physical properties.

During the past few years, more accurate modeling of crosslinked epoxy systems were carried out using elaborated force fields for introducing the atomic interactions and molecular structures. Doherty et al. [6] created PMA networks using a progressive and polymerization reaction scheme.

Yarovsky and Evans [7] developed a static crosslinking technique for construction of a crosslinked low molecular-weight water-soluble epoxy resins cured with different hardening agents. Heine et al.

[8] simulated the dynamically formed end-crosslinked PDMS networks structure and calculate the elastic modulus of structure. Gou et al. [9] built crosslinked epoxy/nanotubes composite system to study interfacial load transfer. Wu and Xu [10] used molecular dynamics simulations for crosslinking of epoxy resins to study the diffusion of water into the crosslinked networks. Fan and Yuen [11] performed MD simulations to predict thermo-mechanical properties of fully cured epoxy network consisting of EPON-862 epoxy resin and TETA (triethylenetetramine) curing agent. Lin and Khare [12] presented a single-step polymerization method for the creation of crosslinked polymers.

Varshney et al. [13] modeled two different systems of EPON-862 resin and DETDA hardener.

They describe different approaches to build highly crosslinked networks and proposed a multistep relaxation procedure. Komarov et al. [14] proposed a new crosslinking method using a coarsegrained level and then mapping into a fully atomistic model. Bandyopadhyay et al. [15] used OPLS force field to construct molecular model of the EPON-862 resin DETDA hardener crosslinked epoxy system to predict the influence of crosslinking density on the thermo-mechanical properties.

Li and Strachan [16] characterized the charge evolution during polymerization of the EPON-862/DETDA system using explicit electronegativity equalization method. Using the DREIDING force field, the obtained thermo-mechanical properties were found to be in excellent agreement with the experiments. Using the same modeling procedure, Li and Strachan [17] recently studied the effects of DGEBA resin and 33DDS hardener epoxy films thickness on the thermomechanical properties. They found that glass transition temperature, elastic modulus and yield stress of the films decrease with decreasing the thickness. Soni et al. [18] recently studied the volumetric and thermal properties of cross-linked epoxy systems consisting of DGEBA resin and POP curing agent with different lengths using molecular dynamics (MD) simulations.

The objective of this work is to propose a simple method in order to construct epoxy molecular systems consisting of diglycidyl ether of bisphenol A (DGEBA) resin crosslinked by diethylenetriamine (DETA). Our aim is not to reproduce the exact chemical kinetics, but to propose an efficient scheme to obtain realistic crosslinked structure. The DGEBA/DETA epoxy is now widely used as a reactive polymer adhesive in the industry. We should also note that the crosslinking of the DGEBA/DETA epoxy system has not been studied yet. To this aim, the OPLS united atom force field is used for describing the atomic interactions and molecular structures as well. A dynamic crosslinking scheme is used to construct the crosslinked network by the molecular dynamics 121 simulation. Then, the predicted properties are compared with the experimental values existing in the literature.

Fig. 1, Molecular structures of DGEBA resin and DETA hardener.

Molecular dynamics modeling

The epoxy system studied in this work is consisting of diglycidyl ether of bisphenol A (DGEBA) resin crosslinked by diethylenetriamine (DETA). The molecular structures of DGEBA and DETA molecules are illustrated in Fig. 1. The OPLS united-atom force field [19][20][21] was used to define the bond, angle and dihedral parameters with Lenard-Jones non-bonded interactions. In the OPLS force field, the total energy of the system is the sum of all the individual energies associated with bond, angle, dihedral, and non-bonded interactions. The OPLS potential is expressed as follows:

the bonding stretching energy, U bond , is expressed by: here, K r is the bond constant. The angle bending energy is expressed as:

here K is the angle energy constant. The dihedral energy, U dih , is expressed as: 122 where, the V factors are the OPLS dihedral constants. Finally, the non-bonded energy has the following relation:

here, the first term in the right hand side is the LJ potential for the non-bonding interactions. The second term is linked to columbic integrations, which is neglected in our modeling. Using the united atom force field, all CH 3 , CH 2 , CH, NH 2 , NH and alkyl groups were modeled as single united atoms with their corresponding masses, except for the C and H atoms in the phenyl rings of resin. Such a simplification reduces the computational cost of simulation. The Lennard-Jones (LJ) potential is used to introduce van der Waals interaction between individual atoms. The equilibrium spacing parameter and potential depth of well were taken to be the arithmetic mean and geometric mean of the respective atoms, respectively. LAMMPS [22], free and open source software is used for all the molecular dynamics simulations in this study. for bonds, angles and dihedrals information were written to a input file in the native LAMMPS format. We should note that the periodic boundary conditions were applied all the all directions.

Therefore, these simulations are representative of bulk epoxy and not the epoxy films. The time increment in all the simulations was set to 0.5 fs. The 5:2 system (consisting of 179 united atoms) is initially energy-minimized using the conjugate gradients method and then equilibrated using Nosé-Hoover thermostat (NVT) method. Next, the system was replicated for 7 times, energy minimized by conjugate gradients method and equilibrated by NVT at temperature of 1000 K and low density of 0.6 g/cm3. In the next step, the system density was increased to 0.9 g/cm 3 by gradually decreasing the periodic simulation box during a NVT run at 1000 K run for 500 ps. The high velocity of atoms at 1000 K and at the low density, help the system to be mixed completely. At this point, the energy of the system is minimized by conjugate gradients method and at the fixed density;

the temperature of the system was decreased slowly from 1000 K to 450 K using the NVT method for 2.5 ns. The final epoxy melt was created by replicating the system for 15 times (consisting of 22,912 united atoms) and performing the equilibration using the NVT method for 100 ps at 450 K and constant density of 0.9 g/cm 3 (The simulation box volume is 7.69×7.69× 7.69 nm 3 ). The obtained epoxy melt is then crosslinked. To this aim, the polymerization is carried out at 450 K and at the fixed density of 0.9 g/cm 3 . The Langevin thermostat [23] with damping parameter of 0.5 fs was used to control the temperature fluctuations during the non-equilibrium crosslinking step.

During the crosslinking process, chemical reactions are simulated in a stepwise manner using a distance-based criterion [16,17]. Bonds are created between reactive atoms if they fall within a cutoff distance of 5 Å which is checked every 10 fs with low probability of 1% to avoid drastic 124 changes in the energy of the system. In order to smooth the energy changes in the crosslinking step, the maximum movements of the atoms was limited to 0.2 Å in every time step. During the crosslinking step, the evolution of crosslinking density at different times was recorded. The crosslinking density is defined as the ratio of the total number of created bonds to the total number of epoxide groups (maximum possible number of bonds between reactive carbon and nitrogen atoms). Fig. 4 illustrates the evolution of crosslinking density as a function of time. As it can be observed, the crosslinking density increases very quickly in the early stage but then slows down while the network grows continuously. We also conducted the crosslinking procedure with higher probability for bond creation between the reactive pairs. We found that in order to avoid the drastic changes in the energy of the system, the probability of for creating the bonds should be less than 5%

during the first 100 ps of the crosslinking simulation. In the current study, the epoxy system is studied at four different crosslinking densities of 84%, 70%, 59% and 46%. All the crosslinked structures are then let to relax using isothermal, isobaric (NPT) simulation with pressure of 1 atm

for 100 ps at 450 K.

Crosslinked epoxy network

The crosslinked epoxy at temperature of 450 K was uniformly cooled down to room temperature using NPT simulations for 0.5 ns. The 84% crosslinked epoxy is shown in the Fig. 5. For this case, the crosslinked network has the density of 1.12 g/cm 3 at low pressure of around 1 atm (The simulation box volume is 7.13×7.13× 7.13 nm 3 ). . The obtained density is in a close agreement with the experimentally reported density of 1.13 g/cm 3 [24] for DGEBA/DETA crosslinked epoxy. We also found that the density of epoxy increases by increasing of the crosslinking density. For instance, the density of 46% crosslinked epoxy at room temperature is calculated to be around 1.08 g/cm 3 . In this way, the negligible underestimation of experimental density by the MD calculations could be explained because of the effect of crosslinking densities.

Fig. 5, The obtained 84% crosslinked epoxy at room temperature.

The spatial density distributions of the obtained crosslinked networks were calculated to investigate the degree of homogeneity of the structure. The homogeneity of the obtained structure is the critical point for the evaluation of the thermo-mechanical properties of the bulk epoxy systems. In this way, it is expected that by increasing the number of the atoms considered in the modeling, the structure would present higher level of homogeneity. Therefore, for a system consisting of a limited number of molecules, one should ensure the homogeneity of the structure before evaluating the properties.

Otherwise, the properties may present anisotropic nature. To address this issue, the simulation box was divided into ten equal slabs along the X, Y and Z axes and the total masses of the atoms in each 126 of those ten slabs were computed. The total masses of atoms in each slab were divided by the volume of the slab to obtain the density distribution of the system. As an instance, the spatial density distributions of the 84% crosslinked DGEBA/DETA epoxy at an arbitrary time are shown in Fig. 6.

As it can be observed, the density distribution of 84% crosslinked system is uniform along all three axes of the simulation box. Such an observation reveals the acceptable homogeneity of the obtained crosslinked epoxy network by the modeling scheme presented in this work. The glass transition temperature (Tg) is one of the important characterized property of the epoxy systems. Since some physical properties of amorphous polymers are linked with this transition point, the understanding of this phenomenon is critical. For instance, the mechanical properties of epoxy decline considerably after the glass transition temperature. Therefore, for a particular load bearing application, one should ensure that the working temperatures are acceptably below the glass transition. As a common method, in this study the glass transition is evaluated by plotting the density versus temperature relation for different crosslinked epoxy networks. To this aim, for each of the four crosslinked epoxy models, the structure was uniformly cooled down from 600 K to 250 K using NPT method for 1.1 ns at pressures of 1 atm. The temperatures and densities of the system 127 were recorded and averaged during every 25 ps intervals. We should also note that we performed the simulations at lower cooling rates for 2 ns and 4 ns and we observed that the obtained densitytemperature relations are convincingly independent of cooling rates. In Fig. 7 the computed densitytemperature relation for four different crosslinked epoxy networks are shown. As the first general finding, we could observe that by increasing the crosslinking density, the density of epoxy increase which is in agreement with previous studies [16][17][18]. This is due to increasing of the number of covalent bonds in the structure by the increasing of the crosslinking density. These covalent bonds make the structure more compressed and rigid as well. above the change in slope. However, in the obtained curves, the change in the slope is usually gradual and we could not observe sharp changes in the slope. As a result, the determination of Tg is highly sensitive to the way that the data points are analyzed. Such a problem was also reported in the similar study of Bandyopadhyay et. al [15] for the evaluation of Tg of EPON-862/DETDA using the OPLS potential function. In this way, it is more accurate to report the Tg as a temperature range rather than a single temperature value. As it is shown in Fig. 7, by fitting a line to the data points at low temperatures, the Tg could be found in the range of 400-430 K. This range is within the experimentally reported Tg of 392 K [24] and 409 K [1] for the cured DGEBA/DETA epoxy systems.

Summary

In summary, we proposed a simple method for the construction of crosslinked epoxy molecular systems. The epoxy system in our work was consisting of diglycidyl ether of bisphenol A (DGEBA) resin crosslinked by diethylenetriamine (DETA). The OPLS united atom force field was employed for introducing the atomic interactions and molecular structures as well. We used dynamic crosslinking scheme, in which chemical reactions are simulated in a stepwise manner using a distance-based criterion. The obtained crosslinked structure at room temperature was found to be considerably homogeneous with a density close to that of the experimentally created epoxy. We also calculated the glass transition temperature of epoxy, which was found to be in agreement with experimental values in the literature. It was found that the density of epoxy increases by increasing the crosslinking density for a wide range of temperatures.

Chapter 5: Multiscale modeling of thermal conduction in graphene epoxy nanocomposites

In this chapter, we present a multiscale scheme on the basis of molecular dynamics (MD) and finite element (FE) methods for evaluating the effective thermal conductivity of graphene epoxy nanocomposites. The proposed hierarchical multiscale approach includes three different scales.

First, we used MD simulations for the investigation of thermal conduction in graphene epoxy assembly at atomic scale. Our results suggest that thermal conductivity of single layer graphene decline by around 30% in epoxy matrix for two different hardener chemicals. Using MD, we also calculated thermal boundary conductance (TBC) between crosslinked epoxy and graphene sheet. In the next step, using the results obtained by the MD method, we developed FE based representative volume elements (RVE) of the nanocomposite in order to evaluate the thermal conductivity at the microscale. Finally, nanocomposite effective thermal conductivity was obtained using FE homogenization of an ensemble of microscale RVEs. The validity of the proposed approach was confirmed by comparing predicted results with experimental results in the literature.

Introduction

As discussed in chapter 3, graphene present unique combination of extraordinary high thermal conductivity and elastic stiffness properties. Having a planar structure with high surface to volume ratio, propose the graphene as an excellent candidate for the reinforcement of thermal, mechanical and barrier properties of polymeric materials. As concluded from the results in chapter 2, finite element calculations confirm the possibility of higher reinforcement in thermal and mechanical properties of polymeric materials by addition of graphene nanoplatelets in comparison with carbon nanotubes. It was recently reported [1] that the thermal conductivity of graphene-filled epoxy could be enhanced by around 23 times higher than that of the pure epoxy, which is an exciting finding for the reinforcement of thermal conductivities of polymer materials. In this way, the more fundamental understanding of graphene epoxy nanocomposite properties is the critical issue for its future applications. Due to difficulties in experimental characterization of nanostructured materials, analytical and numerical simulations are getting more attractive as alternatives.

The objective of this work is to develop an accurate multiscale approach for taking into account the atomistic effects in the evaluation of polymer nanocomposites effective thermal conductivity. The proposed multiscale scheme consists of three different scales. In the first scale, atomistic effects are studied using molecular based calculations by means of molecular dynamics (MD) simulations. In the next scale, finite element modeling was carried out to evaluate thermal conductivity response of nanocomposite representative volume elements (RVE) at microscale on the basis of molecular dynamics results. In the final scale, nanocomposite effective thermal conductivity was calculated by homogenization of results obtained in the second scale using finite element method.

In the molecular dynamics modeling, we used dynamic crosslinking scheme for construction of two different crosslinked epoxy networks. In this study, diglycidyl ether of bisphenol A (DGEBA) resin was crosslinked by diethylenetriamine (DETA) or diaminodiphenyl sulfone (DDS) hardeners with similar crosslinking densities. After obtaining the crosslinked epoxy structures, we modeled graphene epoxy assembly. Using the developed molecular models, we studied the effect of interfacial interactions between epoxy and graphene atoms on the graphene thermal conductivity. We observed that, as a result of phonon scattering, thermal conductivity of graphene may reduce by around 30%.

Using the developed graphene epoxy molecular dynamics model, we also calculated thermal boundary conductance (TBC) between epoxy and graphene. We found that TBC values and reduction in thermal conductivity of graphene are acceptably independent of hardener type and they show limited dependency on interfacial pressure. We also studied the effects of formation of covalent bonds between graphene and epoxy atoms. Molecular dynamics calculations suggest that formation of 5% covalent bonds could enhance TBC by 3 times and on the other hand it results in around 2 times higher reduction in graphene thermal conductivity. In the finite element modeling of nanocomposite RVE, obtained TBCs by molecular dynamics simulations were used to introduce contact interactions between graphene and epoxy. Moreover, in introducing the fillers thermal conductivity in finite element modeling, reduction in graphene thermal conductivity resulting from interaction with epoxy atoms was also considered. Since the aspect ratio of graphene nanoplatelets are so high and their properties are four orders of magnitude higher than that of the pure epoxy, RVE results present anisotropic nature. In this study, we proposed a new finite element based methodology to report effective thermal conductivity of nanocomposites on the basis of anisotropic RVEs results. Finally, we compared our results with recent experimental results of graphene epoxy nanocomposites reported by Shahil and Balandin [1]. Comparison for a wide range of graphene fillers volume fractions reveals a considerable validity of the proposed multiscale scheme in the evaluation of the effective thermal conductivity of polymer nanocomposites. We also note that the finite element modeling methodology proposed in this work could be used to accurately model nanocomposite structures with high volume concentrations and aspect ratios for the fillers.

Molecular dynamics modeling and results

In this study, we used molecular dynamics (MD) simulations for analyzing thermal conduction in graphene epoxy assembly using the LAMMPS [2] molecular dynamics package For the modeling of bonding interactions between carbon atoms in graphene, a wide variety of potential functions exist in the literature. Our recent studies on the thermal conductivity [3] and mechanical properties [4] of graphene have shown considerable validity of the optimized Tersoff potential [5] in the evaluation of graphene properties. Therefore, the bonding interactions between carbon atoms in graphene were modeled using the optimized Tersoff potential [5].

Epoxy is a thermoset polymer which is formed from two different chemicals. These are referred to as the resin and the hardener (crosslinking agent). The resin consists of monomers or short chain polymers with an epoxide group at either end and the hardener consists of polyamine monomers.

When the resin and crosslinking agent are mixed together, the amine groups react with the epoxide groups to form a covalent bond. Each NH group can react with an epoxide group to form a rigid and strong crosslinked structure. As shown in Fig. 1, in this study, the epoxy consists of diglycidyl ether of bisphenol A (DGEBA) resin crosslinked by diethylenetriamine (DETA) or diaminodiphenyl sulfone (DDS) hardeners. Therefore, five DGEBA molecules react with two DETA molecules and two DGEBA molecules react with one DDS molecule to form crosslinked structures. To this aim, OPLS united atom force field [6,7] was used for describing atomic interactions and molecular structures. In the OPLS force field, the total energy of the system is the sum of all the individual energies associated with bond, angle, dihedral, and non-bonded interactions. OPLS potential has been successfully used to construct crosslinked epoxy networks [8]. Using OPLS united atom force field, all CH 3 , CH 2 , CH, NH 2 , NH and alkyl groups were modeled as single united atoms with their corresponding masses, except for the C and H atoms in the phenyl rings. The Lennard-Jones (LJ) potential is used to introduce van der Waals interactions between individual atoms. The equilibrium spacing parameter and potential well depth were taken to be the arithmetic mean and geometric mean of the respective atoms, respectively. The simulation time increment was set to 0.5 fs. The system was initially energy-minimized using the conjugate gradients method and then equilibrated using Nosé-Hoover thermostat (NVT) method at 600 K and low density of 0.9 g/cm 3 .

For creating the epoxy film, two 12-6 LJ potential walls were introduced along the film normal direction. During the crosslinking process, chemical reactions are simulated in a stepwise manner using a distance-based criterion [9]. Bonds are created between reactive atoms if they fall within a cutoff distance of 5 Å being checked every 10 fs with low probability of 3% to avoid dramatic changes in the energy of the system. Evolution of crosslinking density as function of crosslinking time is shown in Fig. 2. Crosslinking density is defined as the ratio of the total number of created bonds to the total number of epoxide groups. As it can be observed, during initial 100ps time of crosslinking procedure, crosslinking density increases to around 60% for both systems. After this point, we could observe that crosslinking density increases only gradually and a saturation limit is approached. By performing the crosslinking procedure for 1 ns, final DGEBA/DETA structure has 137 crosslinking density of 76.5 %. In the case of DGEBA/DDS epoxy, simulation was performed for 2ns to get 77.5% crosslinked structure. The obtained crosslinked films were cooled down from 600 K to 300 K for 500 ps by NVT. The final structure is at low pressure of around 1 atm with the density of ≈1.11 g/cm 3 for DGEBA/DETA epoxy and density of ≈1.14 g/cm 3 for DGEBA/DDS epoxy in close agreements with the experimental densities.

Fig. 2, Evolution of crosslinking densities as functions of crosslinking time for DGEBA/DETA and DGEBA/DDS epoxy networks. At this point, the graphene structure was also assembled on the top of crosslinked epoxy films as illustrated in Fig. 3. The LJ potential was used to introduce non-bonding interactions between graphene and epoxy atoms. At this stage, by applying periodic boundary conditions in all directions, the system was relaxed to atmospheric pressure at room temperature using NPT method. Fig. 3, Crosslinked epoxy graphene assembly for the evaluation of thermal conduction in graphene epoxy nanocomposite. Periodic boundary conditions are applied in all directions which mean that graphene is in contact with epoxy from both sides.

Van der Waals interactions between polymer and fillers are not considered as strong interactions.

Therefore, for enhancement of mechanical properties of polymer nanocomposite materials, nanofillers are commonly surface modified to form covalent bonds with the matrix. We note that pristine graphene cannot form covalent bonds with the epoxy matrix. In the fabrication of graphene polymer nanocomposites, in order to improve load bearing and elastic stiffness properties of nanocomposite, graphene sheets are also modified to form covalent bonds with the polymer matrix [10][11][12]. For this purpose, graphene oxide sheets are usually functionalized to improve graphene epoxy interface strength. In this way, oxide groups on the graphene surface react with epoxide groups of hardener and form covalent bonding [12]. Effects of these covalent bonds on the effective thermal conductivity of graphene polymer nanocomposites have not been studied yet. To provide a general viewpoint, we used a simplified approach for taking into account the formation of these covalent bonds on the thermal conduction in graphene epoxy nanocomposites. In this approach, we introduce bonds between carbon atoms in graphene and carbon atoms in epoxide groups [12]. The bonding ratio is introduced as the number of created bonds with respect to the total of number of carbon atoms in graphene. Moreover, in our modeling we also studied the effect of residual pressure on the thermal conduction in graphene epoxy nanocomposites. To this aim, we also studied the effect of 1 GPa pressure on the thermal conduction of graphene in epoxy matrix.

We used non-equilibrium molecular dynamics (NEMD) method [13] to evaluate the thermal conductivity of pristine and epoxy affected graphene. In this approach, we fix the atoms at the two ends of the structure to avoid interactions through the periodic boundary condition along the heat conduction direction. The graphene sheet (excluding the fixed atoms) along the x direction was divided into 21 slabs. First, the structure was relaxed at 300 K using the NVT method. In the next step, a temperature difference was applied between the first and 21st slabs. In this study, the temperatures at the first (hot reservoir) and 21st (cold reservoir) slabs were set to 310 and 290 K, respectively, using the NVT method, while the remaining slabs are under constant energy (NVE) simulations. In order to keep the temperature differences at two ends, an amount of energy is added to the hot reservoir and at the same time another amount of energy is removed from the cold reservoir by the NVT method.

Fig. 4, Energy values added to hot reservoir and removed from cold reservoir by the NVT method.

Here, graphene is in contact with DGEBA/DDS crosslinked epoxy with 5% covalent bonds.

In Fig. 4, the energy adding and removing to the hot and cold reservoirs in the case of graphene-DGEBA/DDS epoxy with 5% covalent bonds concentration are illustrated. As it is shown, the amount of energy added to the hot reservoir is equal to removing energy from the cold reservoir which accordingly means that the total energy of the system is accurately conserved. Moreover, the slopes of energy curves in both cases are linear which means that the energy is added or removed with a constant rate. These observations imply that the energy of the simulation box remains constant and a constant heat flux is applied throughout the sample. This heat flux is a constant quantity of kinetic energy transferring from the atoms in the hot reservoir to those in the cold reservoir at each simulation time step. The non-equilibrium steady state heat transfer can be achieved after 200 ps of the exchanging process, when a temperature gradient is established along the longitudinal direction. The applied heat flux by the NVT method is calculated based on the slope of energy curves (Fig. 4). The simulations were performed for 5 ns and the averaged temperatures in each slab of the graphene sheets were computed to find a linear temperature profile.

As it is shown in Fig. 5, by neglecting the jumps at the two ends, one could find a linear temperature gradient (dT/dx) in the middle of the specimen. By calculating the established slope, thermal conductivity (K) could be calculated using one-dimensional form of the Fourier law . In order to simulate the interfacial heat transfer between the graphene and epoxy, an initial temperature difference of ∆T(0) was applied. To this aim, the epoxy and graphene atoms temperature were set at T and T+∆T(0), respectively for 100 ps using the NVT method. We should note that in this study, the ∆T(0) was set to 100 K and the temperatures were averaged over multiple runs (six in our case). If ∆T(0) is too small with respect to the substrate temperature, the recorded temperature will be very noisy [14]. To simulate the heat transfer process, the NVT is switched off and the system is allowed to relax by performing the constant energy (NVE) simulations. As shown in Fig. 6, the graphene and epoxy temperatures were recorded during the relaxation process time (t).

∆T(t) decays exponentially with a single relaxation time, τ, such that as: 141

The associated decay time (τ) is related to the graphene surface area, A, heat capacity, C, and interface conductance, TBC, through the following relationship [14,15]: Using 5.6×10 -4 Jm -2 K -1 for the heat capacity per unit area of an atomic layer of graphite [15,16],

thermal boundary conductance between graphene and epoxy could be obtained. In Table 1, obtained results by molecular dynamics calculations are presented. As the first finding we could observe that although the chemical structures of DETA and DDS are different, effect of these hardeners are close to each other. As reported in Table 1, MD simulations suggest that the thermal conductivity of graphene decreases by around 30% as a result of interfacial interactions with epoxy atoms. Although the pressure does not present significant effect on graphene thermal conductivity, we found that thermal boundary conductance increases by around 30% for the case of 1 GPa pressure. On the other hand, in the case of formation of 5% covalent bonds between graphene and epoxy atoms, we could observe that thermal boundary conductance enhanced by around three times of that with no covalent bonds. This ratio is around two times for the case of 2% covalent bonds.

However, as a result of these covalent bonds, we found that thermal conductivity of graphene decreases by around 60% in the case of 5% covalent bonds. We note that since we could not find any direct experimental measurements of thermal boundary conductance between graphene and 142 epoxy matrix, we compared our MD results with experimental measurements in the case of graphene silica [17,18]. As illustrated in Table 1, molecular dynamics results are in the same order of magnitude of reported values for thermal boundary conductance between graphene and silica. Graphene-Silica [17] ---------25 (Experimental) Graphene-Silica [18] ---------50 (Experimental) The mechanism responsible for reduction of graphene thermal conductivity in epoxy could be relatively explained by plotting the graphene and epoxy atoms phonon density of states (DOS). The DOS was computed by calculation of the Fourier transformation of atomic velocity autocorrelation functions during equilibrium, NVE simulations (no heat flux and no fixed atoms). As shown in Fig. 7 (a), the epoxy atoms do not change the relative frequency and magnitude of the maximum number of accessible phonons in the DOS. However, we could observe that at the lower frequencies, the vibrations of graphene atoms surrounded by the epoxy atoms are damped out. This damping is more obvious at frequencies lower than 13 THz. We also calculated the DOS for three epoxy atoms, which are shown in Fig, 7(b). As it can be observed, in the epoxy atoms, the maximum numbers of phonons are available at frequencies lower than 13 THz. However, we could still observe phonons at higher frequencies up to 50 THz. The overlap of epoxy and graphene atom phonons results in the scattering of phonons in graphene in a wide range of frequencies. In the pristine graphene, the thermal transport is ballistic, meaning low collision of phonons. However, as a result of the interactions with epoxy atoms at the interface, phonon scattering leads to reduction of graphene thermal conductivity.

Finite element modeling

In this section, details of the finite element modeling performed in this study are presented. In the multiscale scheme developed in this work, the finite element method was used to evaluate local thermal conductivity constants of nanocomposite RVEs at the micoscale and also for homogenization and evaluation of the nanocomposite effective thermal conductivity at the macroscale. Computational limitations of the finite element method impose limits on the maximum number of elements that are used for introducing geometries in the model. In this way, current simulations of composite structures are limited to modeling of RVEs consisting of limited number of fillers. Finite element modeling in this study has been performed using the ABAQUS (Version 6.10) package along with python scripting. As a common assumption, in our finite element modeling, graphene platelets were modeled by disc geometry. We note that the diameter to thickness ratio of a disc is commonly called aspect ratio. In Fig. 8a, a sample of 3D cubic RVE filled with 10% volume fraction of graphene platelets with aspect ratio of 100 is illustrated. Due to computational limitations and modeling complexities, in our finite element modeling numbers of perfect particles inside the micoscale RVEs were limited to 100. For a particular geometry of the fillers, the size of the RVE was accordingly adjusted based on volume fractions of fillers. In order to accurately model the desired volume concentration of fillers inside RVE, models were constructed in a way to satisfy periodicity criterion. This means that if a filling particle passes one boundary side of the RVE, the remaining part of that particle continues from the opposite side. This way, if one puts these RVE cubes together side by side, no mismatch will be observed and all of the particles will have similar geometry. As shown in Fig. 8a, for a better and more flexible meshing of RVEs, each disc was partitioned into four symmetric parts. Since in this work we aimed to model high volume fraction of graphene platelets, we assume that particles are locally unidirectional. Interaction between graphene discs and epoxy matrix in this work was modeled by introducing thermal conductance contact elements between discs and matrix interfaces. As discussed in section 2.1.2, in the evaluation of thermal conductivity along each direction, we merge two thin auxiliary parts to each RVE box. These thin auxiliary parts have the same section size as that of the RVE and their thermal conductivity are chosen to be one million to that of the matrix. Contacting surfaces of these auxiliary parts with RVE outer surfaces were thermally tied together. For loading conditions, a constant heating surface heat flux (+q) was applied on the external surface of one auxiliary part (hot surface) while on the outer surface of opposite auxiliary part (cold surface) the same magnitude cooling surface heat flux (-q) was applied. As a result of applied heat fluxes, a temperature gradient will form along the RVE length. By the use of this loading condition, the Fourier law conditions would be satisfied more accurately than simply applying the heat fluxes on the outer surfaces of the RVE (see section 2.1.2). Temperature gradient, dT/dz, is evaluated by calculating temperature differences between RVE external hot and cold surfaces. Finally, thermal conductivity of RVE along the desired direction (K z ) can be obtained using one-dimensional form of the Fourier law . In Fig. 9, experimental results are compared with multiscale results for graphene sheets with 0 and 5% covalent bonding with epoxy atoms. In the finite element modeling of nanocomposite RVEs, the thermal conductivity of fillers were modified on the basis of values reported in Table 1. We note that in the experimental study [1], the aspect ratios of fillers were assumed to be 100, thus, in our finite element modeling we also used the same value for the fillers aspect ratio. Thermal boundary conductance values obtained by MD calculations (Table 1) are used to introduce contact conductance between fillers and matrix in the finite element RVE modeling. Thermal conductivity constants of 3 different RVEs were calculated and using the finite homogenization scheme, nanocomposites effective thermal conductivities are reported. As it can be observed, multiscale modeling results are in close agreement with experimental values confirming the validity of the multiscale scheme developed in this work. The reported results in Fig. 9 suggest that formation of covalent bonds between graphene and epoxy atoms slightly decrease (less than 5%) the effective thermal conductivity of nanocomposite. This reduction is because of destructive effect of these 148 covalent bonds on the graphene thermal conductivity. This also suggests that three times enhancement of thermal boundary conductance does not improve effective thermal conductivity of nanocomposite considerably. Here, we note that it is quite well known that the formation of these covalent bonds between epoxy and graphene atoms present critical effects on the mechanical reinforcement of nanocomposites. Therefore, since the reduction in nanocomposite effective thermal conductivity is shown to be insignificant, the formation of covalent bonds between epoxy and graphene is recommended to improve mechanical properties and it is expected that the nanocomposite thermal conductivity response would not be affected considerably.

Fig. 10, A typical developed 3D finite element RVE for 2% concentration of randomly oriented platelet fillers with aspect ratios of 100.

As discussed in section 5.3, in the finite element modeling of RVEs, we assumed that fillers are locally unidirectional. This assumption is mainly due to the fact that in the modeling of a limited number of particles with high aspect ratios and high volume concentrations, the RVEs present anisotropic nature and particles arrange in parallel [20]. To verify the error that might be raised by this assumption, we developed RVEs for 2%concentration of randomly oriented fillers with the same aspect ratios shown in Fig. 10. Although the volume fraction is low, we still observe clusters of particles that are arranged in parallel. Such an observation suggests that the developed RVE also present anisotropic nature. Therefore, we calculated the thermal conductivity constants of 3 different RVEs with randomly oriented fillers and then we used the finite element homogenization scheme in order to calculate effective thermal conductivity of the nanocomposite. Interestingly, we found that the difference in the calculated effective nanocomposite thermal conductivities is less than 3% which shows the accuracy of the developed modeling scheme for higher fillers volume concentrations. Such an observation is interesting and confirms that the finite element modeling methodology proposed in this work could be used accurately for modeling of nanocomposite structures with high volume concentrations and aspect ratios of the fillers. We also note that in the modeling scheme developed in this work we did not take into account the effect of interphase between graphene and epoxy. As discussed in section 2.2, the finite element results show that for the platelet fillers, the interphase effect is negligible for the cases that the aspect ratio is higher than 40.

Therefore, since the aspect ratio of graphene fillers is assumed to be 100, we could accurately neglect the interphase effect. However, we should note that the interphase effect should be considered in the modeling of nanocomposites filled with spherical particles [21].

Summary

In summary, a novel multiscale scheme on the basis of molecular dynamics (MD) and finite element (FE) methods is developed for the evaluation of effective thermal conductivity of graphene epoxy nanocomposites. In the proposed multiscale approach, the molecular dynamics method is used in order to investigate thermal conduction of fillers and matrix at the nanoscale. To this aim, we constructed crosslinked epoxy graphene molecular models. Molecular dynamics results suggest that the thermal conductivity of single layer graphene, in epoxy matrix, decline by around 30%. Using MD method, we also calculated thermal boundary conductance (TBC) between crosslinked epoxy and the graphene sheet. Based on the results acquired by the MD method, we developed finite element based representative volume elements (RVE) in order to evaluate local thermal conductivity constants of the nanocomposite. In this modeling, TBC values obtained by the MD simulations were used for introducing thermal contact interactions between fillers and matrix in the finite element modeling. Moreover, graphene thermal conductivity was modified on the basis of MD results. At the final step, the nanocomposite effective thermal conductivity was calculated using the finite element approach for homogenization of anisotropic thermal conductivity constants of 3 different microscale RVEs. The proposed finite element methodology in this work could be used as

In our studies, we used combination of atomistic and continuum methods for the evaluation of thermal conductivity and elastic stiffness of nanostructured materials. The molecular dynamics method was used to probe and investigate the materials properties at atomic scale. Then, the bulk properties of materials were evaluated using finite element, Mori-Tanaka and strong-contrast, continuum approaches. Finally, predicted properties were compared with existing experimental results and the accuracy of the developed modeling approaches was discussed.

First, we compared three different groups of modeling techniques for the evaluation of thermal conductivity and elastic modulus of random two-phase composite materials. These three methods are: Mori-Tanaka (mean field approach), finite element (numerical method) and strong contrast (statistical continuum approach). Our results show that finite element method is the most accurate method and it could capture more details in the prediction of effective properties of the nanocomposite materials. On the other hand, Mori-Tanaka method is shown to be a fast as well as a valid alternative for the finite element modeling within a limited range of fillers geometries. Our results reveal that the strong contrast method based on statistical two-point correlation functions could not accurately describe the fillers geometry effects. We also developed 3D finite element modeling for the investigation of interphase effects on the elastic modulus and thermal conductivity of nanocomposite materials. Our results show that while the interphase effect is significant for the spherical fillers, it turns to be less effective as the fillers geometry deviate more from spherical shape.

We also proposed a global statistical approach for taking into account the nanocomposite microstructure in the investigation of nanocomposite effective properties filled with spherical particles. The accuracy of the proposed global statistical approach was verified by comparing the modeling results with experimental tests on silica/epoxy nanocomposites. In the continuum modeling of nanocomposites, our works were always limited within a linear response of materials.

For the future work, the developed finite element models should be modified in order to include more complicated constitutive laws for the matrix, fillers and interphase. As an instance, the finite element models could be developed for visco-plastic materials, and compare the results with viscoplastic self consistent approach.

Graphene, the two-dimensional structure of honeycomb carbon atoms is a unique nanostructured material due to its exceptionally high mechanical, thermal and electrical properties. We used classical molecular dynamics simulations for the investigation of the effective parameters on the mechanical response and thermal conductivity of graphene. The validity of the molecular dynamics procedures developed in our works was verified by observing close agreements with the latest experimental results in the literature. For the chemical modification of graphene, boron and nitrogen are the best candidates because of their close atomic size and strong valence bonds with carbon atoms. Thus, we studied the thermal conductivity and tensile response of the boron and nitrogen doped graphene.

Our results show that the thermal conductivity of graphene decreases considerably by small content of doping elements. We observed that the elastic modulus of graphene does not change considerably by addition of boron and nitrogen atoms. Similar to other known materials, experimentally fabricated graphene sheets are not generally perfect and different types of defects may exist. Defects in the graphene are usually in the form of point vacancy, bivacancy, Stone-Wales and grain boundaries. Accordingly, we also studied the thermal conductivity and mechanical response of defective graphene. We found that thermal conductivity of graphene is more sensitive to existence of defects in comparison with chemical doping. On the other hand, we observed that elastic modulus, tensile strength and strain at failure of graphene decrease noticeably by increasing of defect concentrations. Moreover, we also studied the effect of curvature on the thermal conductivity of graphene. Interestingly, we found that the graphene could keep its high thermal conduction properties, even when it goes under high curvature angles. Here, we should note that there are wide varieties of interatomic potentials available for modeling of carbon atoms such as Tersoff, Brenner, AIREBO, REBO, EAM, MEAM and forcefields such as CHARMM, AMBER and ReaxFF.

Therefore, it will be interesting that in the future works compare accuracies of different potential functions in the prediction of thermal and mechanical properties of graphene and graphite. The computational costs of different potential functions should be also compared in order to provide a complete vision for the researchers. Moreover, in order to form covalent bonds between graphene and polymer matrixes, the graphene sheets are commonly oxidized. In the future work it will be also interesting to study the thermal conductivity and mechanical properties of oxidized graphene sheets.

We also proposed a simplified method in order to construct crosslinked epoxy molecular systems.

The OPLS united atoms force field was employed for introducing the atomic interactions. The dynamic crosslinking approach was used to construct the epoxy networks. It was found that the density of epoxy increases by increasing the crosslinking density. In the next step, we developed a combined molecular dynamics-finite element multiscale approach for the investigation of thermal conduction in graphene-epoxy nanocomposites. Our results suggest that the thermal conductivity of single layer graphene may decline by 30% in the epoxy matrix. We also report the thermal boundary conductance between the crosslinked epoxy and single-layer graphene sheet, which is shown to be in the same order of magnitude of experimental estimations. The underlying mechanisms are also

Fig. 1 ,

 1 Fig. 1, exemples d'un volume élémentaire représentatif (VER) cubique avec différentes géométries de charges orientées aléatoirement et recouvertes d'une couche représentant l'interphase. De droite à gauche, les charges sont des disques, des sphères et des cylindres.

Fig. 2 ,

 2 Fig. 2, Illustration schématique de l'approche statistique globale présentée ici (a) une image MEB d'un nanocomposite epoxy renforcé par des particules de silice. (b) digitalisation de cette image pour obtenir une image binaire. (c) fonction de corrélation à deux points (FCDP) calculée sur cette image binaire. (d) Evaluation par la method du "strong contrast" du module élastique local à partir de la fonction FCDP précédente. (e) Modèle par éléments finis utilisé pour moyenner les propriétés obtenues sur chaque image, et finalement d'obtenir les propriétés globales du nanocomposite.

  structure réticulée réaliste. L'époxyde DGEBA/DETA est actuellement beaucoup utilisé comme adhésif réactif dans l'industrie. Nous avons choisi d'utiliser le champ de force « Optimized Potentials for Liquid Simulations -united atom » (OPLS-ua) qui est utilisé pour décrire les interactions atomiques mais aussi les structures moléculaires. Une procédure de réticulation dynamique est développée pour construire le réseau réticulé par la simulation de dynamique moléculaire. La structure réticulée obtenue à température ambiante se caractérise par sa grande homogénéité et une densité proche de celle obtenue expérimentalement. Nous avons également calculé la température de transition vitreuse de l'époxyde, qui se trouve être en bon accord avec les valeurs obtenues expérimentalement dans la littérature. On observe que la densité de l'époxyde augmente en même temps que la densité de réticulation, et ce pour une large gamme de température. La structure réticulée de la DGEBA/DETA obtenue dans cette étude est illustrée dans la figure3.

Fig. 3 ,Fig. 4 ,

 34 Fig. 3, structure réticulée de DGEBA/DETA obtenue à température ambiante.

Fig. 5 ,

 5 Fig.5, Comparaison de la modélisation multi-échelle avec les résultats expérimentaux pour l'amélioration de la conductivité thermique sur un nanocomposite époxyde avec incorporation de graphène ou de graphène multi-couche..

Theory 34 a

 34 specific function of coordinates of the particles, called potential energy. The potential function has the general structure of: where r are position vectors of the particles, and the function Um is called the m-body potential.

Fig. 3 ,

 3 Fig. 3, Short-range potentials: (a) Lennard-Jones and (b) Morse [5].

Fig. 5 ,

 5 Fig.5, Developed 3D cubic periodic RVEs containing100 perfect particles with different geometries, (a) tubes with aspect ratio of 60, (b) spherical inclusions and (c) platelets with aspect ratio of 100. Under the each RVE sample, the meshed specimen is also illustrated.

Fig. 1 ,Fig. 2 ,

 12 Fig.1, Effects of RVE size (number of particles) on the composite effective elastic modulus (Ec/Em) and thermal conductivity (Kc/Km) for spherical (a and b) and for platelets (c and d) inclusions. The reported bands are based on six different simulations. The volume concentration is 1% and contrast in inclusion to matrix properties is 500.

Fig. 5 ,

 5 Fig.5, The effects of fillers geometry on the (a) elastic modulus and (b) thermal conductivity of twophase random composite materials at constant volume fractions of 1% and different fillers to matrix (FMC) contrast in material properties of 20 and 1000. The subscripts c, f and m stand for composite, fillers and matrix, respectively.

Fig. 6 ,

 6 Fig.6, Interphase reinforcement ratio for (a) elastic modulus and (b) thermal conductivity of random nanocomposite materials for different fillers to matrix (FMC) and interphase to matrix (IMC) contrast in material properties. The interphase thickness is 0.4 for all the presented results here.

Fig. 7 ,

 7 Fig.7, The effects of interphase to matrix (IMC) properties contrast, thickness and volume fraction on the elastic modulus and thermal conductivity of the nanocomposite samples filled with spherical particles. The filler to matrix properties contrast is fixed at 20.

Fig. 8 ,

 8 Fig. 8, 3D cubic RVEs for the evaluation of longitudinal elastic modulus and thermal conductivity of unidirectional nanocomposites filled by (a) disc and (b) cylindrical particles.

Fig. 9 , 2 . 4 . 1 10 .

 924110 Fig.9, The effect of fillers geometry on the longitudinal (a) elastic modulus and (b) thermal conductivity of nanocomposite samples filled with unidirectional fillers (the fillers volume fraction is 1%). The effect of interphase on the reinforcement of longitudinal (c) elastic modulus and (d) thermal conductivity of unidirectional filled nanocomposites (the interphase thickness is 0.4).

Fig. 10 ,

 10 Fig. 10, Schematic illustration of the global statistical approach presented in this study. (a) A scanning electron microscopy (SEM) image from 2 wt. % silica+epoxy nanocomposite. (b) Digitization of the SEM image to a binary image. (c) Calculated two-point correlation function (TPCF) of the digitized image. (d) Strong-contrast evaluation of the local elastic modulus using the TPCF information provided from SEM image. (e) The finite element model for averaging of eight different SEM image as the representation of the nanocomposite (each color belongs to one SEM image).

Fig. 12 .Fig. 12 ,

 1212 Fig. 12. Such a RVE was constructed in ABAQUS by developing a Python user subroutine. We also developed a C++ program for creating the random distribution of platelets inside the RVE and obtaining the anisotropic materials constants in the RVE global coordinates. The provided information by the C++ code was used as input data for Python script for the construction of the final models in the ABAQUS.

Fig. 2 ,

 2 Fig. 2, Different types of defects in equilibrated graphene sheets.

Fig. 3 ,

 3 Fig.3, Molecular dynamics simulation setup for prediction of thermal conductivity of the (a) defectfree or chemically doped graphene. (b) Circumferential and (c) longitudinal heat conduction for 180 degree curved graphene sheet. Curved graphene modelings were performed in a way that the desired heat flow is established along the armchair direction and the length of the specimen along the X directions are equal.

Fig. 4 ,

 4 Fig. 4, (a) Energy values added to hot reservoir and removed from cold reservoir by the NVT method. (b) Established steady-state temperature profile along the graphene showing a linear relation in the middle of sample. These results are for the case of 2% concentration of Stone-Wales defects in graphene.

1×10 9 sFig. 5 ,

 95 Fig. 5, Calculated strain energy response for armchair single layer graphene at different loading strain rates by optimized Tersoff (OpT) potential compared with original Tersoff (OrT) potential.

Fig. 6 ,

 6 Fig. 6, (a-d) Ductile deformation behavior of SLG with chirality angle of ≈10º by optimized Tersoff, (e-h) brittle deformation of an armchair SLG obtained by original Tersoff potential.

Fig. 7 ,

 7 Fig. 7, Obtained uniaxial engineering stress-strain relation of SLG with different loading chirality angles at strain rate of 1×10 8 s -1 .

Fig. 10 ,

 10 Fig. 10, Predicted Young's modulus (Y) and ultimate tensile strength (UTS) of armchair single layer graphene with various concentrations of Nitrogen atoms and point vacancies, at strain rate of 1×10 9 s -1 .

Fig. 12 ,

 12 Fig. 12, Tensile deformation of single layer graphene with 6% concentration of nitrogen atoms at various stages of loading, at strain rate of 1×109 s -1 .

Fig. 14 ,

 14 Fig. 14, Elastic modulus, tensile strength and strain at failure of graphene as function of point vacancy, bivacancy and Stone-Wales defects concentrations.

Fig. 20 ,

 20 Fig. 20, Carbon atom longitudinal DOS for SLG for a carbon atom (a) with no doping and (b)directly bonded with a nitrogen atom.

Fig. 23 ,

 23 Fig. 23, Average in-plane and along the thickness normalized heat current auto-correlation function (HCACF) obtained from MD simulations for bulk graphite. The results are averaged for eight different simulations.

Fig. 2 ,

 2 Fig.2, As the amine groups approach, the reactive sites in the epoxy resin are activated by opening the epoxy rings at the ends of the molecule.

Fig. 3 ,

 3 Fig. 3, Different reaction rates for the amines of DETA.

Fig. 4 ,

 4 Fig. 4, Evolution of crosslinking density versus the crosslinking time.

Fig. 6 ,

 6 Fig. 6, Spatial Density distribution for 84% crosslinked DGEBA/DETA epoxy along the length of simulation box (the partition size is 7.13 Å).

Fig. 7 .

 7 Fig. 7. The calculated density temperature relation for 46%, 59%, 70% and 84% crosslinked DGEBA/DETA epoxy systems.

Fig. 1 ,

 1 Fig. 1, Molecular structures of DGEBA, DDS and DETA molecules. The epoxy structure was created by replication of the initial structure. This way, for DGEBA/DETA epoxy, initial 5:2 structure was replicated 96 times (17184 united atoms) and for DGEBA/DDS epoxy, initial 2:1 structure was replicated 216 times (19656 united atoms). The

Fig. 5 ,

 5 Fig. 5, Established steady-state temperature profile along the graphene showing a linear relation in the middle of sample. Here, graphene is in contact with DGEBA/DDS crosslinked epoxy with 5% covalent bonds.

Fig. 6 ,

 6 Fig. 6, Typical averaged temperatures of graphene and epoxy atoms for the evaluation of thermal boundary conductance at T=300K showing an exponential decay in temperature difference, ΔT(t).

Fig. 7 ,

 7 Fig. 7, (a) The phonon density of state (DOS) of pristine and epoxy affected graphene at room temperature, showing the considerable damping of phonons at frequencies lower than 13 THz. (b) Typical DOS for three atoms in the crosslinked DGEBA/DETA epoxy network.

Fig. 8 ,

 8 Fig. 8, (a) Finite element modeling of representative volume (RVE) of graphene epoxy nanocomposite. (b) Finite element modeling for homogenization of anisotropic RVEs thermal conductivity constants.

Fig. 9 ,

 9 Fig. 9, Comparison of multiscale modeling results with experimental results of graphene and multilayer graphene epoxy nanocomposites [1].

  Nous avons également étudié les effets dus à la formation de liaisons covalentes entre les atomes de graphène et d'époxyde. Les simulations de dynamique moléculaire suggèrent que la formation de 5% de liaison covalentes pourrait augmenter d'un facteur 3 la CTI mais pourrait aussi réduire d'un facteur 2 la conductivité thermique du graphène. Lors de la modélisation par éléments finis du volume élémentaire représentatif du nanocomposite, les CTI calculées par DM sont utilisés pour introduire les paramètres d'interactions de contact entre le graphène et l'époxyde. De plus, lors de l'introduction des valeurs de conductivité thermique des charges dans le modèle par éléments finis, la baisse de la conductivité thermique du graphène due à l'interaction avec les atomes de l'époxyde a été prise en compte. En raison de la valeur élevée du facteur de forme des nanofeuillets de graphène et de leurs propriétés 4 fois plus grande que celles de l'époxyde pur, le

(DDS), avec des densités de réticulation similaires dans les deux cas. Après avoir obtenu les structures réticulées, nous avons modélisé les composites graphène-époxyde. Grâce au développement de modèles de dynamique moléculaire, nous avons étudié l'effet des interactions interfaciales entre les atomes de l'époxyde et ceux du graphène afin de calculer comment est perturbé la conductivité thermique du graphène. Nous avons observé que, en raison de la diffusion des phonons, la conductivité thermique du graphène est réduite d'environ 30%. Le mécanisme responsable de cette réduction est étudié en traçant les densités d'état de phonon pour les atomes de graphène et ceux de l'époxyde. Nous avons pu observer que les vibrations des atomes de graphène entourés par ceux d'époxyde sont amorties, ce qui entraîne une réduction de la conductivité thermique du graphène. En utilisant ce modèle de dynamic moléculaire, nous avons aussi calculé la conductivité thermique interfaciale (CTI) entre l'époxyde et le graphène. Nous avons montré que ces valeurs de CTI et que la réduction de la conductivité thermique du graphène sont raisonnablement indépendant du type de durcisseur et ne dépendent que faiblement de la pression interfaciale. VER présente des propriétés anisotropes. Pour préciser ce point, nous avons proposé une nouvelle modélisation par éléments finis pour calculer les conductivité thermiques effectives des nanocomposites sur la base des propriétés anisotropes du VER. Dans la figure 5, Les résultats de la modélisation multi-échelle sont comparés aux résultats expérimentaux du nanocomposite époxyde avec incorporation de graphène ou de multi-couches de graphène, tel que trouvés dans l'article de Shalil et Balandin (Nano Letters 2012, 12, p. 861-867). Comme observé ici, les résultat de la modélisation multi-échelle sont en bon accord avec les tests expérimentaux qui confirment la validité de la stratégie de modélisation adoptée. Les résultats donnés à la figure 5 montrent que la formation de liaisons covalentes entre le graphène et l'époxyde diminue faiblement (moins de 5%) la conductivité thermique effective du nanocomposite. Cette réduction s'explique par l'effet de ces liaisons covalentes qui tend à diminuer la conductivité thermique du graphène. Par conséquent, puisque la réduction de la conductivité thermique effective est faible, la formation de liaisons covalentes entre l'époxyde et le graphène est recommandé pour améliorer les propriétés mécaniques du composite et la baisse associée de la conductivité thermique du nanocomposite est considérée comme négligeable.

Table 1 ,

 1 Comparison between statistical strong contrast results with Finite Element (FE) and Mori-Tanaka (MT) methods for different fillers geometry (aspect ratio, AR)

	E c /E m	K c /K m

  Additional stirring is performed by ultrasonic dispersion 66 during 15 min at 50°C. The aliphatic hardener, DETA, was added in a stoichiometric ratio to epoxy resin under magnetic stirring at 50°C. Curing process was performed for one hour at 60 °C, another one hour at 100 °C and finally 2 hour of post-curing at 120 °C. 2.4.2.2 Dynamic mechanical analyzer Pure epoxy resin samples and nanocomposites samples were cut to target dimensions 35 × 7.5 × 2.5 mm3. Specimens are tested in single cantilever loading in elastic range by a dynamic mechanical analyzer (DMA) (Netzsch 242C). Tests are performed between 25 and 180°C with a heating rate of 3 °C/min and a frequency of 5 Hz. The tests were conducted on both pure epoxy and epoxy nanocomposite with 2 wt. % silica. Average modulus values were calculated at each temperature from two specimens.

2.4.2.3 Scanning electron microscopy (SEM)

Nanocomposite samples are prepared by cutting and polishing to obtain a very smooth surface. Last step is performed with a 1μm diamond polishing suspension. SEM analyses were carried out on a QUANTA 200 FEG from FEI with specimen chamber pressure at 150 Pa in backscattering mode (GAD). A series of pictures is recorded in the center of the sample to avoid influence of sample boundary.

Table 1 ,

 1 Comparison between the elastic constants of graphite obtained by molecular dynamics simulation and experimental studies (units are in GPa).

		C 11	C 12	C 13	C 33	C 44	C 66
	Present MD study	970	95	8	26	---	449
	Ultrasonic and static tests [47]	1060	180	15	36.5	4	440
	Inelastic neutron scattering [48]	1440	---	---	37.1	4.6	460
	Inelastic x-rays [49]	1109	139	0	38.7	5	485

Table 1 -

 1 Calculated thermal boundary conductance and reduction in graphene thermal conductivity for graphene epoxy assembly obtained by molecular dynamics simulations.

	Structure	Pressure Covalent	Reduction in	Thermal boundary
			bond (%)	thermal conductivity	conductance (MW m -2 K -1 )
				(%)
	Graphene-	1 atm	0	29±3	10.8
	DGEBA/DETA
	Graphene-	1 GPa	0	29±3	13.5
	DGEBA/DETA
	Graphene-	1 atm	0	30±3	10.7
	DGEBA/DDS
	Graphene-	1 GPa	0	30±3	14
	DGEBA/DDS
	Graphene-	1 atm	2	50±5	20.2
	DGEBA/DETA
	Graphene-	1 atm	5	57±6	31.7
	DGEBA/DETA
	Graphene-	1 atm	2	47±5	20.5
	DGEBA/DDS
	Graphene-	1 atm	5	65±6	31.8
	DGEBA/DDS

Acknowledgements

an accurate method for the evaluation of composite effective properties filled with high volume fraction and high aspect ratio of fillers. Comparison between our multiscale modeling results with experimental results for a wide range of graphene concentrations confirm the validity of the proposed multiscale scheme in the evaluation of thermal conduction of polymer nanocomposites.

We also studied the effect of formation of covalent bonds between graphene and epoxy atoms on the effective thermal conductivity of nanocomposites. These covalent bonds present critical effects on the stiffness and failure reinforcement of nanocomposite structures. However, the modeling results suggest that formation of covalent bonds between graphene and epoxy atoms reduces slightly (less than 5%) the effective thermal conductivity of nanocomposites.

I am also grateful to the department of Advanced Materials and Structures (AMS

We remind that due to limited number of particles and also unidirectional distribution of graphene platelets, developed RVEs present anisotropic thermal conductivity nature. In order to report effective thermal conductivity of the nanocomposite (macroscale) for a particular volume fraction of graphene fillers, we need to perform homogenization of results obtained from different anisotropic RVEs. To this purpose, we use finite element approach as shown in Fig. 8b. In this method, first we calculate thermal conductivity constants of several (3 samples in this study) RVEs. Then, we created a cubic specimen and meshed this cube by eight-node linear heat transfer brick (DC3D8) elements (represented by 20 3 elements). In the next step, we introduced a new anisotropic material for every created element in the specimen. An anisotropic material has six thermal conductivity constants.

These thermal conductivity constants were calculated by randomly selecting an RVE result and transforming the selected RVE thermal conductivity second order tensor from the global coordinate system to another randomly oriented coordinate system. Since too many elements are used in this modeling, the generated random microstructure presents isotropic nature and could be considered as the representation of nanocomposite effective properties. By applying surface heat fluxes at outer surfaces of the developed cube, effective thermal conductivity of nanocomposite was obtained using the one-dimensional form of the Fourier law. Such a modeling scheme could be also used for homogenization of composite elastic properties. The proposed finite element homogenization scheme could be considered as an accurate solution in response to variation of finite element results in the evaluation of composite effective properties by modeling of a limited number of fillers.

Multiscale modeling results

In Fig. 9, multiscale modeling results are compared with experimental results of graphene and multilayer graphene epoxy nanocomposites [1]. As shown, in the experimental results, distinct difference exists for different types of epoxy matrixes. We note that recent experimental study suggests that the thermal conductivity of graphene decreases by increasing of graphene number of atomic layers [19].

The high quality graphite has in-plane thermal conductivity of 2000 W/m-K and single layer graphene present thermal conductivity of 4100±500 W/m-K [19]. Therefore, we assume the graphene and multi-layer graphene platelets have an average thermal conductivity of 3000 W/m-K.

Conclusions and future works

discussed by comparison of the phonon density of states for pristine and epoxy affected graphene sheet. In the next step, we used finite element modeling for the evaluation of bulk thermal conductivity of nanocomposite on the basis of results obtained by molecular dynamics method. By comparing the modeling results with the latest experimental results in the literature, the validity of the proposed multiscale approach was confirmed. We note that in the developed molecular dynamics models of crosslinked epoxy, for reducing the computational costs of calculations, the columbic interactions were neglected. In the future work, our study could be modified by considering the columbic interactions in the modeling of epoxy.

Publications:

B. Mortazavi Abstract:

Nanostructured materials are gaining an ongoing demand because of their exceptional chemical and physical properties. Due to complexities and costs of experimental studies at nanoscale, computer simulations are getting more attractive as experimental alternatives. In this PhD work, we tried to use combination of atomistic simulations and continuum modeling for the evaluation of thermal conductivity and elastic stiffness of nanostructured materials. We used molecular dynamics simulations to probe and investigate the thermal and mechanical response of materials at nanoscale.

The finite element and micromechanics methods that are on the basis of continuum mechanics theories were used to evaluate the bulk properties of materials. The predicted properties are then compared with existing experimental results.