Introduction

Polymer fibers have been widely explored in many engineering applications due to the advantages of their anti-corrosion, high-toughness, and low-density features. However, bulk polymers, either (semi)-crystalline or amorphous, exhibit low thermal conductivities (TCs) of only 0.1–0.3 W m−1 K−1 at room temperature.1 Such poor performance in heat conduction limits their applications in electronic packaging,2 aerospace materials,3 and automotive parts.4 The low TC in bulk polymers have been attributed to various fundamental mechanisms, including the disordered arrangement of polymer chains in the amorphous phase, the presence of grain boundaries between polycrystalline domains, inter-chain entanglements, nanovoids, and other defects and impurities, all of which result in significant phonon scatterings that hinder heat conduction.5,6,7,8 Nevertheless, fully extended (with all the dihedrals in the trans state) single polyethylene (PE) chain has been proved to possess extremely high TC in the order of 100 W m−1 K−1 through all-atomistic molecular dynamics (AA-MD) simulations.9 This highlights the potential of highly efficient heat conduction in bulk polymers, if all the chains are extended and aligned along the heat conducting direction.

Numerous experimental studies have shown that tensile elongation in the thermal-drawing process can effectively enhance the TC in bulk polymer fibers. Shen et al.10 fabricated high-quality ultra-drawn polyethylene (PE) nanofibers with room-temperature TC up to 104 W m−1 K−1 at a draw ratio of 400 due to the restructuring of polymer chains towards alignment and crystallization. Single crystalline PE mats with a draw ratio of λ = 350 were also found to exhibit excellent TC in the order of 40 W m−1 K−1, attributed to the presence of ~85% of long needle-shaped crystals.11 Recently, room-temperature PE microfibers12 with a high TC of 51 W m−1 K−1, as well as PE nanofibers13 with both ultra-high strength of 11 GPa and TC of around 50 W m−1 K−1 were reported. The current strategies to align individual PE chains and increase local crystallinity have been proven successful in enhancing the TC of polymeric materials; however, most studies rely on trial-and-error experimentations with few theoretical guidance. In particular, the optimization of the thermal-drawing protocols, including the annealing profile, draw ratio (tensile strain), and draw rate (strain rate), remains largely unexplored at the molecular level.

To complement experimental studies, AA-MD simulations have been widely applied to study chain alignment and heat conduction in many polymer systems.14,15,16,17,18,19,20,21 Using AA-MD simulations, Liu and Yang22 show that the TC of bulk PE increases to 1.3 W m−1 K−1 at a strain of 2, and interestingly, such enhancement in TC is more significant when the strain rate is smaller. They found that the alignment of polymer chains, measured by the Herman’s orientational order parameter, also increases with strain, which correlates well with the increase in TC. By stretching amorphous polymer fibers at 500 K, followed by quenching at 300 K, Zhang and Luo23 achieved a 176-fold increase in TC at a strain of ~9 in AA-MD simulations. However, the TC vs. strain predicted that there does not reach a plateau as observed in experiments,11,24 but instead exhibits singularity at strains larger than 10. This is probably due to the fact that the chains are covalently bonded across the periodic boundary of the simulation cell, which would undergo unphysical backbone stretching at high strains and lose the opportunity for inter-chain sliding or detachment to release the residue stress. Zhang and Luo25 also investigated the morphology-dependent TC in amorphous PE by artificially tuning the chain stiffness. They found that the enhanced crystallinity and chain alignment along the heat conducting direction can lead to an increase in TC by more than one order of magnitude compared to amorphous PE.25

In spite of the above experimental and theoretical advances, the impacts from the fiber thermal-drawing protocols on polymer chain alignment, crystallinity, and the resulting phonon transport properties remain obscure. In particular, there are great interests for a scalable modeling framework to cover both thermal phonon transport properties of individual polymer chains at the nanoscale and morphological transitions during the thermal-drawing process at the mesoscale. For example, Liu and Yang22 simulated individual PE chains consisting of up to 400 carbon atoms, corresponding to a chain molecular weight of 5602 g mol−1, while the typical molecular weight of high-density PE is up to 106 g mol−1,26 which would have a critical influence on the TC of polymers.27,28 From experiments29,30 and simulations,28 it is well known that longer polymer chains could lead to larger TCs, if the polymers are well aligned. This is probably due to the resulting longer thermal pathways through the chain, along which phonons can transport with less scattering and longer mean free path (MFP).31 In addition to the length scale above, the time scale involved in the fiber thermal-drawing process is also beyond the scope of conventional AA-MD simulations. Such time-scale limitation in AA-MD requires very high strain rates to be applied, which would inhibit relaxations of residue stresses in polymers. A scalable modeling framework, enabled by coarse-grained (CG) MD simulations, could greatly extend the length and time scales of conventional AA-MD simulations by two orders of magnitude.32,33,34,35 However, to the best of our knowledge, the application of CG-MD simulations to precisely model heat conduction in thermally drawn bulk polymer fibers has not been reported yet.

In this work we develop a predictive CG force field (Fig. 1a) for bulk PE (up to a relatively high molecular weight of 10,500 g mol−1), capable of reproducing the correct thermal transport and phase transition features. This is achieved through target property matching and the multistate-iterative Boltzmann inversion (MS-IBI) method36 with our AA-MD simulation results serving as the benchmark. Utilizing this CG force field, we conduct thermal-drawing MD simulations and thermal-conductivity calculations on bulk PE fibers to study the effects of thermal-drawing protocols, including the stretching (drawing) temperature, tensile strain level, and strain rate, on the observed enhancement in room-temperature TCs. PE chain alignment and crystallinity (global and local) are quantified to discover the correlation between PE morphology (Fig. 1b) and TC.

Fig. 1: Simulation model and procedure.
figure 1

a Schematic of the coarse-grained (CG) model for PE, in which every five methylene groups (-CH2-) are represented by one bead. b CG-MD simulation snapshots of PE fibers thermally strained at 400 K under a strain rate of 5 × 106 s−1 under different strains after stress relaxation. c Flowchart of the simulated fiber thermal-drawing protocol with details discussed in Methods.

Full size image

More interestingly, we discover that a combination of moderate drawing temperature and strain rate is optimal in enhancing TCs by matching the temperature-dependent Rouse viscoelastic relaxation time and the applied strain rate. Such optimal thermal-drawing protocol is condensed to similar Deborah numbers (only differ across the melting temperature) over different drawing temperature and strain rate combinations considered here. Furthermore, phonon transport properties of thermally drawn PE fibers are computed to reveal the molecular mechanism behind the enhanced thermal transport. Using the effective-medium-theory model developed here based on serial-parallel thermal resistors, we show that both chain alignment and crystallinity are important to TCs. Finally, we find that the greatly enhanced TC along the fiber-drawing direction is mainly attributed to the increased intrinsic phonon MFP, as well as the increased longitudinal phonon group velocities resulting from the enlarged elastic modulus along that direction.

Results and discussion

Chain alignment and degree of crystallinity of thermally drawn PE fibers

To reveal the impacts from fiber thermal-drawing conditions on polymer morphology and to develop a process–structure relationship, we investigated the chain alignment and crystallinity after thermal strains (see Methods). Figure 2 shows the orientational order parameters P2z over all the beads after completing the thermal-drawing process outlined in Fig. 1c. At a drawing temperature of 500 K and a strain rate of 5 × 107 s−1, P2z increases with increasing strain, reflecting more chain alignments (see the probability distribution of P2z in Supplementary Fig. 1), and then reaches a plateau of ~0.8. However, under other conditions (e.g., a drawing temperature of 400 K and a strain rate of 5 × 107 s−1), P2z first increases and then decreases with increasing strain. This is different from the trend right after thermal drawing and before removing the external stresses, where P2z increases with increasing the strain, and then reaches a plateau (Supplementary Fig. 2c), indicating partial reversibility of the thermal-drawing process upon removing the external stresses. The decrease in P2z is probably attributed to excessive stresses built up under high strains (Supplementary Fig. 3b), the release of which results in structural shrinking and chain alignment destruction after external stresses (or strains) are removed. This also results from the increased volumes of nanovoid-type vacancies (Supplementary Fig. 4d), which provides room for chain folding after removing the strains from high strain levels.

Fig. 2
figure 2

The change of orientational order parameter of bulk PE as a function of strain obtained from different drawing temperatures and a strain rate of 5 × 107 s−1 a and a drawing temperature of 400 K b and 500 K c and different strain rates followed by a 20-ns relaxation at 300 K.

Full size image

In addition, P2z also depends nonlinearly on the drawing temperature under the same strain rate. Figure 2a shows that P2z reaches the highest values at intermediate drawing temperatures considered here (i.e., 400 and 500 K). The tensile stress decreases with increasing the drawing temperature under the same strain level and strain rate (Supplementary Fig. 3a), reaching a steady value above 500 K. On one hand, excessive stress concentration and vacancy formation during the tensile loading at lower drawing temperatures (e.g., 300 K) are both harmful to chain alignment and would diminish P2z (Supplementary Fig. 4b). On the other hand, higher tensile forces (stress) are still required to stretch polymer chains and compete against entropic chain randomization at higher temperatures.37 Therefore, lower stresses at higher drawing temperatures (i.e., 600 K) fails to provide enough tensile forces to align the chains well, resulting in lower P2z.

Furthermore, strain rate also exhibits strong nonlinear influence on P2z at the same drawing temperature. For example, under a drawing temperature of 400 K, P2z increases monotonically with decreasing the strain rate down to 5 × 106 s−1 (the smallest strain rate considered here due to the computational cost, see Fig. 2b). Under a lower drawing temperature of 300 K, the PE fiber suffers from early fracture at small tensile strains of ~2, and therefore fails to align the polymer chains well (see the P2z in Supplementary Fig. 5a). Interestingly, under a drawing temperature of 500 K, P2z reaches a maximum under an intermediate strain rate of 2 × 108 s−1 (see Fig. 2c). A similar trend is observed under a higher drawing temperature of 600 K (Supplementary Fig. 5b), showing maximum P2z under an intermediate strain rate of 5 × 108 s−1. It is noteworthy that the optimal strain rate to achieve the maximum P2z increases from 5 × 106 to 5 × 108 s−1 when increasing the drawing temperature from 400 to 600 K.

We further examined the impacts from fiber-drawing conditions on the resulting polymer crystallinity, which, interestingly, do not match exactly the trend observed in chain alignment P2z, especially under larger strains. Similar to P2z, the system-averaged degree of crystallinity Ψ firstly increases with increasing the strain with local crystallites distributed uniformly in the bulk polymer (see Fig. 3a). Similar to P2z, Ψ reaches a maximum of 77% at an intermediate drawing temperature of 400 K, as well as an intermediate strain rate of 5 × 107 s−1 (see Fig. 3b). Specifically, the degree of crystallinity increases with decreasing the strain rate down to 5 × 106 s−1 at a drawing temperature of 400 K (Fig. 3c). Similarly, the maximum degree of crystallinity also increases with decreasing the strain rate down to 5 × 106 s−1 under a drawing temperature of 300 K (Supplementary Fig. 6a). Note that it is challenging to further decrease the strain rate due to the high computational cost, and therefore we are not able to determine the optimal strain rate to maximize the degree of crystallinity under a drawing temperature of 300 K. The degrees of crystallinity reach maxima under strain rates of 2 × 108 and 5 × 108 s−1, respectively, at drawing temperatures of 500 and 600 K (Fig. 3d, e).

Fig. 3: Local crystallinity and system-averaged degrees of crystallinity Ψ in thermally drawn PE fiber.
figure 3

a CG-MD simulation snapshots showing local crystallinity in thermally drawn PE fiber at strains ε = 0, 2, 4, 6, and 10. Samples are obtained from a strain rate of 5 × 106 s−1 and a drawing temperature of 400 K. The zoomed-in snapshots in a illustrate the morphology difference between the amorphous and crystalline domains. Each polymer bead is regarded as a crystalline bead when the local degree of crystallinity is >0.7. For strained samples, only the central regions of the simulation boxes are displayed. b System-averaged degrees of crystallinity Ψ with increasing the strain under a fixed strain rate of 5 × 107 s−1 and various drawing temperatures. c–e System-averaged degrees of crystallinity Ψ with increasing the strain under various strain rates and constant drawing temperatures of 400, 500, and 600 K.

Full size image

The optimal strain rate to achieve the maximum crystallinity Ψ increases from 5 × 106 to 5 × 108 s−1 when increasing the drawing temperature from 400 to 600 K, which agrees with the trend observed in P2z (see Fig. 4). This is due to the fact that the alignment of polymer chains also enhances thread-like precursors responsible for the formation of the shish-kebab crystals,38 reflected from the further increase in Ψ after annealed stress relaxation, compared with the minor decrease in Ψ before annealed stress relaxation (Supplementary Fig. 7). Therefore, a combination of moderate drawing temperature and strain rate is essential to achieve both highest degrees of chain alignment and crystallinity. However, high polymer crystallinity does not guarantee perfect chain alignment (i.e., polymer chains in crystalline domains do not perfectly align) along the fiber-drawing direction, especially under high strain rates. This is reflected from the fact that Ψ does not drop off as much as P2z under high strain rates. The structural shrinkage and chain alignment destruction after removing excessive stresses under high strains, which results in the drop off of P2z, however, does not affect Ψ. Higher degrees of polymer crystallinity are expected to enhance their TC due to the higher TC in the crystalline phase (however, only along the chain direction) than that in the amorphous phase.25,39

Fig. 4
figure 4

The maximum strain-induced P2z a and degree of crystallinity Ψ b achieved after annealed stress relaxation as a function of fiber-drawing temperature and strain rate. The optimal strain rates to achieve maximum P2z and degree of crystallinity at 400, 500, and 600 K are 5 × 106, 2 × 108, and 5 × 108 s−1, respectively.

Full size image

Viscoelastic relaxation modulated chain alignment and crystallinity

From the viscoelastic mechanics aspect, the tensile stress increases with increasing the strain rate under the same strain, as expected from the limited extent of stress relaxation during fast fiber drawing (Supplementary Fig. 3). In addition to the higher local stresses, higher strain rate will also lead to larger volume of vacancy in the form of nanovoids (Supplementary Fig. 4), which will further destroy the chain alignment. On the other hand, lower strain rate suffers from extra stress relaxation and unwanted configurational randomization of polymer chains, and therefore cannot provide enough tensile elongation forces to untie the entangled chains and align them further. In general, thermal energy transport along strong covalent bonds in well-aligned individual chains is more effective than that via weak van der Waals (vdW) interactions between different chains. Therefore, an optimal combination of fiber-drawing temperature and strain rate is essential to best align polymer chains and achieve the maximum TC along the fiber-drawing direction.

To unfold the viscoelastic mechanism behind the observed optimal drawing temperature and strain rate, we predicted the viscoelastic relaxation times for undrawn bulk PE at the above drawing temperatures (see Methods). The computed stress relaxation modulus G(t) decays as a function of time and decays faster when increasing the temperature, as shown in Fig. 5. The plateau modulus of G(t), defined after the linear Rouse dynamics ends and within the entanglement dynamics region in G(t), is ~2.5 MPa at 400 K (see Supplementary Fig. 8) in this work, which matches the experimental value of 2.6 MPa at 413 K.40 The Rouse relaxation time τR,41 which measures the diffusion time needed for chains shorter than the entanglement length to release stresses, is defined here when G(t) decreases to ~107 Pa. This cutoff in G(t) is used to define τR because of the Rouse dynamics region, G(t) ~ t−0.5,37 and therefore the log(G(t)) vs. log(t) curves remain linear with a constant slope for each temperature case (see Fig. 5), until G(t) drops to ~10 MPa where the Rouse dynamics ends. Accordingly, τR for PE fibers drawn at 400, 500, and 600 K are determined to be 5 × 10−10, 1 × 10−10, and 5 × 10−11 s, respectively.

Fig. 5
figure 5

Time evolutions of stress relaxation moduli, G(t), for undrawn amorphous PE at various temperatures. The various stages of chain relaxation are denoted in the figure.

Full size image

To verify the computed Rouse relaxation time τR, we have converted the dimensionless τR of polymers, predicted by Likhtman et al.41 using MD simulations and dimensionless force field units (in which the dimensionless time is normalized by \(t^ \ast = \sigma \sqrt {\frac{m}{{k_{\mathrm{B}}T}}}\)), to the physical τR values presented using real force field units here. A dimensionless τR of 1000 for the 100-bead polymer melt at 365 K by Likhtman et al.41 is converted into a physical τR of 2.4 × 10−9 s based on the σ, m, and T values used here. This is within an order of magnitude from our computed τR values of 5 × 10−10 s at 400 K for the 150-bead PE melt in this work. The differences could be due to the differences in the temperatures applied, the molecular weights of the polymers simulated, and the functional forms used in the two force fields. Therefore, our τR predictions could be considered reasonably accurate when compared to other theoretical studies.

The decrease in τR correlates perfectly well with the increase in the optimal strain rates when increasing the drawing temperature from 400 to 600 K. This is due to the fact that when the stress relaxes faster, shorter processing time (i.e., higher strain rate) is required to lock stresses for aligning polymer chains under higher temperatures. Excessively low strain rates lead to a liquid-like polymer, which diffuses too fast into disordered structures, while excessively high strain rates lead to a solid-like polymer, which is easy to break and form nanovoids. These two competing factors indicate that the optimal combination of drawing temperature and strain rate, a delicate balance between stress localization and chain diffusion (for viscoelastic relaxation), is critical for the polymer thermal-drawing process.

With the above in mind, we utilize the Deborah number De,42 a dimensionless number measuring the competing effects between chain diffusion (which leads to viscoelastic relaxation) and stress localization (which helps chain alignments) rates, to condense the optimal drawing temperature (which is reflected by the corresponding viscoelastic relaxation time) and strain rate into a single value. De is defined here as the ratio between the temperature-dependent Rouse relaxation time τR and the thermal-drawing processing time τP. Note that τP, estimated by the simulation time when TC reaches the maximum under tensile loading, is a good measure of the stress localization rate and inversely correlated with the strain rate. To achieve optimal strain rates at drawing temperatures of 400, 500, and 600 K, τP (see Fig. 7b–d) are 2 × 10−6, 1 × 10−7, and 3.8 × 10−8 s, respectively. The corresponding values of De are 2.5 × 10−4, 1 × 10−3 and 1.3 × 10−3, respectively. Although we have not exhaustively tested other drawing temperature and strain rate combinations here, we could imply from the current data that a De in the order of 1 × 10−3 would lead to the maximum degrees of chain alignment and crystallinity for thermal-drawing above the melting temperature (>409 K for PE here, see Supplementary Discussion and Supplementary Fig. 19). For the thermal-drawing process below the melting temperature, a much smaller De in the order of 1 × 10−4 would be optimal to process a glassy polymer fiber due to significantly slower chain diffusions (measured by τR).

To validate the magnitude of the De measured in this work, we estimated and compared to De in Shrestha et al.’s experimental work.13 In that work, the PE microfiber was locally drawn under a strain rate of 1400 s−1 and a drawing temperature of 450 K. The strain (or draw-down ratio) of PE was estimated to be 225 because the diameter of the PE fiber deceases from 1.5 to 100 nm. Therefore, the processing time should be 225/1400 = 0.16 s. The molecular weight M of PE in Shrestha et al.’s work is 3–6 × 106 g mol−1. The Rouse relaxation time τR follows a power scaling law of τR ~ M2.37 Considering that our computed τR values are reasonably accurate, we then used the above-scaling law to estimate the τR value for PE with a molecular weight of M = 4.5 × 106 g mol−1 studied in Shrestha et al.’s work. The τR in Shrestha et al.’s work is estimated to be ~6 × 10−5 s, based on τR values of 5 × 10−10 and 1 × 10−10 s at 400 and 500 K, respectively, with M = 10,500 g mol−1 used here. As a result, De in Shrestha et al.’s work is estimated to be 6 × 10−5 s/0.16 s = 3.75 × 10−4, which matches reasonably well with our computed De of 2.5 × 10−4 and 1 × 10−3 at 400 and 500 K, respectively, to achieve the optimal fiber morphology and TC. These interesting findings suggest that it would take around 103 and 104 times of the Rouse relaxation time of polymers of different molecular weights, respectively, to optimally draw a polymer fiber in the melt and the glassy state.

Impacts from fiber-drawing conditions on TC

In general, the TC of bulk PE along the drawing direction increases with increasing the strain (a few exceptions are discussed below), as shown in Fig. 6a. Contributions to the total TC include those from the virial (i.e., bonded, angular, and nonbonded interactions here) and the convection (kinetic energy) terms (see Methods). Figure 6a shows that the contribution from bonded interactions to the total TC is dominant and also increases with increasing the strain, which is expected, as high-frequency (energy) phonons are carried by bond vibrations. The contribution from nonbonded interactions to the total TC ranks second, but is almost independent of strain. Contributions from angular interactions and convections to the total TC are both negligible. As the strain increases, the contribution from bonded interactions to the total TC increases from 69 to 97%, while that from nonbonded interactions decreases from 25 to only 2% (Fig. 6b). This finding is reasonable because the backbones of polymer chains reorient towards the drawing direction with increasing the strain, such that heat can transfer in parallel to the backbones more efficiently through stronger bonded interactions rather than weaker vdW interactions.

Fig. 6
figure 6

TC contributions from the virial (bonded, angular, and nonbonded interactions) and the convection terms a and their ratios with respect to the total TC b as a function of strain under a drawing temperature of 500 K and a strain rate of 5 × 107 s−1.

Full size image

In addition, we observe that the TC of bulk PE along the stretching direction correlates well with both the chain alignment P2z and the degree of crystallinity Ψ, since both morphological features are enhanced with increasing the strain, except for the cases under high strain rates. For example, as shown in Fig. 7a, under an intermediate strain rate of 5 × 107 s−1, TC reaches maximum at an intermediate drawing temperature of 500 K. A nonlinear trend between TC and the strain level is observed under high strain rates, when both P2z and Ψ drop off due to limited stress relaxation and the resulting internal structural defects. For example, TC decreases with further elongation at a drawing temperature of 400 K and all the strain rates considered here, as shown in Fig. 7b.

Fig. 7: TC of bulk PE along the drawing direction.
figure 7

a TCs of bulk PE as a function of strain under various drawing temperatures and a fixed strain rate of 5 × 107 s−1. TCs as a function of strain under various strain rates and fixed drawing temperatures of b 400 K, c 500 K, and d 600 K. The half black and white circles in a denote the model-predicted TC values based on the effective-medium theory developed here.

Full size image

Similar to the trend observed in chain alignment P2z and degree of crystallinity Ψ, moderate strain rate and drawing temperature are both essential to the improvement of TC, as shown in Fig. 7. Lower optimal strain rates for TC are achieved under lower drawing temperatures, a trend that can be explained by the delicate balance between viscoelastic relaxation and stress localization, where such optimal combination can be condensed into similar Deborah numbers (see above). Among all the drawing conditions considered here, TC reaches a maximum of 1.57 W m−1 K−1 at a strain level of 20, a drawing temperature of 500 K, and a strain rate of 2 × 108 s−1, which is 14 times higher than that of undrawn amorphous PE. Such TC achieved through thermal drawing is outstanding, reaching 80% of the theoretical limit of bulk PE, that is, its pure crystalline counterpart with TC κc∥ = 1.953 W m−1 K−1 parallel to the chains in a PE crystal (see below for details). It is important to note that achieving the theoretical TC limit of a single PE chain in bulk PE could be challenging due to excessive inter-chain phonon scatterings due to vdW-type weak nonbonded interactions.

However, neither P2z nor Ψ can correlate perfectly well with TC to yield a predictive model using only one of these two structure ordering parameters. For example, TC data scatter around with respect to P2z, although they do show some exponential increase with increasing P2z, which agree with the model proposed by Liu and Yang,22 but with many high TC outliers (Fig. 8a). Similarly, TC data also scatter around with respect to Ψ, showing weak increasing trend with increasing Ψ (Fig. 8b). Since both factors have significant impacts on the reported TC here along PE fibers, we need to develop a new model to correlate TC with both P2z and Ψ.

Fig. 8
figure 8

Variations of calculated TC as a function of a the orientational order parameter P2z and b the degree of crystallinity ψ. The model-predicted TC profile by Liu and Yang22 is shown as a dashed line in a for comparison.

Full size image

Based on the effective-medium theory,43 we developed a new model here for the TC of semicrystalline polymers under tensile deformation. Note that the Choy–Young’s model39 and the series model44,45 underestimate the TC at high degrees of crystallinity.46 The reason is that both models consider the TC of the amorphous domain as a constant. However, the TC of amorphous domain should change with strain due to the corresponding change in the orientational order parameters of PE chains in the amorphous domain (see Supplementary Fig. 9). In addition, the orientational distribution of PE chains along the drawing direction is not uniform due to the local residue stresses, and therefore the thermal resistors in a model should vary along the drawing direction. With the above in mind, an accurate effective-medium theory should include the above two factors. On the other hand, measuring the bridge numbers and the TC of each bridge are challenging in the modified Takayanagi’s model with varying thermal resistors,46 so the PE morphology information obtained from CG-MD simulations here would be very useful in developing such a model.

In our effective-medium theory model, the simulated polymer system is divided into Ns serial slabs along the drawing direction to address the issue of non-uniform chain orientations, with the thickness of each slab equal to the diameter of a CG bead (5.8 Å). Each slab is composed of both crystalline and amorphous domains, mimicking thermal resistors in parallel, due to the fact that they blend well with each other, as shown from the cross-sectional area of the PE fibers (Fig. 3a). This new model combines serial and parallel thermal resistors systematically, as illustrated in Fig. 9a. Finally, the effective TC κeff along the drawing direction is estimated by the following equation:

$$\frac{1}{{\kappa _{\mathrm{eff}}}} = \mathop {\sum}\limits_{i = 1}^{N_{\mathrm{S}}} {\left[ {{\mathrm{\Psi }}_i\left( {\kappa _{{\mathrm{c}} \bot }\frac{{2 - 2P_{2z{\mathrm{c}}}^i}}{3} + \kappa _{{\mathrm{c}}||}\frac{{2P_{2z{\mathrm{c}}}^i + 1}}{3}} \right) + \left( {1 - {\mathrm{\Psi }}_i} \right)\kappa _0{\mathrm{exp}}\left( {aP_{2z{\mathrm{a}}}^i} \right)} \right]^{ - 1}},$$
(1)

where Ψi is the degree of crystallinity of the ith slab; κc∥ and κc⊥ are the TCs of the crystalline domain in parallel and perpendicular to the fiber axis, respectively; \(P_{2z{\mathrm{c}}}^i\) and \(P_{2z{\mathrm{a}}}^i\) are the orientational order parameters over crystalline and amorphous beads in the ith slab, respectively; κ0 is the TC of undrawn isotropic amorphous PE; and a is the only fitting parameter obtained by minimizing the relative deviations between CG-MD simulated and model-predicted TCs under all the conditions here. Note that the exponential correlation between the TC contribution from the amorphous domain and \(P_{2z{\mathrm{a}}}^i\) on the right of Eq. 1 is based on the model developed by Liu and Yang,22 where pre-factor a was introduced as a fitting parameter.

Fig. 9: Theoretical model of TC and comparison between simulated and model-predicted TCs.
figure 9

a Illustration of the effective-medium theory developed here based on a serial-parallel thermal resistor model for the TC of semicrystalline polymers under tensile deformations. The Ns slabs, each composed of crystalline and amorphous domains in parallel, are serially connected to each other. b 3D plots of CG-MD simulated and model-predicted TCs under different Ψ and P2z. The color of the triangular symbols represents the relative deviations (in absolute values) between CG-MD simulated and model-predicted TCs.

Full size image

To determine κc∥ and κc⊥, we conducted equilibrium molecular dynamics (EMD) simulations (see Methods) using the Green–Kubo method.47 A perfect PE crystal structure consisting of 36 150-bead polymer chains is used here and shown in Supplementary Fig. 10a. To obtain reliable predictions of κc∥ and κc⊥, EMD simulations were repeated five times using different random seeds for generating the initial bead velocities. The size dependence of κc∥ is also tested by studying a simulation cell with double the length along the chain axis, and no discernible size dependence is found. In addition, due to the low energy carried by long wavelength phonons in polymers, we do not expect the long wavelength-cutout effect in our EMD simulations with a large scale of ~100 nm to greatly affect the TC results. The resulting ensemble-averaged κc∥ = 1.953 W m−1 K−1, κc⊥ = 0.131 W m−1 K−1, κ0 = 0.105 W m−1 K−1, and the fitting parameter a = 1.5. Note that κ0 of 0.105 W m−1 K−1 for unstrained amorphous PE is lower than those of semicrystalline PE measured in experiments due to both the absence of crystallinity and the relatively low molecular weight considered here. Our predicted κ0 value does match the experimental value of ~0.11 W m−1 K−1 for fully amorphous polymers well.48

The model-predicted TC profile under a drawing temperature of 500 K, strain rate of 5 × 107 s−1, and different strain levels (see Fig. 7a) agrees well with the CG-MD simulated TC profile. For all thermal-drawing conditions and strain levels, relative deviations between the simulated and the model-predicted TCs (see Fig. 9b) are mostly within 20%, except for a few cases under excessively high strain rates and strain levels. For example, under a drawing temperature of 600 K, strain rate of 5 × 109, and strain levels of 7 and 8, the relative deviations are ~50%. Excluding these outliers, the average relative deviation is ~11.5% over all the thermal-drawing conditions.

Enhanced phonon transport properties from crystallization

To elucidate the underlying mechanism behind thermal transport in thermally drawn polymer fibers, we investigate the phonon properties of PE samples at 300 K. Without loss generality, here we consider samples under different strain levels prepared at a representative strain rate of 5 × 107 s−1 and a drawing temperature of 500 K. According to the kinetic theory of phonon transport, the TC can be estimated as:47

$$\kappa = \frac{1}{V}\mathop {\sum }\limits_{p,k} C_{p,k}v_{p,k}l_{p,k} = \frac{1}{{3VN_0}}\mathop {\sum }\limits_{p,k} \hbar \omega _{p,k}D_{p,k}\frac{{\partial f_{p,k}}}{{\partial T}}v_{p,k}l_{p,k},$$
(2)

where Cp,k is the phonon specific heat, vp,k is the group velocity, lp,k is the phonon MFP, V is the volume of the sample, N0 is the total number of beads in the system, Dp,k is the phonon vibrational density of states (VDOS), and fp,k is the Bose–Einstein distribution function for phonon excitations. The subscripts “p” and “k” indicate phonon polarization and wave vector, respectively. Equation 2 unveils that the mode-specific phonon MFP, VDOS, and group velocity are the primary factors dictating phonon thermal transport. Note that lp,k could be decomposed into contributions from the intrinsic phonon–phonon scattering lph–ph and the phonon–boundary scattering lph–b due to the finite simulation box dimension.

Figure 10a shows the calculated intrinsic phonon MFP lph–ph and the intrinsic TC for a macroscopic sample κ∞, obtained from length-dependent TC results where the impact from the finite system size-induced phonon–boundary scattering is removed from NEMD simulations (see Methods and Supplementary Fig. 11). Through tensile deformation, lph–ph increases from around 10 Å at zero strain to 3367 Å at a strain of 18. This is due to the fact that strain-induced chain alignment and crystallization reduce structural disorders, weaken three-phonon Umklapp scattering, and thus, enhance lph–ph.

Fig. 10: Intrinsic phonon transport properties of bulk PE at 300 K with increasing strain under a strain rate of 5 × 107 s−1 and a drawing temperature of 500 K.
figure 10

a Phonon intrinsic mean free path lph–ph and intrinsic TC κ∞. The error bars represent standard errors from linearly fitting the 1/κ(L) vs. 1/lph–b curve (see Supplementary Fig. 11). b Phonon VDOS (top) and the corresponding spectral specific heat (bottom). c Longitudinal (vL) and transverse (vT) group velocities.

Full size image

The top of Fig. 10b shows the phonon VDOS of amorphous PE under different strain levels and that of a perfect PE crystal consisting of 2025 150-bead chains (structure shown in Supplementary Fig. 10b). When increasing the strain, the two major peaks, originally located at 5 and 27 THz, in the VDOS shift slightly to the right. In addition, when increasing the strain, the height of the first peak increases, while that of the second peak decreases. The VDOS peak at 27 THz is attributed to the vibrations of covalent bond, ωbond, based on the simple normal mode equation: \(\omega _{\mathrm{bond}} = 2\sqrt {k_e/m}\), where ke is the force constant of the bond stretching potential and m is the mass of the bead. The VDOS peak at 5 THz thus comes from angular bending and nonbonded interactions. The above trends in VDOS are expected, as the VDOS of amorphous PE under larger strain converges to that of a perfect PE crystal. According to Eq. 2, the frequency-dependent specific heat approximately correlates with the product of frequency and VDOS. The bottom of Fig. 10b shows that phonon specific heat in low-frequency domain increases, while that in high-frequency domain decreases when increasing the strain. These differences in frequency-dependent VDOS and specific heat are expected to have minor impacts on TCs of thermally drawn PE fibers.

Figure 10c shows that the longitudinal group velocity vL increases when increasing the strain, while the transverse group velocity vT is not greatly affected by the strain (see Methods). Through tensile deformation, vL at a strain of 18 is three times larger than that at zero strain due to the enhanced elastic modulus along the fiber-drawing direction. Therefore, the increase in TC is mainly attributed to the improvements in the intrinsic phonon MFP and the longitudinal group velocity, while the phonon VDOS has a minor contribution to the enhancement in TC.

In summary, a predictive CG-MD model for simulating the thermal-drawing process and the TC of PE fibers was developed here. This carefully trained model, facilitated by target property matching and the MS-IBI method, could successfully describe the structural, phase transition, and TC properties of PE fibers. CG-MD simulations are performed to investigate the impacts of annealed tensile deformation (mimicking fiber-drawing conditions, such as drawing temperatures, strain levels, and strain rates) on chain alignment, crystallinity, phonon transport, and the resulting TCs in bulk PE. Our simulation results show that P2z, a quantitative measure of chain alignment, increases and then reaches a plateau with increasing the strain, except for the cases under excessively high strain rates, where P2z increases and then drops off sharply due to extreme stress concentrations and the resulting structural damage. Likewise, the degree of crystallinity, Ψ, shows similar tendency to P2z with respect to strain, but differs from P2z under high strain levels (i.e., Ψ only drops slightly). Both P2z and Ψ exhibit nonlinear dependences on the drawing temperature and strain rate. Lower drawing temperatures will lead to slower stress relaxation, and possibly destroy the chain alignment and crystallinity after releasing the residual stress, while excessively high drawing temperatures cannot provide enough tensile forces (stress dissipates too fast due to chain diffusions) to align polymer chains well or promote crystallization. Meanwhile, higher strain rates will lead to extreme stress localization and structural damage, similar to the cases at lower temperatures, while lower strain rates suffer from extra stress relaxation and chain diffusions, similar to the cases at higher temperatures.

An optimal combination of moderate drawing temperature and strain rate is essential to achieve highest degrees of chain alignment P2z, crystallinity Ψ, and the resulting TCs. From the viscoelastic relaxation aspect, we found that such an optimal combination requires a delicate balance between stress localization (to align polymer chains) and chain diffusion (for viscoelastic relaxation). After examining the Deborah number, we found that it would take around 103 and 104 times of the Rouse relaxation time of melt and glassy polymers (with the molecular weight considered here), respectively, to achieve the optimal fiber performance. The Deborah number could serve as an important and general dimensionless metric to design the best thermal-drawing protocols for a wide range of polymer types under various combinations of drawing temperature and strain rate.

For phonon transport, the contribution from bonded interactions to the total TC is more significant than contributions from nonbonded interactions, angular interactions and convections. The influences of strain level, drawing temperature, and strain rate on TC are similar to those on P2z and Ψ. However, neither P2z nor Ψ can correlate very well with the computed TC due to the interplay between P2z and Ψ. An effective-medium-theory model, considering impacts from both P2z and Ψ, as well as the serial-parallel nature of semicrystalline polymers, has been successfully applied here to predict the TCs of thermally drawn PE fibers. Furthermore, the enhancement in TC by tensile deformation is mainly attributed to the increases in the intrinsic phonon MFP and the longitudinal group velocity. We believe that this large-scale molecular modeling work could provide fundamental insights into the thermal-drawing process of semicrystalline polymer fibers and provide a complete process–structure–property correlation. This work could also guide future experimental efforts to optimize thermal-drawing protocols to promote chain alignment, degree of crystallinity, and the resulting phonon transport performance.

Methods

Modeling the annealed fiber-drawing process via thermal strains

All the MD simulations are conducted using the LAMMPS software package49 with periodic boundary conditions in three dimensions. The CG force field for PE is developed based on AA-MD simulation results using the AIREBO potential.50 For the CG force field training approach and validation results, please refer to Supplementary Method and Supplementary Discussion. To investigate the polymer fiber-drawing process modeled via thermal straining, a system consisting of 2025 150-bead PE chains was relaxed at 500 K using the double-bridging51 algorithm for 200 ns. The randomized amorphous PE system was first quenched to 300 K at a rate of 2 K n s−1, and then heated to 600 K at the same rate, during which samples at different temperatures were collected. After preparing amorphous PE samples under 300, 400, 500, and 600 K, CG-MD tensile deformation simulations were performed at the above temperatures under various strain rates (see Table 1). Note that the strain rates here are higher than experimental values, for example, 1400 s−1 in Shrestha et al.’s experimental work.13 The reason is that the PE’s Rouse relaxation time here is much shorter than those in experiments due to the lower molecular weight used here. Therefore, high strain rates are required to match the short Rouse relaxation time here to yield a suitable Deborah number (relaxation time × strain rates) that can achieve the optimal morphology and TC of PE fibers. Uniaxial deformation was applied along the z-axis of the simulation box, until reaching an engineering strain of 20. The pressures in the lateral dimensions were kept constant at 1 bar during the deformation. When the desired strain value was achieved, each strained sample was quenched to 300 K instantly and stress relaxed under the NPT ensemble (P = 1 bar and T = 300 K) for 20 ns to release all the residue stresses, mimicking the status of post-processed PE fibers. Figure 1b shows the simulation snapshots of the PE fiber morphology under different strain levels after the quenching and the stress relaxation processes at 300 K. Figure 1c depicts the flowchart of the simulated fiber thermal-drawing process.

Table 1 Thermal-drawing processing parameters for bulk PE fibers considered here.
Full size table

Stress relaxation modulus calculations

The stress relaxation modulus G(t) for undrawn bulk amorphous PE were measured using EMD simulations through the stress autocorrelation function:

$$G\left( t \right) = \frac{V}{{k_{\mathrm{B}}T}}\left\langle {\sigma _{ab}\left( 0 \right)\sigma _{ab}\left( t \right)} \right\rangle,$$
(3)

where kB is the Boltzmann constant, σab are the off-diagonal xz and yz shear components of the stress tensor, and z is the fiber-drawing direction. In EMD simulations, the system was first relaxed under the NPT ensemble at 1 bar (anisotropic pressure coupling, where each dimension is able to change) and corresponding temperatures (300, 400, 500, and 600 K) to remove the residual stresses. Then, the shear stresses were collected every 20 fs over a long period of 1 μs under the NVT ensemble to calculate G(t) using the multiple-tau correlator method.52

TC calculations

The TCs of single PE chain and bulk PE system were calculated using the NEMD simulations through the Müller-Plathe algorithm.53 A schematic of the simulation system used to compute the TCs of bulk PE systems is shown in the inset of Supplementary Fig. 16b. In NEMD, we induce a heat flux Jtotal by exchanging the velocity of the hottest atom in the cold slab with that of the coldest atom in the hot slab. For the systems to reach steady states in which the temperature profile is stable, it typically takes 1 ns for AA-MD simulations of single PE chains, and 20 ns for CG-MD simulations of both single PE chains and bulk PE systems. After reaching steady states, the temperature of each slab was averaged over another 1 ns for AA-MD simulations and 60 ns for CG-MD simulations. Supplementary Fig. 16b shows a typical temperature profile of the simulation system along the z-axis. The temperature gradient \(\partial T/\partial z\) was then obtained by fitting the linear regions of the temperature profile using the least-squares method. Note that the artificial selection of linear regions could result in deviations in the computed TC. For a single PE chain, the deviation is quite small (see Supplementary Fig. 12), while the deviation for bulk PE systems could be negligible due to the highly linear temperature profile obtained from the NEMD simulations (see Supplementary Fig. 16b). Sequentially, the total TC can be computed from the Fourier’s law:

$$\kappa _{\mathrm{total}} = - \frac{{J_{\mathrm{total}}}}{{\left( {\partial T/\partial z} \right)A}},$$
(4)

where A is the cross-sectional area for heat flux. For single PE chain systems, A = 18 Å2, based on the unit cell dimensions of a single chain in an idealized bulk PE crystal structure.9 For bulk PE systems, A equals to the cross-sectional area of the simulation box. All TC calculations were performed at 300 K using NEMD simulations.

TC contributions from various mechanisms are revealed by decomposing the total heat flux in NEMD simulations, which is given by the following equation:54,55

$$J_{\mathrm{total}} = \frac{1}{V}\left( {\mathop {\sum }\limits_i e_iv_i - \mathop {\sum }\limits_i S_iv_i} \right),$$
(5)

where ei and vi are the internal (potential and kinetic) energy and velocity of bead i, respectively. Si for bead i is calculated using the following equation:

$$S_{ab} = - \left[ {\begin{array}{*{20}{c}} {\frac{1}{2}\mathop {\sum }\limits_{n = 1}^{N_{\mathrm{p}}} \left( {r_{1a}F_{1b} + r_{2a}F_{2b}} \right) + \frac{1}{2}\mathop {\sum }\limits_{n = 1}^{N_{\mathrm{b}}} \left( {r_{1a}F_{1b} + r_{2a}F_{2b}} \right)} \\ { + \frac{1}{3}\mathop {\sum }\limits_{n = 1}^{N_{\mathrm{a}}} \left( {r_{1a}F_{1b} + r_{2a}F_{2b} + r_{3a}F_{3b}} \right)} \end{array}} \right],$$
(6)

where Sab is the ab components of Si; r1, r2, and r3 are positions of interacting beads 1, 2, and 3, respectively; F1, F2, and F3 are the forces acting on beads 1, 2, and 3, respectively; and Np, Nb, and Na are the numbers of nonbonded (vdW here), covalent bonded and angular interactions acting on bead i, respectively. The first and second terms on the right of Eq. 5 are the convection and virial contributions due to different interactions, respectively. Therefore, the total heat flux, Jtotal, is decomposed into contributions from convection and virial contributions due to different interactions:

$$J_{\mathrm{total}} = J_{\mathrm{convection}} + J_{\mathrm{nonbonded}} + J_{\mathrm{bond}} + J_{\mathrm{angle}}$$
(7)

The total TC, κtotal, is then decomposed into contributions from convection, nonbonded, bonded, and angular interactions by the corresponding heat flux term using the Fourier’s law.56

To calculate TCs of perfect PE crystals composed of chains with finite length, we need to consider the interfacial thermal resistance when two chains are joint head-to-tail via nonbonded interactions. As a result, we expect discontinuities in the temperature profile for a PE crystal, if we calculated TC using the Müller-Plathe algorithm. Therefore, the TCs of perfect PE crystals at 300 K were determined using the EMD simulation via the Green–Kubo method.47 In EMD, after energy minimization, the system was relaxed under the NPT ensemble with constant pressure of 1 bar for 20 ns, followed by an extensive sampling for 1.6 μs to integrate the heat flux autocorrelation function for computing κi along dimension i (x, y or z). Specifically,

$$\kappa _i = \frac{V}{{k_{\mathrm{B}}T^2}}\mathop {\int }\nolimits_0^\infty \left\langle {J_i\left( 0 \right)J_i\left( t \right)} \right\rangle {\mathrm{d}}t,$$
(8)

where Ji is the heat flux along the ith dimension.

Orientational order parameter and degree of crystallinity

Orientational order parameter22,57,58 is a very useful measure for the alignment of polymer chains along the applied strain direction. The local chain orientation at each bead i can be computed by the chord vectors connected to the bead: ei = (ri + 1 − ri − 1)/|ri + 1 − ri − 1|. The total orientational order parameter P2z along the z direction was calculated as:

$$P_{2z} = \frac{{\mathop {\sum }\nolimits_{i = 1}^N 3\left\langle {\left( {{\mathbf{e}}_i \cdot {\mathbf{e}}_z} \right)^2} \right\rangle - 1}}{{2N}},$$
(9)

where ez is the unit vector along the direction of the applied strain and N is the total number of beads.

To characterize each bead i as either crystalline or amorphous, the local orientational order parameter P2,i,59 also a measure of local degree of crystallinity, is calculated as:

$$P_{2,i} = \frac{{\mathop {\sum }\nolimits_{j = 1}^n 3\left\langle {\left( {{\mathbf{e}}_i \cdot {\mathbf{e}}_j} \right)^2} \right\rangle - 1}}{{2n}},$$
(10)

where ej is the unit chord vector at bead j and n is the number of beads j that are within a cutoff distance of 10 Å from bead i. If P2,i is larger than a threshold value of 0.7, bead i is regarded as a crystalline bead, otherwise it is regarded as an amorphous bead. Then, the degree of crystallinity Ψ of the system is calculated as Ψ = Nc/N, where Nc is the number of crystalline beads. We also estimated Ψ by the density method used in regular experiments. Specifically, for a PE system with an averaged density ρ, the degree of crystallinity can be defined as \(\Psi = \frac{{\rho \, - \, \rho _0}}{{\rho _{\mathrm{cry}} \, - \, \rho _0}}\), where ρ0 and ρcry are the densities of amorphous and perfect crystalline PE systems, respectively. The Ψ values determined by the density method are lower than those by the order parameter method in this work, but they share the same trend with increasing the strain (see Supplementary Fig. 13). Therefore, we expect that the differences derived from these two methods would not change the main conclusions in this work.

Intrinsic phonon property calculations

The velocity autocorrelation function (VACF) averaged over all CG beads was computed from 40 ps of CG-MD simulation under the NPT (P = 1 bar and T = 300 K) ensemble with velocities sampled every 20 fs. The normalized phonon VDOS is obtained by taking the discrete Fourier transform of the VACF. Specifically, the VDOS, D(ω), can be expressed by

$$D\left( \omega \right) = \frac{1}{{2\pi }}\left[ {\mathop {\int }\nolimits_0^\infty \frac{{\left\langle {v\left( 0 \right)v\left( t \right)} \right\rangle }}{{\left\langle {v\left( 0 \right)v\left( 0 \right)} \right\rangle }}\exp \left( { - {\mathrm{i}}\omega t} \right){\mathrm{d}}t} \right]^2,$$
(11)

where ω is the phonon angular frequency, and v is the velocity of each CG bead.

According to the kinetic theory of phonon transport;60 in general, the TC κ is proportional to the effective phonon MFP leff through:

$$\kappa \! \sim \mathop {\sum }\limits_i C_iv_il_i = l_{\mathrm{eff}}\mathop {\sum }\limits_i C_iv_i,$$
(12)

where Ci, vi, and li denote the phonon specific heat, group velocity, and MFP of mode-i phonon, respectively. For a fully condensed bulk polymer system without voids/vacancies, the intrinsic Umklapp phonon–phonon scattering and the phonon–boundary scattering are considered as the primary phonon scattering mechanisms in NEMD simulations. Assuming that phonons scatter through the Matthiessen’s rule based on the phonon relaxation time approximation,61,62,63 the κ value can be correlated to different phonon MFPs as:

$$\frac{1}{\kappa }\sim \frac{1}{{l_{\mathrm{eff}}}} = \frac{1}{{l_{\mathrm{ph - ph}}}} + \frac{1}{{l_{\mathrm{ph - b}}}},$$
(13)

where lph–ph is the intrinsic phonon MFP due to three-phonon Umklapp scattering and lph–b is the phonon MFP due to phonon–boundary scattering. lph–b can be estimated by the half lengths of the simulation cell along the z-axis here due to the extra heat source and sink introduced in the Müller-Plathe algorithm.53 To study the intrinsic phonon–phonon scattering and estimate TCs for macroscopic polymer fibers, samples at strains of 0, 2, 4, 6, 12, and 18 were replicated two and four times along the strain direction (z-axis) to obtain different cell lengths. For these simulation cells at different lengths, the NEMD simulations were performed for 80–800 ns depending on the size of the system in order to obtain converged temperature profile and TC. The intrinsic MFP lph–ph was then evaluated by linearly fitting the 1/κ(L) vs. 1/lph–b curve (Supplementary Fig. 11).

The primary phonon group velocities can be estimated by the elastic moduli of the system using the following equation:60

$$v_{\mathrm{L}} \approx \sqrt {E/\rho } \;{\mathrm{and}}\;v_{\mathrm{T}} \approx \sqrt {G/\rho },$$
(14)

where vL and vT are the longitudinal and transverse group velocities, respectively; E and G are the Young’s modulus and the shear modulus, respectively; and ρ is the mass density. CG-MD tensile and shear deformation simulations were both conducted to obtain the stress–strain curves, from which the elastic moduli E and G (see Supplementary Fig. 14) could be estimated by linearly fitting the linear elastic region. Note that the stress–strain curves were obtained from post-drawn samples (after stress relaxation and quenched to 300 K) collected under various strain levels; therefore, we can obtain strain-dependent elastic moduli and the corresponding primary group velocities.