Recent research has highlighted the potential to achieve high-thermal-conductivity polymers by aligning their molecular chains. Combined with other merits, such as low-cost, corrosion resistance, and light weight, such polymers are attractive for heat transfer applications. Due to their quasi-one-dimensional structural nature, the understanding on the thermal transport in those ultra-drawn semicrystalline polymer fibers or films is still lacking. In this paper, we built the ideal repeating units of semicrystalline polyethylene and studied their dependence of thermal conductivity on different crystallinity and interlamellar topology using the molecular dynamics simulations. We found that the conventional models, such as the Choy-Young's model, the series model, and Takayanagi's model, cannot accurately predict the thermal conductivity of the quasi-one-dimensional semicrystalline polyethylene. A modified Takayanagi's model was proposed to explain the dependence of thermal conductivity on the bridge number at intermediate and high crystallinity. We also analyzed the heat transfer pathways and demonstrated the substantial role of interlamellar bridges in the thermal transport in the semicrystalline polyethylene. Our work could contribute to the understanding of the structure–property relationship in semicrystalline polymers and shed some light on the development of plastic heat sinks and thermal management in flexible electronics.

Polymers have been widely used in a variety of energy-related applications, such as solar cells,1,2 thermoelectric devices,3 and batteries4 because they are lightweight, corrosion-resistive, stable, nontoxic, and inexpensive. However, the low thermal conductivity of pure amorphous polymers, which is typically 0.10.5W/(mK),5–7 to a large extent limits their applications in the heat dissipation area. Efforts to improve the thermal conductivity of polymers-based materials can be classified into two categories: (1) blending the polymer with materials with higher thermal conductivity, such as manufacturing polymer composites with hybrid fillers,8,9 or metal powders;10 (2) mechanically aligning the polymer chains and taking advantage of the high thermal conductivity along the molecular chain backbone.11–13 Various theoretical calculations and simulations have proved that the thermal conductivity of a single polymer chain or a crystalline polymer unit is considerably high along the chain direction, which can even exceed 100 W/(m K).12,14 In this spirit, experimentally ultra-drawn polyethylene (PE) nanofibers were reported to have as extremely high thermal conductivity as ∼100 W/(m K); thermal conductivity of the commercial ultra-drawn single PE microfibers (e.g., spectra 2000) has been measured as ∼10–16 W/(m K).15–17 Zhu et al. showed that the thermal conductivity of those commercial microfibers could increase ∼80%–150% with a moderate draw ratio 2–6.16 Because fibers provide limited practical benefits, researchers start to propose methods to massively fabricate large-area polyethylene films with high thermal conductivity. Ghasemi et al. proposed to manufacture ultra-drawn highly aligned PE polymer films through a sol-gel and ultra-drawn process.18 The measured thermal conductivity of these preliminary high-crystallinity highly aligned PE polymer films can be comparable to a single commercial fiber, 16W/(mK). Recently, an even higher thermal conductivity, ∼60 W/(m K), was achieved in scalably fabricated PE films by a roll-to-roll system.19 In all of these nanofiber, microfiber, and film samples, the PE is in a semicrystalline state with high structural anisotropy while the thermal transport in those materials is not well understood.

Semicrystalline polymer is a two-phase mixture system with the ordered crystalline lamellar phase and the disordered amorphous phase. Isotropic semicrystalline polymer is generally treated as a spherulite with the folded chains forming an isotropic spheroidal structure; anisotropic semicrystalline polymer can be considered as an alternative stack of the two phases to simplify the physical models. An ideal repeating unit of semicrystalline polymers is a sandwich structure with two crystalline regions at the sides called lamella and the amorphous region in the middle called interlamellar region.20,21 Theoretically, Choy22 proposed a thermal conductivity model based on the Maxwell effective medium equation to predict how the thermal conductivity of semicrystalline polymers depends on crystallinity. For highly oriented semicrystalline polymers with an orientation angle θ between the lamellar axis and the drawn direction, Choy-Young's model is

KeKaKe+2Ka=X[Kc/Ka1Kc/Ka+2sin2θ+Kc/Ka1Kc/Ka+2cos2θ],
(1)

where Ke, Ka, Kc, Kc, and X are the thermal conductivity of the semicrystalline polymer, the amorphous counterpart, crystalline counterpart (along the chain), the crystalline counterpart (perpendicular to the chain), and crystallinity, respectively. This model fits well with the experimental data of isotropically distributed lamellas;22 nevertheless, it fails in highly drawn semicrystalline polymers. Besides, the series model23,24 is widely used to apply the thermal resistance network concept in such a two-phase mixture system. However, it underestimates the thermal conductivity of semicrystalline polymers at high crystallinity as well. The reason is that both models treat the entire interlamellar region as one amorphous region and underestimates the role of the chain topologies. Specifically, interlamellar bridges that connect two adjacent crystalline lamellas increase their fractions and change their chain alignment dramatically with stretching, which play an important role in the heat conduction. Takayanagi25 proposed a model that specifically treats those bridges the same as the crystalline stems, which are aligned in the extrusion direction and have the same thermal conductivity  Kc. The thermal conductivity of the semicrystalline region depends on a fitting parameter b, which was roughly assumed to be the width ratio of the region formed by these highly oriented bridges (defined as the ratio of the width of the bridge region to the width of the whole semicrystalline region). However, the Takayanagi model fails to predict the temperature dependence of thermal conductivity due to overestimating the contribution from bridges. Similar to the Takayanagi's model, Halpin-Tsai equation26 was adopted based on the Cox model27 that considered the bridges as a fiber phase connecting at least two lamellas and forming a continuous crystal region, while the thermal conductivity was fitted with the aspect ratio of the fibers (defined as the ratio of the diameter to the length of the fibers) as a free parameter.28 Both Takayanagi's and Cox's models treat bridges as tie chains (straight ultra-aligned chains), which could be true at ultrahigh draw ratios, however, would lose their predictive power at low and intermediate draw ratios when bridges are in a much coiled state. Moreover, all the models above could not reveal the heat transfer pathways among different chain topologies in the interlamellar region.

Recently, the enhanced Monte Carlo (EMC)21,29 method developed by in't Veld et al. was widely adopted to simulate the mechanical and structural properties of semicrystalline polymers, especially under tensile deformation.30,31 This method is a powerful tool that performs random Monte Carlo (MC) moves to mimic the formation of interlamellar region of semicrystalline systems and allows the investigation into different interlamellar topology. In this work, we built the ideal one-dimensional repeating unit of semicrystalline PE systems with different crystallinities and interlamellar topology using the EMC method. Their thermal conductivities were calculated by the nonequilibrium molecular dynamics (NEMD) simulations. Results showed that besides the dependence on crystallinity, thermal conductivity highly depends on bridge numbers at intermediate and high crystallinity. Our study explained why the Choy-Young's model and series model lost its predictive power at high crystallinity. Moreover, we modified the Takayanagi's model by incorporating the thermal conductivity of bridges and analyzing the heat transfer pathways along bridges, loops, and tails in the interlamellar region, providing a guide on evaluating the thermal conductivity of high-crystallinity polymers.

Figure 1 shows the schematics of a typical semicrystalline polyethylene model in our simulations. We chose the minimum representative repeating unit, including two crystalline regions and an interlamellar region, as our model. We used the crystalline PE unit cell to represent the crystalline region, which is a different but simplified version from the switchboard model of a polymer crystalline lamella. After quenching, the interlamellar region evolves into an amorphous region and a transitional region (i.e., inter-phase regions). Typically, we take into account three types of chain topologies in the interlamellar region, i.e., bridges, loops, and tails. Among these, bridge is defined as the segments connecting the crystalline stems throughout the interlamellar region; loop is defined to connect crystalline stems at the same side; and for tail, one side of it is connected to the crystalline stems while the other side terminates.

FIG. 1.

Schematic diagram of semicrystalline structure, which contains two crystalline regions (white regions on both sides), two interphase regions (gray shaded regions), and an amorphous region (intermediate white region). The inter-phase regions and the amorphous region sandwiched in the middle are defined as an interlamellar region. These regions were distinguished according to the density profiles shown in  Appendix B. Crystalline stems are denoted by purple slashes, bridges are denoted by green solid lines, loops are denoted by blue dashed lines, and tails are denoted by orange dotted lines.

FIG. 1.

Schematic diagram of semicrystalline structure, which contains two crystalline regions (white regions on both sides), two interphase regions (gray shaded regions), and an amorphous region (intermediate white region). The inter-phase regions and the amorphous region sandwiched in the middle are defined as an interlamellar region. These regions were distinguished according to the density profiles shown in  Appendix B. Crystalline stems are denoted by purple slashes, bridges are denoted by green solid lines, loops are denoted by blue dashed lines, and tails are denoted by orange dotted lines.

Close modal

The EMC method is adopted to build those models following the procedures described in Refs. 30–33. We briefly describe the major procedures as follows: (1) initially, a crystalline PE structure comprising 4×6×76 unit cells was generated with the lattice constant a=7.72Ǻ,b=4.45Ǻ,c=2.53Ǻ, resulting in 7296 atoms and 48 chains in the simulation box. We used the united-atom (UA) force field34 (see Table III in  Appendix A for parameters in our simulation) to describe the interatomic potential. The chains were then tilted to make sure that the {201} plane is normal to the z-direction, following the procedure in Refs. 21, 31, and 35, because such an orientation in semicrystalline PE has been suggested as one of the most common configurations due to its minimum interfacial energy.36,37 We selected and fixed 10, 15, 24, 30, 33, and 38 layers of atoms on each side in the z-direction as the crystalline region in 0%-, 18%-, 50%-, 72%-, 83%-, and 100%-effective crystallinity configurations, respectively. Each layer contains 4×6×1 unit cells. Here, the effective crystallinity X is defined as X=Lc/Ls, with Lc and Ls as the effective length of crystalline and the whole semicrystalline regions used for fitting in the following NEMD thermal conductivity calculations, respectively. We note here that even for the 0%-crystallinity configurations, there are still 20 layers crystalline region on each side because we need to reserve the heat source and heat sink regions for NEMD thermal conductivity calculations where the linear fitting region is entirely amorphous. (2) Next, 16 of the 48 chains were cut randomly and 1201, 1000, 600, 343, and 214 united atoms in these chains were removed in 0%-, 18%-, 50%-, 72%-, and 83%-effective crystallinity configurations, respectively, to match the experimental amorphous density of 0.855g/cc(Ref.38) at T =300K and P=1atm. This procedure generated 32 tails and 32 bridges in the interlamellar region. (3) After that, all the united atoms in the interlamellar region underwent random Monte Carlo (MC) moves, including single-site displacement, end-rotation, re-bridging, end-reptation, and end-bridging, to form the interphase and amorphous regions. Among these moves, bond formation and breakage were involved. Specifically, single-site displacement displaces a random site slightly to change the bond lengths and angles; end-rotation rotates the sites at a tail end with one to three sites involved in the rotation; re-bridging cuts three sites, which are consecutively bonded, from a bridge and then reinserts them in the bridge by searching for the possible insertion positions without altering the lengths or angles; end-reptation excises and appends a segment of one tail end to another tail end with random torsion angles, forming new tails; end-bridging move introduces the conversion between bridge-tail pairs and loop-tail pairs. During this procedure, the number of tails and the total number of bridges and loops remained the same. Detailed information can be found in Refs. 21 and 31. These random MC moves were performed at 10 000 K in order to promote a rapid randomization of the interlamellar configuration, generating the amorphous phase.31 (4) Finally, these configurations were equilibrated at 10 000 K for 20 000 Monte Carlo cycles, followed by quenching from 10000 K to 300 K stepwise during 600 000 MC cycles. The quenching was performed in the sequence of 5000 K, 2000 K, 1000 K, 750 K, 500 K, 400 K, 350 K, and 300 K in 20 000, 40 000, 40 000, 40 000, 60 000, 80 000, 160 000, and 160 000 MC cycles, respectively. A merit in building configurations by the EMC method is that the atoms in the interphase region are covalently bonded to the crystalline stems and amorphous regions instead of simply stacked through the van der Waals interactions. We prepared 8 to 10 independent initial configurations for each crystallinity system. A configuration with X =50% built by the EMC method is shown in Fig. 2(a) for illustration.

FIG. 2.

(a) Schematic diagram for NEMD simulations. Atoms in crystalline regions are in purple, while atoms in the interlamellar regions are distinguished according to different chain topologies: bridges (green), loops (cyan), and tails (orange). Vacuum regions were added on each side along the z-direction. Atoms in the boundary regions were fixed. Heat source and sink were applied in the following regions. Fitting region was chosen far away enough from the thermostats. (b) A typical temperature profile in the NEMD simulation and the fitting parameter for effective thermal conductivity in semicrystalline PE. (c) A typical temperature profile in the NEMD and the fitting line for amorphous thermal conductivity in 0%-crystallinity PE.

FIG. 2.

(a) Schematic diagram for NEMD simulations. Atoms in crystalline regions are in purple, while atoms in the interlamellar regions are distinguished according to different chain topologies: bridges (green), loops (cyan), and tails (orange). Vacuum regions were added on each side along the z-direction. Atoms in the boundary regions were fixed. Heat source and sink were applied in the following regions. Fitting region was chosen far away enough from the thermostats. (b) A typical temperature profile in the NEMD simulation and the fitting parameter for effective thermal conductivity in semicrystalline PE. (c) A typical temperature profile in the NEMD and the fitting line for amorphous thermal conductivity in 0%-crystallinity PE.

Close modal

MD simulations using the LAMMPS package39 were performed to equilibrate the ensembles generated by the EMC methodology and to study their thermal properties. The same UA force field and periodic boundary conditions in all three directions were applied. Time step equals 1 fs in all simulations because there are no fast-vibrating hydrogens in the model. Independent configurations for each crystallinity were equilibrated in the NPT (constant total number of particles N, pressure P, and temperature T) ensemble for 8 ns, during which the pressure P =1 atm and temperature T =300 K were maintained. The pressure and temperature damping parameters were 1000 fs and 100 fs, respectively. After 8 ns of NPT equilibration, the fluctuations of the system volume are within 0.3% and energy converges within 0.7%. Then, 0.2 ns NVE ensemble (constant total number of particles N, the volume V, and the total energy E of the system) was performed for relaxation. During MD simulations, the numbers of bridges, loops, and tails do not change.

We summarized the number of bridges, loops, and tails and the simulation box parameters of each configuration after MD equilibrations in Table I. The length in the z-direction in different systems was roughly the same (163Ǻ). The fitting length Ls for the effective thermal conductivity is fixed as 120Ǻ, and thus, the lengths of crystalline (amorphous) region for fitting Lc (Li) are 0 (120), 22 (98), 60 (60), 86.5 (33.5), and 99.2 (20.8) Ǻ at 0%-, 18%-, 50%-, 72%- and 83%-crystallinity, correspondingly. These crystalline and interlamellar (inter-phase and amorphous) regions were distinguished by the density profiles, with details shown in the  Appendix B. The three interlamellar topologies (bridges, loops, and tails) of independent configurations are also summarized in Table I. Note that the absolute numbers of bridges and loops depend on the random MC moves in EMC and they will reach an equilibrium number after sufficient long simulation time. We listed the numbers of bridges and loops for each crystallinity. We also note that some of the independent configurations at the same crystallinity have the same bridge number; however, the alignment and the atom number of each chain are different due to the random MC moves. Among all the configurations, the tail number is fixed as 32, and the average number of bridge (loop) is 2 (30), 4.2 (27.8), 11.7 (20.3), and 20.6 (11.4) at 18%-, 50%-, 72%- and 83%-crystallinity, respectively. Obviously, the length of interlamellar region decreases with increasing crystallinity, and thus, more bridges tend to be formed in higher crystallinity systems, which is at the cost of loop number.

TABLE I.

Parameters in each independent configuration at different crystallinities (X) after MD equilibration. The bridge, loop, and tail numbers are itemized as Bridge No., Loop No., and Tail No., respectively. Lx, Ly, and Lz are the average lengths in the x, y, and z directions of the simulation boxes, respectively. Lc (Li) is the length of crystalline (interlamellar) region when fitting for the effective thermal conductivity.

CrystallinityX =0%X =18%X =50%X =72%X =83%X =100%
Bridge No. … 1, 2, 4 2, 3, 4, 6, 7 9, 10, 11, 12, 14, 15 16, 18, 19, 20, 21, 22, 23, 27 … 
Loop No. … 31, 30, 28 30, 29, 28, 26, 25 23, 22, 21, 20, 18, 17 16, 14, 13, 12, 11, 10, 9, 5 … 
Tail No. … 32 32 32 32 … 
Lx(Ǻ) 33.9±0.9 34.1±0.9 35.2±0.4 35.9±0.5 36.3±0.5 38.5 
Ly(Ǻ) 27.8±0.3 27.8±0.2 27.8±0.1 28.1±0.2 28.3±0.2 27.1 
Lz(Ǻ) 165.7±3.5 166.5±4.3 165.8±1.7 163.0±3.0 161.5±1.6 149.4 
Lc(Ǻ) 0.0 22.0 60.0 86.5 99.2 120.0 
Li(Ǻ) 120.0 98.0 60.0 33.5 20.8 0.0 
CrystallinityX =0%X =18%X =50%X =72%X =83%X =100%
Bridge No. … 1, 2, 4 2, 3, 4, 6, 7 9, 10, 11, 12, 14, 15 16, 18, 19, 20, 21, 22, 23, 27 … 
Loop No. … 31, 30, 28 30, 29, 28, 26, 25 23, 22, 21, 20, 18, 17 16, 14, 13, 12, 11, 10, 9, 5 … 
Tail No. … 32 32 32 32 … 
Lx(Ǻ) 33.9±0.9 34.1±0.9 35.2±0.4 35.9±0.5 36.3±0.5 38.5 
Ly(Ǻ) 27.8±0.3 27.8±0.2 27.8±0.1 28.1±0.2 28.3±0.2 27.1 
Lz(Ǻ) 165.7±3.5 166.5±4.3 165.8±1.7 163.0±3.0 161.5±1.6 149.4 
Lc(Ǻ) 0.0 22.0 60.0 86.5 99.2 120.0 
Li(Ǻ) 120.0 98.0 60.0 33.5 20.8 0.0 

The thermal conductivity was calculated by the NEMD simulations.40,41 In order to establish a temperature gradient along the z-direction without introducing the coupling between thermostats, vacuum regions were added on each side, and the atoms in the boundary regions along the z-direction were fixed. Then, the Langevin thermostats were applied to the heat source and sink regions as shown in Fig. 2(a).

The non-equilibrium simulations were performed for 5 ns to make sure that the system is in steady sate. In steady state, the heat flux Jz along the z-direction is

Jz=12A(|ΔEsourceΔt|+|ΔEsinkΔt|).
(2)

Here, A is the cross-sectional area in the x-y plane, Δt is the time step, ΔEsource and ΔEsink are the amounts of heat added to the heat source and subtracted from the heat sink, respectively. Results show that ΔEsourceΔt and ΔEsinkΔt are the same, confirming the symmetry of our system.

The effective thermal conductivity Ke of the whole semicrystalline region is defined as

Jz=KeΔTΔL,
(3)

where ΔT is the temperature difference across the fitted region, ΔL is the length of the whole fitted region. All the parameters were calculated in the z-direction. Regions used to fit the thermal conductivity are far away enough from the heat bath regions, as shown in Fig. 2(b). The amorphous thermal conductivity Ka was fitted in the 0%-crystallinity configurations, as shown in Fig. 2(c), and the middle region was adopted to fit Ka.

Figure 3(a) shows the dependence of the effective thermal conductivity of all the model polymers on their crystallinity X. Thermal conductivities of independent configurations with the same X were averaged statistically. To compare with the Choy-Young's model, we calculated the intrinsic thermal conductivity along the chain Kc∥ and perpendicular to the chain Kc of a crystalline PE system and plotted the prediction with the consideration of the tilted angle between the crystalline chain direction and the z-direction. To calculate Kc and Kc, we built configurations with the {001} plane normal to the z direction so that the chain direction is along the z-direction. Since the NEMD method has a strong size effect on calculating the thermal conductivity in the chain direction, we first calculated the thermal conductivity of this configuration with an effective length 12 nm in the chain direction, Kc = 35 ±3 W/(m-K), and we assume the intrinsic thermal conductivity along the chain direction Kc∥ to be larger than this value. We also found that the Choy-Young's model is not sensitive to Kc when Kc > 10 W/(m-K), so we used Kc = 35 ±3 W/(m-K) in the model. Kc was obtained using the NEMD method by applying the heat baths on the two sides in the x-direction with an effective length 10 nm. Kc = 0.10±0.02 W/(m-K), which is comparable to that of the amorphous counterpart. For the series model, we used the calculated averaged thermal conductivity of 100%- and 0%- crystallinity systems as Kc and Ka and plot the series model predictions Ke in the shaded area. Using the UA force field and NEMD, Kc is 21±2 W/(m-K), and Ka is 0.12±0.02 W/(m K). Those calculated thermal conductivity values agree reasonably well with other simulation results42,43 and are consistently lower than the experimental results15,20,44 due to the simplified UA force field45,46 and size effect.47 The UA force field has been used to study the thermal conductivity of amorphous polymers, such as PE48 and polyamide-6,6.49 Both simulations agree well with experimental values.5–7,50 Besides, the UA force field has also been used to simulate the thermal conductivity of single crystalline PE chains,51 which gave an underestimated value compared with those using AIREBO potential.14 As pointed out in Ref. 45, using the UA force field, the specific heat of single PE chain at low temperature could be overestimated, and an inaccurately predicted phonon dispersion could lead to an underestimation in the acoustic sound velocity at low reduced wave vectors. We note here that the absolute thermal conductivity values from our simulations are not critical as long as those values that we used to testify the models are consistently from the same force field and simulation method. The use of UA model enables us to study the thermal properties in a much larger system with the limited computational resources. For polymers with different crystallinity, even though a length-dependent thermal conductivity might apply for the crystalline region, both the Choy-Young's model and the series model are not sensitive to the accurate value of Kc that is two orders of magnitude higher than Ka.

FIG. 3.

(a) Effective thermal conductivity (Ke) of the semicrystalline PE as a function of crystallinity. The orange dots represent the average value of MD simulation results at each X. Estimations from the Choy-Young's model and the series model with uncertainties from thermal conductivity inputs are shown as the blue and green shaded area, respectively. (b) Effective thermal conductivity of the semicrystalline PE as a function of bridge number (n). Three regions were observed according to the behavior of thermal conductivity as a function of n. Region I includes configurations with 18%, where the thermal conductivity of the interlamellar region can be roughly described by the pure amorphous thermal conductivity Ka. Region II includes configurations with X =50% and 72%, where Ke increases linearly with n. The orange and blue shaded areas are the predictions from our modified Takayanagi's model with Kb=0.9±0.1W/(mK). Region III includes configurations with X =83%, where Ke increases linearly with n in a much larger slope. The prediction results using Kb=0.53+0.076nW/(mK) with ±10% uncertainty are shown in the green shaded area. The intrinsic thermal conductivity of bridges Kb is estimated from a separate simulation that only allows heat flowing in the bridges.

FIG. 3.

(a) Effective thermal conductivity (Ke) of the semicrystalline PE as a function of crystallinity. The orange dots represent the average value of MD simulation results at each X. Estimations from the Choy-Young's model and the series model with uncertainties from thermal conductivity inputs are shown as the blue and green shaded area, respectively. (b) Effective thermal conductivity of the semicrystalline PE as a function of bridge number (n). Three regions were observed according to the behavior of thermal conductivity as a function of n. Region I includes configurations with 18%, where the thermal conductivity of the interlamellar region can be roughly described by the pure amorphous thermal conductivity Ka. Region II includes configurations with X =50% and 72%, where Ke increases linearly with n. The orange and blue shaded areas are the predictions from our modified Takayanagi's model with Kb=0.9±0.1W/(mK). Region III includes configurations with X =83%, where Ke increases linearly with n in a much larger slope. The prediction results using Kb=0.53+0.076nW/(mK) with ±10% uncertainty are shown in the green shaded area. The intrinsic thermal conductivity of bridges Kb is estimated from a separate simulation that only allows heat flowing in the bridges.

Close modal

Our NEMD calculation results agree well with the Choy-Young's model and the series model in their general trends, i.e., the effective thermal conductivity increases with increasing crystallinity. Quantitatively, the Choy-Young's model provides an excellent estimation of the average value at X =18% and 50%, but loses the predictive power at intermediate and high crystallinity, i.e., 72% and 83%, especially at X =83%. The discrepancy between the calculated value and the maximum Choy-Young's model prediction increases from ∼35% at 72%-crystallinity to ∼68% at 83%-crystallinity. Our results confirm the failure of Choy-Young's model at high crystallinity in the highly anisotropic semicrystalline polymer. Similarly, using the series model, the discrepancy between the calculated value and the maximum prediction from the model increases from ∼10% at 18%-crystallinity to ∼72% at 83%-crystallinity. We attribute such a large deviation of both models to the increased amount of bridges so that a single thermal conductivity of the amorphous counterpart is not able to describe the thermal conductivity in the interlamellar region. This statement is supported by the calculated “local effective thermal conductivity,” which is defined as the ratio of the heat flux to the local temperature gradient (see  Appendix C for the plots). The local effective thermal conductivity of the amorphous region is neither the thermal conductivity of the amorphous counterpart nor the crystalline counterpart.

Figure 3(b) shows the thermal conductivity of semicrystalline PE as a function of their bridge number n. We group the calculated results into three regions: (I) at a low crystallinity, for example, X =18%, their bridge numbers are mainly in the range of 1–4. In these configurations, the effective thermal conductivity of the entire amorphous regions can be roughly treated the same as Ka; (II) at an intermediate crystallinity, i.e., X =50% and 72%, and their bridge numbers are in the range of 3–15. The role of bridges and their interaction with atoms in loops and tails becomes important, and the effective thermal conductivity increases linearly with n, where predictions from the series model start to deviate; (III) at a high crystallinity (X =83%), the bridge number can even achieve 27, thus heat transfer along bridges dominates so that the effective thermal conductivity increases linearly with n in a much larger slope. Those results indicate that not only the crystallinity, but also the number of bridges should be considered to estimate the effective thermal conductivity of semicrystalline polymer systems with high crystallinity.

Even though Takayanagi's model takes into consideration the role of bridges, the failure of this model lies in the treatment of interlamellar bridges the same as crystalline stems. In order to investigate the role of bridges at high crystallinity, we seek to modify the Takayanagi's model. Instead of treating bridges as a continuous region of crystalline stems, we treated them separated from crystalline region, with a thermal conductivity Kb, and we built the model based on a thermal resistance network concept. Figure 4(a) shows our modified model. In the interlamellar region, bridges are in parallel with the amorphous part formed by the rest of the chain topologies, i.e., loops and tails. And the interlamellar region is in series with the crystalline region. The modified effective thermal conductivity is

Ke=[XKc+1XbKb+(1b)Ka]1.
(4)

Here, b is the width ratio of the bridge region, which assumes to be b =n/N, where n is the bridge number, and N is the total chain number (In our model, N is 48).

FIG. 4.

(a) Modified Takayanagi's model. The paralleled amorphous and bridge regions are in series with the crystalline region. The amorphous region is assumed to be formed by the loops and tails. Therefore, the modified effective thermal conductivity of semicrystalline system is Ke=[XKc+1XbKb+(1b)Ka]1. b is the ratio of width of bridge region, and X is the crystallinity. (b) Schematic diagram for possible pathways (P1, P2, and P3) among groups u (blue dashed lines), l (orange dashed lines), and b (green solid lines). The symbols u, l, and b represents the upper region, the lower region, and the bridges, respectively.

FIG. 4.

(a) Modified Takayanagi's model. The paralleled amorphous and bridge regions are in series with the crystalline region. The amorphous region is assumed to be formed by the loops and tails. Therefore, the modified effective thermal conductivity of semicrystalline system is Ke=[XKc+1XbKb+(1b)Ka]1. b is the ratio of width of bridge region, and X is the crystallinity. (b) Schematic diagram for possible pathways (P1, P2, and P3) among groups u (blue dashed lines), l (orange dashed lines), and b (green solid lines). The symbols u, l, and b represents the upper region, the lower region, and the bridges, respectively.

Close modal

The interfacial thermal resistance between the crystalline and interlamellar region was negligible at 300 K in our configurations because (1) the covalent bonding between the two regions is much stronger than the van der Waals interactions. The interfacial thermal conductance has been proved to increase dramatically with stronger bonding strength.52–54 (2) The thermal resistance of the interface region is much smaller compared with that in the crystalline region and amorphous region. Since thermal conductivity of the amorphous region is two orders of magnitude smaller than that of the crystalline region, the effective thermal conductivity is dominated by the interlamellar thermal conductivity. Therefore, the interfacial thermal resistance between the crystalline and interlamellar region can be reasonably neglected. This negligible interface thermal resistance can be seen from the slope of temperature profile at the interface region, which does not show a noticeable deviation from that of the crystalline region in the NEMD simulation [Fig. 2(b)], proving the estimation in Ref. 20.

To further understand the transition from region I to III, we divided the atoms in the interlamellar region into three groups, which include: (1) the upper group, where tails and loops are connected to the upper crystalline chain stem (group u); (2) the lower group, where tails and loops are connected to the lower crystalline chain stems (group l); (3) bridges, which connect the upper and lower crystalline chain stems (group b). The contribution from each group to the total heat flux was then studied. The possible pathways of heat transfer among these groups due to non-bonded interactions are simplified in Fig. 4(b).

In LAMMPS, the heat flux vector J in a specific group can be calculated by summing the contributions from all the atoms in this group39 

J=1V[ieiviiSivi],
(5)

where V is the volume of the specific group, i is the atom index, ei is the energy, i.e., potential and kinetic energy, of atom i, vector vi is the velocity of atom i, and Si is the stress tensor for atom i.

Since heat flows from heat source to heat sink in our simulations and atoms in group l are not bonded to other groups, we reasonably assume the heat obtained by this group l is mainly due to the non-bonded interactions with atoms in the other two groups. Therefore, to study the role of non-bonded interactions in our model, we excluded specific pathways between different groups, relaxed the system in 0.2 ns NVT ensemble at 300 K (temperature damping parameter is 100 fs) followed by 0.2 ns NVE ensemble, and calculated the contributions from each group to heat flux and the thermal conductance of the interlamellar region after performing 5 ns non-equilibrium simulations.

Seven representative configurations were considered as summarized in Table II: (i) 18%-crystallinity with 2 bridges (Conf. i), (ii) 50%-crystallinity with 3 bridges (Conf. ii), (iii) 50%-crystallinity with 7 bridges (Conf. iii), (iv) 72%-crystallinity with 12 bridges (Conf. iv), (v) 83%-crystallinity with 20 bridges (Conf. v), (vi) 83%-crystallinity with 23 bridges (Conf. vi), and (vii) 83%-crystallinity with 27 bridges (Conf. vii).

TABLE II.

Representative configurations for the heat flux contribution calculations.

ConfigurationsConf. iConf. iiConf. iiiConf. ivConf. vConf. viConf. vii
Crystallinity (%) 18 50 50 72 83 83 83 
Bridge No. 12 20 23 27 
ConfigurationsConf. iConf. iiConf. iiiConf. ivConf. vConf. viConf. vii
Crystallinity (%) 18 50 50 72 83 83 83 
Bridge No. 12 20 23 27 

Figure 5 shows the comparison among the contributions from different groups to heat flux. When we excluded pathways P2 and P3, i.e., non-bonded interactions between groups u and l, and that between groups l and b, heat in group u can only be transferred by non-bonded interactions within group b. Figure 5(a) shows the corresponding contributions to the total heat flux from groups b and u in different configurations. When n = 2, X =18% (region I), the contribution from loops and tails is around 30%, demonstrating the importance of pathway P1. Therefore, the interlamellar region can be treated roughly as a pure amorphous region. Although group u is not bonded to group b, the large number of atoms within the non-bonded interaction cut-off distance ensures this pathway competitive with those along bridges, which is consistent with the results in Ref. 55. To prove this idea visually, snapshots with each chain individualized from Visual Molecule Dynamics (VMD)56 are shown in Fig. 6(a). In contrast, with the increase of n, contributions from bridges become dominant, especially when n20, where contributions from bridges are even larger than 97%. Thus, the modified Takayanagi's model has to be introduced in these configurations.

FIG. 5.

Comparison among the contributions from different channels to heat flux. (a) Contributions from groups b and u to the total heat flux when excluding pathways P2 and P3 in seven configurations listed in Table II. (b) Contributions from groups b and l to the total heat flux when excluding pathways P1 and P2 in seven configurations listed in Table II. (c) Contributions from groups b and u + l to the total heat flux when excluding pathways P1 and P3 in seven configurations listed in Table II. (d) Estimated Kb when excluding all the pathways in representative configurations at X = 50%, 72%, and 83%.

FIG. 5.

Comparison among the contributions from different channels to heat flux. (a) Contributions from groups b and u to the total heat flux when excluding pathways P2 and P3 in seven configurations listed in Table II. (b) Contributions from groups b and l to the total heat flux when excluding pathways P1 and P2 in seven configurations listed in Table II. (c) Contributions from groups b and u + l to the total heat flux when excluding pathways P1 and P3 in seven configurations listed in Table II. (d) Estimated Kb when excluding all the pathways in representative configurations at X = 50%, 72%, and 83%.

Close modal
FIG. 6.

VMD snap shots. (a) Representative snapshots from VMD to show the distance among atoms in groups u, l, and b. (Left, black atoms are in a bridge from group b, blue atoms are in a loop from group l; Middle, black atoms are in a bridge from group b, purple atoms are in a tail from group u; Right, blue atoms are in a loop from group l, purple atoms are in a tail from group u.) (b) Representative snapshot from VMD to show the configuration of bridges at X = 50% (left), 72% (middle), and 83% (right). Apparently, bridges at X = 83% are much aligned.

FIG. 6.

VMD snap shots. (a) Representative snapshots from VMD to show the distance among atoms in groups u, l, and b. (Left, black atoms are in a bridge from group b, blue atoms are in a loop from group l; Middle, black atoms are in a bridge from group b, purple atoms are in a tail from group u; Right, blue atoms are in a loop from group l, purple atoms are in a tail from group u.) (b) Representative snapshot from VMD to show the configuration of bridges at X = 50% (left), 72% (middle), and 83% (right). Apparently, bridges at X = 83% are much aligned.

Close modal

Similarly, if we excluded the pathways P1 and P2—the non-bonded interactions between groups u and l, and that between groups u and b—heat from group l is totally transferred by non-bonded interactions with group b. As shown in Fig. 5(b), similar results were observed, proving the effect of the non-bonded interactions between groups u and b and that between groups u and l are almost the same. We also simulated the case when the pathways P1 and P3 (non-bonded interactions between groups u and b, and that between groups l and b) were excluded. In this case, heat from group u can be transferred to group l through non-bonded interactions, but there is no non-bonded interaction between group b and the other two groups. As shown in Fig. 5(c), when n is small, the contribution to heat flux from groups u and l is even larger than that from group b. However, the contribution from group b increases with n and becomes dominant at high n number. Especially, when n20, the contribution from group b is larger than 96%. Such a dominant role of bridges at high crystallinity is consistent with the results in Ref. 48, where they found that the strong covalent bonds dominate in the heat transfer compared with the nonbonded interactions.

In the following, we focus on the role of bridges of configurations in regions II and III. When all the non-bonded pathways among groups l, u, and b were excluded, heat is mainly transferred along bridges, and therefore, the thermal conductivity of bridges can be roughly estimated. One thing to note here is, at X =50% and 72%, their bridges are much coiled, excluding non-bonded interactions could introduce a slight difference in the chain alignment, and thus our results only show a rough estimation of Kb [see Fig. 5(d)]. Here, Kb is the intrinsic thermal conductivity of group b by substituting A in Eq. (2) with the effective cross-sectional area (n48A). When n is in region II, Kb fluctuates around 0.9 W/(m K), while for n in region III, Kb increases linearly with a slope of 0.076. Thus, we calculated the effective thermal conductivity using our modified model with Kb=0.9±0.1W/(mK) for region II and Kb=0.53+0.076n for region III.

We attribute such a transition of Kb from region II to region III to two reasons: (i) bridges at X =83% are much aligned compared with those at X =50% and 72%. Figure 6(b) shows the VMD snap shot. Namely, coiled bridges at X =50% and 72% offer more phonon scattering centers; the increased number of scattering centers and the increased bridge number affect Kb in an opposite and competitive manner, resulting in Kb fluctuating around a constant value; (ii) at X =82%, with the large increase of bridge number, the role of non-bonded interaction becomes negligible, and most of the heat is transferred along chains through bond interactions. VMD snapshot also shows the more aligned bridge chains at high X, which is expected since the bridge chains tend to align when the interlamellar region is short (high X) and more bridge chains (large n) are present there. As simulations57 have pointed out that thermal conductivity along polymer chains exponentially scales with the chain orientation, such a less-coiled configuration could greatly enhance the heat transfer capability.

Finally, we test our modified model by predicting the thermal conductivity of PE with a crystallinity X > 20%, which is plotted in Fig. 3(b) to directly compare with the NEMD calculated data. The thermal conductivity values used in the model were from separate simulations: Ka=0.12±0.02W/(m-K), Kc=21±2W/(m-K),Kb=0.9±0.1W/(m-K) for X =50% and 72%, Kb varies with n as Kb=0.53+0.076n with ±10% uncertainty for X =83%. We found that the predictions from our modified model fit excellently well with the effective thermal conductivity calculated from NEMD, which proves that this modified model captures the heat transfer pathways in the semicrystalline PE and the important role of bridges in the interlamellar region. In all the previous experimental studies of semicrystalline polymers, due to the challenges of defining and measuring the parameter b, only crystallinity X is considered for the data fitting and explanation. Our study shows that at high crystallinity, such fittings could give large errors (e.g., ∼70%) and how to treat the interlamellar bridges is the key for the modeling.

In summary, we investigated the effective thermal conductivity of semicrystalline PE with crystallinity varying in a large range, from 18% to 83%. Our results clearly show that the predictions from Choy-Young's model and the series model for the thermal conductivity value of semicrystalline polyethylene fail, especially at high crystallinity. We attributed this deviation to the substantial role of bridges in carrying the heat by studying their contributions to the heat flux in the interlamellar region. We modify Takayanagi's model to take into account the effective thermal conductivity of bridges, and this modified model successfully predicts the effective thermal conductivity of semicrystalline PE at intermediate and high crystallinity. Our results prove that there is a transition from the intermediate crystallinity (X =50% and 72%) to high crystallinity (X =83%) due to the non-bonded interactions and the alignment of bridges, and suggest that both crystallinity and the bridge number in the interlamellar region need to be considered to explain the thermal conductivity of polymers with high crystallinity. As a preliminary study, the simple one-dimensional model in this study could contribute to the understanding of thermal transport in semicrystalline polymers. We expect a structure-property guide based on a more sophisticated model that can describe the true structure of high-crystallinity polymer fibers or films in the future work. Our work could shed some light on the development of plastic heat sinks and the thermal management in flexible electronics.

The authors would like to thank Gregory C. Rutledge and Raghavan Ranganathan for their helpful discussions on building the semicrystalline PE model. K.K., T.L., and J.L. acknowledge the financial support from the NC State Faculty Research and Professional Development Funds. T.L. also acknowledges the funding support from China Scholarship Council. X.L. acknowledges the National Natural Science Foundation of China (NSFC) Grant No. 51506062. J.Z. acknowledges the program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning Grant No. TP2014012. G.C. acknowledges DOE-BES Award No. DE-FG02-02ER45977. Computational resources were provided by the High Performance Computing Center at North Carolina State University.

The united atom (UA) Force field was used in the EMC modeling and MD simulations consistently. The parameters are listed in Table III.

TABLE III.

Force field parameters used in the EMC modeling and MD simulations.

Non-bondEn=4ϵσr12σr6,r<rc
TypeϵKcalmole1σÅrcÅ
CH2-CH2 0.09339913958 4.009 10.0225 
CH3-CH3 0.226542543 4.009 10.0225 
Non-bondEn=4ϵσr12σr6,r<rc
TypeϵKcalmole1σÅrcÅ
CH2-CH2 0.09339913958 4.009 10.0225 
CH3-CH3 0.226542543 4.009 10.0225 
BondEb=Krr02
TypeKKcal·mole1·Å2r0Å
CH2-CH2 449.5 1.53  
CH2-CH3 449.5 1.53  
CH3-CH3 449.5 1.53  
BondEb=Krr02
TypeKKcal·mole1·Å2r0Å
CH2-CH2 449.5 1.53  
CH2-CH3 449.5 1.53  
CH3-CH3 449.5 1.53  
AngleEa=Kθθ02
TypeKKcal·mole1·rad2θ0deg
CH2-CH2-CH2 60 112  
CH2-CH2-CH3 60 112  
CH3-CH2-CH3 60 112  
AngleEa=Kθθ02
TypeKKcal·mole1·rad2θ0deg
CH2-CH2-CH2 60 112  
CH2-CH2-CH3 60 112  
CH3-CH2-CH3 60 112  
DihedralEd=12K11+cosϕ+12K21cos2ϕ+12K31+cos3ϕ
TypeK1Kcal·mole1K2Kcal·mole1K3Kcal·mole1
CH2-CH2-CH2-CH2 1.6 −0.867 3.24 
CH2-CH2-CH2-CH3 1.6 −0.867 3.24 
CH3-CH2-CH2-CH3 1.6 −0.867 3.24 
DihedralEd=12K11+cosϕ+12K21cos2ϕ+12K31+cos3ϕ
TypeK1Kcal·mole1K2Kcal·mole1K3Kcal·mole1
CH2-CH2-CH2-CH2 1.6 −0.867 3.24 
CH2-CH2-CH2-CH3 1.6 −0.867 3.24 
CH3-CH2-CH2-CH3 1.6 −0.867 3.24 

The regions of crystalline, interphase, and amorphous were distinguished according to the density profile illustrated in Fig. 7. The average density of the crystalline region is 1g/cc and that of the amorphous region is 0.87g/cc, which are reasonable compared with the experimental results. The drop in the density from crystalline to amorphous occurs in a transition region which is defined as the interphase region. The thickness of the interphase is estimated as 1nm.58 

FIG. 7.

A representative density profile (density as a function of z) at X =50%.

FIG. 7.

A representative density profile (density as a function of z) at X =50%.

Close modal

The local effective thermal conductivity is defined as the ratio of heat flux to the local temperature gradient. The local effective thermal conductivity KL of representative samples was fitted. Four adjacent data points on each side of the fitted point were adopted for fitting. While for the fitted point close to the interface, one adjacent data point on each side was used for fitting. Figure 8 shows the KL in configurations with X =72% and 83%. The KL in the interphase region behaves as a transitional region between the crystalline and the amorphous region. Moreover, comparing Figs. 8(a) and 8(b), the amorphous region has a higher KL with higher crystallinity because of the increasing bridge number and better alignment.

FIG. 8.

Representative local effective thermal conductivity profiles at (a) X =72% and (b) X =83%. The heat baths were employed on both sides, shown as the shaded areas. The red shaded area denotes the heat source, and the blue shaded area denotes the heat sink. The symbols c, I, and a represent crystalline, interphase, and amorphous regions that are distinguished according to the density profiles shown in Fig. 7.

FIG. 8.

Representative local effective thermal conductivity profiles at (a) X =72% and (b) X =83%. The heat baths were employed on both sides, shown as the shaded areas. The red shaded area denotes the heat source, and the blue shaded area denotes the heat sink. The symbols c, I, and a represent crystalline, interphase, and amorphous regions that are distinguished according to the density profiles shown in Fig. 7.

Close modal
1.
W. U.
Huynh
,
J. J.
Dittmer
, and
A. P.
Alivisatos
,
Science
295
(
5564
),
2425
2427
(
2002
).
2.
J. Y.
Kim
,
K.
Lee
,
N. E.
Coates
,
D.
Moses
,
T.-Q.
Nguyen
,
M.
Dante
, and
A. J.
Heeger
,
Science
317
(
5835
),
222
225
(
2007
).
3.
S. N.
Patel
,
A. M.
Glaudell
,
K. A.
Peterson
,
E. M.
Thomas
,
K. A.
O'Hara
,
E.
Lim
, and
M. L.
Chabinyc
,
Sci. Adv.
3
(
6
),
e1700434
(
2017
).
4.
T.
Janoschka
,
N.
Martin
,
U.
Martin
,
C.
Friebe
,
S.
Morgenstern
,
H.
Hiller
,
M. D.
Hager
, and
U. S.
Schubert
,
Nature
527
(
7576
),
78
(
2015
).
5.
D. J.
David
and
A.
Misra
,
Relating Materials Properties to Structure with MATPROP Software: Handbook and Software for Polymer Calculations and Materials Properties
(
CRC Press
,
2001
).
6.
X.
Xie
,
D.
Li
,
T.-H.
Tsai
,
J.
Liu
,
P. V.
Braun
, and
D. G.
Cahill
,
Macromolecules
49
(
3
),
972
978
(
2016
).
7.
X.
Xie
,
K.
Yang
,
D.
Li
,
T.-H.
Tsai
,
J.
Shin
,
P. V.
Braun
, and
D. G.
Cahill
,
Phys. Rev. B
95
(
3
),
035406
(
2017
).
8.
G.-W.
Lee
,
M.
Park
,
J.
Kim
,
J. I.
Lee
, and
H. G.
Yoon
,
Composites, Part A
37
(
5
),
727
734
(
2006
).
9.
K. I.
Winey
,
T.
Kashiwagi
, and
M.
Mu
,
MRS Bull.
32
(
4
),
348
353
(
2007
).
10.
Y. P.
Mamunya
,
V.
Davydenko
,
P.
Pissis
, and
E.
Lebedev
,
Eur. Polym. J.
38
(
9
),
1887
1897
(
2002
).
11.
Z.
Zhong
,
M. C.
Wingert
,
J.
Strzalka
,
H. H.
Wang
,
T.
Sun
,
J.
Wang
,
R.
Chen
, and
Z.
Jiang
,
Nanoscale
6
(
14
),
8283
8291
(
2014
).
12.
J.
Liu
and
R.
Yang
,
Phys. Rev. B
86
(
10
),
104307
(
2012
).
13.
S.
Shen
,
A.
Henry
,
J.
Tong
,
R.
Zheng
, and
G.
Chen
,
Nat. Nanotechnol.
5
(
4
),
251
255
(
2010
).
14.
A.
Henry
and
G.
Chen
,
Phys. Rev. Lett.
101
(
23
),
235502
(
2008
).
15.
X.
Wang
,
V.
Ho
,
R. A.
Segalman
, and
D. G.
Cahill
,
Macromolecules
46
(
12
),
4937
4943
(
2013
).
16.
B.
Zhu
,
J.
Liu
,
T.
Wang
,
M.
Han
,
S.
Valloppilly
,
S.
Xu
, and
X.
Wang
,
ACS Omega
2
(
7
),
3931
3944
(
2017
).
17.
See http://www.goodfellow.com/E/Polyethylene-UHMW-Fibre.html for the thermal conductivity of commercial ultra-drawn single PE microfibers (spectra 2000).
18.
H.
Ghasemi
,
N.
Thoppey
,
H.
Xiaopeng
,
J.
Loomis
,
L.
Xiaobo
,
J.
Tong
,
W.
Jianjian
, and
G.
Chen
, in
2014 IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm)
(IEEE,
2014
), pp.
235
239
.
19.
D. K. Y.
Xu
,
B.
Song
,
Z.
Jiang
,
J.
Zhou
,
J.
Loomis
,
J.
Wang
,
M.
Li
,
H.
Ghasemi
,
X.
Huang
,
X.
Li
, and
G.
Chen
, “
Nanostructured polymer films with metal-like thermal conductivity
,” preprint arXiv:1708.06416 (
2017
).
20.
21.
P. J.
in't Veld
and
G. C.
Rutledge
,
Macromolecules
36
(
19
),
7358
7365
(
2003
).
22.
C.
Choy
and
K.
Young
,
Polymer
18
(
8
),
769
776
(
1977
).
23.
J.
Hennig
,
J. Polym. Sci., Part C: Polym. Symp.
16
,
2751
(
1967
).
24.
I. M.
Ward
and
J.
Sweeney
,
Mechanical Properties of Solid Polymers
(
John Wiley & Sons
,
2012
).
25.
M.
Takayanagi
,
K.
Imada
, and
T.
Kajiyama
,
J. Polym. Sci. C
15
,
263
(
1967
).
26.
J.
Affdl
and
J.
Kardos
,
Polym. Eng. Sci.
16
(
5
),
344
352
(
1976
).
27.
H.
Cox
,
Br. J. Appl. Phys.
3
(
3
),
72
(
1952
).
28.
A.
Gibson
,
G.
Davies
, and
I.
Ward
,
Polym. Eng. Sci.
20
(
14
),
941
948
(
1980
).
30.
S.
Lee
and
G. C.
Rutledge
,
Macromolecules
44
(
8
),
3096
3108
(
2011
).
31.
P. J.
in 't Veld
,
M.
Hütter
, and
G. C.
Rutledge
,
Macromolecules
39
(
1
),
439
447
(
2006
).
32.
J. M.
Kim
,
R.
Locker
, and
G. C.
Rutledge
,
Macromolecules
47
(
7
),
2515
2528
(
2014
).
33.
I.-C.
Yeh
,
J. L.
Lenhart
,
G. C.
Rutledge
, and
J. W.
Andzelm
,
Macromolecules
50
(
4
),
1700
1712
(
2017
).
34.
W.
Paul
,
D. Y.
Yoon
, and
G. D.
Smith
,
J. Chem. Phys.
103
(
4
),
1702
1709
(
1995
).
35.
I.-C.
Yeh
,
J. W.
Andzelm
, and
G. C.
Rutledge
,
Macromolecules
48
(
12
),
4228
4239
(
2015
).
36.
S.
Gautam
,
S.
Balijepalli
, and
G.
Rutledge
,
Macromolecules
33
(
24
),
9136
9145
(
2000
).
37.
D.
Bassett
and
A.
Hodge
,
Proc. R. Soc. A
377
,
25
37
(
1981
).
38.
G.
Dee
,
T.
Ougizawa
, and
D.
Walsh
,
Polymer
33
(
16
),
3462
3469
(
1992
).
39.
S.
Plimpton
,
J. Comput. Phys.
117
(
1
),
1
19
(
1995
).
40.
C.
Oligschleger
and
J.
Schön
,
Phys. Rev. B
59
(
6
),
4125
(
1999
).
41.
P.
Schelling
,
S.
Phillpot
, and
P.
Keblinski
,
J. Appl. Phys.
95
(
11
),
6082
6091
(
2004
).
42.
A.
Henry
,
G.
Chen
,
S. J.
Plimpton
, and
A.
Thompson
,
Phys. Rev. B
82
(
14
),
144308
(
2010
).
43.
X.
Wei
,
T.
Zhang
, and
T.
Luo
,
Phys. Chem. Chem. Phys.
18
(
47
),
32146
32154
(
2016
).
44.
Y.
Lu
,
J.
Liu
,
X.
Xie
, and
D. G.
Cahill
,
ACS Macro Lett.
5
(
6
),
646
650
(
2016
).
45.
A.
Henry
and
G.
Chen
,
Nanoscale Microscale Thermophys. Eng.
13
(
2
),
99
108
(
2009
).
46.
D.
Luo
,
C.
Huang
, and
Z.
Huang
,
J. Heat Transfer
140
(
3
),
031302
(
2017
).
47.
J.
Shiomi
,
Annu. Rev. Heat Transfer
17
,
177
203
(
2014
).
48.
T.
Zhang
and
T.
Luo
,
J. Phys. Chem. B
120
(
4
),
803
812
(
2016
).
49.
E.
Lussetti
,
T.
Terao
, and
F.
Müller-Plathe
,
J. Phys. Chem. B
111
(
39
),
11516
11523
(
2007
).
50.
R. J.
Crawford
,
Plastics Engineering
(
Butterworth-Heinemann
,
1998
).
51.
H.
Guo-Jie
,
C.
Bing-Yang
, and
L.
Yuan-Wei
,
Chin. Phys. Lett.
31
(
8
),
086501
(
2014
).
52.
M. D.
Losego
,
M. E.
Grady
,
N. R.
Sottos
,
D. G.
Cahill
, and
P. V.
Braun
,
Nat. Mater.
11
(
6
),
502
506
(
2012
).
53.
M.
Hu
,
P.
Keblinski
, and
P. K.
Schelling
,
Phys. Rev. B
79
(
10
),
104305
(
2009
).
54.
P. J.
O'Brien
,
S.
Shenogin
,
J.
Liu
,
P. K.
Chow
,
D.
Laurencin
,
P. H.
Mutin
,
M.
Yamaguchi
,
P.
Keblinski
, and
G.
Ramanath
,
Nat. Mater.
12
(
2
),
118
122
(
2013
).
55.
V.
Rashidi
,
E. J.
Coyle
,
K.
Sebeck
,
J.
Kieffer
, and
K. P.
Pipe
,
J. Phys. Chem. B
121
(
17
),
4600
4609
(
2017
).
56.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
,
J. Mol. Graphics
14
(
1
),
33
38
(
1996
).
57.
J.
Liu
and
R.
Yang
,
Phys. Rev. B
81
(
17
),
174122
(
2010
).
58.
A.
Ghazavizadeh
,
G. C.
Rutledge
,
A. A.
Atai
,
S.
Ahzi
,
Y.
Rémond
, and
N.
Soltani
,
J. Polym. Sci., Part B: Polym. Phys.
51
(
23
),
1692
1704
(
2013
).