Polymer-BN Composites as Thermal Interface Materials for Lithium-Ion Battery Modules: Experimental and Simulation Insights

1. Introduction

As lithium-ion batteries (LIBs) become increasingly indispensable in electric vehicles (EVs), portable electronics, and stationary energy storage systems, their thermal management has emerged as a central engineering challenge [1,2,3,4,5,6,7,8]. Battery operation generates substantial heat, particularly under fast-charging or high-discharge conditions, leading to temperature gradients, accelerated degradation, and, in severe cases, thermal runaway [9,10]. Maintaining thermal uniformity and dissipating heat efficiently are critical for battery performance, safety, and lifespan.

In typical battery modules, heat is dissipated through passive and active cooling strategies, including airflow, liquid cooling manifolds, and heat pipes [2,3,4,5,11,12,13]. However, a significant thermal bottleneck persists at the interface between battery cells and their surrounding cooling infrastructure. Imperfect contact, surface roughness, and the low thermal conductivity of adhesives or structural materials contribute to significant thermal contact resistance [9,11,14,15,16,17,18]. This makes developing effective thermal interface materials (TIMs) crucial to advancing battery thermal management systems (BTMS).

TIMs are designed to bridge the thermal resistance gap between solid interfaces, such as battery casings and heat sinks. Polymer-based TIMs are often favored due to their mechanical flexibility, electrical insulation, and ease of processing [9,11,15,16,18,19]. Yet, the low intrinsic cross-plane thermal conductivity of most polymers (typically < 0.5 W m⁻¹ K⁻¹) severely limits their ability to conduct heat. To overcome this, high-conductivity fillers—such as graphite, metal oxides, and carbon nanostructures like graphene—are incorporated into the polymer matrix to form thermally conductive composites [20,21,22].

Graphene, with its exceptional in-plane thermal conductivity (up to 5000 W m⁻¹ K⁻¹), high aspect ratio, and tunable surface chemistry, has emerged as one of the most promising TIM fillers [18,23,24,25,26,27]. Its incorporation into polymer matrices can dramatically enhance thermal performance, particularly when oriented or networked to form continuous heat conduction paths [15,28,29,30,31,32,33]. However, graphene composites are also electrically conducting, which is a disadvantage in preventing undesired short circuits within the battery module.

Hexagonal boron nitride (h-BN), often dubbed “white graphene,” offers similarly high thermal conductivity along with electrical insulation, making it attractive for TIM applications where dielectric performance is important [17]. BN-polymer composites have demonstrated thermal conductivities up to 10 W m⁻¹ K⁻¹ depending on filler loading, morphology, and processing method [16,18,34,35,36,37]. A table summarizing prior results of thermal conductivity for various TIMs is presented in Table S1 [38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55]. Recent advances in filler surface functionalization, hybrid filler strategies, and anisotropic alignment have further improved the thermal transport efficiency of such composites .

While the thermal performance of individual materials is important, their integration into real-world systems poses additional design constraints. For example, TIMs must conform to irregular geometries, maintain thermal performance under cyclic loading, and avoid adding excessive mass or volume to the battery pack [11]. A comprehensive evaluation of TIMs requires multiscale assessment—from intrinsic thermal properties such as thermal diffusivity to system-level metrics like inter-cell temperature gradients and transient cooling behavior [9]. Finite-element simulations, particularly those implemented in platforms such as COMSOL Multiphysics, have become essential in linking material properties with thermal field evolution in complex battery assemblies [35,56].

In this study, we explored the fabrication and thermal performance of two classes of polymer composite TIMs: graphene-PLA (GPLA) produced via 3D printing, and BN-loaded thermoplastic polyurethane (TPU), cast with 20 and 40 wt.% BN content. Building on prior demonstrations of enhanced thermal conductivity in similar systems [11,15,15,17,18,28,35,56], we performed a comparative analysis using thermogravimetric analysis (TGA) and thermal diffusivity measurements to characterize material stability and heat transport properties. To evaluate real-world applicability, we implemented the TIMs in a compact battery module comprising six cylindrical 18650 LIB cells with an active liquid cooling manifold. We conducted thermal cycling experiments by heating the module to ${45}^{\circ }$C under various TIM configurations and measuring the decay rates and inter-cell temperature differentials ($\Delta T$). We employed finite-element simulations in COMSOL using experimentally derived parameters to model heat transfer within the battery module. We also studied the changes in the cell temperature with different TIMS at different C-rates as a function of charge-discharge cycles. Our simulations strongly agree with experimental data, particularly in reproducing the trends observed T and $\Delta T$ values for different TIMs. In addition to finite-element simulations performed in COMSOL using experimentally derived parameters, we developed an analytical framework linking the cooling time constant directly to the thermal diffusivity of the TIM. Furthermore, our study reveals that effective thermal management in battery modules depends on balancing two key processes: vertical heat sinking into the cooling infrastructure (cross-plane) and lateral thermal diffusion (in-plane) across adjacent cells. An ideal TIM must therefore achieve both a high heat decay rate and strong lateral thermal spreading to minimize inter-cell temperature gradients. Our findings contribute to the growing body of work demonstrating the feasibility of additive-manufactured, thermally-engineered polymer composites for thermal management in LIB systems.

2. Materials and Methods

2.1. Materials

Polylactic acid or PLA filament was purchased from Hacthbox while graphene-loaded PLA or GPLA filament was purchased from Blackmagic3D.com. Thermoplastic polyurethane or TPU filament and TPU pellets were purchased from TCPoly. Hexagonal boron nitride (h-BN) micropowder and N, N-dimethylformamide (DMF) were purchased from Sigma Aldrich.

2.2. Fabrication of Thermal Interface Materials

For 3D printing, a 70 x 30 x 1.4 mm slab was designed using the geometry module of COMSOL Multiphysics. Upon finalizing the geometry, the design was exported and sliced in PrusaSlicer software to generate the code for 3D printing. The exported code was then printed using a Prusa i3 MKS3+ 3D printer. We used commercially available polylactic acid or PLA, graphene-loaded PLA filament or GPLA, and thermoplastic polyurethane or TPU filaments for comparison with our polymer-nano composites.

For preparing TPU-BN composites, TPU pellets were dried in a vacuum oven at 100 °C for 4 hours. Subsequently, exfoliated micro h-BN micropowder was dispersed in N, N-dimethylformamide (DMF) using a tip sonicator at 10 watts for 30 minutes. After sonication, the mixture was placed on a hot plate at 80 °C and magnetically stirred at 150 rpm. Dried TPU pellets were added to the mixture slowly. Subsequently, the mixture was stirred at 80 °C and 150 rpm for 4 hours. Finally, the viscous mixture was cast onto a thin kitchen aluminum foil using an MTI doctor blade casting system to maintain a uniform thickness of 1.4 mm. The films were dried overnight at room temperature inside a chemical hood (120 cfm air flow) to remove all of the solvent slowly. Once dried, the TIMs were peeled from the foil and cut into 70 x 30 x 1.4 mm slabs.

2.3. Characterization Techniques

All TIMs were characterized using thermogravimetric (TGA) analysis (TA instrument SDT Q600) under nitrogen gas flow at 100 ml/min from room temperature to 600 °C with a ramp of 20 °C/min. Cross-plane thermal diffusivity was measured using the laser flash technique (Linseis LZT meter). Before diffusivity measurements, all samples were coated with graphite 33 spray to ensure uniform contact for heat flow. To analyze the performance of the TIMs, we made a prototype battery module with six 18650 cells (Panasonic NCR18650GA 3300 mAh); an aluminum manifold with internal tubes (similar to the manifold used in Tesla model 3) was used to flow liquid water coolant around the battery module at 40 ml/min using an electrical pump. TIMs were placed between the cooling tube and battery interfaces on both sides of the battery module. Six k-type thermocouples were inserted into the battery module to monitor the temperature of each battery using Pico software. FLIR E60 camera was used to capture IR images at different temperatures during heating and cooling. The battery module was heated to around 45 °C and then cooled down inside a hood with controlled air flow at 120 cfm. All the measurements were performed in triplicates (n=3) to ensure statistical robustness. Scanning electron microscope images, presented in Figures S1 and S2 of the supporting information, were obtained using a Hitachi 6600 SEM. High resolution X-ray diffraction (HR-XRD) is performed using a RIGAKU Ultima IV diffractometer, employing Cu K$\alpha$ radiation, on powder samples that were held by a standard Al sample holder. Quantitative analysis using Rietveld refinement is performed on the XRD peaks using PDXL software. Atomic force microscopy (AFM) measurements were performed in a non-contact mode using AIST-NT SPM Smart system and cantilevers (HQ: NSC14/Al BS-50) from Micromasch. AIST-NT image analysis and processing (Version 3.2.14) software was used for AFM image analysis. Malvern Zetasizer 90 was used for dynamic light scattering measurements. Based on the DLS measurements, the average lateral size of BN particles is 745 ± 92 nm, which is consistent with the average size deduced from AFM measurements that showed an average lateral size of 600-800 nm with a thickness of 70-100 nm (see Figure S3 in the supporting information). Differential scanning calorimetry was performed to extract the specific heat values (Figure S4) using TA instruments Discovery DSC from 25- 600 °C.

2.4. COMSOL Modeling

The computational modeling part of the heat transfer mechanism was achieved using COMSOL Multiphysics 6.2 employing finite element method (FEM). The 18650 cells were modeled using a lithium-ion battery module (liion) with LMO-${\mathrm{LiMnO}}_4$ Spinel, graphite (${\mathrm{LixC}}_6$, MCMB) and LiPF₆ in 3:7 EC: EMC as respective cathode, anode, and electrolyte. Electrochemical reactions at the electrode/electrolyte interfaces follow Butler-Volmer kinetics given by

$$i=i_0\left(exp\left({\displaystyle \frac{{\alpha }_{a}F\eta }{RT}}\right)-exp\left(-{\displaystyle \frac{{\alpha }_{c}F\eta }{RT}}\right)\right)$$

(1)

where ${\alpha }_a$ (/${\alpha }_c$) is the anodic (/cathodic) charge-transfer coefficient, F is the Faraday constant, $\eta$ is the overpotential, R is the universal gas constant, and T is the temperature in Kelvin. The electrolyte transport through porous electrodes was modeled via the Bruggeman relation.

A 3D heat transfer in solids and fluids (ht) module was further coupled with the Li-ion module using the average values of heat generated ($Q_{h\left(avg\right)}$) and temeprature ($T_{avg}$). The thermal behavior of the module was then studied with a maintained coolant flow (water) through aluminum tubes similar to commercial battery packs. TIMs of size 70 x 30 x 1.4 mm (same as the ones used in experiments) were placed in between the batteries and the liquid cooling manifold. For the thermal cycling part, the module heating was achieved by a boundary heat source at the bottom of the cell with a 10 W heat source with a step function that triggers heat decay at 2400 s (cooling). A non-slip boundary condition was used for the fluid flow with flow rate limited to 1.08 x ${10}^{-7}$ m²/s using the laminar flow (spf) module. COMSOL further uses the nonisothermal flow (nitf) module to couple the heat transfer spf and ht interfaces. Normal meshing (max. element size 11.4 nm) was used for all domains, with coolant tube boundaries defined using fine (max. element size 5.3 nm). For the charge-discharge studies, the six cells were assumed to be charging independency at the set c-rate simultaneously. An alternating charge-discharge cycle of current flow was employed using a square wave function with the current and cycle time as the variable to control the C-rate.

3. Results

Figure 1a presents the TGA curves for five materials: neat PLA, GPLA, neat TPU, and TPU loaded with 20 wt.% and 40 wt.% BN). All samples were heated under nitrogen atmosphere up to 600°C, and weight loss was monitored to evaluate thermal degradation behavior and filler content. Neat PLA and GPLA both exhibit a sharp one-step degradation between 300°C to 400°C, consistent with known depolymerization and chain scission mechanisms of PLA [11]. However, the GPLA sample demonstrates a slightly earlier onset of decomposition and retains significantly more mass at high temperature due to the presence of graphene fillers. Specifically, at 450 °C, PLA retains less than 2% of its initial mass, whereas GPLA retains approximately 15%. This difference indicates that there is ∼ 15 wt. % of graphene incorporated into GPLA. Neat TPU shows a broader and more gradual weight loss profile with a single step beginning near 300 °C, and a residual mass of around 45% at 450 °C. Such a single-step decomposition process is characteristic of the breakdown of both hard and soft segments in the polymer matrix. The addition of BN significantly alters this behavior. The 20%BN–80%TPU composite exhibits a two-step decomposition with a second decomposition onset at 400 °C and retains approximately 52% of its initial mass at 450 °C. The 40%BN–60%TPU sample demonstrates the highest thermal stability among all samples, retaining nearly 70% of its mass at the same temperature. The two-step decomposition is more evident in 40%BN–60%TPU sample. The progressive increase in residue with h-BN loading is indicative of the high thermal stability of h-BN, which acts as an inert filler resisting oxidative or thermal degradation [16]. The emergence of a second decomposition step in BN–TPU composites at 400 °C can be attributed to interfacial interactions between the BN fillers and the TPU chains. The high thermal stability and inertness of BN not only provide a thermal shielding effect but also restrict the mobility of TPU chains through non-covalent interactions such as van der Waals forces or hydrogen bonding (e.g., N–H⋯BN). These interfacial effects delay the decomposition of segments closely associated with BN surfaces, resulting in the observed higher-temperature degradation step. Therefore, the two-step decomposition is indicative of filler-mediated stabilization and heterogeneous breakdown of the polymer matrix in the presence of thermally robust BN domains.

Figure 1b shows the thermal diffusivity of the same five samples as a function of temperature up to 70°C. The diffusivity data reveals a clear hierarchy in thermal transport performance, strongly dependent on the type and concentration of thermally conductive fillers. The highest thermal diffusivity is observed in the 40%BN–60%TPU composite, with values consistently above $2.5\times {10}^{-2}$ cm² s^–1 across the entire temperature range. This value is more than an order of magnitude higher than that of the base TPU, which maintains a diffusivity near $3.5\times {10}^{-3}$ cm² s^–1. The substantial enhancement can be attributed to the formation of percolated thermal pathways within the TPU matrix due to the high BN content. BN’s high intrinsic thermal conductivity, good dispersion and high aspect ratio, likely contribute to efficient phonon transport even at modest temperature gradients [17]. The 20%BN–80%TPU composite shows intermediate performance, with a diffusivity of approximately $8.5\times {10}^{-3}$ cm² s^–1. While this is a notable improvement over neat TPU, the sub-percolation regime of filler content may limit the formation of continuous conduction networks. This aligns with existing studies suggesting that filler loadings above the percolation threshold are essential for achieving significant enhancements in thermal diffusivity [18]. The GPLA sample exhibits only a marginal improvement over neat PLA, with diffusivity values of about $3.0\times {10}^{-3}$ cm² s^–1 compared to $2.1\times {10}^{-3}$ cm² s^–1 for PLA. This modest increase suggests graphene in PLA matrix is suboptimal for effective thermal transport.

To assess the effectiveness of the TIMs under realistic conditions, a module comprising six cylindrical Li-ion cells was assembled with the TIM inserted between the cells and the liquid cooling manifold (Figure 2a,b). The module was subjected to controlled heating inside a forced-air convection hood with an airflow rate of 120 CFM. Heating was continued until the average temperature of the battery module reached approximately 45 °C, at which point the heat source was turned off, initiating the cooling phase. Each cell was attached with an independent thermocouple to monitor transient temperature profiles throughout the heating and cooling cycles. The aqueous liquid coolant was pumped through the Al manifold throughout the heating and cooling cycles. Typical IR images of the cells at different stages of heating are shown in Figure 2c,d.

Figure 3a shows the representative thermal response for all six cells when PLA was used as the TIM. Although all cells were exposed to the same environmental conditions, the temperature evolution was non-uniform, with distinct differences in heating rates and peak temperatures. Notably, Battery-1 reached the highest temperature while Battery-6 remained significantly cooler, exhibiting a temperature difference of nearly 4°C at the peak. These differences reflect the anisotropic thermal conductance within the module, asymmetries in cell-to-manifold contact, and intrinsic differences in cell self-heating behavior. Figure 3b shows the average temperature ($\overline{T}$) of the entire module, obtained from the temperature readings from all six thermocouples. The average curve shows a smooth and symmetric profile with a clear demarcation between the heating and cooling phases. To evaluate thermal uniformity within the module, the temperature difference $\Delta T_i\left(t\right)=|T_i\left(t\right)-\overline{T}\left(t\right)|$ was computed for each cell i with respect to the instantaneous average temperature $\overline{T}\left(t\right)$. The results are plotted in Figure 3c. The $\Delta T$ curves show that thermal nonuniformity peaks near the end of the heating phase, where some cells lag in temperature rise compared to others. During the cooling phase, $\Delta T$ gradually decreases as heat is more uniformly extracted via the manifold, but notable differences persist over the full transient period.

The transient cooling behavior of the battery module, following the removal of the external heat source, was analyzed using the following thermal conduction model. TIM serves as the dominant path for heat dissipation from the battery cells to the cooling manifold. Assuming negligible heat generation during the cooling phase and lumped thermal capacity of the cells, the rate of temperature decay can be modeled by equating the rate of internal energy loss to the conductive heat flux through the TIM. The rate of heat loss from the battery can be expressed as:

$$-\dot{Q}=C_{p,\mathrm{dev}}{M}_{\mathrm{dev}}{\displaystyle \frac{dT}{dt}}$$

(2)

where $C_{p,\mathrm{dev}}$ is the specific heat capacity and $M_{\mathrm{dev}}$ is the effective thermal mass of the battery module. The heat flux through the TIM is modeled as:

$$\dot{Q}={\displaystyle \frac{\kappa A_{\mathrm{TIM}}}{d_{\mathrm{TIM}}}}(T-T_{\infty })$$

(3)

where $\kappa$ is the thermal conductivity of the TIM, $A_{\mathrm{TIM}}$ is the contact area, $d_{\mathrm{TIM}}$ is the TIM thickness, T is the instantaneous temperature of the battery, and $T_{\infty }$ is the steady-state temperature.

Equating these expressions yields the governing differential equation for temperature decay:

$$C_{p,\mathrm{dev}}{M}_{\mathrm{dev}}{\displaystyle \frac{dT}{dt}}=-{\displaystyle \frac{\kappa A_{\mathrm{TIM}}}{d_{\mathrm{TIM}}}}(T-T_{\infty })$$

(4)

which simplifies to:

$${\displaystyle \frac{dT}{dt}}=-{\displaystyle \frac{\kappa A_{\mathrm{TIM}}}{C_{p,\mathrm{dev}}{M}_{\mathrm{dev}}{d}_{\mathrm{TIM}}}}(T-T_{\infty })$$

(5)

Introducing the thermal diffusivity $\alpha =\kappa /\left(\rho C_{p,TIM}\right)$, the time constant $\tau$ governing exponential decay is given by:

$$\tau ={\displaystyle \frac{C_{p,\mathrm{dev}}{M}_{\mathrm{dev}}{d}_{\mathrm{TIM}}}{{\rho }_{\mathrm{TIM}}{C}_{p,\mathrm{TIM}}{\alpha }_{\mathrm{TIM}}{A}_{\mathrm{TIM}}}}$$

(6)

showing that $\tau$ is inversely proportional to the thermal diffusivity ${\alpha }_{\mathrm{TIM}}$ of the interface material. Thus, the solution to the temperature decay equation is an exponential approach to $T_{\infty }$:

$$T\left(t\right)=T_{\infty }+(T_0-T_{\infty })e^{-t/\tau }$$

(7)

where $T_0$ is the initial temperature at the start of the cooling phase. Equation (7) was used to fit the experimental cooling data for each TIM configuration. The extracted decay constant $\tau$ provides a quantitative measure of the rate at which the battery module sheds heat. All the fit values are provided in Table S2 in the supporting information. Lower values of $\tau$ correspond to faster cooling and more thermally effective TIMs. The fitted results were subsequently compared across different TIM materials and correlated with their independently measured thermal diffusivities to verify the predicted inverse relationship $\tau \propto 1/{\alpha }_{\mathrm{TIM}}$. A representative plot (see Figure 3d) confirms this trend.

Figure 4 presents the time-resolved temperature deviations $\Delta T_i\left(t\right)=|T_i\left(t\right)-\overline{T}\left(t\right)|$ for each of the six cells in the attery module using five different TIMs: (a) PLA, (b) GPLA, (c) TPU, (d) 20%BN–80%TPU, and (e) 40%BN–60%TPU. The metric $\Delta T$ quantifies the deviation of each cell’s temperature from the instantaneous module average $\overline{T}\left(t\right)$, providing a dynamic measure of inter-cell thermal uniformity. Interestingly, the PLA TIM—despite having the lowest thermal diffusivity among the tested materials—exhibits the narrowest spread in $\Delta T$, with all six cells maintaining relatively close thermal trajectories. The peak $\Delta T$ remains below 3 °C throughout the thermal cycle. This result appears counterintuitive given PLA’s poor thermal conductivity. However, it can be rationalized by considering the relatively uniform thermal impedance across all cell–manifold interfaces: in the absence of a highly conductive path, each cell dissipates heat at a similar, albeit slow, rate (as indicated by a high $\tau$ for PLA in Figure 2d). The system behaves almost adiabatically on short time scales, and thermal gradients develop primarily due to intrinsic differences in cell positioning or contact pressure rather than material-mediated heat spreading.

In contrast, the introduction of higher-diffusivity TIMs such as GPLA and TPU increases the magnitude and spread of $\Delta T$. Notably, the GPLA TIM (Figure 4b) exhibits the largest inter-cell thermal heterogeneity, with Battery-6 showing a pronounced deviation exceeding 5 °C. This suggests that while the thermal conductivity of GPLA is higher than PLA, it remains insufficient to rapidly redistribute heat, and may accentuate local differences arising from subtle variations in thermal contact resistance. Similarly, the TPU-only TIM (Figure 4c) produces a broader $\Delta T$ spread than PLA, albeit less pronounced than GPLA. These observations underscore the sensitivity of thermal equilibration to both filler dispersion and the quality of interface coupling between cells and the cooling manifold. For the 20%BN–80%TPU system (Figure 4d), peak $\Delta T$ values slightly increased, with Battery-1 deviating significantly from the average. In the 40%BN–60%TPU composite (Figure 4e), the inter-cell temperature differences become markedly lower, and the spread among all six cells narrows considerably. This suggests that only beyond a critical threshold of thermal diffusivity does the TIM effectively serve its homogenizing function, enabling lateral heat spreading across the module to suppress local temperature peaks.

It should be noted that an ideal TIM should simultaneously exhibit both a high rate of cooling (see Figure 2d) and a low spread of $\Delta T$, which is achieved in our case through 40%BN–60%TPU composites. The physical mechanism underlying these trends can be understood in terms of the competition between two thermal processes: (1) vertical heat sinking into the cooling manifold, and (2) lateral thermal diffusion across adjacent cell interfaces via the TIM. At low diffusivity (e.g., PLA), both mechanisms are weak, resulting in symmetric but slow heat loss across the module. As diffusivity increases without reaching the percolation threshold (e.g., GPLA, 20%BN–TPU), lateral diffusion becomes anisotropic and non-uniform, which can amplify local hotspots, especially near the geometric or thermal boundaries of the module. Once the TIM achieves sufficiently high diffusivity (e.g., 40%BN–TPU), lateral heat transfer dominates, allowing efficient redistribution of thermal energy from warmer to cooler regions and thus reducing $\Delta T$. Edge effects are also evident in these data due to the finite size of the battery module and the small number of monitored cells (N = 6). Cells located at the periphery (e.g., Battery-1 or Battery-6) often show the largest deviation from the average temperature. This behavior is consistent with the expectation that edge cells experience asymmetric thermal environments—differing airflow exposure, partial insulation, or varying proximity to the cooling manifold boundaries. To further elucidate the relationship between TIM thermal properties and inter-cell thermal uniformity, we extracted the peak values of $\Delta T$ for each configuration and plotted them against the measured thermal diffusivity of the corresponding materials, as shown in Figure 4f. Error bars represent the standard deviation of the peak $\Delta T$ values across three replicates for each TIM.

To evaluate the dynamic thermal response of the battery module under practical operating conditions, COMSOL Multiphysics simulations were conducted across a range of charge–discharge rates (1C, 2C, 4C, and 8C) [2,57,58,59,60,61,62]. A six-cell Li-ion battery module was subjected to 10 continuous charge–discharge cycles at varying C-rates (1C, 2C, 4C, and 8C), followed by a 2.5-hour cooling period with various TIMs placed between the cells, as shown in Figure 5. Figure 5a presents a representative image of the temperature profile at the end of the first cycle using PLA as the TIM. The thermal distribution among the cells is relatively symmetric. However, we notice an increase in the cell temperature, as evidenced by different colors of the liquid-cooling manifold (hot pink) compared to the cells (golden yellow). Figure 5b–f depict the temperature profiles at the end of the tenth cycle for the different TIMs. For PLA (Figure 5b), a significant temperature increase is evident, with pronounced hotspots in TIM toward the left edge of the cell. The GPLA and TPU TIMs (Figure 5c and d) demonstrate improved thermal performance relative to PLA, but still exhibit notably asymmetry as noticed from different colors for the liquid-cooling manifold (purple) and the cells (orange for GPLA and hot pink for TPU). These results are consistent with their moderate thermal diffusivities, which enable some degree of vertical heat sinking but are insufficient for complete lateral thermal equilibration. The 20% BN–TPU composite (Figure 5e) provides a further reduction in maximum temperature and a modest improvement in uniformity. The most striking improvement is observed with the 40% BN–TPU composite (Figure 5f), where the temperature field is both lower in magnitude and more homogeneous with the liquid-cooling manifold and the cells both showing similar color.

To analyze more quantitatively, we present the average temperature profiles in Figure 6. Each subplot compares the temperature rise for different TIM configurations: PLA, GPLA, TPU, 20% BN–TPU, and 40% BN–TPU. Across all conditions, a clear and systematic increase in temperature with increasing C-rate is observed, consistent with the enhanced internal heat generation governed by the relation $Q\propto I^{2}R$, where I is the current and R is the internal resistance of the cell. At low C-rates (e.g., 1C in Figure 6a), the thermal load is moderate leading to maximum temperature ≤ 22 °C. At 1C, all TIMs appear reasonably effective in preventing excessive temperature rise, although subtle differences in performance are already discernible with 40% BN–TPU showing the best performance, as expected. As the C-rate increases (particularly at 4C and 8C) in Figure 6b-d, the distinctions among TIMs become significantly more pronounced. Notably, the PLA TIM, which exhibits the lowest thermal diffusivity among the tested materials, results in the highest average temperatures at every C-rate. In contrast, the 40% BN–TPU composite consistently demonstrates the lowest temperature profiles across all C-rates, highlighting the role of enhanced thermal diffusivity in facilitating vertical heat sinking and efficient thermal management. Furthermore, all materials exhibit a thermal staircase in peak temperature with increasing C-rate, indicating cumulative heating effects due to incomplete thermal recovery between successive cycles.

At the highest C-rate of 8C, PLA exhibited the largest $T_{avg}$, reaching slightly above 50 °C over the course of the cycles (Figure 7a). In contrast, the 40%BN–TPU TIM, which demonstrated the highest experimental thermal diffusivity, limited the thermal rise to below 50 °C, despite identical electrochemical and convective boundary conditions. Figure 7b-e present the average inter-cell temperature deviation, denoted as $\Delta T_{\mathrm{avg}}\left(t\right)$, for each TIM configuration across charge–discharge cycles at different C-rates. Across all C-rates, a clear hierarchy in $\Delta T_{\mathrm{avg}}$ emerges. Low-diffusivity TIMs such as PLA and GPLA exhibit the lowest inter-cell deviations in contrast to high-diffusivity TIMs (similar to results discussed in Figure 4). As discussed earlier, this can be rationalized in terms of the balance between vertical heat transfer to the manifold vs. the lateral heat spreading. While high diffusivity TIMs enable faster heat decay through vertical heat transfer (cf. Figure 3d), there is asymmetric lateral heat diffusion resulting in higher inter-cell deviations.

4. Conclusions

We systematically evaluated the thermal performance of graphene-PLA (GPLA) and boron nitride (BN)-loaded TPU composites as thermal interface materials (TIMs) for lithium-ion battery modules. Material characterization confirmed that the incorporation of high thermal conductivity fillers significantly enhanced both thermal stability and thermal diffusivity. Through a combination of analytical modeling, finite-element simulations in COMSOL Multiphysics, and experimental thermal cycling studies, we established a direct inverse relationship between the TIM thermal diffusivity and the heat decay rate of the battery module. Our analysis further revealed that optimal TIMs must balance vertical heat sinking into the cooling infrastructure with lateral thermal diffusion across adjacent cells to minimize inter-cell temperature gradients. Among the materials studied, the 40%BN–60%TPU composite demonstrated superior performance, achieving both rapid cooling and thermal uniformity. These findings underscore the potential of filler-engineered, additive-manufactured polymer composites for scalable, efficient thermal management in high-power lithium-ion battery applications.