Nanoscale thermal properties of carbon nanotubes/epoxy composites by atomistic simulations

Carbon nanotubes/epoxy composites are increasingly employed in several industrial fields, because of the enhanced material properties provided by the nanofillers. In particular, the thermal conductivity of these nanocomposites is determined by heat transfer mechanisms occurring over multiple scales, thus causing a complex relation between effective response and microscopic characteristics of the material. Here, the thermal properties of epoxy composites reinforced by carbon nanotubes are investigated using atomistic simulations. For a better understanding of how the effective thermal conductivity arises from the characteristics of the composite at the nanoscale, the thermal properties of its constituents are studied separately according to different geometrical, physical and chemical characteristics. The thermal conductivity of carbon nanotubes and epoxy resin alone is first investigated by molecular dynamics; then, the Kapitza resistance at the nanotube–nanotube and nanotube–epoxy interfaces is studied as well. The effective thermal conductivity of the carbon nanotubes/epoxy composite is finally computed and the observed behavior interpreted on the basis of the properties of the nanofillers, matrix and interfaces alone. Results – verified against effective medium theory predictions – show that, for the considered configurations, the effective thermal conductivity of the nanocomposite increases with the nanotube length and volume fraction, with the curing degree of the epoxy and system temperature. In perspective, the presented approach could be employed to investigate other constitutive materials or properties of nanocomposites.


Introduction
Polymers are currently employed in a multitude of engineering and industrial applications because of their lightweight, low cost, chemical and mechanical stability, and many other application-specific properties. Due to their good mechanical properties and ability to withstand environmental degradation, epoxy resins are particularly suitable for various applications in the aerospace and automotive industry. However, epoxies typically show thermal conductivities ( ) between 0.1 and 0.2 W m −1 K −1 [1], thus making them ill-suited for applications where efficient heat transfer is required instead [2,3]. While the introduction of fibers into the polymer matrix provides a substantial enhancement of its mechanical properties [4], higher filler loadings are usually required to observe an appreciable increase in the effective thermal conductivity, which lead to significant processing issues [5]. Hence, the recent research has focused on the possibility to reinforce epoxy matrices with ultra-high thermal conductive nanofillers as well, in order to improve the thermal performance of the resulting composite [2,[5][6][7][8][9][10]. Among different reinforcements for polymer nanocomposites (PNCs), in ; whereas, Huxtable et al. [19] a 15% enhancement at 1 wt%. Surprisingly, in some cases it has been observed that introducing carbon fillers reduces the overall of nanocomposites with respect to the matrix alone [20].
The moderate increase in the effective thermal conductivity of PNCs reinforced by highly conductive nanofillers may have different reasons. First, carbon nanofillers tend to agglomerate in bundles when introduced in the polymer matrix, due to the Van der Waals interactions between them [21]. This segregation of fillers could have a detrimental effect on the of the composite [22]. Second, Kapitza resistance ( ) could be a severe bottleneck in the heat transfer through PNCs [23]. Because in these composites heat transfer is mainly ruled by phonon transfer, such interfacial thermal resistance is caused by the phonon scattering at filler-matrix and filler-filler interfaces [24,25], which hinders the heat flux through the composite thus reducing its effective . Different factors affect the magnitude of , for instance the concentration of defects, vacancies or functionalizations in the fillers or matrix, as well as the reciprocal orientation between the dispersed fillers [26][27][28][29].
Since the length scale of the interfacial effects responsible for Kapitza resistance is in the order of few angstroms, the conventional experimental techniques and continuum modeling are unable to characterize them. Instead, atomistic simulations such as Molecular Dynamics (MD) can be used to unveil how different interfaces regulate the thermal transport across nanocomposites [30][31][32]. Here, the most challenging part is to generate a realistic molecular model of the interfacial regions, which can be then employed to quantify with different geometrical, physical and chemical characteristics of the PNC [33][34][35].
The present study, for the first time, investigates by atomistic simulations the thermal conduction properties of composites based on DGEBA-DETA (diglycidyl ether of bisphenol A-diethylenetriamine) epoxy reinforced by carbon nanotubes. For a better understanding of how the effective thermal conductivity arises from the characteristics of these composites at the nanoscale, we have divided our research into different parts (see a scheme of this approach in Fig. 1). First, the thermal conductivity of carbon nanotubes ( ) and epoxy ( ) alone are studied and validated against the literature, according to different length and curing degree, respectively. Second, the Kapitza resistance at the nanotube-nanotube ( , ) and nanotube-epoxy ( , ) interface is computed with different nanotube orientations, length and curing degree, and validated by literature results. The influence of the polymer matrix on the thermal conductivity of the immersed nanofillers is assessed as well. Third, the effective thermal conductivity of the resulting nanocomposite ( ) is assessed and compared with Effective Medium Theory (EMT) [36] predictions. Finally, a sensitivity analysis is carried out to quantify how the thermal properties of the composite could be adjusted by changing some geometrical, physical and chemical characteristics of its constituents.

Theoretical background
Molecular dynamics. Reverse Non-Equilibrium Molecular Dynamics (RNEMD) based on the Müller-Plathe's approach [37] has been employed to compute the thermal conductivity of carbon nanotubes and epoxy alone. In this approach, the simulation box is divided into bins, say, along the axis, where the terminal bins '1' and ' ' are specified as the regions at lower temperature (blue regions in Fig. 2) and the middle bin '( ∕2) + 1' the one at higher temperature (red region in Fig. 2). Heat flux is then imposed on the system as a primary perturbation. In the RNEMD method, energy is continuously transferred, in an artificial manner, from the cold regions to the hot one by exchanging the velocity of the hottest atom ( ℎ ) in the cold region with the one of the coldest atom ( ) in the hot region [37]. This implies a more rapid convergence of the simulation with respect to the non-equilibrium MD (NEMD) approach, where the thermal gradient is imposed and the heat flux measured instead. Indeed, the simulated heat flux typically presents large fluctuations and thus slower convergence; whereas, the thermal gradient, being averaged on time and space, tends to converge in a faster way [37]. The RNEMD process eventually imposes an artificial heat flux ( ) from the cold to the hot region [38], which can be computed as: where is the simulation time, is the cross-sectional area perpendicular to the heat flux direction and is the mass of atoms, with the summation done over all the velocity exchange processes. Since the total energy is conserved during the exchange events, the amount of energy per unit time and cross-sectional area exchanged from the hot to the cold region via heat conduction is the same as . As a result, generates a temperature gradient ⟨ ∕ ⟩ throughout the system, whose can be then calculated using the Fourier's law [37]: Clearly, the heat flux follows the binning direction and, thus, the computed refers to that specific direction. The RNEMD approach is also adopted for computing the Kapitza resistance at the filler-filler interfaces. At steady state condition, the computed temperature drop ( ) and heat flux ( ) across the CNT-CNT interface allows to obtain [33,[39][40][41] = . ( Instead, the Kapitza resistance between CNT and epoxy is here measured using a transient MD simulation. This transient method represents a good solution when one phase (polymer matrix) has a considerably lower thermal conductivity with respect to the other (CNT). In these cases, the Kapitza resistance is directly related to the matrix-filler interface area, the filler heat capacity, and the relaxation time of the temperature decay [31,42,43]. In this approach [31], the whole computational box is first thermalized to an equilibrium temperature (for example, = 300 K). Once the system is fully equilibrated, the CNT is heated up for 200 ps to a temperature equal to + (0), where is the previous equilibration temperature and (0) = 100 K. During  the CNT heating, the temperature of epoxy is kept constant at the equilibration using a second thermostat. After that, the whole system is allowed to relax in a micro-canonical ensemble (NVE). The filler and epoxy temperatures are recorded during the relaxation process. After a characteristic relaxation time, thermal equilibration is achieved in the whole simulation box. The difference between the average CNT temperature and the epoxy one ( ) can be fitted by the Newton's law of cooling: ( ) = (0) exp(− ∕ ), being the relaxation time constant [31]. Since the interfacial is expected to be the predominant resistance in the simulated heat transfer process [42], the fitted over the MD trajectories is related to the filler surface area ( ), its heat capacity ( ) and the , at the CNT-epoxy interface as , = . (4) Finally, the effective thermal conductivity of nanocomposite has been evaluated by a Non-Equilibrium Molecular Dynamics approach (NEMD). Indeed, in heterogeneous systems such as nanocomposites, the interchange of atoms' velocity between the hot and cold regions could involve atoms of different phases (e.g. CNTs and epoxy), which are characterized by different vibrational modes and may lead to an inhomogeneous heat flux throughout the computational domain and thus biased results. Then, to generate a homogeneous heat flux, we used a NEMD approach in the micro-canonical ensemble, with Langevin thermostats coupled separately to the epoxy and CNTs atoms of each cold and hot region [31]. In this method, the system is firstly equilibrated in the isothermal-isobaric ensemble (NPT) at 300 K and 1 bar. Then, the simulation box is divided into bins (see Fig. 2): the central region is heated up to a higher temperature (e.g. 400 K); the two extreme regions are cooled down to a lower temperature (e.g. 200 K); the remaining bins are kept in the micro-canonical ensemble (NVE). The thermostated regions are coupled to Langevin thermostats, which allow a direct evaluation of the heat exchange across the domain.
Effective medium theory. Theoretical models based on the effective medium theory [36,44] allow estimating the thermal conductivity of composites as a function of the properties of the matrix and fillers. However, most of these models consider perfect filler-matrix interfaces (i.e. no interfacial thermal resistance) and idealized morphology (i.e. uniform distribution of fillers). Therefore, EMT models provide a quick approximation of the properties of nanocomposites under idealized conditions. In this work, the Maxwell-Garnett model [36] adapted by Nan et al. [45] is used to verify simulation results. This model includes the effects of interfacial thermal resistance at filler-matrix interfaces ( , ), as well as of filler shape, size and distribution. In the case of aligned continuous fibers, the values of are evaluated depending on the direction of the heat flux with respect to the main axis of the fibers. Considering a transverse heat flux with respect to the fillers orientation, the thermal conductivity is evaluated as: where and are the thermal conductivity of matrix and filler, respectively, and being = 3 ∕ 1 the filler aspect ratio, 3 and 1 the radii of two perpendicular axes of the filler, and = , the Kapitza radius. Instead, considering a heat flux parallel to the fillers, the effective thermal conductivity of the composite along that direction is estimated as Evidence from the literature has proved that this EMT model can predict accurately the of composites at low ; whereas, predictions are no more accurate (i.e. more than 50% deviation from experimental results) at high volume fractions (e.g. > 0.3) [46].

Computational details
In the present work, the considered nanofillers are single-walled carbon nanotubes ( Fig. 3(a)) with different length and diameter (and thus volume fraction of nanocomposite), while the epoxy resin and curing agent are DGEBA and DETA, respectively ( Fig. 3(b)). DGEBA is a bi-functional reactant with two epoxide groups at the ends, while DETA has five reactive sites including both primary and secondary amine groups, hence it is a multifunctional (fivefold-functional) reactant. Therefore, DGEBA and DETA are able to produce 3D cross-linked epoxy polymers: each DETA molecule can react at most with five DGEBA molecules, each of which is capable of being connected to another DETA molecule through its second epoxide head. Thus, the normal composition ratio of DGEBA/DETA in the blend is 5:2 [47,48], and this has been considered in the simulations. The curing reaction takes place between the epoxy groups of DGEBA and the amine groups of DETA to form a cross-linked structure, whose curing degree determines the physical properties of the polymer [49]. Several cross-linking procedures have been proposed in the literature to simulate epoxy resin structures characterized by physical properties adherent with the experimental ones [50][51][52]. Here, a multi-step cross-linking approach has been employed [53], in which every reactive pair within a suitable spatial criteria is bonded iteratively [54]. The partially cross-linked epoxy system is relaxed by energy minimization followed by a short MD simulation in the canonical ensemble. This cross-linking process is repeated until the desired curing degree (or maximum reaction cutoff) is reached. If not differently specified, the systems simulated in this work present a cross-linking degree of 68%. It is well known that the physical properties of epoxy depend on the degree of cross-linking, with better mechanical and heat transfer performance obtainable at higher cross-linking degrees. However, from a numerical standpoint, it is computationally expensive to perform multi-step cross-linking beyond 70% cross-linking degrees. Also, high degrees of cross-linking result in excessively strained epoxy, which may be difficult to equilibrate in numerical computations [47,55]. Therefore, we have chosen a 68% cross-linking degree for most of the computations as a compromise between the accuracy of material properties and the computational cost.
All the presented simulations have been performed using the LAMMPS molecular dynamics package [56,57]. The COMPASS [58] and Tersoff [59] force fields have been used for the epoxy and CNTs, respectively. Lorentz-Berthelot mixing rules have been considered for the Lennard-Jones interaction potentials between epoxy and nanotubes. All the setups have been initially energy minimized and then equilibrated in the canonical ensemble at 300 K for up to 300 ps (depending on the system size), with a time step of 1 fs unless specified. Afterward, thermal properties are computed by means of RNEMD, NEMD or transient MD, as discussed previously. To guarantee reliable statistics, production runs are continued up to 2 ns after that steady temperatures and heat flux are achieved in the simulated systems. Further computational details and simulation methods are discussed in each specific subsection of Results. The input and output files of all the presented simulation types are archived in the Zenodo repository [60].

Results and discussion
Molecular dynamics has been employed to study the thermal properties of the components of carbon nanotubes/epoxy composites, namely carbon nanotube fillers, epoxy matrix and their interface. Then, the resulting nanocomposite is simulated as well. The obtained results are presented and discussed separately in the following sub-sections.

Thermal properties of carbon nanotubes
Determining how the thermal properties of CNTs are influenced by their structural characteristics is important for predicting the effective thermal conductivity of nanocomposites. Since previous works have extensively investigated the thermal behavior of CNTs [38,40,61], here the relation between CNT length and thermal conductivity is studied by RNEMD approach with the only aim to verify our simulation setup against literature evidence. A schematic about the employed configuration is reported in Fig. 4(a), where the nanotube is placed in an empty box. In detail, the thermal conductivity of single-walled carbon nanotubes with length ( ) ranging from 5 to 20 nm and fixed chirality (5,5) is computed. Results in Fig. 5(a) show that increases with due to the presence of phonons with longer wavelengths, following a power-like trend has been previously investigated in the literature [61], and it is not studied here again.
Further, the Kapitza resistance between adjacent carbon nanotubes ( , ) without any surrounding media is computed using a triple CNTs setup, as shown in Fig. 4(b). In this setup, the three CNTs belong to the left ( ), right ( ) and middle ( ) part of a CNT network, which is considered as a representative building block of the more complex dispersion of nanotubes in PNCs. In Fig. 4(b), and ℎ represent the horizontal overlap and vertical normal distance between between contiguous nanotubes. The distance between overlapped CNTs (ℎ) has been chosen as 4 Å, which is closely related to the equilibrium distance between graphene layers in graphite or multi-walled CNTs [62,63] and short enough to allow non-bonded interactions between the CNTs and thus heat transfer.
each pair of CNTs. The calculated Kapitza resistance at the CNT-CNT interface strongly depends on the distance between the overlapping CNTs (parameter ℎ in Fig. 4(b)). The chosen value for ℎ, i.e. ℎ = 4 Å, is closely related to the equilibrium distance between graphene layers in graphite or CNT walls in multi-walled CNTs (which is about 3.5 Å) [62,63]. The considered short distance between the fillers could be achieved in case of either CNT structures percolating throughout the epoxy matrix, or nanocomposites with high CNT loadings. The impact of horizontal overlap on

Thermal properties of epoxy matrix
The thermal properties of the pristine epoxy resin have been computed for a cubic simulation box (7.8 nm x 7.8 nm x 7.8 nm) containing 820 DGEBA and 328 DETA molecules and different curing degrees ( = 60-90%). Note that, in experiments, the curing degree could be controlled by manipulating the conditions at which the polymerization process occurs, such as temperature, time, pressure and presence of catalyst. Here, we first investigate the effect of on the density of polymer ( ). To this purpose, the cured epoxy resin is simulated for 1 ns at ambient temperature (300 K) and pressure (1 bar) in the isothermal-isobaric ensemble. As shown in Fig. 6(a), the density of polymer increases with the cross-linking percentage, due to the higher entanglement of epoxy chains. Simulations show that ranges from 1.13 to 1.15 g cm −3 at different cross-linking degree, and this is comparable with the available experimental evidence for DGEBA-DETA epoxy (i.e. 1.16 g cm −3 [64]). Results in Fig. 6(a) are best fitted (R 2 = 0.93) by = 5 6 , where 5 = 0.92 g cm −3 and 6 = 0.05. Further, the RNEMD method is employed to determine the thermal conductivity of DGEBA-DETA epoxy ( ) for different curing degree. As depicted in the schematic in Fig. 2, a heat flux is imposed by  Müller-Plathe algorithm between the terminal and mid bins of the simulated epoxy domain. The computed increases with the crosslinking degree (see Fig. 6(b)) due to the denser epoxy matrix, and such behavior can be fitted by = 7 8 + 9 with 7 = 3.23 × 10 −7 W m −1 K −1 , 8 = 2.67 and 9 = 0.14 W m −1 K −1 optimized coefficients (R 2 = 0.95). The measured values of are found to be in the range 0.165-0.205 W m −1 K −1 for = 60-90%, which are in good agreement with the literature evidence, where = 0.150-0.250 W m −1 K −1 [65][66][67]. In particular, Licari [68] reported that the experimental value of for DGEBA cured with DETA hardener is 0.2 W m −1 K −1 .

Thermal properties of carbon nanotubes/epoxy interface
The Kapitza resistance between the nanotube fillers and epoxy matrix ( , ) has been computed by transient MD. In Fig. 7 we report the temperature profile at the CNT-epoxy interface at different time steps, and the average temperature evolution in time of both the matrix and filler.
Firstly, we evaluate the interface thermal resistance at different cross-linking density, as shown in Fig. 8(a). In these simulations, a 7.8 nm long (5,5) CNT filler is immersed in a 3. with the curing degree can be interpreted by the better affinity between the two phases due to polymerization, as well as by the higher density of epoxy matrix at larger (see Fig. 6(a)), which could improve the CNT-epoxy contact at their interface. Higher curing degrees lead to improved structural properties of the resin in terms of stiffness and density. Eslami et al. [69,70] showed that such higher density allows a better organization of the polymer interface layers, thus reducing the interface thermal resistance. Varshney et al. [32] reported a 20% increase in the interface thermal conductance with the curing degree due to the enhancement of the matrix stiffness, since it reduced the differences in the vibrational modes between the two phases.
Then, we assess the effect of nanofiller length on , , considering (5,5) CNTs with = 4, 6 and 7.8 nm length immersed in an epoxy domain (3.4 nm x 3.4 nm x 7.8 nm). An increase of , with the CNT length has been observed (see Fig. 8(b)), in agreement with the trend reported by Gardea et al. [71]. These results can be fitted by the exponential function The influence of temperature on , is also explored and reported in Fig. 9(a). Again, these simulations consider a 7.8 nm long (5,5) CNT filler immersed in a 3.4 nm x 3.4 nm x 7.8 nm box of epoxy resin. Results show that , decreases linearly with higher temperatures, possibly due to a better matching between the vibrational characteristics of the two phases [72,73]. Such behavior follows a trend , = 16 + 17 , being 16 = −8.93 × 10 −10 m 2 W −1 and 17 = 3.95 × 10 7 m 2 K W −1 (R 2 = 0.98). Similar relations between , and temperature has been also reported in previous works [74,75].
Finally, we investigate whether the presence of epoxy matrix has an effect on the thermal properties of CNT fillers or not. To quantify the magnitude of this effect, the thermal conductivity of a (5,5) CNT placed in an empty box ( ) or immersed in an epoxy matrix ( − ) has been evaluated by RNEMD simulation. The normalized thermal conductivity of the CNT ( * = − ∕ ) is plotted in Fig. 9(b) for . It can be observed that * = 0.7 in average, i.e. the thermal conductivity of CNTs surrounded by epoxy matrix decreases by 30% with respect to the relative value in vacuum, as a result of interfacial interactions between nanotube and polymer atoms. As previously seen in literature, one of the most important issues that hinders the thermal performance of carbon fibers in composite materials is in fact the damping effect of the surrounding polymer matrix on carbon fibers, which limits the thermal conductivity enhancement provided by highly conductive fillers in nanocomposites. Mortazavi et al. [31] performed combined molecular dynamics-finite element modeling of thermal conduction in graphene-epoxy nanocomposites, reporting that the thermal conductivity of single layer graphene declines by around 30% because of interfacial interactions with epoxy atoms. Zhang et al. [28] also concluded that simulations indicate that the polymer molecules wrapped around CNT walls have a negative influence on the thermal conductivity of the nanotube along the longitudinal direction.

Effective thermal conductivity of nanocomposites
Once the thermal properties of each component (and their interfaces) of carbon nanotubes/epoxy composites have been explored, the effective thermal conductivity of the resulting PNC is finally investigated. To this aim, a cubic simulation box with 7.8 nm edge filled with DGEBA-DETA and reinforced by five (5,5) SWCNTs is initially considered (Fig. 10(a)), therefore achieving a volume fraction of fillers equal to = 0.03. In this configuration, the total number of atoms is 48750, the length of the CNT fillers is 7.8 nm and periodic boundary conditions are applied along all three Cartesian directions. A time step of 0.5 fs is used for all simulations, to guarantee a better numerical convergence. A specific computational protocol has been followed for simulating PNCs. The five CNTs are first introduced in the simulation domain, which is then filled by DGEBA and DETA molecules at the considered ratio. The domain is energy minimized and relaxed for 1 ns in the isothermal-isobaric ensemble with Nosé-Hoover thermostat (300 K, 1 bar). Once the steady state is eventually achieved, the crosslinking procedure is performed through successive equilibration runs in the isothermal-isobaric ensemble. When the target density of epoxy is reached, the system is finally relaxed in a micro-canonical ensemble for 100 ps. NEMD simulations can be then performed on this fully relaxed configuration to evaluate the effective thermal conductivity of the nanocomposite.
Both the along the direction of CNT axis ( direction in Fig. 10(a)) and the one perpendicular to it ( direction in Fig. 10(a)) have been evaluated. The mean temperature profile throughout the binned computational domain and the measured during the run along either or direction are reported in Fig. 11. These computations have been done at = 68%. Figs. 11(b) and 11(d) show that the measured thermal conductivity of nanocomposite in the parallel and perpendicular direction with respect to the CNT axis are equal to , = 1.67 W m −1 K −1 and , = 0.32 W m −1 K −1 , respectively. On the one side, , is nine times larger than , demonstrating the potential thermal effectiveness of percolating CNTs at low concentration ( = 0.03) already. On the other side, , shows only a 1.5-fold enhancement with respect to the thermal conductivity of the pristine epoxy matrix. In fact, comparing the temperature profiles in   the and directions (Figs. 11(a) and 11(c)), pronounced temperature discontinuities can be noticed when the heat flows along the direction perpendicular to the CNT: such temperature jumps are provided by the Kapitza resistance at the filler-matrix interfaces, which leads to the lower value of , . Hence, a significant anisotropy is evident for the thermal conduction properties of PNCs with continuous aligned fibers: in the case of , , the positive contribution of CNTs to heat conduction is predominant; while, in the case of , , the Kapitza resistance at the filler-matrix interfaces affects dramatically the effective heat transfer.
Simulation results have been compared with the thermal conductivity of nanocomposite predicted by the EMT-based model for continuous aligned fibers in Eqs. (5) and (7), considering the damping effect of epoxy on the thermal conductivity of nanotubes (see Fig. 9(b)). Results in Figs. 11(b) and 11(d) show that both the , and the , measured by MD are verified against the EMT predictions, since discrepancies lower than 7% can be noticed.
To achieve a better understanding on the main factors determining the effective thermal conductivity of carbon nanotubes/epoxy composites, some geometrical, chemical or physical parameters of the PNC have been varied and the response evaluated. First, the possible effect of simulation box size on the computed , is checked. The smaller computational domain made of epoxy resin reinforced by one (5,5) carbon nanotube shown in Fig. 10(b) is built to keep the same volume fraction as the previous case (i.e. = 0.03). This reduced computational domain has 3.4 nm x 3.4 nm x 7.8 nm size, being the nanotube 7.8 nm long. A NEMD simulation with = 68% and = 300 K is carried out. The measured , for this smaller system is equal to 1.60 W m −1 K −1 , that is less than 1% lower than the value computed for the larger system depicted in Fig. 10(a). This evidence confirms that the effective thermal conductivity along the direction is not significantly affected by the size of the simulation box in the and directions. Hence, smaller simulation domains along and directions could be safely employed when computing , . Therefore, to save computational time, the effect of different crossliking degrees of the epoxy matrix, namely = 0%, 68% and 80%, on the thermal conductivity of the PNC along the main axis of the CNTs ( , ) is evaluated on the smaller computational domain ( Fig. 10(b)). Results in Fig. 12(a) show that , increases with . This positive correlation is due to both the enhanced thermal conductivity of polymer matrix (see Fig. 6(b)) and reduced Kapitza resistance at the filler-matrix interface (see Fig. 8(a)) at increasing curing degree. However, all the systems present higher values of thermal conductivity if compared to the pristine epoxy, demonstrating that the CNT fillers play a dominant role in the heat transfer mechanism with respect to the epoxy resin.
The effect of temperature on , has been studied for = 250, 300 and 350 K, considering the smaller computational domain in Fig. 10(b) with = 68%. As depicted in Fig. 12(b), , tends to increase with temperature. A similar behavior has been reported also in the literature [74], and it originates from the direct proportionality between the thermal conductivity of CNTs [76] and epoxy resin [77] with temperature (at least in the considered range), as well as by the reduced Kapitza resistance at the CNT-epoxy interface (see Fig. 9(a)).   The effect of CNT length on , has been then studied considering the smaller computational domain in Fig. 10(b) with = 68%. Simulations have been done varying from 4 to 7.8 nm. Fig. 13(a) shows the enhancement in , of the reinforced systems with respect to pure epoxy, being * = , ∕ . For both the CNTs with 4 nm and 6 nm length, , is higher than the of pristine epoxy, but sensibly lower than , of the system with continuous fillers (i.e. = 7.8 nm). In fact, while , tends to increase with the length of fillers ( Fig. 8(b)), is positively affected by (Fig. 5(a)) and -most importantly -fillers shorter than the box size present two additional , at the CNTepoxy interfaces that, in the series of thermal resistances throughout the PNC, affect dramatically the overall thermal properties.
Finally, the influence of the volume fraction of CNTs on the overall thermal conductivity of nanocomposite has been explored considering a cubic computational domain (7.8 nm x 7.8 nm x 7.8 nm) made of epoxy resin ( = 68%) reinforced with two SWCNTs with = 7.8 nm. Three systems are simulated with different chirality (i.e. (5,5), (10,10) and (12,12)), which correspond to = 1.2% 4.9% and 6.5% volume fraction, respectively. Notice that the SWCNTs are not filled with epoxy. Fig. 13(b) shows that , linearly increases with , in agreement with Eq. (7) and previous study [78].

Conclusion
Epoxy-based composite materials are gaining increasing attention because of their mechanical strength, lightweight and ease of process, especially in industrial sectors such as automotive, aerospace, energy and electronics. Enhanced thermal conduction properties are of crucial importance in some applications, e.g. in thermal storage devices or heat exchangers. To develop such materials, the epoxy matrix can be reinforced by carbon-based nanofillers with superior thermal conductivity, for instance carbon nanotubes or graphene nanoribbons. The effective thermal conductivity of the resulting nanocomposite, however, is affected by the properties of matrix, fillers and their interface that occur over scales spanning from the nano to the macro ones. Such multi-scale nature of nanocomposites makes difficult to predict the macroscopic thermal response of the material starting from the molecular characteristics of its constituents.
In the present work, we have developed a original computational approach to investigate how the atomistic characteristics of nanofillers, polymer matrix and their interfaces determine the effective thermal conductivity of the resulting nanocomposite. As a representative case study, carbon nanotubes/epoxy composites are considered. Atomistic simulations have been used to assess how the geometrical, chemical and physical features of the constituents of the nanocomposite affect its thermal response. On the one side, the thermal conductivity of carbon nanotubes increases with their length; whereas, the one of DGEBA-DETA epoxy with the curing degree. Moreover, the epoxy matrix surrounding the nanotubes has a damping effect on their thermal conductivity. On the other side, the thermal boundary resistance at the nanotube-nanotube interface reduces with the overlap between nanofillers; whereas, the one at the nanotube-epoxy interface decreases with the curing degree and increases with the nanotube length and temperature. These numerical evidence have been verified against the previous literature. The effective thermal conductivity of the carbon nanotubes/epoxy composite is finally computed through molecular dynamics, and the observed behavior interpreted on the basis of the thermal properties of the nanofillers, matrix and interfaces alone.
These atomistic results may show only a partial overview of the overall thermal transfer through nanocomposites, due to the limited size of the computational domains that can be simulated by molecular dynamics in a reasonable amount of time. In fact, heat transfer phenomena occurring at higher scales (e.g. thermal percolation) cannot be studied by atomistic models alone. Nevertheless, the correlations between the molecular characteristics of the nanocomposite and its thermal behavior presented in this work are essential to develop comprehensive multi-scale models able to predict the thermal response of components made of carbon nanotubes/epoxy. As an example, the Authors of the present work have employed the reported correlations as input parameters for mesoscopic and continuum thermal models of macroscopic samples of nanocomposites, which could be finally compared against experimental results [79]. In perspective, the presented numerical approach could be employed to investigate other constitutive materials or properties of composite materials.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.