Multiscale simulation of fullerene reinforced composite structures: From molecular dynamics to finite element continuum mechanics

The aim of the present study is to propose a multiscale computational technique for the prediction of the elastic mechanical properties of nanoreinforced composites. The proposed method utilizes a molecular dynamics (MD) based numerical scheme to capture the mechanical behaviour of the nanocomposite at nanoscale and then a classical continuum mechanics (CM) analysis based on the finite element method (FEM) to characterise the microscopic performance of the nanofilled composite material. The material under investigation is polyamide 12 (PA 12) randomly reinforced with fullerenes C60. At the first stage of the analysis, in order to capture the atomistic interfacial effects between C60 and PA 12, a very small cubic unit cell containing a C60 molecule, centrally positioned and surrounded by PA 12 molecular chains, is simulated via MD. Interand intra-molecular atomic interactions are described by using the Condensed Phase Optimized Molecular Potential for Atomistic Simulation Studies (COMPASS). According to the elastic properties data arisen by the MD simulations, an equivalent finite element volume with the same size as the unit cell is developed. At the second stage, a CM micromechanical representative volume element (RVE) of the C60 reinforced PA 12 is developed via FEM. The matrix phase of the RVE is discretised by using solid finite elements which represent the PA 12 mechanical behaviour while each C60 location is meshed with equivalent solid finite elements. Several multiscale simulations are performed to study the effect of the nanofiller volume fraction on the mechanical properties of the C60 reinforced PA 12 composite. Comparisons with other corresponding experimental results are attempted, where possible, to test the performance of the proposed method. * Corresponding author: nanif@mech.upatras.gr © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). MATEC Web of Conferences 188, 01013 (2018) https://doi.org/10.1051/matecconf/201818801013 ICEAF-V 2018


Abstract.
The aim of the present study is to propose a multiscale computational technique for the prediction of the elastic mechanical properties of nanoreinforced composites. The proposed method utilizes a molecular dynamics (MD) based numerical scheme to capture the mechanical behaviour of the nanocomposite at nanoscale and then a classical continuum mechanics (CM) analysis based on the finite element method (FEM) to characterise the microscopic performance of the nanofilled composite material.
The material under investigation is polyamide 12 (PA 12) randomly reinforced with fullerenes C60. At the first stage of the analysis, in order to capture the atomistic interfacial effects between C60 and PA 12, a very small cubic unit cell containing a C60 molecule, centrally positioned and surrounded by PA 12 molecular chains, is simulated via MD. Inter-and intra-molecular atomic interactions are described by using the Condensed Phase Optimized Molecular Potential for Atomistic Simulation Studies (COMPASS). According to the elastic properties data arisen by the MD simulations, an equivalent finite element volume with the same size as the unit cell is developed.
At the second stage, a CM micromechanical representative volume element (RVE) of the C60 reinforced PA 12 is developed via FEM. The matrix phase of the RVE is discretised by using solid finite elements which represent the PA 12 mechanical behaviour while each C60 location is meshed with equivalent solid finite elements.
Several multiscale simulations are performed to study the effect of the nanofiller volume fraction on the mechanical properties of the C60 reinforced PA 12 composite. Comparisons with other corresponding experimental results are attempted, where possible, to test the performance of the proposed method.

Introduction
Polymer composites are increasingly required for numerous industrial applications such as electrical and thermal insulators, strong but lightweight components in automotive and aeronautical applications and many others. However, the use of polymeric materials, generally, presents specific limitations which need to be overcome and it seems that the discovery of carbon based nanostructured materials showed the appropriate path for eliminating them.
Carbon is capable of forming many nanostructured allotropes such as graphene, carbon nanotubes (CNTs) and fullerenes which present excellent structural, electronic, optical, chemical and mechanical properties. Evidently, these carbon allotropes due to their unique material characteristics have a strong potential to be used as fillers in polymeric materials with the main purpose of enhancing their mechanical, thermal, chemical and electrical performance. Nowadays, despite the high cost of the relevant manufacturing processes, carbon nanomaterial based nanocomposites have already illustrated a significant growth potential in various industry sectors such as manufacturing, automotive, aerospace, electronics, green technologies and medical applications [1]. Nanocomposites which are reinforced with carbon allotropes belong to an emerging class of organic hybrid materials that exhibit improved mechanical properties at very low loading levels compared with conventional microcomposites [2,3]. Possibly, graphene and CNTs play the most promising role because of their superior structural and functional properties such as high aspect ratio, high mechanical strength and high electrical properties [4]. Both graphene and CNT based nanocomposites usually possess the properties of fracture toughness, lightweight, fatigue strength and cost-effectiveness. Nanocomposites reinforced with graphene and CNTs are also utilized in energy storage and conversion devices as well as in biomedical applications such as gene delivery and nano-medicine. However, despite the fact that CNTs presents a similar mechanical behavior with graphene, graphene nanocomposites have gained most of the research attention due to their better electrical and thermal performance.
Although graphene and CNT based nanocomposites have attracted tremendous scientific and industrial interest, the reinforcement of conventional materials with fullerenes has not be studied to the same extent, despite the fact the specific nanomaterials have gained one of the first places in modern material science related to carbon since their discovery in 1985 [5]. The C60 molecule [6] which is the most known and representative of the fullerene family, has 12 pentagons located around the vertices of an icosahedron and 20 hexagon rings placed at the centers of icosahedral faces. This spherical molecular structure leads to unique mechanical property characteristics [4,[7][8] which may be proved valuable for the development of new promising isotropic-like nanocomposites. In order to utilize the outstanding properties of fullerenes, they may be incorporated into polymeric matrices in the form of filler materials.
Considering the above, the present study aims to investigate the mechanical behaviour of a polymeric material reinforced with fullerenes C60. In view of the fact that experimental procedures [9][10][11][12] are excessively costly for the material characterization of polymerfullerene nanocomposite systems, numerical techniques such as MD [13,14] seem to be very attractive for predicting and investigating their performance. Therefore, in the present study a hybrid numerical multiscale technique based on MD, FEM and unit cell methods [15] has been implemented to characterise the elastic mechanical properties of a polymeric material reinforced uniformly with molecules C60. Among numerous candidate polymer materials, the PA 12 has been chosen as the matrix material of the investigated nanocomposite since it has already been experimentally demonstrated that the addition of graphene nanoplatelets or CNTs, which are analogous to fullerenes, at even small concentration up to 2 wt% within pure PA 12, may lead to significant improvements of its mechanical and electrical behaviour [12,16].

Multiscale Numerical Approach
The composite material under investigation consists of C60 molecules as reinforcements and PA 12 as a polymeric matrix. The molecular structure of C60 and PA 12 polymeric chain are depicted in Figure 1a and b, respectively. The uniform dispersion of carbon nanofillers in the polymer matrix is a general requirement for accomplishing the desired mechanical and physical properties. Thus, in the present analysis, the reinforcements C60 are assumed to be periodically and homogenously distributed in the PA 12 matrix material. Given the previous assumption, it becomes evident that the nanocomposite may be analysed according to the ABCDEFGH cubic RVE which is illustrated in Figure 2. Note, that small amounts of fullerenes are commonly added into a polymeric matrix in order to prevent agglomerate phenomena. Moreover, it is proved that carbon based nanofillers in the range of 3%-5% by weight promise substantial improvement of mechanical properties and may achieve the same reinforcement as 20%-30% of microsized fillers [3]. Thus, rather small loadings of C60 reinforcement up to 8% by weight are studied here. However, small concentrations of nanofillers inevitably lead to large RVEs, which require substantial computational effort when atomistic methods such as MD are to be followed exclusively. On the other hand, the mechanical properties of the nanocomposites are significantly influenced by interfacial interactions between the nanofillers and the polymer matrix. Therefore, the accurate representation of interfacial effects is essential and may only be achieved via a numerical analysis at an atomistic level.
Given the above, a multiscale approach is proposed here, which simulates the nanocomposite in two distinct stages. In the first one, the small central unit cell abcdefgh depicted in Figure 2 is formulated via MD. The specific unit cell is constructed in such a way so that it contains a centrally positioned single C60 crystal at a high concentration of 25 wt%. The arisen MD model created by using "Materials Studio 2017" [17] is illustrated in Figure 3. The MD model is constructed by using the "Amorphous Cell Calculation" module and by adopting a cubic type of lattice at a room temperature of 298 K. The potential energy of the system is assumed to be described by the COMPASS force field [18], the first and only ab initio force field that enables an accurate and simultaneous prediction of various gasphase and condensed-phase properties of organic and inorganic materials.
The lattice parameter a of the cubic unit cell may be estimated each time via the following equation: where the VC60 and VP12 are the volumes of the C60 reinforcement molecule and the PA 12 matrix material given, respectively, by: where the DC60 is the van der Waals diameter of the C60 molecule taken equal to 1.1nm [19], mC60 is the mass of the C60, i.e., 1.1956×10 -21 kg, ρPA12 is the density of the PA 12 polymer typically taken equal to 900kg/m 3 and, finally, wC60 denotes the weight percentage (wt%) of the C60 in the unit cell.
To determine the size of the cubic unit cell simulated by MD, the C60 weight percentage (wt%) is set to wC60 = 25. In this way, a lattice parameter equal to a = 1.673nm is arisen via Eqs. (1) to (3). The stable configuration of the cubic MD model (Figure 3) is achieved by using the specific side length and by performing energy minimization via the steepest descent algorithm available. Next, the calculation of the necessary mechanical properties is accomplished by applying small strains via the "Forcite" module. The simulation is performed by utilizing a canonical ensemble (constant temperature, number of atoms, and volume). Finally, the mechanical properties of the unit cell are estimated via its lamé parameters λ, μ which are computed by the MD simulations and by assuming that the nanocomposite unit cell behaves like an isotropic medium. Then, the Young's modulus EUC and the Poisson's ratio νUC of the unit cell are calculated, respectively, via the following well-known equations: The above described numerical procedure leads to a Young's modulus equal to Euc = 4.186GPa and a Poisson's ratio equal to νuc = 0.21 for the PA 12/C60 unit cell with wC60 = 25 of Figure 3.
In the second stage of the analysis, the arisen mechanical property MD based data (Euc = 4.186GPa and νuc = 0.21) regarding the abcdefgh unit cell with mC60 = 25 are utilized to develop a continuum same sized cubic volume within the whole ABCDEFGH RVE of the nanocomposite with wC60 = 2, 4, 6, 8. Note that Eq. (1) is appropriate for the estimation of the size of the whole RVE as well. According to the followed micromechanical analysis and due to the symmetry, only the one quarter of the problem domain may be formulated, which is denoted by LMNOKFIJ in Figure 2 and 4a, given that appropriate symmetry conditions are adopted. An orthogonal Cartesian coordinate system is used as reference with x, y and z axes aligned with the main dimensions of the one quarter of the RVE. Both the central PA 12/ C60 unit cell volume as well as the outer PA 12 volume (Figure 4a) are discretised with isoparametric, hexahedral, quadratic, twenty-noded finite elements having three degrees of freedom per node (displacements ux, uy, uz) [20] as Figure 4b illustrates.
Since the interfacial phenomena are already incorporated in the first stage atomistic analysis, the connection between the central unit cell volume and the outer PA 12 volume is obviously considered perfect. Regarding the material properties of the outer polymeric-only volume, it is assumed that the PA 12 is isotropic with a Young's modulus equal to EPA12 = 2.420GPa and a Poisson's ratio equal to νPA12 = 0.39 [12]. In order to determine the elastic properties of the nanocomposite, appropriate boundary conditions must be applied. For the calculation of the Young's modulus ERVE and Poisson's ratio νRVE, a small uniform strain εy = 0.01 is applied on the boundary LKFM. In addition, the constraints ux = 0, uy = 0, and uz = 0 are applied on the boundaries OLMN, OJIN, and OJKL whereas the boundaries JKFI and NIFM are kept parallel to their original shape by nodal coupling since shear stresses on these boundaries must be zero due to the symmetry. The Young's modulus of the nanocomposite ERVE is computed via the ratio of average stress, obtained from the sum of reactions in the ground y = 0, to the applied strain εy = 0.01. The Poisson's ratio νRVE is computed through the ratio of the arisen average transverse strain, to the applied normal strain.

Numerical Results
Small filler loadings of 2, 4, 6, and 8 wt% are numerically tested. Figure 5a and b depict the contours of the equivalent von Mises stresses and strains for the filler loading case of 2% by weight. The stress concentrations within the reinforced unit cell become obvious while strains tend to increase within the outer PA 12 volume and especially on the deformed boundary above the reinforcement. The variation of the nanocomposite Young's modulus ERVE versus the C60 weight percentage (%) is given in Figure 6a. The numerical outcome shows that the addition of C60 nanoparticles within PA 12 by 8% by weight may lead to a stiffness improvement up to 17%. Some corresponding experimental results are as well illustrated in Figure 6a. Very good agreement between the numerical results and those measured via appropriate experimental procedures [12] may be observed only for the filler loading cases 2 and 6 wt%. Note that the experimental measurements [12] show a distinct stiffness loss when the reinforcement concentration is above 4 wt%. This is perhaps due to the non-bonded complicated interactions between neighboring C60 molecules which get very close to each other under high concentrations. Unfortunately, since the adopted technique is based on the use of a unit cell containing a single isolated C60, it seems that such phenomena may not be incorporated by the proposed FEM-only models (Figure 4). The dependence of the nanocomposite Poisson's ratio νRVE on the percentage mass fraction of the C60 reinforcement is illustrated in Figure 6b. Generally, the stiffness of the nanocomposite increases linearly as the percentage mass fraction of the nanofiller increases. On the other hand, the Poisson's ratio of the nanocomposite is reduced as the filler mass fraction gets higher values due to the small Poisson' ratio of C60. It has to be noted here that a denser mesh than that depicted in Figure 4b, has a negligible effect on the FEM numerical solutions.

Conclusions
A three dimensional, multiscale micromechanical numerical method has been formulated in order to predict the elastic mechanical behaviour of PA 12 polymer which is uniformly reinforced with C60 nanoparticles at small loading levels.
The complicated interfacial phenomena between the nanofiller and the polymeric matrix require analysis at an atomistic level while small nano-reinforcement mass fractions demand continuum numerical treatment to reduce the enormous computational cost.
Therefore, a two stage analysis is proposed according which an MD simulation is firstly performed to extract the mechanical behaviour of a PA 12/C60 unit cell of a very high filler mass fraction, and secondly an appropriate volume element, which contains a same sized unit cell and, thus, represents the whole nanocomposite, is simulated in a continuum manner via FEM. Elastic material parameters such as Young's modulus and Poisson's ratio have been computed for several loading levels up to 8% by weight. The results have shown that the increase of C60 mass fraction leads to a significant increase of the nanocomposite stiffness but also to a decrease of its Poisson's ratio.