Review of Recent Developments on Using an Off-Lattice Monte Carlo Approach to Predict the Effective Thermal Conductivity of Composite Systems with Complex Structures

Here, we present a review of recent developments for an off-lattice Monte Carlo approach used to investigate the thermal transport properties of multiphase composites with complex structure. The thermal energy was quantified by a large number of randomly moving thermal walkers. Different modes of heat conduction were modeled in appropriate ways. The diffusive heat conduction in the polymer matrix was modeled with random Brownian motion of thermal walkers within the polymer, and the ballistic heat transfer within the carbon nanotubes (CNTs) was modeled by assigning infinite speed of thermal walkers in the CNTs. Three case studies were conducted to validate the developed approach, including three-phase single-walled CNTs/tungsten disulfide (WS2)/(poly(ether ether ketone) (PEEK) composites, single-walled CNT/WS2/PEEK composites with the CNTs clustered in bundles, and complex graphene/poly(methyl methacrylate) (PMMA) composites. In all cases, resistance to heat transfer due to nanoscale phenomena was also modeled. By quantitatively studying the influencing factors on the thermal transport properties of the multiphase composites, it was found that the orientation, aggregation and morphology of fillers, as well as the interfacial thermal resistance at filler-matrix interfaces would limit the transfer of heat in the composites. These quantitative findings may be applied in the design and synthesis of multiphase composites with specific thermal transport properties.


Introduction
Polymer composites can combine the merits of fillers and polymer matrix in a hybrid system. They can even replace the traditional materials (e.g., glass-based materials, metals and alloys) in a number of applications, due to their unique optical, mechanical, thermal and electrical properties [1]. During recent years, carbon-based two-phase polymer composites, such as carbon nanotube (CNT)-polymer composites and graphene-polymer composites, have attracted much attention from both academy and industry, owing to their superior strength, high thermal conductivity and good electrical conductivity [2][3][4]. In general, the thermal conductivity and electrical conductivity of carbon additives (e.g., CNTs, graphene, carbon black, carbon fiber) are both much higher than those of polymer matrices.
With the additions of carbon fillers into polymers, both the thermal and electrical conductivity of composites can be significantly enhanced, which may limit the applications of composites in some specific fields. For instance, in electronic packaging, the packaging materials require high thermal conductivity but low electrical conductivity. Two-phase composites with only one type of carbon filler are difficult to achieve this requirement. Multi-types of fillers are required to enhance thermal conductivity but retain electrical insulation for composites. For example, Dai et al. developed polyimide (PI) composites containing one-dimensional (1D) SiC nanowires and two-dimensional (2D) graphene nanosheets as hybrid fillers [5]. As the 1D SiC nanowires prevent the graphene nanosheets forming networks, the three-phase composite achieves a low electrical conductivity (~1.32ˆ10´1 0 S/cm) but a high thermal conductivity (~0.577 W/m¨K) [5].
On the other hand, as electronic devices are developing to the integrated and micro-scale level, thermal interface materials (TIMs) with high thermal conductivity are necessary to effectively dissipate the heat and then to protect electronic devices [6,7]. Traditional TIMs with one type of filler require high loading of fillers to achieve the expected thermal conductivity of 0.5-5 W/m¨K, which dramatically increases the cost and the viscosity of the composite. Recently, it has been found that combining multiple fillers in polymer may significantly enhance the thermal conductivity of the composites owing to the synergistic effect of fillers. Zhao et al. obtained a three-phase composite by adding 0.2 vol % graphene foam (GF) and 2.7 vol % multilayer graphene flakes (MGFs) into polydimethylsiloxane (PDMS). Due to the synergistic effect of GF and MGFs, the MGF/GF/PDMS composites achieved a thermal conductivity of up to 1.08 W/m¨K, which is 440%, 184% and 80% higher than those of pure PDMS, 0.2 vol % GF/PDMS and 2.7 vol % MGF/PDMS composites, respectively [8].
In addition to the above multiphase composites, other diverse types of multiphase composites have been developed for their specific applications, such as CNT/inorganic nanoparticle/polymer [9,10], CNT/graphene/polymer [11][12][13], graphene/inorganic nanoparticle/polymer [14,15], and polymer blends [16,17]. Such multiphase composites can combine the merits of fillers and matrix to achieve advanced properties for diverse applications. Experimental measurements have shown that the effective thermal conduction properties of multiphase composites can be significantly enhanced due to the synergistic effects of fillers. However, computational studies in this field are still limited due to the complex structure and composition of multiphase composites. The effective medium theory has been used to estimate the effective thermal conductivity (K eff ) of multiphase composites [18]. The K eff is calculated in separated steps, which neglects the synergistic effect of fillers. An effective modeling of heat transfer in multiphase composites can reveal the thermal transport phenomena and limitations in them, which may help to design multiphase composites with specific thermal properties.
In this article, we review an effective approach to model the heat transfer in multiphase composites with complex structures based on our work about the computational studies on the effective thermal properties of multiphase composites. The paper is organized as following: In Simulation Methods, an off-lattice Monte Carlo approach for modeling heat transfer in a three-phase composite is presented in detail. Three case studies are discussed in Results and Discussion to demonstrate the novelty, accuracy and superiority of the developed approach. The case-study systems include three-phase single-walled CNT (SWNT)/tungsten disulfide (WS 2 )/poly(ether ether ketone) (PEEK) composite, four-phase SWNT/SWNT bundle/WS 2 /PEEK composite, and complex graphene/poly(methyl methacrylate) (PMMA) composite. We quantitatively studied the factors influencing the thermal transport in multiphase composites, such as the morphology, dispersion, alignment and concentration of fillers, as well as the effects of the interfacial thermal resistance between any pair of components. Such quantitative findings may guide researchers to design and synthesize multiphase composites with specific thermal transport properties.

Simulation Methods
In this section, we chose a three-phase SWNT/WS 2 /PEEK composite to illustrate the off-lattice Monte Carlo approach. The off-lattice Monte Carlo approach regards the heat flow as the result of the movement of discrete thermal walkers with random motion [19]. The Monte Carlo approach has been successfully applied to model heat or mass transport in flow through porous materials [20,21] and in convective flows [22,23]. A three-dimensional (3D) model was built based on the SWNT/WS 2 /PEEK composite fabricated by Naffakh et al. [9]. As shown in Figure 1a, SWNTs (500 nm length and 2 nm diameter) and WS 2 spherical nanoparticles (110 nm diameter) were randomly and uniformly dispersed in PEEK matrix (925 nm side length). The mass fractions of SWNTs, WS 2 nanoparticles and PEEK were 0.5 wt. %, 0.5 wt. % and 99.0 wt. %, respectively. The 3D model is a representative volume element (RVE) of the three-phase composite, which can be repeated to accurately replicate the realistic composite samples. In this section, we chose a three-phase SWNT/WS2/PEEK composite to illustrate the off-lattice Monte Carlo approach. The off-lattice Monte Carlo approach regards the heat flow as the result of the movement of discrete thermal walkers with random motion [19]. The Monte Carlo approach has been successfully applied to model heat or mass transport in flow through porous materials [20,21] and in convective flows [22,23]. A three-dimensional (3D) model was built based on the SWNT/WS2/PEEK composite fabricated by Naffakh et al. [9]. As shown in Figure 1a, SWNTs (500 nm length and 2 nm diameter) and WS2 spherical nanoparticles (110 nm diameter) were randomly and uniformly dispersed in PEEK matrix (925 nm side length). The mass fractions of SWNTs, WS2 nanoparticles and PEEK were 0.5 wt. %, 0.5wt. % and 99.0 wt. %, respectively. The 3D model is a representative volume element (RVE) of the three-phase composite, which can be repeated to accurately replicate the realistic composite samples. is placed in the center of a PEEK cube with a side length of 925 nm, while 317 SWNTs (2 nm diameter and 500 nm length) are randomly distributed in the PEEK cube. The WS2 particle is painted red and the nanotubes are black in the figure. Constant heat flux is applied along the x direction by creating a hot surface and a cooled surface (Reproduced with permission from [18]. Copyright Elsevier, 2015); (b) a contour plot of thermal walker distribution in the center xy plane of the SWNT/WS2/PEEK model at the thermal steady state (Reproduced with permission from [24]. Copyright American Chemical Society, 2015).
To model the heat transfer in the composite, a constant heat flux was applied along x-direction by continuously releasing a large quantity (e.g., 40,000) of thermal walkers from both sides in each time step (i.e., hot walkers from x = 0, and cold walkers from the other side). The two sides can be treated as a hot and a cold surface, as presented in Figure 1a. The hot walkers and cold walkers (with negative energy) have same absolute value of energy, so the energy of the whole system is conserved. A thermal walker jumps randomly following a Brownian motion after released from the surfaces. The Brownian motion is described by position changes of thermal walkers in each direction. The position changes take values from a normal distribution with a zero mean and a standard deviation, σ, expressed as [25]: where m D is the thermal diffusivity of the matrix and t  is the time step duration of the simulation.
Thermal walkers travel randomly from the initial release surface. Once a thermal walker jumps to the interface between PEEK and a SWNT, it is allowed to either jump into the SWNT with a probability of fm-SWNT, or still remains in the PEEK matrix with a probability of (1 − fm-SWNT). The fm-SWNT is related to the interfacial thermal resistance (known as Kapitza resistance) between the PEEK matrix and SWNTs, which can be estimated from the acoustic mismatch theory, as follows [26]: is placed in the center of a PEEK cube with a side length of 925 nm, while 317 SWNTs (2 nm diameter and 500 nm length) are randomly distributed in the PEEK cube. The WS 2 particle is painted red and the nanotubes are black in the figure. Constant heat flux is applied along the x direction by creating a hot surface and a cooled surface (Reproduced with permission from [18]. Copyright Elsevier, 2015); (b) a contour plot of thermal walker distribution in the center xy plane of the SWNT/WS 2 /PEEK model at the thermal steady state (Reproduced with permission from [24]. Copyright American Chemical Society, 2015).
To model the heat transfer in the composite, a constant heat flux was applied along x-direction by continuously releasing a large quantity (e.g., 40,000) of thermal walkers from both sides in each time step (i.e., hot walkers from x = 0, and cold walkers from the other side). The two sides can be treated as a hot and a cold surface, as presented in Figure 1a. The hot walkers and cold walkers (with negative energy) have same absolute value of energy, so the energy of the whole system is conserved. A thermal walker jumps randomly following a Brownian motion after released from the surfaces. The Brownian motion is described by position changes of thermal walkers in each direction. The position changes take values from a normal distribution with a zero mean and a standard deviation, σ, expressed as [25]: where D m is the thermal diffusivity of the matrix and ∆t is the time step duration of the simulation. Thermal walkers travel randomly from the initial release surface. Once a thermal walker jumps to the interface between PEEK and a SWNT, it is allowed to either jump into the SWNT with a probability of f m-SWNT , or still remains in the PEEK matrix with a probability of (1´f m-SWNT ). The f m-SWNT is related to the interfacial thermal resistance (known as Kapitza resistance) between the PEEK matrix and SWNTs, which can be estimated from the acoustic mismatch theory, as follows [26]: where ρ m , C P m , ν m and R bd are the density of the PEEK matrix, the specific heat capacity of PEEK, the speed of sound in PEEK, and the interfacial thermal resistance between PEEK and SWNT, respectively. When a thermal walker jumps into a SWNT, it is assumed to travel with an infinite speed, due to the ballistic phonon transport and ultrahigh thermal conductivity of SWNTs [27]. The implementation of this assumption occurs by randomly placing the walker anywhere within the SWNT. The random placement is based on a uniform distribution function and it occurs in a single time step. The diffusivity inside the CNTs would need to be considered explicitly (using a second Gaussian random motion) if the CNTs were on the same order of magnitude or longer than the wavelength of a phonon. However, this case would present itself for much longer CNTs than the ones consider presently. When inside a SWNT, thermal walkers may exit the SWNT based on another probability, designated as f SWNT-m . This probability is determined from f m-SWNT , as [19]: where V SWNT and A SWNT are the volume and surface area of SWNTs, and C f´SWNT is a thermal equilibrium factor at the PEEK-SWNT interface which depends on the geometry of SWNTs and the interfacial area between PEEK and a SWNT. The thermal equilibrium factor is introduced in order to preserve the second law of thermodynamics at thermal equilibrium. The above relation between f SWNT-m and f m-SWNT can be explained as follow: When reaching the thermal steady state, the heat flux exiting a SWNT should be equal to that entering the SWNT. In a time step, all thermal walkers inside a SWNT may travel to the PEEK matrix owing to their infinite speed. However, among the thermal walkers in PEEK, only those around the SWNT surface may jump into the SWNT due to the Brownian motion. In order to maintain a balanced heat flux exiting and entering a SWNT, the relation as described in Equation (3) should be satisfied. When inside a WS 2 nanoparticle, a thermal walker jumps randomly, similar to the way it travels in the PEEK matrix, but with a different thermal diffusivity (that of WS 2 ) used in Equation (1). At the PEEK-WS 2 nanoparticle interface, thermal walkers from either the PEEK or the WS 2 behave similarly to the walkers crossing the PEEK-SWNT interface from the PEEK side. However, the two probabilities ( f m-WS 2 and f WS 2´m ) have different relation from that described in Equation (3), due to the different motion of walkers in WS 2 nanoparticles and SWNTs. As thermal walkers have Brownian motion in WS 2 nanoparticles, f m-WS 2 and f WS 2´m are related as: where f WS 2´m and f m-WS 2 are the walker travelling probabilities from the WS 2 nanoparticles to the PEEK matrix and the reverse. Variables r, σ m and σ WS 2 are the radius of the WS 2 nanoparticles, the standard deviation of Brownian motion in the PEEK matrix, and the standard deviation of Brownian motion in the WS 2 nanoparticles, respectively. The parameter C f´WS 2 is the thermal equilibrium factor at the PEEK-WS 2 nanoparticle interface, which can be numerically determined in the same manner as C f´SWNT [27][28][29]. Thermal walkers will be bounced back when they jump outside of the model to maintain a constant heat flux. Periodic boundary conditions are applied in the y and z directions. The computational domain was divided into 300ˆ300ˆ300 grids to calculate the temperature profile, which can be obtained by counting the number of hot walkers and then subtracting the number of cold walkers in a grid cell. Figure 1b shows a typical contour plot of thermal walker distribution in the thermal steady state. With constant heat flux applied along the x direction, the temperature along the x direction should be a straight line which has a slope inversely proportional to the thermal conductivity of the SWNT/WS 2 /PEEK composite [30]. A reference model of pure PEEK matrix was also built to estimate the effective thermal conductivity (K eff ) of the composite. With the same heat flux and boundary conditions, the temperature profiles along x direction in the composite model and the pure matrix model are related as: where q", T c and T m are the applied constant heat flux, the temperature in the composite and the temperature in the pure PEEK matrix, respectively. K m is the thermal conductivity of the pure PEEK matrix, which is known to be 0.23 W/m¨K from the literature [31]. Thus, as expressed in Equation (5), the effective thermal conductivity of the composite (K eff ) can be calculated based on the temperature profiles in the composite and the pure PEEK.

Results and Discussion
In this section, three case studies are presented and discussed to demonstrate the capabilities of the developed off-lattice Monte Carlo approach. The case study systems include three-phase SWNT/WS 2 / PEEK composites, SWNT/WS 2 /PEEK composites with SWNT bundles and graphene/PMMA composites with different sized graphene sheets. The typical parameters influencing or limiting the thermal transport properties of the composites were quantitatively investigated. The quantitative findings may provide an overview of the thermal transport phenomena and limitation mechanisms in multiphase composites.

Validation of the Developed Off-Lattice Monte Carlo Approach
The developed approach was validated by comparing the simulation results with the measured K eff of SWNT/WS 2 /PEEK composites with different compositions [9], as presented in Figure 2. The different compositions (e.g., 0.5/0.5/99.0) were marked with the mass fractions of the SWNTs (0.5%), the WS 2 nanoparticles (0.5%) and the PEEK matrix (99.0%). In the models for different compositions, the dimensions of the SWNTs and the diameter of the WS 2 nanoparticle were maintained, while the side length of the PEEK cube was varied. As shown in Figure 2, the simulation results from our approach are in good agreement with the experimental data, which validates the developed approach.
Keff of SWNT/WS2/PEEK composites with different compositions [9], as presented in Figure 2. The different compositions (e.g., 0.5/0.5/99.0) were marked with the mass fractions of the SWNTs (0.5%), the WS2 nanoparticles (0.5%) and the PEEK matrix (99.0%). In the models for different compositions, the dimensions of the SWNTs and the diameter of the WS2 nanoparticle were maintained, while the side length of the PEEK cube was varied. As shown in Figure 2, the simulation results from our approach are in good agreement with the experimental data, which validates the developed approach. Effective medium theory approaches (EMT) are commonly applied to predict the K eff of composites [33][34][35][36]. For comparison with our developed approach, two widely-used EMTs (i.e., the Maxwell-Garnett model (MG) [35,37] and Nan et al.'s model [38]) were utilized to predict the K eff of the SWNT/WS 2 /PEEK composites. The comparison among different models is presented in Figure 2. The MG-EMT overestimated the K eff of the SWNT/WS 2 /PEEK composites, which is likely due to the neglect of the interfacial thermal resistance at the SWNT-PEEK interface. On the contrary, Nan et al.'s EMT model underestimated the K eff of the SWNT/WS 2 /PEEK composites. This is because that model does not account for the synergistic effects of SWNTs and WS 2 nanoparticles in the composites.

Effects of Interfacial Thermal Resistances on the K eff of SWNT/WS 2 /PEEK Composites
The interfacial thermal resistance (R bd ) at the nanofiller-matrix interface (i.e., SWNT-PEEK, WS 2 -PEEK) is caused by the difference in the vibrational phonon spectra of each side, and by interfacial defects [39,40]. The R bd is found to greatly limit the K eff of SWNT/WS 2 /PEEK composites [10]. The R bd at the SWNT-PEEK (R SWNT-PEEK ) and WS 2 -PEEK (R WS 2´P EEK ) interfaces was varied to investigate the effects on the K eff of SWNT/WS 2 /PEEK composites, and the results are presented in Figure 3. The R bd is interpreted by an average phonon transmission probability at the interfaces in our approach, as described already above. It has been reported in the literature that the R bd in the nanoscale falls into the range of 1.0ˆ10´9-1.0ˆ10´6 m 2¨K /W [41][42][43][44], so the R SWNT-PEEK was varied between 1.158ˆ10´9 and 1.158ˆ10´6 m 2¨K /W, which corresponds to an average phonon transmission probability of 0.001-1.0. Similarly, the R WS 2´P EEK ranged from 2.32ˆ10´1 0 to 2.32ˆ10´8 m 2¨K /W, corresponding to an average phonon transmission probability from 0.005 to 0.5.
is interpreted by an average phonon transmission probability at the interfaces in our approach, as described already above. It has been reported in the literature that the Rbd in the nanoscale falls into the range of 1.0 × 10 −9 -1.0 × 10 −6 m 2 ·K/W [41][42][43][44], so the PEEK -SWNT R was varied between 1.158 × 10 −9 and 1.158 × 10 −6 m 2 ·K/W, which corresponds to an average phonon transmission probability of 0.001-1.0. Similarly, the PEEK -WS 2 R ranged from 2.32 × 10 −10 to 2.32 × 10 −8 m 2 ·K/W, corresponding to an average phonon transmission probability from 0.005 to 0.5.  As shown in Figure 3, the K eff of SWNT/WS 2 /PEEK composites decreases with both R SWNT-PEEK and R WS 2´P EEK . In Figure 3b, compared with that in the parallel SWNT case, the effect of R WS 2´P EEK on the K eff of SWNT/WS 2 /PEEK composites is much weaker in the random and perpendicular SWNT cases. R SWNT-PEEK significantly impedes the transfer of heat between SWNTs and the PEEK matrix, weakening the enhancement of SWNTs on the K eff . At a high R SWNT-PEEK , SWNTs with three different orientations lead to a similar K eff , which is close to the K PEEK , indicating that SWNTs do not enhance the heat conduction in the composite. A low R SWNT-PEEK is desired to obtain the composite with high K eff , which may be achieved by proper functionalization of SWNTs to couple the phonon spectra of SWNTs and PEEK [45][46][47]. The influence of R WS 2´P EEK on the K eff of composites is much weaker than that of R SWNT-PEEK . This may be because: (i) the ultrahigh thermal conductivity of SWNTs allows them to dominate the heat transfer through composites; (ii) the much larger interfacial area of SWNTs make the influence of the SWNT-PEEK interface more significant than that of the WS 2 -PEEK interface for heat transfer; (iii) the long cylindrical SWNTs are more effective than the spherical WS 2 nanoparticles to enhance the K eff of the composites [48], offering a larger characteristic length scale over which heat can be transferred.

Effects of the Morphology of SWNTs on the K eff of SWNT/WS 2 /PEEK Composites
Since SWNTs dominate the heat transfer in the three-phase composites, the morphology of SWNTs (i.e., length and diameter) was varied to study its effect on the K eff of SWNT/WS 2 /PEEK composites. The length was varied from 100 to 900 nm while keeping a constant diameter of 2 nm. The diameter was varied from 2 to 8 nm while the length was kept at 500 nm. As presented in Figure 4, when SWNTs are parallel or random to the heat flux, longer SWNTs induce higher K eff , which is consistent with previous studies [39,49]. Due to the ballistic phonon transport in SWNTs, longer SWNTs are more effective than short ones to transport heat through the composites, leading to a higher K eff [50]. For the study of SWNT diameter, a composite with randomly oriented SWNTs was chosen for investigation as it represents most realistic composites. As shown in Figure 4b, the K eff of SWNT/WS 2 /PEEK composites increases with the decrease of SWNT diameter, which may be ascribed to the interfacial area of SWNTs. With same volume fraction of SWNTs, a larger interfacial area between SWNTs and PEEK is obtained for SWNTs having smaller diameter. The larger interfacial area can afford more effective heat transfer channels between SWNTs and PEEK, thus inducing a larger enhancement of the K eff . The diameter used here covers the range of SWNTs (ď2 nm), double-walled CNTs (2-4 nm), and multi-walled CNTs (ě4 nm). Therefore, it can be speculated that SWNTs are more efficient than other CNTs to enhance the K eff of SWNT/WS 2 /PEEK composites.
SWNTs are more effective than short ones to transport heat through the composites, leading to a higher Keff [50]. For the study of SWNT diameter, a composite with randomly oriented SWNTs was chosen for investigation as it represents most realistic composites. As shown in Figure 4b, the Keff of SWNT/WS2/PEEK composites increases with the decrease of SWNT diameter, which may be ascribed to the interfacial area of SWNTs. With same volume fraction of SWNTs, a larger interfacial area between SWNTs and PEEK is obtained for SWNTs having smaller diameter. The larger interfacial area can afford more effective heat transfer channels between SWNTs and PEEK, thus inducing a larger enhancement of the Keff. The diameter used here covers the range of SWNTs (≤2 nm), doublewalled CNTs (2-4 nm), and multi-walled CNTs (≥4 nm). Therefore, it can be speculated that SWNTs are more efficient than other CNTs to enhance the Keff of SWNT/WS2/PEEK composites. The length was varied from 100 to 900 nm, corresponding to an aspect ratio from 50 to 450. The diameter was varied from 2 to 8 nm while the length was kept as 500 nm. Different compositions with randomly orientated SWNTs were chosen to study the effect of SWNT diameter. The error bars represent the standard deviation of the results obtained from 3 separate simulations with different distribution of SWNTs. Reproduced with permission from [18]. Copyright Elsevier, 2014. Figure 4. Effects of (a) length and (b) diameter of SWNTs on the K eff of SWNT/WS 2 /PEEK composites. The length was varied from 100 to 900 nm, corresponding to an aspect ratio from 50 to 450. The diameter was varied from 2 to 8 nm while the length was kept as 500 nm. Different compositions with randomly orientated SWNTs were chosen to study the effect of SWNT diameter. The error bars represent the standard deviation of the results obtained from 3 separate simulations with different distribution of SWNTs. Reproduced with permission from [18]. Copyright Elsevier, 2014.

Model of SWNT/WS 2 /PEEK Composites with SWNT Bundles
SWNTs tend to aggregate into bundles during the composite synthesis process due to strong van der Waals forces. The thermal conductivity of SWNT bundles is generally lower than that of an individual SWNT in the bundle, due to the interfacial thermal resistance among adjacent SWNTs [51][52][53]. Thus, the SWNT bundles may limit heat transfer in the composites [39]. A SWNT/WS 2 /PEEK composite model with SWNT bundles was built to shed some light on the influence of SWNT bundles on the thermal transport properties of the composites. As shown in Figure 5, there were individual SWNTs, SWNT bundles, WS 2 nanoparticles, and PEEK matrix in the system. In the present work, only the straight bundles with line contacts were taken into account. The CNT bundles with complicated interconnected networks were out of the scope of this article. SWNTs tend to aggregate into bundles during the composite synthesis process due to strong van der Waals forces. The thermal conductivity of SWNT bundles is generally lower than that of an individual SWNT in the bundle, due to the interfacial thermal resistance among adjacent SWNTs [51][52][53]. Thus, the SWNT bundles may limit heat transfer in the composites [39]. A SWNT/WS2/PEEK composite model with SWNT bundles was built to shed some light on the influence of SWNT bundles on the thermal transport properties of the composites. As shown in Figure 5, there were individual SWNTs, SWNT bundles, WS2 nanoparticles, and PEEK matrix in the system. In the present work, only the straight bundles with line contacts were taken into account. The CNT bundles with complicated interconnected networks were out of the scope of this article.  Figure 6a. The models with different SWNT orientations (parallel, random and perpendicular to the heat flux) were all built, as illustrated in Figure 5. As shown in Figure 6a, when SWNT bundles are parallel or randomly oriented relative to the direction of the heat flux, the Keff of

Effects of the Morphology of SWNT Bundles on the K eff of SWNT/WS 2 /PEEK Composites
The number of SWNT bundles was changed from 0 to 48 (no individual SWNTs) to investigate the effect of SWNT dispersion state on the K eff of the SWNT/WS 2 /PEEK composites, and the results are presented in Figure 6a. The models with different SWNT orientations (parallel, random and perpendicular to the heat flux) were all built, as illustrated in Figure 5. As shown in Figure 6a, when SWNT bundles are parallel or randomly oriented relative to the direction of the heat flux, the K eff of the composites slightly decreases with an increase of SWNT bundles. This may be caused by the non-uniform distribution of SWNTs with the presence of SWNT bundles, as well as the SWNT-SWNT thermal resistance within the bundles, which prevents heat from conducting along the heat flux direction [54]. Different from the above two SWNT orientations, more SWNT bundles perpendicular to the heat flux can induce a higher K eff of the SWNT/WS 2 /PEEK composites. The K eff slightly increases by 40% when the bundle number increases to 48. This is likely due to the bigger diameter of SWNT bundles compared with individual SWNTs, which accelerates the heat transfer in the radial direction (also the heat flux direction).
The number of individual SWNTs in each bundle was varied from 10 to 25, at an increment of 5, to study its influence on the K eff of the composites. The total number of SWNTs (960) and the number of SWNT bundles (36) were kept constant. As shown in Figure 5b, when SWNT bundles are parallel or randomly oriented relative to the heat flux, more SWNTs in each bundle induce a lower K eff . With the same number of bundles, more SWNTs per bundle lead to a worse distribution of SWNTs, which may reduce the heat transfer along the heat flux. On the other hand, more SWNTs in a bundle increase the stiffness of the bundle, which weakens the phonon coupling between SWNTs and PEEK via low-frequency vibrations [55], thus leading to a larger SWNT-PEEK thermal resistance [56]. When SWNT bundles are perpendicular to the heat flux, the number of SWNTs in each bundle has no apparent effect on the K eff of the composites. will exist among SWNTs within the bundles. Previous models cannot take into account this resistance for Keff prediction, whereas the developed approach can investigate this resistance by controlling the motion of thermal walkers at the SWNT-SWNT interface [24]. The SWNT SWNT R was varied from 6.153 × 10 −10 to 6.153 × 10 −7 m 2 ·K/W, corresponding to an average phonon transmission probability of 0.001-1.0 [57]. As presented in Figure 7, the Keff of the SWNT/WS2/PEEK composites decreases with the increase of the  3.2.2. Effects of the SWNT-SWNT Thermal Resistance on the K eff of SWNT/WS 2 /PEEK Composites When SWNTs aggregate into bundle structures, SWNT-SWNT thermal resistance (R SWNT-SWNT ) will exist among SWNTs within the bundles. Previous models cannot take into account this resistance for K eff prediction, whereas the developed approach can investigate this resistance by controlling the motion of thermal walkers at the SWNT-SWNT interface [24]. The R SWNT-SWNT was varied from 6.153ˆ10´1 0 to 6.153ˆ10´7 m 2¨K /W, corresponding to an average phonon transmission probability of 0.001-1.0 [57]. As presented in Figure 7, the K eff of the SWNT/WS 2 /PEEK composites decreases with the increase of the R SWNT-SWNT . A higher R SWNT-SWNT more greatly prevents heat from transferring among the bundled SWNTs, leading to a lower K eff of the composites.
A critical SWNT-SWNT thermal resistance, R c , was found to dominate the heat transfer in the composites, which was estimated to be 0.155ˆ10´8 m 2¨K /W. When R SWNT-SWNT < R c , more SWNT bundles may induce a higher K eff of the composite. This is because at low R SWNT-SWNT , heat prefers to transfer through SWNT-SWNT contacts, thus the SWNT bundles would be more effective heat transfer channels than the mono-dispersed SWNTs. It can be inferred that the detrimental effect of the SWNT bundles could be reduced by decreasing the R SWNT-SWNT to be less than R c , which may be achieved by covalent functionalization to enhance the phonon coupling and weaken the phonon scattering at a SWNT-SWNT interface [58,59].
will exist among SWNTs within the bundles. Previous models cannot take into account this resistance for Keff prediction, whereas the developed approach can investigate this resistance by controlling the motion of thermal walkers at the SWNT-SWNT interface [24]. The SWNT SWNT R was varied from 6.153 × 10 −10 to 6.153 × 10 −7 m 2 ·K/W, corresponding to an average phonon transmission probability of 0.001-1.0 [57]. As presented in Figure 7, the Keff of the SWNT/WS2/PEEK composites decreases with the increase of the

R
, heat prefers to transfer through SWNT-SWNT contacts, thus the SWNT bundles would be more effective heat transfer channels than the mono-dispersed SWNTs. It can be inferred that the detrimental effect of the SWNT bundles could be reduced by decreasing the SWNT SWNT R to be less than c R , which may be achieved by covalent functionalization to enhance the phonon coupling and weaken the phonon scattering at a SWNT-SWNT interface [58,59].

Model of Graphene/PMMA with Complex Structure
In the graphene-based polymer composites, graphene sheets have a large size distribution in the x-y plane (50-500 nm) due to the oxidation and reduction of graphene sheets during the composite fabrication [3]. The existing models commonly ignore the size distribution of graphene sheets, resulting in an inaccurate prediction of the thermal conductivity. In the developed approach, one can take into account the size distribution (various length, width and thickness), volume fraction and orientation of graphene sheets, as well as the interfacial thermal resistance at the graphene-polymer interface.
In this subsection, a graphene/PMMA composite model was built to validate that the develop approach can be applied to study graphene/polymer composites. Graphene sheets with various length (50-500 nm), width (50-500 nm), thickness (2.4-9.0 nm) and different orientations (parallel, random and perpendicular to the heat flux) were generated in the model.
The developed graphene/PMMA model was validated by comparing the predicted K eff with the measured results, as shown in Figure 8a. As the interfacial thermal resistance (R bd ) at the graphene-PMMA interface was an input value, the model was validated as follows: Estimate the R bd by matching the simulated K eff with the measured value for one composite, and then use this estimated R bd as input to calculate the K eff of other composites. The predicted K eff showed a good agreement with the experimental data, validating the developed model. The R bd was estimated to be 1.906ˆ10´8 m 2¨K /W. For comparison, a modified effective medium theory (EMT) was used to calculate the K eff of the graphene/PMMA composites, and the results are presented in Figure 8b. The estimated K eff from the modified EMT is much higher than the experimental results, and even higher than the predicted K eff of the model with parallel graphene. This is likely because the modified EMT cannot take into account the size distribution of graphene sheets. In the modified EMT, the length and width were treated as infinite when compared with the thickness of graphene sheets, which failed to take into account the graphene sheet with relatively short length and width (e.g.,~50 nm) [60][61][62]. It should be noted that if given the K eff of a specific composite, the interfacial thermal resistance between the nanofillers and the matrix can be estimated by using the developed approach and an inverse calculation procedure.
Rbd as input to calculate the Keff of other composites. The predicted Keff showed a good agreement with the experimental data, validating the developed model. The Rbd was estimated to be 1.906 × 10 −8 m 2 ·K/W. For comparison, a modified effective medium theory (EMT) was used to calculate the Keff of the graphene/PMMA composites, and the results are presented in Figure 8b. The estimated Keff from the modified EMT is much higher than the experimental results, and even higher than the predicted Keff of the model with parallel graphene. This is likely because the modified EMT cannot take into account the size distribution of graphene sheets. In the modified EMT, the length and width were treated as infinite when compared with the thickness of graphene sheets, which failed to take into account the graphene sheet with relatively short length and width (e.g., ~50 nm) [60][61][62]. It should be noted that if given the Keff of a specific composite, the interfacial thermal resistance between the nanofillers and the matrix can be estimated by using the developed approach and an inverse calculation procedure.  More details of the experimental set-up can be found in Reference [4]; (b) the thermal conductivity of GA-PMMA composites as a function of graphene volume fraction. The interfacial thermal resistance between graphene sheets and PMMA was estimated to be R bd = 1.906ˆ10´8 m 2¨K /W in the developed model. The same value was utilized in the composites with parallel-and perpendicular-oriented graphene. In the modified EMT, the utilized thermal conductivity of graphene and the R bd of graphene-PMMA were 100 W/m¨K and 1.0ˆ10´8 m 2¨K /W, respectively. The error bars represent the standard deviation of the results from separate measurements of thermal conductivity of graphene/PMMA composites. Reproduced with permission from [4]. Copyright Elsevier, 2015.

Conclusions
In summary, the developed off-lattice Monte Carlo approach has proved to be accurate as a computational model for heat transfer phenomena and heat transfer mechanisms for multiphase composites with complex structures. The developed approach not only provides a more accurate method to predict the K eff of the multiphase composites than existing EMT models, but also offers an effective computational approach to estimate the interfacial thermal resistance between the nanofillers and the matrix. The quantitative findings presented herein showed that multiphase composites with higher K eff can be obtained by (a) reducing the interfacial thermal resistances at filler-matrix interfaces; (b) aligning the fillers along the heat flux direction; (c) using fillers with larger interfacial area; and (d) improving the dispersion of fillers to be more uniform rather than forming bundles. Through proper modifications of the geometry and thermal properties of the components in the model, the developed approach may be applied to study the thermal transport properties of other multiphase systems, such as CNT/graphene/polymer composites, graphene stabilized polymer blends, CNT stabilized emulsions and other multiphase organic or inorganic composites.