Heat Transfer at the Interface of Graphene Nanoribbons with Different Relative Orientations and Gaps

: Because of their high thermal conductivity, graphene nanoribbons (GNRs) can be employed as ﬁllers to enhance the thermal transfer properties of composite materials, such as polymer-based ones. However, when the ﬁller loading is higher than the geometric percolation threshold, the interfacial thermal resistance between adjacent GNRs may signiﬁcantly limit the overall thermal transfer through a network of ﬁllers. In this article, reverse non-equilibrium molecular dynamics is used to investigate the impact of the relative orientation (i.e., horizontal and vertical overlap, interplanar spacing and angular displacement) of couples of GNRs on their interfacial thermal resistance. Based on the simulation results, we propose an empirical correlation between the thermal resistance at the interface of adjacent GNRs and their main geometrical parameters, namely the normalized projected overlap and average interplanar spacing. The reported correlation can be beneﬁcial for speeding up bottom-up approaches to the multiscale analysis of the thermal properties of composite materials, particularly when thermally conductive ﬁllers create percolating pathways.

Since the thermal conductivity (λ) of pristine polymers is generally low (0.11-0.44 W/m·K [20]), it would be expected that adding carbon nanofillers with high λ (100-5000 W/m·K [21,22]) would greatly enhance the effective thermal conductivity of the resulting mixture. However, empirical and numerical studies have shown that PNCs with carbon nanofillers may present λ well below the value predicted by effective medium approximations [23][24][25][26]. In fact, the presence of defects (e.g., atom substitutions, atomic vacancies, Stone-Wales dislocations) [27], edge chirality [28], dumping due to the surrounding polymer [29], and interfacial thermal resistance (R k , also known as Kapitza resistance) at the filler-matrix and filler-filler interfaces significantly hinder the heat transfer through the network of fillers [30,31]. In particular, R k originates from the mismatch of phonon spectra at the filler-filler and filler-matrix interfaces, which causes phonon scattering and thus a bottleneck for the heat transfer across the interfaces [32][33][34]. The resulting additional interfacial thermal resistances may thus lead to a dramatic reduction in the effective λ of PNCs. The considerable impact of R k on the effective λ of PNCs has stimulated the investigation of new approaches to reducing thermal resistances at interfaces. For example, the presence of edge-grafted molecular junctions or cross-linkers between graphene nanofillers was experimentally and numerically found to reduce the thermal resistance at the filler-filler interface [35,36]. Furthermore, the controlled chemical functionalization of carbon nanofillers improved the thermal transport across filler-matrix interfaces [37][38][39][40], while reducing the thermal conductivity of the nanofillers at the same time [41].
On the one hand, when the filler volume fraction is below percolation threshold, fillers are unlikely to create large interconnected chains with high thermal conductivity [42]; therefore, the R k at the filler-matrix interface is one of the primary factors determining the thermal performance of the composite. On the other hand, with filler concentrations larger than the percolation threshold, an interconnected and highly conductive network of fillers is typically formed [42]. In this case, it is important to know how relative orientation of neighbouring fillers affects heat transfer at their interface. In fact, above the percolation threshold, R k at filler-filler interfaces has a fundamental role in determining the effective thermal conductivity of PNCs, since direct contact between contiguous fillers takes place [43]. Fasano et al. have studied the effect of geometry and filler functionalization on the overall thermal transmittance of networks of carbon nanotube fillers [44]. In particular, they noticed that the thermal resistance at the filler-filler interface strongly decays with both the overlapping area and concentration of covalent cross-linkers between adjacent carbon nanotubes. A similar behaviour has recently been observed also in the case of networks of Graphene Nanoribbons (GNRs) [36].
In the present study, a novel correlation between the relative orientation of adjacent isolated graphene nanoribbons and the interfacial thermal resistance at their interface has been found by molecular dynamics simulations. This correlation allows to compute the interfacial thermal resistance between two graphene nanoribbons at different relative orientations, namely overlap along the x, y, z directions and rotation around the x, y, z axes. The correlation obtained could be a valuable tool for a quick prediction of the R k between graphene nanofillers, therefore simplifying the multiscale simulation of the effective λ of polymeric composites including percolating networks of graphene nanoribbons.

Methods
The effect of the relative orientation of contiguous GNRs on their interfacial thermal resistance has been studied in this work by classical molecular dynamics simulations. A setup including three GNRs (see Figure 1a) was considered as a minimal building block of more complex networks of nanofillers. Each GNR had width W = 24 Å and length L = 200 Å, and contained 1968 atoms. Carbon-carbon bonded interactions were modelled by the adaptive intermolecular reactive empirical bond order (AIREBO) potential [45] implemented in the LAMMPS molecular dynamics package [46], whereas, non bonded van der Waals interactions between the GNRs were modelled by 12-6 Lennard-Jones potential. All simulation setups were initially energy minimized and equilibrated for 200 ps in the canonical ensemble (300 K, 0.5 fs time step). Afterwards, Muller-Plathe's algorithm [47] was adopted to impose a fixed heat flux through the system, from the central (red area in Figure 1a) to the lateral (blue areas in Figure 1a) regions of the setup, while the remaining atoms were kept in the micro-canonical ensemble [48]. This reverse non-equilibrium molecular dynamics simulation was performed for 2 ns, to achieve a steady state temperature profile and heat flux. Further details on the simulation protocol are reported elsewhere [36].
As depicted in Figure 1, the horizontal overlap (a, along x axis), vertical overlap (b, along y axis), interplanar spacing (h, along z axis), and angular displacement ( θ = (θ x , θ y , θ z )) between adjacent GNRs were varied to assess the correlation between R k and the network geometry. The configuration with 3 GNRs depicted in Figure 1a was employed to explore the effect of a, b, h, θ x , and θ z variations on R k , while a setup consisting of 2 GNRs (the central and the rightmost one in Figure 1a) was considered to test different θ y , since θ y = 0 • breaks the symmetry of the 3 GNR setup. The positive value of the angular displacements considered is highlighted by the green arrows in Figure 1b. Figure 1. Networks of graphene nanoribbons (GNRs) studied by reverse non-equilibrium molecular dynamics simulations. (a) Configuration adopted to investigate the effect of horizontal overlap (a, along x axis), vertical overlap (b, along y axis) and interplanar spacing (h, along z axis) on the R k between GNRs (W width, L length). In the inset, the definition of b is highlighted. Muller-Plathe's algorithm is used to assure a constant heat flux from the region of the setup with higher temperature to the lower ones. (b) Configurations adopted to study the effect of angular displacement (θ x , θ y , θ z ) on the average R k between GNRs. The tested configurations have been chosen according to the one-at-a-time sensitivity approach [49].

Results
The nanoscale gap in the overlapping area between GNRs limits the phononic heat transfer through the network, therefore reducing its overall thermal transmittance. In fact, in the configuration in Figure 1, the weak van der Waals forces between the GNRs are the sole cause of their interaction; as a consequence, heat transfer is limited by an additional thermal resistance at the interface between the GNRs.
The interfacial thermal resistance between a pair of neighbouring GNRs can be computed as: where ∆T is the average temperature jump across the interface (see Figure 2) caused by the specific heat flux (q x ) through the system. Note that, due to the configuration considered, the heat flow can be approximated as one-directional along the x axis. The effect of each geometrical parameter has been explored one at a time (one-at-a-time sensitivity approach [49]), namely by varying one parameter in a specific range while keeping the other ones fixed at the reference value. A detailed list of configurations tested is reported in Table A1. Note that the configurations explored in the sensitivity analyses have been selected with the aim of avoiding both thermal contact (e.g., R k 10 −10 m 2 K/W) and insulation (e.g., R k 10 −7 m 2 K/W) between the GNRs (see the configurations N. 1-20 in Table A1). Consequently, the configurations presenting no horizontal and vertical overlap have not been further considered, since they typically show R k > 10 −7 m 2 K/W (see, for instance, the configurations N. 21 and 22 in Table A1). Firstly, the horizontal and vertical overlaps between the GNRs have been varied. The impact of horizontal overlap, a, is simulated in the range 20-80 Å, with b fixed at 0 Å, h at 4 Å and angular displacement set to θ = (0,0,0). As depicted in Figure 3a and previously reported [36], R k decreases with larger horizontal overlap following a power law, because this allows enhanced van der Waals interactions between GNRs and thus better phononic heat transfer [44]. A similar behaviour is observed for varying values of vertical overlap (b tested in the range 0-18 Å, with a fixed at 40 Å, h at 4 Å and θ = (0,0,0)), since less overlapping area (i.e., higher b) leads to a power increase of R k (see Figure 3b). The impact of interplanar spacing between the GNRs has been instead studied with h in the range 2.5-6 Å, while a is fixed at 40 Å, b at 0 Å and θ = (0,0,0). Since van der Waals interactions tend to drastically decrease with distance (12-6 Lennard-Jones potential), the results in Figure 3c confirm a power-law-like increase of R k with h, in good accordance with previous works [36].
Secondly, the effect of angular displacement between graphene nanoribbons on R k has been studied around the x, y and z axes (see Figure 1b). The central GNR is rotated around the x axis with θ x ranging in the interval 0-7 • , while the other geometrical parameters are kept constant (a = 40 Å, b = 0 Å, h = 4 Å, θ y = θ z = 0 • ). The results in Figure 3d show no significant variation of R k with θ x . This could be due to the fact that, in the reference system considered, rotations around the x axis cause half of the central GNR to approach the lateral GNRs (i.e., increased heat transfer), and the other half to move away from it (i.e., reduced heat transfer). These competing phenomena may balance each other, therefore leading to an approximately constant R k with θ x . Conversely, rotations around the y axis (θ y from 0 • to 2 • ; a = 40 Å, b = 0 Å, h = 4 Å, θ x = θ z = 0 • ) induce larger R k (Figure 3e), because of the increasing average distance between GNRs and thus the reduced phononic heat transfer through van der Waals interactions. Finally, the central sheet is rotated around the z axis in the 0-45 • range, while keeping the other geometrical parameters fixed as a = 40 Å, b = 0 Å, h = 4 Å and θ x = θ y = 0 • . The results in Figure 3f show a strong increase of R k with θ z , because of the progressively reduced overlapping area between the GNRs and thus phononic heat transfer at the interface. a b c d e f Inspired by the approach typically employed for classic heat transfer phenomena, the main geometrical parameters influencing the R k at the GNR-GNR interface have been grouped into meaningful dimensionless quantities. In particular, we suggest the following set of possible dimensionless quantities: the normalized interfacial thermal resistance namely the ratio between interfacial (GNR-GNR) and internal (single GNR) thermal resistance; the normalized projected overlap that is the ratio between the projected overlapping area and the GNR size; the normalized average interplanar spacing namely the ratio between the average interplanar spacing between the adjacent GNRs and the distance of their closest possible approach at equilibrium conditions (van der Waals radius). Because of their definitions, A n ∈ [0, 1] and h n ∈ [0, ∞). In the simulated setups, λ = 96 W/m·K and σ = 3.4 Å [36]. In the reference system considered, A is estimated by projecting along the z axis the overlapping area between two adjacent GNRs (see, for example, the green surfaces in Figure 4a). The average interplanar spacing, instead, is computed by averaging the interplanar spacing (h) between two adjacent GNRs in the overlapping area (see, for example, the green lines in Figure 4b). The complete list of R * k , A n and h n values for the simulation setups considered is reported in Table A1. To achieve a quantitative correlation between R * k and the dimensionless quantities suggested, the Artificial Intelligence-powered modelling engine by Eureqa software [50,51] has been employed. In detail, Eureqa employs genetic algorithms to identify the mathematical equation f that best match a set of data, that is R * k = f (A n , h n ) in this case. Similarly to the typical empirical expressions adopted for heat transfer, only basic formula building-blocks (i.e., constant, multiplication, division, power) have been considered for the Eureqa fitting. As a result, the simulation results in Table A1 (configurations N. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] have been found to be best fitted (R 2 = 0.99, see Figure 4b) by the following expression: being k 0 = 4.460 × 10 −1 , k 1 = 2.646 and k 2 = 1.491. Coherently with the simulation evidence depicted in Figure 3, the interfacial thermal resistance increases with the average interplanar spacing, while larger overlapping areas lead to enhanced heat transfer between adjacent GNRs (i.e., lower R * k ). Notably, R * k is more sensitive to h n variations compared to A n ones. The full list of R * k values obtained by Molecular Dynamics (MD) simulations (red dots) and predicted by Equation (5) (black lines) is reported graphically in Figure 3 and numerically in Table A1. It is worth noting that Equation (5) is valid for A n > 0. However, configurations characterized by no overlapping area (A n = 0) typically show R k > 10 −7 m 2 K/W and thus the GNRs can be considered as thermally insulated from each other (see, for instance, the configurations N. 21 and 22 in Table A1): this would eventually avoid the creation of a percolating network through the polymer matrix in which they are included.

Discussion
If GNRs are employed as fillers to improve the properties of polymeric matrix, the estimation of R k at the interface between the GNRs is helpful to compute the effective thermal conductivity (λ e f f ) of the resulting composite. Considering, for example, the modified Maxwell-Garnett effective medium approximation by Shahil and Balandin [52], the λ e f f can be computed as where: λ p and λ m are the thermal conductivities of filler particles and polymeric matrix; φ is the volume fraction of fillers; R B is the interfacial thermal resistance at the filler-matrix interface; H is the thickness of fillers. Since graphene fillers tend to disperse non-homogeneously within the polymeric matrix and to create aggregates or networks with different sizes [53,54], λ p should be computed for the graphene aggregates rather than for the graphene nanoparticles alone. For instance, considering the setup in Figure 1a as a representative aggregate of GNRs, the equivalent λ p of the network can be computed by a one-dimensional approximation of heat conduction along the main axis (x, in this case), that is where a series of lumped thermal resistances due to the N GNRs and M GNR-GNR interfaces is taken into account. In Equation (7), the following notations are adopted: L p , overall length of the fillers network along x axis; L i and λ i , respectively the length and thermal conductivity of the i-th GNR filler along x axis; R k,j , the interfacial thermal resistance at the j-th filler-filler interface, which could be readily computed by Equation (5) given the relative orientations of GNRs. Clearly, more accurate estimates of λ p could be achieved by three-dimensional models of heat conduction through the aggregate, for example employing effective medium approximations [55], shape factors [56], mesoscopic [42] or continuum [57] simulation models, which would anyway require an accurate estimation of R k,j .

Conclusions
In this study, the effect of relative orientation of isolated adjacent graphene nanoribbons on the thermal resistance at their interface has been systematically studied by reverse non-equilibrium molecular dynamics. A simulation setup consisting of three GNRs was considered as a minimal and relevant building block of more complex networks of nanofillers, such as those that can be found in polymer nanocomposites above the percolation threshold. The analysis included different geometrical characteristics of the network, such as the GNR-GNR overlap along the x, y, z directions and rotation around the x, y, z axes. Simulation results allowed an empirical and comprehensive correlation between the aforementioned geometrical parameters and R k to be determined for thermally percolating networks of GNRs, namely for R k 10 −7 m 2 K/W. This correlation was found to accurately predict the simulated interfacial thermal resistance between contiguous GNRs, at least in the geometrical ranges considered. The proposed correlation could be utilized for the bottom-up prediction and optimization of the effective thermal transmittance of graphene nanofiller networks [58]. Future studies should focus on assessing the possible effect of irregular GNR shapes on the R k at the filler-filler interface [59], and determining similar correlations for the R k at the filler-matrix interface taking into account the effect of matrix and filler functionalizations. In a broader context, this study may facilitate the design of carbon reinforced polymer nanocomposites with adjustable/enhanced thermal properties.  Acknowledgments: Computational resources were provided by HPC@POLITO (http://www.hpc.polito.it). We thank Eliodoro Chiavazzo and Pietro Asinari for useful discussion on data interpretation.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

GNRs
Graphene Nanoribbons PNCs Polymer Nanocomposites MD Molecular Dynamics AIREBO Adaptive Intermolecular Reactive Empirical Bond Order LAMMPS Large-scale Atomic/Molecular Massively Parallel Simulator Appendix A Table A1. List of Molecular Dynamics (MD) and fitting (Equation (5)) results of R k at the GNR-GNR interface, as illustrated in the setup in Figure 1. Notice that the correlation in Equation (5)