Simulation of Elastic and Fatigue Properties of Epoxy/SiO<sub>2</sub> Particle Composites through Molecular Dynamics

Gaoge Zhao; Jianzhong Chen; Yong Lv; Xiaoyu Zhang; Li Huang

doi:10.32604/cmes.2021.015388

Simulation of Elastic and Fatigue Properties of Epoxy/SiO2 Particle Composites through Molecular Dynamics

Hubei Key Laboratory of Theory and Application of Advanced Materials Mechanics, School of Science, Wuhan University of Technology, Wuhan, 430070, China
*Corresponding Author: Yong Lv. Email: lvyonghl@163.com
Received: 15 December 2020; Accepted: 08 April 2021

Abstract: The influence of different nanoparticle sizes on the elastic modulus and the fatigue properties of epoxy/SiO2 nanocomposite is studied in this paper. Here, the cross-linked epoxy resins formed by the combination of DGEBA and 1,3-phenylenediamine are used as the matrix phase, and spherical SiO2 particles are used as the reinforcement phase. In order to simulate the elastic modulus and long-term performance of the composite material at room temperature, the simulated temperature is set to 298 K and the mass fraction of SiO2 particles is set to 28.9%. The applied strain rate is 109/s during the simulation of the elastic modulus. The results show that the elastic modulus of the material increases with the increase in particle size. Furthermore, fatigue simulation under strain control is performed on the model with SiO2 nanoparticle radius of 12 Å. The results indicate that the influence trend of variable frequencies on the fatigue mechanical response is similar, and the mean stress decreases with the increase in number of cycles. In addition, the smaller the loading period and the more the number of cycles, the greater the mean stress reduction. Finally, the change in energy and free volume fraction are evaluated under fatigue loading condition.

Keywords: SiO2 particle; elastic modulus; fatigue; free volume; molecular dynamics

Composite materials are used widely in the fields of aerospace, transportation, construction, machinery and industry [1–3] because of their advantages of higher strength-to-weight ratio, longer durability, and better environmental resistance compared with metal materials. Particle reinforced composite material is characterized by discontinuous particles. Its stiffness and strength are slightly lower while the material isotropy is better compared with continuous fiber reinforced material. Therefore, particle reinforced composite is still widely valued. As is known, the elastic modulus and fatigue performance are crucial indicators for particle reinforced composite.

In previous studies, the evaluation method for elastic modulus and fatigue was mainly experimental method. Some scholars have studied the effect of particle reinforcement on material elastic modulus. Su et al. [4] investigated the elastic modulus of bioinspired 3D printed two-phase composites, and found that elastic modulus generally decreases as unit cell size decreases for discontinuous phase composite. Li et al. [5] investigated the interfacial property of carbon fiber/epoxy composites and found that the interfacial property is enhanced by addition of silica nanoparticles. Also, some researchers have investigated the effects of particle reinforcement on fatigue property of composite materials. Zeltmann et al. [6] studied the fatigue property of glass fiber reinforced epoxy enhanced by nanoparticles and found that composites with small amounts of nanoparticles added into the matrix have bending fatigue strength similar to that obtained for the neat glass fiber reinforced epoxy matrix composite. As is known, the properties of materials are also closely related to the molecular structure. With the continuous improvement of computer performance, the emergence of various supercomputing platforms has greatly promoted the development of molecular dynamics simulation. Molecular dynamics is a general simulation method based on the Newtonian equations and force functions, which provides the opportunity to analyze the mechanical properties of materials on the atomic scale. Some mature theoretical systems of molecular dynamics simulation have been established. Molecular dynamics simulation has been widely used to calculate and analyze the stress-strain and failure behaviors of materials [7–9]. The development of molecular dynamics allows easy comparison of the elastic modulus of different types of materials. Adnan et al. [10] studied the influence of filler size on elastic properties of nanoparticle reinforced polymer composites using molecular dynamics simulations. Their results showed that the elastic properties of nanocomposites are significantly enhanced with the reduction of inclusion size. Similarly, the fatigue performance of materials can be simulated on the atomic scale. Some fatigue simulations [11,12] have been performed on pure polymers. In these previous studies, the creep-fatigue relationship of polymers was researched using tensile-tensile fatigue under stress-controlled conditions, and it was found that fatigue and creep cause significant changes in the van der Waals and dihedral potential energies. Besides, the fatigue behavior of amorphous polyethylene under strain-controlled conditions has been studied as well and the results show that stress relaxation is enhanced by cyclic loading.

Most of the previous works have mainly studied the fatigue property of polymers or the mechanical properties (such as elastic modulus, Poisson’s ratio, etc.) of single inclusion composite materials. There is limited research on the fatigue mechanism of composite materials and mechanical properties of randomly distributed nanoparticle materials. In this article, the mechanical properties of randomly distributed nanoparticle materials were investigated. Therefore, three models were established with different particle radii with the same mass fraction. The particle radii were 5, 8, and 12 Å, respectively. Furthermore, the fatigue performance simulation of the model with the nanoparticle radius equal of 12 Å under strain control was carried out through molecular dynamics method. Specifically, cross-linked epoxy resin was used as the composite material matrix phase and spherical SiO2 was used as the reinforcement phase.

There are some well-established force fields that can be selected to simulate the molecular system. For example, condensed-phase optimized molecular potentials for atomistic simulation studies (compass) force field, polymer consistent force field (pcff) and consistent valence force field (cvff) [13]. Yaphary et al. [14] studied the adhesion of epoxy-silica interface in salt environment using cvff and the results showed that sodium chloride solution weakens the adhesion significantly. In this work, cvff was chosen to simulate the elastic modulus and fatigue properties of nano reinforcement composite materials.

The initial model was constructed through the commercial software Material Studio (MS) 7.0, as it has a self-avoidance walking algorithm which is useful for construct a better model than others. After the initial model was constructed, the Large-scale Atomic/Molecular Massively Parallel Simulator (lammps) [15] to run molecular dynamics simulation due to the costlier computation time of MS software. 1,3-Benzenediamine was used as the curing agent and DGEBA was used as the epoxy resin in this paper, and the molecular models are shown in Figs. 1a and 1b, respectively. The reaction atoms of curing agent and epoxy resin are shown in red, named as active atom. Each curing agent can connect four epoxy resins in theory. One curing agent is employed to connect three epoxy resins as representative units of the cross-link model considering that the degree of cross-linking cannot reach 100%. The model is shown in Fig. 1c. This modeling method was applied in Yu et al. [16] and Guru et al. [17] paper and the results obtained are basically consistent with the conventional cross-link methods. Also, this method is convenient for the construction of nano reinforcement model. Figs. 1d–1f show the spherical particles with various radii. A model with low initial density was established first. However, after relaxation, the equilibrium model obtained by this model has large voids, and the structure is not reasonable. This is possibly because the initial density is not high enough. After some tests, it was found that initial density of 1.7 g/cm3 could solve the problem of voids in the equilibrium model. Then, the initial model density was set to 1.7 g/cm3 which could establish a more compact structure. The SiO2 mass fraction was set to 28.9% in this paper in order to compare the performance of different sizes of SiO2 with the same mass fraction. This mass fraction may be slightly higher than the proportion of SiO2 added in practical applications. However, it is suitable for the study of the size effect, it was worthwhile to apply. If the mass fraction of SiO2 particles is too small, it may result in only one SiO2 particle being placed when building a model with larger particles in a certain size box, which is not good for comparing the performance of random particle distribution models. Here, three models were built with 28.9% SiO2 mass fraction and nanoparticle radius of 5, 8, and 12 Å, respectively. Further increase in the particle size will result in fewer particles or even a single inclusion model, so it is not discussed herein. Figs. 2–4 present the initial model of SiO2 particles with radius of 12, 8, and 5 Å, respectively. The specific components of each model are provided in Tab. 1. The mass fraction is controlled by adjusting the number of particles and cross-linked epoxy. After the initial models were built in MS 7.0, they were converted to the lammps executable file under cvff for the molecular dynamics simulations.

Figure 1: The molecular models of (a) 1,3-benzenediamine (b) DGEBA (c) The representative unit of cross-linked epoxy (d) Spherical SiO2 (radius of 5 Å) (e) Spherical SiO2 (radius of 8 Å) (f) Spherical SiO2 (radius of 12 Å)

The simulation was divided into two parts: the first part was the simulation of elastic modulus of nanoparticles with different radii, including three models with particle radius of 5, 8 and 12 Å, respectively; and the second part was the simulation of the fatigue performance of the model with the particle radius equal to 12 Å. The second part was performed to compare the effects of different loading frequencies on the fatigue performance of Nano-reinforcement composite materials. Therefore, only one of these three models was selected for fatigue simulation and analysis.

A Nose–Hoover thermostat and barostat were used to control the temperature and pressure. Cut-off distance for long-range interactions was 12 Å. The entire simulation process usually includes three steps: minimizing energy, balancing the system and performance simulation. The detailed description of the three steps is given below.

(1) Energy minimizations were performed on the initial structure before running the MD simulations by applying the steepest descent algorithm which could optimize the structure to the greatest extent. Then, the conjugate gradient algorithm was employed for further optimization;

(2) The isothermal isobaric ensemble, which is abbreviated as NPT (Constant Number of particles, Pressure and Temperature) ensemble, was used to improve the structure to reach a balanced state, which is suitable for starting simulations.

During the balance process, the pressure in the three directions (x, y, z) was set to zero atmosphere. First, Gaussian distribution was used to give all particles a random speed to make the temperature of the model 100 K. Then, the system was linearly heated from 100 to 500 K through 100 picosecond (ps) simulation and simulated at 500 K for 300 ps with a time step equal to 0.5 fs. Finally, the model was cooled down to 298 K with the cooling rate equal to 1 K/ps, which was aimed to simulate temperature and relaxation for 1 ns with a time step equal to 1.0 fs at the temperature to reach the equilibrium state. After 1 ns relaxation, the densities of the models 1, 2 and 3 tend to be linear and equal 1.27, 1.28 and 1.28 g/cm3, respectively.

The final construction of equilibrium state of particle radius equal 12 Å has shown in Fig. 5 which is output from Open Visualization Tool (Ovito) [18] software, and the other models are not shown here.

The stress tensor for atom is given by the following formula, where a and b take on values x, y, z to generate the components of the tensor in different directions:

where the first term is a kinetic energy contribution for atom I, and the second term is the virial contribution which is shown in Eq. (2).

Wab=12∑n=1Np(r1aF1b+r2aF2b)+12∑n=1Nb(r1aF1b+r2aF2b)+13∑n=1Na(r1aF1b+r2aF2b+r3aF3b)+14∑n=1Nd(r1aF1b+r2aF2b+r3aF3b+r4aF4b)+14∑n=1Ni(r1aF1b+r2aF2b+r3aF3b+r4aF4b)+Kspace(ria,Fib) (2)

Figure 5: The construction of equilibrium state after 1 ns relaxation with SiO2 particle radius of 12 Å

The first term is a pairwise energy contribution where n loops over the Np neighbors of atom I, r1 and are the positions of the 2 atoms in the pairwise interaction, F1 and F2 are the forces on the 2 atoms resulting from the pairwise interaction. The next four terms are the contributions of similar form for the bond, angle, dihedral, and improper interactions, respectively. The last term is the K-Space contribution from long-range Coulombic interactions.

The calculation process of elastic modulus of three different nanoparticle size models are described below. A 5 ps equilibrium simulation was added to obtain another equilibrium conformation after the model was in equilibrium, and a further 5 ps simulation was added to obtain another equilibrium conformation in order to avoid accidental errors caused by one simulation. Finally, three slightly different equilibrium conformations were used for simulation and we take the average of the results was taken for the elastic modulus in the loading directions. The other two directions maintained a state of zero atmospheric pressure when the strain loading was applied in x direction (similarly when the loading direction is y or z). The applied strain rate was 109/s during the uniaxial tensile simulation, and take the first 3% deformation results were taken for linear fitting. The slope of fitting results represents the elastic modulus in the tensile direction.

In general, fatigue can be divided into tension-compression fatigue, compression-compression fatigue and tension-tension fatigue. The tensile-tensile fatigue was applied in the simulations.

In this simulation process, the mean strain was set to 5%, and the fatigue simulation was carried out with a period equal to 1, 2, and 4 ps with a strain range of 4.5%–5.5%. The model was subjected to 5% tensile strain in x direction at a strain rate 109/s before cyclic loading.

One 400 ps simulation was performed on three different frequencies loading conditions, which have mentioned above in order to facilitate the comparison of the relationship between different loading frequencies. Therefore, that 400 cycles were simulated with a period of 1 ps, 200 cycles were simulated with a period of 2 ps, and 100 cycles were simulated with a period of 4 ps.

The loading diagram is shown in Fig. 6. Black represents period of 1 ps of the cyclic loading. Similarly, red represents period of 2 ps of the cyclic loading, and blue represents 4 ps of cyclic loading. The initial linear loading is not shown here.

As shown in Tab. 2, the elastic modulus values of the individual are the average of three results. The model should be isotropic in theory, so the average of the three-direction modulus was taken as the elastic modulus of the entire model, considering that there will be deviations in the molecular models calculated in the respective directions. The results show that as the radius of the nanoparticles increases, the elastic modulus of the composite material also increases when the mass fraction is 28.9%.

The diameter of the largest particle packing was only 2.4 nm in this simulation and the added nanoparticles could reach 100 nm and even larger in practical applications. It is reasonable to predict that the elastic modulus of the composite material gradually increase as the particle size increases within a certain range. It will be possible to use molecular simulations to calculate the size of the actually added nanoparticles on the atomic scale with the further development and application of computer computing power. The simulations can be regarded as virtual experiments to guide the design of materials for a more accurate prediction of the performance while saving time and resources.

In the fatigue loading, the strain loading was a standard sinusoidal curve as shown in Fig. 4. Only one model (with SiO2 radius of 12 Å) was selected for fatigue simulation and analysis as mentioned in beginning of Section 2.2. Here, the mean stress and mean strain were used to study the trend between various frequencies and the average value in one period was used as the mean stress and mean strain [9].

The relationship between the mean stress of various frequencies and the number of cycles is shown in Fig. 7. Three different loading frequencies were compared in the same period considering that time is also an important reason for the stiffness reduction of fatigue cycle. Here, we take the cycle process with a maximum period of 4 ps as a benchmark, and compare the two cycles with a period of 2 ps as one cycle. Similarly, the four cycles with a period of 1 ps are compared to one cycle. A logarithmic model, which is defined in Eq. (3), is used to fit the results.

The stress level is the same stress level at the beginning of fatigue loading duo to these three loading methods use the same model. Here, the stress value is 198 MPa as the ‘a’ value of fitting equation when tensile strain is equal to 5%. The specific fitting equation is shown in Tab. 3.

It can be qualitatively inferred from the fitted curve that as the cycle period decreases, the mean stress reduces faster. This can also be quantitatively seen from the value of ‘b’ in the equation, as the smallest cycle period has the largest ‘b’ value.

The variation of van der Waals energy is shown in Fig. 8. The van der Waals energy with a period of 1 ps is significantly higher than the other two results, and the van der Waals energy with a period of 2 ps has a slight increase compared to the period of 4 ps. The magnitude of van der Waals energy is related to the distance between molecules. The increase of van der Waals force implies an increase in the distance between molecules in the entire simulation system. The results also show that the molecular motion is faster when the period is smaller, and it is more difficult for the molecular structure to reach a relative equilibrium state between molecules.

The mean stress reduction is greater with smaller period as expected. The smaller period has the larger value, which is the opposite of the stress reduction considering only the van der Waals energy. Hence, there should be more significant contributions from other energies. The contribution of bonded energy is mainly derived from bond tensile energy and bond angle energy. In addition, there are also contributions from dihedral angle and improper dihedral angle energy. However, these two terms have very small contributions compared with the bond tensile and bond angle energy. Therefore, only bond tensile energy and bond angle energy were analyzed here.

The result of bond tensile energy is shown in Fig. 9. It can be seen that the first several periods lead to greatly reduced mean stress when the period is 1 ps, and over all, the bond energy is also significantly smaller than other two results. This is the main reason for the faster reduction of the mean stress. It can also be seen that the bond energy of period of 2 ps is slightly smaller than that of 4 ps. The result of bond angle energy is shown in Fig. 10. It can be seen that the period of 2 ps is lower than 4 ps, and it shows a decreasing trend as the number of cycles increases. The difference in bond angle energy is the main reason why cause the mean stress of period equal to 2 ps is lower than 4 ps.

3.2.3 The Free Volume Analysis under Fatigue State of Nano-Reinforcement Material

In general, the change in molecular energy is closely related to the molecular structure and distribution. The free volume change during the entire simulation process is analyzed in below.

Free volume is an index for calculating the volume occupied by particles and it is based on cutting the entire box into small square grids with a size of 0.3 Å, which is enough to obtain accurate results. This grid is considered to be occupied within the van der Waals radius of the specified particle, or it is regarded as grid free and not occupied. The free volume results were processed by VMD [19] software to output higher resolution graphics after calculation by open access software Multiwfn [20]. The free volume in the case of loading with a period of 1 ps is shown here. The left image in Fig. 11 shows the distribution of particles and free volume. The blue color represents free volume and the other colors represent different types of atoms. The periodic display of x, y, and z directions is added on this basis since there is also a lot of free volume inside the model that cannot be displayed. The distribution of free volume is clearly shown in the right image of Fig. 11. The terms as same as Fig. 11 after fatigue loading at period of 1 ps is shown in Fig. 12. Although the result graphics are given with sufficient accuracy, the model is a porous structure with a large free volume and the free volume after cyclic loading will not change significantly. It is expressed in a specific numerical form as shown in the Tab. 4, the free volume when linearly loaded to 5% is 35.95%, and the free volume after 400 ps cycles is 35.89% when the period is 1ps. When the period is 2 and 4 ps, the final free volume is 35.64%.

Figure 11: The free volume at initial strain of 5% (the free space is shown by blue, the left is the simulation model, the right has applied the periodic at three directions based on left)

The free volume is reduced by 0.32% when the period is 2 and 4 ps, and the free volume is only reduced by 0.07% when the period is 1 ps. The free volume reduction rate for periods of 2 and 4 ps is 4.57 times compared with periods of 1 ps. The result of this free volume change is in agreement with the result of Van der Waals energy, indicating that the smaller the period, the more difficult it is to recover the relative distance between molecules under fatigue loading. At the same time, the distance between molecules will decrease due to the relaxation effect in the loading. This result implies that the fatigue effect is higher when the period is 1ps, and the relaxation effect is not obvious. However, when the period is 2 and 4 ps, the free volume is the same and there just a slight difference in van der Waals energy between them, which indicates that the relaxation effect plays a leading role. We will continue to study this issue in depth to determine their inner connection.

Figure 12: The free volume after fatigue loading at period of 1 ps (same as Fig. 11)

In this paper, the mechanical properties of epoxy/SiO2 particle composites with different sizes of spherical SiO2 and fatigue properties with SiO2 particle radius of 12 Å were studied using molecular dynamics, and the following conclusions were obtained:

1. As the size of the nanoparticles increases, the elastic modulus of the material increases when the mass fraction of SiO2 particle is 28.9%. This result is different from that of the single inclusion [13] model. This is because smaller SiO2 particles can agglomerate more easily, resulting in a weaker interface binding capacity with epoxy matrix, while there is no particle agglomeration phenomenon in the single inclusion model.

2. Throughout the fatigue simulation process, the van der Waals energy shows a clear downward trend, which means that the molecular gap becomes smaller as expected. At the same time, the shorter the cycle period, the higher the van der Waals energy, which means that the distance between the molecules is larger.

3. When the frequency in the fatigue loading process is higher, it is harder to achieve relative equilibrium between the molecules, and the mean stress decreases faster, which is mainly due to the lower bond energy. The main reason for the lower mean stress with a cycle period of 2 ps compared to 4 ps is the change in bond angle energy.

4. There is a relaxation effect in the fatigue process of strain control, which is similar to the creep effect of fatigue under stress control. This can also be seen from the change in free volume. The smaller the period, the less the free volume decreases. Hence, the effect of relaxation is less obvious.

Funding Statement: The project supported by the Fundamental Research Funds for the Central Universities (WUT: 2018IB001) and the Fundamental Research Funds for the Central Universities (WUT: 2019III130CG).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Rana, S., Fangueiro, R. (2016). Advanced composite materials for aerospace engineering. UK: Woodhead Publishing. [Google Scholar]

2. Ma, X. C., He, G. Q., He, D. H., Chen, C. S., Hu, Z. F. (2008). Sliding wear behavior of copper-graphite composite material for use in maglev transportation system. Wear, 265(7–8), 1087–1092. DOI 10.1016/j.wear.2008.02.015. [Google Scholar] [CrossRef]

3. Komus, A., Beley, N. (2018). 3.15 Composite applications for ground transportation. Comprehensive Composite Materials II, 3, 420–438. DOI 10.1016/B978-0-12-803581-8.09950-1. [Google Scholar] [CrossRef]

4. Su, F. Y., Sabet, F. A., Tang, K., Garner, S., Pang, S. et al. (2020). Scale and size effects on the mechanical properties of bioinspired 3D printed two-phase composites. Journal of Materials Research and Technology, 9(6), 14944–14960. DOI 10.1016/j.jmrt.2020.10.052. [Google Scholar] [CrossRef]

5. Li, L. Y., Liu, L. Q., Tang, L. C., Gao, Y., Zhang, Z. (2012). Using silica nanoparticles reinforce interfacial property of carbon fiber/epoxy composites. Journal of Materials Engineering, 2(6), 32–36. DOI 10.5923/j.ijme.20120204.01. [Google Scholar] [CrossRef]

6. Zeltmann, S. E., Prakash, K. A., Doddamani, M., Gupta, N. (2017). Prediction of modulus at various strain rates from dynamic mechanical analysis data for polymer matrix composites. Composites Part B: Engineering, 120(1), 27–34. DOI 10.1016/j.compositesb.2017.03.062. [Google Scholar] [CrossRef]

7. Pedone, A., Malavasi, G., Menziani, M. C., Segre, U., Cormack, A. N. (2008). Molecular dynamics studies of stress–strain behavior of silica glass under a tensile load. Chemistry of Materials, 20(13), 4356–4366. DOI 10.1021/cm800413v. [Google Scholar] [CrossRef]

8. Liu, C., Ning, W., Tam, L. H., Yu, Z. (2020). Understanding fracture behavior of epoxy-based polymer using molecular dynamics simulation. Journal of Molecular Graphics and Modelling, 101(14–15), 107757. DOI 10.1016/j.jmgm.2020.107757. [Google Scholar] [CrossRef]

9. Linga, G., Ballone, P., Hansen, A. (2015). Creep rupture of fiber bundles: A molecular dynamics investigation. Physical Review E, 92(2), 22405. DOI 10.1103/PhysRevE.92.022405. [Google Scholar] [CrossRef]

10. Adnan, A., Sun, C. T., Mahfuz, H. (2007). A molecular dynamics simulation study to investigate the effect of filler size on elastic properties of polymer nanocomposites. Composites Science and Technology, 67(3–4), 348–356. DOI 10.1016/j.compscitech.2006.09.015. [Google Scholar] [CrossRef]

11. Sahputra, I. H., Echtermeyer, A. T. (2015). Creep-fatigue relationship in polymer: Molecular dynamics simulations approach. Macromolecular Theory and Simulations, 24(1), 65–73. DOI 10.1002/mats.201400041. [Google Scholar] [CrossRef]

12. Sahputra, I. H., Echtermeyer, A. T. (2014). Molecular dynamics simulations of strain-controlled fatigue behaviour of amorphous polyethylene. Journal of Polymer Research, 21(11), 1–13. DOI 10.1007/s10965-014-0577-2. [Google Scholar] [CrossRef]

13. Gaedt, K., Hltje, H. (1998). Consistent valence force-field parameterization of bond lengths and angles with quantum chemical ab initio methods applied to some heterocyclic dopamine D3-receptor agonists. Journal of Computational Chemistry, 19(8), 935–946. DOI 10.1002/(SICI)1096-987X(199806)19:8%3C935::AID-JCC12%3E3.0.CO;2-6. [Google Scholar] [CrossRef]

14. Yaphary, Y. L., Yu, Z., Lam, R. H., Hui, D., Lau, D. (2017). Molecular dynamics simulations on adhesion of epoxy-silica interface in salt environment. Composites Part B: Engineering, 131(1), 165–172. DOI 10.1016/j.compositesb.2017.07.038. [Google Scholar] [CrossRef]

15. Plimpton, S. (1995). Fast parallel algorithms for short-range molecular dynamics. Journal of Computational Physics, 117(1), 1–19. DOI 10.1006/jcph.1995.1039. [Google Scholar] [CrossRef]

16. Yu, S., Yang, S., Cho, M. (2009). Multi-scale modeling of cross-linked epoxy nanocomposites. Polymer, 50(3), 945–952. DOI 10.1016/j.polymer.2008.11.054. [Google Scholar] [CrossRef]

17. Guru, K., Sharma, T., Shukla, K. K., Mishra, S. B. (2016). Effect of interface on the elastic modulus of CNT nanocomposites. Journal of Nanomechanics and Micromechanics, 6(3), 4016004. DOI 10.1061/(ASCE)NM.2153-5477.0000109. [Google Scholar] [CrossRef]

18. Stukowski, A. (2009). Visualization and analysis of atomistic simulation data with OVITO-the open visualization tool. Modelling and Simulation in Materials Science and Engineering, 18(1), 15012. DOI 10.1088/0965-0393/18/1/015012. [Google Scholar] [CrossRef]

19. Humphrey, W., Dalke, A., Schulten, K. (1996). VMD: Visual molecular dynamics. Journal of Molecular Graphics, 14(1), 33–38. DOI 10.1016/0263-7855(96)00018-5. [Google Scholar] [CrossRef]

20. Lu, T., Chen, F. (2012). Multiwfn: A multifunctional wavefunction analyzer. Journal of Computational Chemistry, 33(5), 580–592. DOI 10.1002/jcc.22885. [Google Scholar] [CrossRef]