Computational investigation of an adaptive PID-controlled PCM-enhanced hybrid thermoelectric-air cooling system for lithium-ion battery packs in electric vehicles

1 Introduction

The global transition toward sustainable transportation has accelerated the adoption of Electric Vehicles (EVs) and Hybrid Electric Vehicles (HEVs), driven by the urgent need to reduce greenhouse gas emissions, urban air pollution, and dependence on fossil fuels (Grebtsov and Kubasov, 2024; Nandakishora et al., 2025). Lithium-ion batteries have emerged as the dominant energy storage technology for EVs owing to their high energy density, favorable power-to-weight ratio, long cycle life, and low self-discharge rate (Panmuang et al., 2024; Dan et al., 2023). A lithium-ion battery generates electricity through the reversible movement of lithium ions between an anode and a cathode separated by an electrolyte and a separator (Mukane et al., 2021; Nitesh and Ravichandra, 2021). These characteristics make lithium-ion batteries well suited for vehicle propulsion systems, enabling longer driving ranges and faster charging capabilities (Hamednia et al., 2023).

However, the thermal behavior of lithium-ion battery packs during high-rate charging and discharging cycles poses significant challenges to vehicle performance, safety, and battery longevity (Weragoda et al., 2023; Zhao et al., 2022). Excessive heat generation within battery packs originates from three primary mechanisms: (i) irreversible Joule heating due to the internal resistance of cells, tabs, and interconnects; (ii) electrochemical polarization at the electrode–electrolyte interfaces; and (iii) reversible entropic heat associated with the entropy change during intercalation and deintercalation reactions (Rahman et al., 2025; Bernardi et al., 1985). If this heat is not effectively dissipated, the cell temperature can rise beyond the safe operating window (typically 25 °C–60 °C), leading to accelerated capacity fade, increased internal resistance, electrolyte decomposition, separator shrinkage, and in extreme cases, thermal runaway with the risk of fire or explosion (Togun and Basem, 2025; Das et al., 2025). Furthermore, non-uniform temperature distribution across cells within a pack exacerbates cell-to-cell state-of-charge imbalance, reducing the effective pack capacity and service life (Awasthi et al., 2025).

Numerous Battery Thermal Management System (BTMS) strategies have been investigated to address these thermal challenges. Air cooling remains the simplest and most cost-effective approach. Wankhede et al. (2025) conducted experimental investigations on air-cooled BTMS for formula student EVs and demonstrated effective temperature control at low to moderate discharge rates; however, the heat removal capacity became insufficient at elevated currents. Zhao et al. (2022) proposed a trapezoid air-cooling configuration that improved temperature uniformity but required high fan power at increased airflow velocities.

Liquid cooling systems offer superior heat transfer coefficients compared to air cooling. Chang et al. (2025) integrated battery liquid cooling with HVAC dehumidification and heat recovery, achieving effective thermal regulation; however, the absence of a proper airflow delivery system resulted in non-uniform temperature distribution. Rahmani et al. (2025) simulated cold plate cooling with water–ethylene glycol mixtures under diverse driving cycles in MATLAB Simulink, demonstrating adaptive BTMS strategies; however, higher charge rates significantly increased the refrigerant power consumption. Saxena et al. (2025) combined PCM with liquid cooling using serpentine and zig-zag cold plates coupled with Computational Fluid Dynamics (CFD) modeling, achieving an 83.9% reduction in pumping power; however, low coolant temperatures created local temperature gradients leading to uneven heat distribution. While effective, liquid cooling introduces system complexity, additional weight, and potential leakage risks (Zhao et al., 2023).

Phase Change Material (PCM)-based systems provide passive thermal buffering through latent heat absorption during the solid–liquid phase transition. Werlich et al. (2025) demonstrated multi-objective optimization of a prismatic LFP battery pack with PCM, identifying Pareto-optimal configurations; however, suboptimal thermal conductivity led to increased temperatures at high power levels. Awasthi et al. (2025) reviewed PCM-based hybrid BTMS for prismatic cells, concluding that pure PCM systems require supplementary active cooling for sustained high-load operation. Togun and Basem (2025) investigated hybrid nanofluids combined with PCM in cylindrical pack cooling and demonstrated improved heat transfer; however, the complexity of nanofluid preparation and long-term stability remain practical challenges.

Thermoelectric cooling using Peltier modules has emerged as a compact, solid-state alternative that operates without moving parts or refrigerants. Bhattacharyya et al. (2025) optimized thermoelectric coolers for light EVs using a coupled finite-difference axial cell model, balancing efficiency and mechanical stability; however, the model lacked transient thermal adaptability under dynamic operating conditions. Panmuang et al. (2024) experimentally investigated water cooling and thermoelectric cooling techniques for batteries and demonstrated that TECs provide precise temperature control at low heat loads. Nandakishora et al. (2025) modeled a Peltier-hybrid cooling system for EV batteries and achieved improved thermal regulation; however, the energy consumption of the TEC unit was not quantified relative to the overall system efficiency.

Hybrid approaches combining multiple cooling methods have shown promising results. Karimi et al. (2022) combined paraffin PCM with aluminum heat sinks, copper heat pipes, and air cooling for a high-power lithium-ion module; however, insufficient cell spacing led to hot spots and uneven temperature distribution. Narimani et al. (2024) presented a unified thermal management framework integrating HVAC and battery cooling loops using a gain-scheduling PID controller; however, the battery experienced overheating as the circulating air could not remove sufficient heat. Qin et al. (2025) implemented a cooperative thermal management system with a dual-evaporator air conditioning design for EVs; however, the system struggled to cool the cabin and battery independently due to their different temperature requirements. Dai and Long (2023) designed a fast-charging BTMS with direct refrigerant cooling, achieving a battery temperature of 45 °C; however, large coolant temperature differences between inlet and outlet affected temperature uniformity.

Control-based approaches have also been explored for optimizing cooling system operation. Hajidavalloo et al. (2023) demonstrated non-linear model predictive control-based thermal management for cold-weather conditions; however, the high computational effort increased internal heat generation. Zhuang et al. (2021) developed an intelligent thermal management system that achieved 25%–40% energy savings through optimized airflow; however, the system did not investigate multi-mode switching among fundamentally different cooling technologies or PCM integration. Guo and Jiang (2021) proposed a novel EV thermal management system based on refrigerant cooling and heating; however, localized hotspots and temperature gradients persisted across the battery pack.

Despite these advances, several critical research gaps remain. First, the combined integration of PCM, thermoelectric cooling, and forced air convection as a triple-hybrid system has not been systematically investigated for cylindrical 18,650 cell packs. Second, adaptive multi-mode control strategies that dynamically switch among five or more cooling modes based on real-time thermal conditions and rate-of-temperature-change are largely unexplored. Third, comprehensive energy efficiency analysis comparing the Coefficient of Performance (COP) across multiple cooling configurations under identical operating conditions is lacking. Fourth, parametric sensitivity analysis identifying the single most influential design parameter for hybrid BTMS optimization has received limited attention. Fifth, the coupled effect of battery aging on thermal management performance has not been addressed in the context of hybrid PCM-TEC systems.

To address these gaps, the present study proposes a computational framework with the following specific objectives:

• Develop and validate a calibrated lumped-parameter thermal model for a 13S8P lithium-ion battery pack against published experimental data.
• Evaluate five distinct cooling modes—natural convection, forced convection, thermoelectric convection, hybrid (forced + TEC), and a novel PCM-enhanced hybrid—across discharge rates from 5A to 30A and airflow speeds from 2 to 15 m/s.
• Design an adaptive PID controller that dynamically selects the optimal cooling mode and fan speed based on real-time thermal conditions, and compare its performance against fair baselines including a threshold-triggered strategy.
• Conduct a comprehensive energy efficiency (COP) analysis comparing all active cooling configurations under identical operating conditions.
• Perform parametric sensitivity analysis across four key design parameters to identify the most influential factor for BTMS optimization, and evaluate the impact of battery aging and driving cycle variations on thermal performance.

The significance of this study lies in providing a comprehensive computational framework that simultaneously addresses thermal performance, energy efficiency, adaptive control, and design sensitivity for EV battery packs. The paper is organized as follows: Section 2 describes the mathematical model and methodology, Section 3 presents the results and discussion, and Section 4 concludes the work with practical implications and future recommendations.

2 Mathematical model and methodology

2.1 Battery pack configuration

The thermal model is based on a 13-series 8-parallel (13S8P) battery pack comprising 104 Panasonic NCR18650BD lithium-ion cells arranged in a zig-zag pattern, consistent with the experimental setup reported by Krishnam Raju and Vedaraj (2025). The pack provides a nominal voltage of 13 × 3.6 = 46.8 V and a total capacity of 8 × 3.2 = 25.6 Ah. The detailed cell specifications are presented in Supplementary Table S1. Seven thermocouple locations (T1–T7) are defined at critical positions within the pack to capture the spatial temperature distribution resulting from the zig-zag arrangement and non-uniform airflow access. The cells are interconnected using nickel-plated tabs and a pulse spot-welding machine, with a Heltec 13S 40A/50A Battery Management System (BMS) for cell balancing, overcharge protection (4.25 V/cell), and over-discharge protection (2.8 V/cell).

2.2 Heat generation model

The heat generated by each cell during operation is computed using the Bernardi equation (Bernardi et al., 1985), as given in Equation 1:where I cell = I pack / n p is the per-cell current (n p = 8 for the 13S8P configuration), R eff T is the temperature-dependent effective internal resistance, T is the absolute temperature in Kelvin, and d V OC / d T is the entropic coefficient. The first term represents irreversible Joule heating (the dominant contribution at high C-rates), while the second term accounts for the reversible entropic heat that can be either exothermic or endothermic depending on the state of charge (SOC). The entropic coefficient is approximated as given in Equation 2:

Q gen , cell = I cell 2 · R eff T + I cell · T · d V OC d T (1)

d V OC d T = − 0.00025 · 1 − 2 · SOC V / K (2)

The effective internal resistance under real pack operating conditions is significantly higher than the single-cell datasheet value (35 mΩ) due to several factors: (i) spot-weld contact resistance at each terminal joint (+5–15 mΩ per cell), (ii) nickel strip resistance along the current path (+2–5 mΩ), (iii) rate-dependent resistance increase at high C-rates (+20–40%), and (iv) interfacial resistance at BMS connections. This phenomenon is well documented; Gao et al. (2021) and Wu et al. (2023) reported that measured pack-level resistance is typically 2–3× the single-cell datasheet value. The calibrated value is consistent with an approximate additive budget per cell: ∼35 mΩ intrinsic DC resistance (datasheet, 1 kHz AC value is lower), ∼10–15 mΩ from two spot-welded nickel-strip terminations, ∼3–5 mΩ distributed nickel-strip and busbar resistance, and a ∼20%–30% rate- and state-of-charge-dependent increment at the C-rates examined here, together yielding an effective 70–85 mΩ range that brackets the calibrated 73.5 mΩ. The calibration therefore recovers a physically attributable quantity rather than an arbitrary fitting factor. The calibrated effective resistance is expressed as given in Equation 3:where R base = 73.5 mΩ is the calibrated effective resistance at 25 °C (corresponding to a calibration factor of 2.10× the datasheet value), and α = 0.001 °C⁻¹ is the temperature coefficient. The total pack-level heat generation additionally includes 8% for tab and interconnect losses and 0.8 W for BMS quiescent power dissipation, as given in Equation 4:

R eff T = R base · 1 + α · ⁡ max 0 , T C − 25 (3)

Q gen , pack = n cells · Q gen , cell · 1.08 + 0.8 W (4)

2.3 Thermal energy balance

The lumped-parameter energy balance for the battery pack is expressed as:where m eff = n cells · m cell + m housing is the effective thermal mass (including the 0.8 kg enclosure, tabs, and BMS PCB), C p , eff = 750 J/(kg·K) is the effective specific heat capacity of the pack assembly, and Q cool is the total heat removed by the active cooling system. Equation 5 is solved using explicit Euler time integration with a 1-s timestep.

m eff · C p , eff · d T d t = Q gen − Q cool (5)

It is emphasised that the present work employs a lumped-parameter nodal network, not a computational fluid dynamics (CFD) formulation. No spatial mesh, geometry discretisation, or Navier–Stokes solution is involved; consequently mesh-independence statistics are not applicable. The computational domain consists of seven thermal nodes corresponding to the seven monitored cell locations (T1–T7), each governed by the per-location temperature update presented in Section 2.6. The equivalent boundary and initial conditions are: (i) a uniform initial condition T j 0 = T amb = 25 °C for all nodes; (ii) a convective (Robin-type) surface condition at each node, q j ″ = h eff f c , j T j − T amb, with h eff supplied by the active cooling model of Section 2.4; (iii) a volumetric source term Qgen fh,j from the heat generation model of Section 2.2; (iv) an adiabatic external enclosure boundary, so that all generated heat leaves only through the modelled cooling path; and (v) a fixed far-field ambient reservoir at 25 °C. Convergence is governed by the temporal step alone, which was verified by step halving. Spatially resolved airflow within the enclosure is outside the scope of this formulation and is identified as future CFD work in Section 3.12.

2.4 Cooling system models

Five cooling modes are evaluated in this study, ranging from passive natural convection to the proposed novel PCM-enhanced hybrid system. The models are described below.

2.4.1 Natural convection

In the absence of active cooling, heat dissipation occurs through natural convection following Newton’s law of cooling, as given in Equation 6:where h nat = 5.5 W/(m²·K) is the natural convection heat transfer coefficient for horizontal cylinders in still air (typical range: 3–10 W/(m²·K) depending on orientation and surface temperature), and A exp = 0.085 m² is the effective exposed surface area of the pack enclosure. No external power is consumed in this mode.

Q nat = h nat · A exp · T pack − T amb (6)

2.4.2 Forced convection

A 12V axial fan circulates air over the battery pack at velocities ranging from 2 to 15 m/s. The heat transfer coefficient is computed using the Churchill–Bernstein correlation (Churchill and Bernstein, 1977) for crossflow over a cylinder:where Re = ρ v air d / μ is the Reynolds number based on cell diameter, and Pr = 0.71 for air. The air properties are evaluated at the film temperature: density ρ = 1.225 kg/m³, dynamic viscosity μ = 1.81 × 10 − 5 Pa·s, and thermal conductivity k air = 0.026 W/(m·K). To account for duct losses, bypass flow around the zig-zag cell arrangement, and non-uniform air distribution, a duct effectiveness factor is applied as given in Equation 8:

Nu = 0.3 + 0.62 · Re 1 / 2 · ⁡ Pr 1 / 3 1 + 0.4 / ⁡ Pr 2 / 3 1 / 4 · 1 + Re 282000 5 / 8 4 / 5 (7)

η duct = min 0.25 + 0.025 · v air , 0.65 (8)

This yields an effective heat transfer coefficient h eff = Nu · k air / d · η duct within the published range of 20–65 W/(m²·K) for cylinder banks. The forced convection heat removal is given by Equation 9:

Q forced = h eff · A exp · T pack − T amb (9)

The fan power consumption is P fan = 3.0 W.

2.4.3 Thermoelectric (Peltier) convection

The active cooling unit employs a TEC1-12706 Peltier module (40 mm × 40 mm). The cold-side cooling capacity is modeled using the Peltier effect:where S = 0.0527 V/K is the Seebeck coefficient, R TEC = 1.98 Ω is the electrical resistance, K TEC = 0.55 W/K is the thermal conductance (all from the TEC1-12706 datasheet), I TEC = 4.0 A is the operating current, and T c and T h are the cold-side and hot-side absolute temperatures, respectively. The cold side is attached to a 40 mm × 40 mm × 15 mm aluminum heat sink, while the hot side is cooled by an 80 mm × 80 mm fan attached to a 90 mm × 60 mm × 25 mm heat sink. A delivery efficiency of 50% accounts for thermal losses in the insulated duct between the TEC cold side and the battery pack enclosure. The TEC power consumption is P TEC = V TEC · I TEC = 12 × 4 = 48 W.

Q c = S · I TEC · T c − 1 / 2 I 2 R TEC · R TEC − K TEC · T h − T c (10)

2.4.4 Hybrid (forced + TEC) convection

The hybrid mode combines both forced convection and thermoelectric cooling. When both systems operate simultaneously, the pre-cooled air from the TEC enhances the forced convection effectiveness. However, the combined heat removal is subject to diminishing returns because both systems act on the same thermal boundary. A synergy factor of 0.82 accounts for this overlap, as given in Equation 11:

Q hybrid = 0.82 · Q forced + Q TEC (11)

The total power consumption is P hybrid = P fan + P TEC = 51 W.

2.4.5 PCM-enhanced hybrid convection (novel)

The proposed PCM-enhanced hybrid mode integrates RT42 paraffin wax (Rubitherm GmbH) with the hybrid cooling system. The RT42 has a solidus temperature of 38 °C, a liquidus temperature of 44 °C, and a latent heat of fusion of L = 165 kJ/kg. The total PCM mass is 1.25 kg (12 g per cell encapsulated in thin aluminum shells surrounding each cell group). The Stefan number, which characterizes the ratio of sensible to latent heat, is given by Equation 12:

Ste = C p , PCM · Δ T L = 2000 · 6 165000 = 0.073 (12)

The Stefan number expresses the ratio of sensible to latent thermal energy stored during phase change. A value of Ste = 0.073 ≪ 1 carries three physical consequences. First, latent storage exceeds sensible storage by roughly a factor of fourteen, so almost all absorbed energy is consumed by melting rather than by raising the PCM temperature. Second, the melt front advances slowly relative to thermal diffusion, so the quasi-steady Stefan approximation used in Equation 13 is valid and sensible-heat transients within the PCM may be neglected. Third, and most importantly for thermal management, the PCM behaves as a near-isothermal heat sink: the pack temperature is buffered within the 38 °C–44 °C melting band while the latent reservoir is being consumed, which is precisely the plateau behaviour observed in Section 3.3. A material with Ste ≳ 1 would instead heat continuously and provide no such plateau. During the mushy zone (T solidus ≤ T ≤ T liquidus + 3 K and melt fraction ϕ < 1), the PCM absorbs latent heat at a rate:where the melt fraction rate d ϕ / d t depends on the local temperature and is capped to ensure numerical stability. Below the solidus temperature, the PCM remains solid and contributes only as additional thermal mass. Above the liquidus, the fully melted PCM provides sensible heat storage at C p , liquid = 2400 J/(kg·K). The total latent heat capacity is 1.25 × 165 = 206 kJ. The total PCM-enhanced hybrid cooling capacity is given by Equation 14:

Q PCM = m PCM · L · d ϕ d t (13)

Q PCM ‐ hybrid = Q hybrid + Q PCM (14)

Critically, the PCM is a passive component, consuming no additional electrical energy beyond the hybrid mode. This distinguishes the PCM-enhanced hybrid from approaches that require active refrigeration to re-solidify the PCM.

2.5 Adaptive PID controller

An adaptive PID controller is designed to dynamically select the optimal cooling mode and fan speed based on the real-time thermal state of the battery pack. The control output is computed using Equation 15:where e t = T avg t − T target is the temperature error relative to the target of 35 °C, and K p = 1.5, K i = 0.05, K d = 0.8 are the proportional, integral, and derivative gains, respectively. The derivative term uses the actual rate of temperature change rather than the error derivative, providing anticipatory control that responds to thermal transients. Anti-windup protection is applied by clamping the integral term to − 100 , 100. The control output determines the active cooling mode according to the thresholds in Table 1. Within the forced convection and hybrid modes, the fan speed is also adaptively computed as v air = v min + v max − v min · ⁡ min u / 30 , 1, providing continuous modulation between 2 and 15 m/s.

u t = K p · e t + K i ∫ 0 t e τ d τ + K d · d T avg d t (15)

TABLE 1. Adaptive PID controller mode selection thresholds and parameters.

2.6 Spatial temperature variation model

To capture the non-uniform temperature distribution within the zig-zag cell arrangement, spatial heat generation factors (f h , j) and cooling effectiveness factors (f c , j) are assigned to each of the seven thermocouple locations (j = 1 , … , 7). Center cells (e.g., T3) experience higher heat accumulation (f h , 3 = 1.22) due to reduced air access and thermal coupling with neighboring cells, while peripheral cells (e.g., T7) exhibit lower heating (f h , 7 = 0.85) and enhanced cooling (f c , 7 = 1.18). Additionally, center cells with f h , j > 1.10 have their cooling effectiveness reduced by a factor of 0.65 to model the limited airflow penetration into the pack interior. The per-location temperature update is:

T j n + 1 = T j n + Q gen · f h , j − Q cool · f c , j · Δ t m cell , eff · C p , eff (16)

2.7 Model calibration and validation

The model was calibrated against the experimental data reported by Krishnam Raju and Vedaraj (2025), specifically the maximum temperature progression of the 13S8P battery pack under natural convection at 30A discharge over 60 min. The experimental reference data (from Table 5b of the cited work) provides temperature values at 10-min intervals: 28 °C (0 min), 56 °C (10 min), 78 °C (20 min), 96 °C (30 min), 112 °C (40 min), 126 °C (50 min), and 138 °C (60 min). The calibration parameters (R base and α) were optimized using the Nelder–Mead simplex algorithm to minimize the Root Mean Square Error (RMSE) between simulated and experimental maximum temperatures. The optimization yielded R base = 73.5 mΩ (calibration factor 2.10) and α = 0.001 °C⁻¹.

2.8 Simulation parameters and procedure

All simulations were implemented in Python 3.11 (NumPy, SciPy, Matplotlib) and employed explicit Euler time integration with a 1-s timestep over a 3600-s (60-min) duration; timestep independence was confirmed by halving the step with no appreciable change in peak temperature. Constant-current cases are fully deterministic. The UDDS and aggressive cycles contain stochastic micro-trip variability and are therefore generated from a seeded pseudo-random number generator; all results derived from these cycles are reported as the mean ± standard deviation of 20 repetitions (seeds 1000–1019), with every cooling strategy evaluated on the identical seed sequence so that comparisons are strictly paired. The ambient temperature was set to 25 °C with an initial pack temperature of 25 °C. Four driving cycle profiles were implemented based on United States Environmental Protection Agency (EPA) driving pattern characteristics: (i) constant current at 20A, (ii) an Urban Dynamometer Driving Schedule (UDDS) pattern urban cycle with an average current of 11.2A and peak of 30A featuring frequent acceleration–cruise–deceleration–idle micro-trips, (iii) a Highway Fuel Economy Test (HWFET) pattern highway cycle with an average of 17.8A and peak of 26A characterized by steady-state cruising with gentle variations, and (iv) an aggressive cycle with an average of 16.6A and peak of 30A representing high-power demand driving. Four battery pack configurations (13S8P, 14S6P, 10S10P, 16S6P) and five aging factors (1.0, 1.2, 1.5, 1.8, 2.0, corresponding to fresh through approximately 2000 charge–discharge cycles) were evaluated. The complete simulation parameters are summarized in Table 2.

TABLE 2. Summary of simulation parameters.

3 Results and discussion

3.1 Model validation

Figure 1 presents the comparison between the simulated and experimental maximum temperatures under natural convection at 30A. So that the validation may be reproduced and checked independently of the source publication, the complete calibration dataset is reproduced here in full.

Table 3.

These seven points constitute the entire basis of the calibration; the experimental column is taken from the maximum-temperature progression reported by Krishnam Raju and Vedaraj (2025) and the simulated column is deterministic and exactly reproducible from the accompanying source code. From these values the model achieves an RMSE of 8.0 °C, a mean absolute error of 6.2 °C, and an R 2 of 0.9515. The mean absolute percentage error is 6.2% when computed over all seven points and 7.3% when the initial condition is excluded; we state both because the t = 0 point carries zero error by construction, both curves being initialised at the same ambient temperature, and its inclusion therefore flatters the percentage metric. The stricter 7.3% figure is the appropriate one for judging predictive accuracy. The prediction errors range from −6.8 °C at 20 min to +17.1 °C at 60 min. The slight underestimation at 10–30 min is attributed to the simplified treatment of thermal mass redistribution during the initial transient phase, where the internal cells heat faster than the pack housing. The overestimation at 60 min arises from the linear resistance–temperature model, which does not capture the saturation of resistance increase at very high temperatures. Despite these localized discrepancies, the overall validation demonstrates that the calibrated model is suitable for comparative analysis across cooling modes under the same operating conditions.

3.2 Comparative analysis of five cooling modes

Figure 2 illustrates the maximum temperature and temperature non-uniformity profiles for all five cooling modes at a constant 20A discharge rate. Table 3 summarizes the quantitative results at 60 min. Natural convection reaches 80.1 °C, significantly exceeding the 60 °C safety limit, with a temperature non-uniformity of 21.64 K. This confirms that passive cooling alone is inadequate for the 13S8P pack at moderate discharge rates, consistent with the experimental findings of Krishnam Raju and Vedaraj (2025) and the simulations of Wankhede et al. (2025).

TABLE 3. Thermal performance comparison of five cooling modes at constant 20A discharge (60 min).

Forced convection at 10 m/s reduces the peak temperature to 56.3 °C, just below the safety limit, while consuming only 3.0 Wh of energy over the 60-min period. The Thermoelectric (TEC) convection operating independently reaches 73.3 °C despite consuming 48.0 Wh—a 16× higher energy cost than forced convection for inferior thermal performance. This confirms the well-documented limitation of standalone TEC modules at sustained moderate heat loads (Bhattacharyya et al., 2025; Panmuang et al., 2024).

The hybrid (Forced + TEC) mode achieves 54.4 °C with 51.0 Wh energy consumption. The PCM-enhanced hybrid mode achieves the lowest peak temperature of 40.2 °C, representing a 26.1% reduction compared to the standard hybrid and a 49.8% reduction compared to natural convection. The temperature non-uniformity is also significantly reduced: 8.87 K for the PCM-enhanced hybrid versus 21.45 K for the standard hybrid, representing a 58.6% improvement in thermal uniformity. This improvement is attributed to the latent heat absorption of the RT42 PCM during the mushy zone (38 °C–47 °C), which buffers the temperature rise precisely in the critical operating range.

The uniformity gain is a direct consequence of where the pack is hottest. As described in Section 2.6, interior nodes such as T3 carry a heat-generation factor f h , 3 = 1.22 while receiving a cooling factor of only f c , 3 = 0.75, further reduced by the 0.65 airflow-penetration penalty applied to interior locations; peripheral node T7 carries f h , 7 = 0.85 with f c , 7 = 1.18. Interior cells are hotter because they are geometrically shielded—each is surrounded on all sides by neighbouring cells that both obstruct the airflow path and act as thermal sources in their own right, whereas peripheral cells have a free convective face and fewer heated neighbours. Purely convective cooling cannot correct this imbalance: since Q ∝ h T − T amb with h set by the local flow, raising the fan speed scales the cooling of every node in proportion and leaves the gradient between them essentially intact. This is why forced convection and the standard hybrid retain Δ T values of 20.09 K and 21.45 K respectively despite very different absolute temperatures. The PCM operates on a different principle. Because melting proceeds isothermally, the latent sink engages preferentially and most vigorously wherever the local temperature first enters the 38 °C–44 °C band—that is, at the hottest interior nodes—and remains dormant at cooler peripheral nodes. The absorption is thus self-targeting, clamping the hot nodes toward the melting plateau while the cool nodes are unaffected, which compresses Δ T from 21.45 K to 8.87 K, a 58.6% reduction obtained with no change in fan power.

3.3 Discharge rate analysis

Figure 3 presents the peak temperature at 60 min across discharge rates from 5A to 30A for all five cooling modes. Table 4 provides the complete numerical results. At 5A, all active cooling modes maintain temperatures near ambient (24 °C–29 °C) with negligible differences, as the heat generation rate (∼2.5 W) is easily managed by any cooling method. The PCM provides no additional benefit at 5A (26.1 °C for both hybrid and PCM-hybrid) because the temperature remains below the solidus of 38 °C.

TABLE 4. Peak temperature (°C) at 60 min for varying discharge rates.

As the discharge rate increases, the differences among cooling modes become pronounced. At 15A, the PCM-enhanced hybrid (38.5 °C) begins to separate from the standard hybrid (40.7 °C) as the temperature enters the phase-change window. At 20A, the PCM-enhanced hybrid maintains 40.2 °C while the standard hybrid reaches 54.4 °C. At 30A, the advantage is dramatic: the PCM-enhanced hybrid achieves 48.9 °C compared to 97.5 °C for the standard hybrid, representing a 49.8% reduction. The nonlinearity follows directly from a latent-energy budget. The installed PCM stores m PCM L = 1.25 × 165 = 206 kJ of latent capacity. At 20A the pack dissipates approximately 26 W, so the reservoir represents about 132 min of buffering—longer than the 60 min test—and the pack is held near the melting plateau throughout, which is why the PCM curve is nearly flat in Figure 3. At 30A the dissipation rises to roughly 58 W (Joule heating scales as I 2), reducing the buffering interval to about 59 min: the reservoir is very nearly exhausted at exactly the point the test ends, which is why 30A is where the PCM advantage is largest yet visibly beginning to erode. Below 15A the pack never reaches the 38 °C solidus, no melting occurs, and the PCM contributes nothing beyond a small sensible thermal mass—accounting for the identical hybrid and PCM-hybrid values at 5A and 10A in Table 4. The benefit is therefore bounded by installed latent capacity and by the residence of pack temperature inside the melting band, not by a general improvement in heat transfer. At 30A, natural convection (151.7 °C) and TEC alone (143.7 °C) reach extreme temperatures that would cause thermal runaway in practice. Notably, only the PCM-enhanced hybrid maintains the temperature below the 60 °C safety limit at all discharge rates up to 30A.

3.4 Airflow speed sensitivity

Figure 4 presents the effect of airflow speed (2–15 m/s) on peak temperature for all five cooling modes at 10A, 20A, and 30A; the corresponding numerical values are given in Supplementary Table S2. Two of the five modes are airflow-invariant by construction and are plotted as dashed reference lines. Natural convection involves no forced flow, so its heat transfer coefficient is fixed at h nat = 5.5 W/(m²·K) irrespective of the abscissa. Standalone thermoelectric convection delivers pre-cooled air through a fixed-duty delivery fan and therefore also does not vary with the swept fan speed. Explicit computation confirms this: the peak-temperature range of both modes across 2–15 m/s is 0.00 °C at every discharge rate (Supplementary Table S2). Their inclusion nevertheless establishes the upper thermal bound against which the three airflow-dependent modes are assessed.

The three airflow-dependent modes follow the scaling anticipated by the Churchill–Bernstein correlation of Equation 7. Because Nu ∝ Re 1 / 2 and Re ∝ v air, the convective coefficient scales as h ∝ v air 1 / 2, so equal increments of fan speed yield progressively smaller thermal returns. The computed data reproduce this square-root behaviour: for forced convection at 30A the peak temperature falls from 139.6 °C at 2 m/s to 98.0 °C at 10 m/s (41.6 °C for the first 8 m/s of increase) but only to 80.6 °C at 15 m/s (a further 17.4 °C for the next 5 m/s). The hybrid mode exhibits the same exponent with a parallel downward offset of 8 °C–11 °C contributed by the Peltier stage.

The PCM-enhanced hybrid departs from this scaling in a manner that is diagnostic of the latent-heat mechanism. Its sensitivity to airflow collapses precisely where the pack enters the melting band: the peak-temperature range across the full speed sweep is 5.13 °C at 10A, 1.32 °C at 20A, and 16.62 °C at 30A, compared with 8.45, 26.10, and 59.04 °C respectively for forced convection. At 10A the PCM curve is numerically identical to the standard hybrid because the pack never reaches the 38 °C solidus and the latent reservoir is never engaged. At 20A the pack sits inside the melting band across the entire speed range, so the near-isothermal latent sink dominates and airflow becomes almost irrelevant. At 30A the latent capacity is partially exhausted within the hour, so airflow dependence partially returns. This progression, invisible in the three-mode presentation, constitutes direct evidence that the observed benefit originates from phase change rather than from additional convective area.

From a design standpoint, diminishing returns become pronounced above approximately 12 m/s for all airflow-dependent modes, so further fan speed increases offer marginal thermal benefit at the cost of increased fan power, acoustic noise, and vibration.

3.5 Energy efficiency (COP) analysis

Figure 5 presents the instantaneous cooling power, cumulative energy consumption, and average COP for the four active cooling modes at 20A. The COP is defined as the ratio of heat removed to electrical energy consumed: COP = Q removed / P input. Forced convection exhibits the highest COP of 22.18, attributed to the low fan power consumption (3.0 W) relative to the substantial convective heat removal. This suggests that forced air cooling is the most energy-efficient approach for moderate thermal loads. The TEC operating independently has the lowest COP of 0.61. The disparity is thermodynamic in origin rather than a matter of component quality. Forced convection is not a heat pump: it transports heat down an existing temperature gradient, and the fan performs work only against fluid friction, not against the thermal gradient itself. Its COP is therefore unbounded by the Carnot limit and rises without penalty as the pack–ambient difference grows. A Peltier module, by contrast, pumps heat up a gradient and is bounded above by COP Carnot = T c / T h − T c, further degraded by two irreversibilities intrinsic to the device: Joule dissipation 1 2 I 2 R TEC distributed through the junction, and parasitic Fourier back-conduction K TEC T h − T c across the same elements that carry the Peltier current (Equation 10). Because both loss terms grow while the useful term S I T c does not, the achievable COP falls steeply as Δ T widens, which is why the standalone TEC consumes 48.0 Wh to achieve a poorer thermal outcome than a 3.0 W fan. The computed value lies within the 0.5–2.0 range reported by Su et al. (2021).

The improvement of the hybrid COP from 1.40 to 1.89 upon adding RT42 paraffin follows directly from the definition COP = Q removed / P input. The PCM is passive: it contributes a latent absorption term Q PCM to the numerator while adding nothing to the denominator, since the electrical draw remains the 51 W of the fan and Peltier stage. The 35% COP gain is thus obtained without any additional electrical expenditure, and is the thermodynamic counterpart of the temperature reduction reported in Section 3.3. The gain is, however, capacity-limited: once the 206 kJ latent reservoir is exhausted the numerator reverts to its hybrid value, so the enhancement is a transient advantage bounded by PCM mass rather than a permanent efficiency shift. The hybrid mode achieves a COP of 1.40, while the PCM-enhanced hybrid improves the COP to 1.89, representing a 35% enhancement. This improvement occurs because the PCM absorbs heat passively, effectively increasing the total cooling capacity without proportional energy cost.

3.6 Adaptive PID controller performance

Figure 6 presents the performance of the adaptive PID controller under the UDDS driving cycle compared to two baseline strategies. Table 5 summarizes the comparison. The always-on hybrid maintains the lowest peak temperature (38.67 °C ± 0.06 °C) but consumes 50.99 ± 0.00 Wh over the 60-min cycle by running continuously regardless of the thermal state. The threshold-triggered baseline, which activates the hybrid mode when T > 33 °C and deactivates when T < 30 °C, reduces energy consumption to 38.54 ± 0.03 Wh (24.4% ± 0.06% savings) with a modest temperature increase to 39.50 °C ± 0.06 °C.

TABLE 5. Adaptive PID controller comparison under the UDDS driving cycle. All values are mean ± standard deviation over 20 seeded repetitions (seeds 1000–1019); all three strategies were evaluated on the identical seed sequence.

The adaptive PID controller achieves the greatest energy savings at 27.07 ± 0.08 Wh (46.90% ± 0.14% versus always-on, 29.75% ± 0.09% versus the fair threshold-triggered baseline) while maintaining the peak temperature at 49.21 °C ± 0.08 °C, safely below the 60 °C limit with an 10.8 °C margin. The mechanism of the saving is a shift from time-based to state-based actuation. The threshold controller is binary and hysteretic: once triggered it commits the full 51 W of the fan and Peltier stage until the lower threshold is recovered, so it repeatedly overshoots the cooling actually required. The PID controller instead maps the instantaneous error, its integral, and the temperature rate onto a graded actuation ladder (Table 1), spending 0 W in natural convection whenever the thermal margin permits and escalating only to the least costly mode that arrests the observed rate of rise. Because the UDDS cycle spends a substantial fraction of its duration at low current with intermittent idle periods, the controller can remain in the zero-power state for extended intervals. Safety is preserved not by conservative always-on operation but by the derivative term, which responds to the rate of temperature rise before the absolute limit is approached; the resulting 10.8 °C margin is therefore achieved with 47% less energy than continuous operation. The mode selection timeline (Figure 6C) shows that the PID controller initially operates in natural convection mode during the low-temperature phase (first ∼25 min), progressively engaging forced convection as heat accumulates, and dynamically switching between natural and forced modes during the stop-and-go portions of the UDDS cycle. This intelligent thermal-aware energy management demonstrates that significant energy savings are achievable without compromising safety, which is particularly important for extending the driving range of EVs.

3.7 Battery pack configuration analysis

Figure 7 and Table 6 compare four pack configurations under the UDDS cycle with PCM-enhanced hybrid cooling. The 10S10P configuration (100 cells, 36.0 V, 32.0 Ah) achieves the lowest peak temperature, 33.3 °C ± 0.04 °C, against 38.1 °C ± 0.01 °C for the 13S8P baseline, 38.6 °C ± 0.02 °C for 14S6P and 38.9 °C ± 0.02 °C for 16S6P—a spread of 5.6 °C across the topologies examined.

TABLE 6. Battery pack configuration comparison under the UDDS cycle with PCM-enhanced hybrid cooling. Values are mean ± standard deviation over 20 seeded repetitions.

The ordering is governed entirely by the parallel cell count n p through the quadratic current dependence of Joule heating. For a fixed pack-level demand I pack, each cell carries I cell = I pack / n p, and the per-cell dissipation is Q cell = I cell 2 R = I pack 2 R / n p 2. Summing over the n s n p cells of the pack gives a total Joule dissipation Q pack = n s I pack 2 R / n p, so pack heat generation scales as n s / n p rather than with total cell count. Raising n p from 8 to 10 lowers the per-cell current at a 20A demand from 2.50A to 2.00A, a 20% reduction that yields a 36% reduction in per-cell heating (0.8 2 = 0.64) and a 20% reduction in total pack heating. Conversely, 14S6P and 16S6P both operate at n p = 6 and therefore at 3.33A per cell, which explains why they cluster together near 38.6 °C–38.9 °C despite differing in series count and total cell number; the small residual difference between them reflects their differing n s, consistent with the n s / n p scaling. Temperature non-uniformity follows the same ordering, falling from 11.24 K (13S8P) to 8.24 K (10S10P).

The design trade-off is that raising n p at fixed cell count lowers the pack voltage (36.0 V for 10S10P versus 46.8 V for 13S8P), which shifts burden onto the DC–DC converter and raises conductor currents downstream. These results highlight that pack topology significantly influences thermal performance: maximizing the parallel cell count reduces thermal stress at the expense of lower pack voltage, which may require a higher-ratio DC–DC converter.

3.8 Battery aging impact

Figure 8 presents the aging impact on thermal performance at 25A with PCM-enhanced hybrid cooling. As the effective internal resistance increases from the fresh value (1.0×, 73.5 mΩ) through intermediate aging states to the end-of-life value (2.0×, 147.0 mΩ), the peak temperature rises from 42.4 °C to 69.0 °C. The 60 °C safety limit is crossed at approximately 1500 cycles (60.1 °C, 1.8× resistance), indicating that the PCM-enhanced hybrid cooling designed for fresh batteries becomes insufficient approximately three-quarters through the battery’s expected lifespan. The temperature non-uniformity increases from 11.7 K (fresh) to 22.1 K (2000 cycles), exacerbating cell-to-cell imbalance. These results have important practical implications: (i) the cooling system must be designed for end-of-life resistance values, not fresh conditions; (ii) the PID controller thresholds should be recalibrated periodically as the battery ages; and (iii) increasing the PCM mass (the dominant parameter identified in Section 3.10) can extend the effective cooling lifespan.

3.9 Driving cycle analysis

Supplementary Figure S1 compares the thermal response under four driving cycles with PCM-enhanced hybrid cooling. The constant 20A discharge produces the highest peak temperature (40.2 °C) due to sustained continuous heat generation without any idle periods for cooling recovery. The UDDS urban cycle produces the lowest peak temperature (38.1 °C) owing to its lower average current (11.2A) and frequent idle periods (lasting 20 s every 180 s) that allow partial heat dissipation. The HWFET highway cycle (39.6 °C) exhibits steady-state thermal behavior with gentle temperature oscillations reflecting the highway speed variations. The aggressive cycle (39.4 °C) shows rapid temperature spikes during the 30A peak demand phases but recovers during the lower-current intervals. All four cycles remain below the 60 °C safety limit under PCM-enhanced hybrid cooling, demonstrating the robustness of the proposed system across diverse driving conditions.

3.10 Parametric sensitivity analysis

Figure 9 presents the sensitivity of the peak temperature to four key design parameters at 25A with PCM-enhanced hybrid cooling. Ranked by absolute influence on peak temperature, the four parameters order as: (1) PCM mass (−14.1%), (2) effective internal resistance (+4.8%), (3) airflow speed (−1.7%), and (4) TEC current (−0.5%). PCM mass is therefore roughly three times more influential than the next parameter and an order of magnitude more influential than the two actuator settings combined. The PCM mass exhibits the highest absolute sensitivity index of −14.1%, indicating that a 50% increase in PCM mass (from 1.25 kg to 1.875 kg) reduces the peak temperature by approximately 3.4 °C. Conversely, a 50% reduction in PCM mass increases the temperature from 42.4 °C to 47.0 °C. The effective internal resistance shows moderate sensitivity (+4.8%), confirming that accurate resistance calibration is important for model fidelity. The airflow speed (−1.7%) and TEC current (−0.5%) have comparatively minor influence, suggesting that the forced convection and TEC parameters are already near their effective operating points. This analysis identifies PCM mass as the primary design lever for thermal optimization, recommending that increasing the PCM loading from 12 g to 18 g per cell could yield the most cost-effective thermal improvement with a modest weight penalty of 0.625 kg at the pack level.

3.11 Comparison with existing literature

The proposed PCM-enhanced hybrid system demonstrates thermal performance consistent with or superior to published computational BTMS studies. Saxena et al. (2025) reported 8%–15% temperature reduction using PCM with liquid cooling; the present study achieves up to 49.8% reduction through the triple-hybrid integration of PCM, TEC, and forced air. Bhattacharyya et al. (2025) reported TEC COP values of 0.5–2.0 for light EVs, consistent with the present TEC COP of 0.61. Zhuang et al. (2021) achieved 25%–40% energy savings through intelligent thermal management; the present adaptive PID achieves 47.2% savings versus always-on operation and 30.1% versus a fair threshold-triggered baseline, both within a comparable range. Karimi et al. (2022) identified temperature non-uniformity as a critical challenge; the present PCM-enhanced hybrid reduces Δ T by 58.6% compared to the standard hybrid, addressing this limitation.

3.12 Limitations and future work

The present study has several limitations: (1) the lumped-parameter model does not resolve intra-cell temperature gradients; (2) the PCM-enhanced hybrid results are computational predictions requiring experimental validation; (3) the driving cycles are pattern-based approximations of EPA profiles; (4) a single ambient temperature (25 °C) was evaluated; (5) the PCM model uses a simplified mushy-zone formulation; (6) cell-to-cell thermal coupling is modeled through spatial factors; (7) the TEC COP is assumed constant; (8) no capacity degradation feedback was included; (9) the airflow pattern was not validated through CFD; (10) PID gains were not formally optimized; (11) the model was calibrated and validated against a single experimental condition (natural convection at 30A) drawn from a companion study by the present authors that was under review at the time of writing; the seven calibration points are therefore reproduced in Section 3.1 so that readers may verify the validation independently, but performance of the active cooling modes nonetheless constitutes extrapolation beyond the validated envelope and no independent experimental corroboration of those modes yet exists; and (12) the 12 g/cell PCM loading adds 1.25 kg to a 4.78 kg cell mass—an approximately 26% gravimetric penalty, excluding encapsulation—which is not carried through into a vehicle-level energy or range accounting, and PCM containment, cycle-life stability and leakage were not modelled. Future work should address experimental validation of the PCM-enhanced hybrid mode, multi-objective optimization of PID gains, hot and cold climate analysis, and CFD simulation of internal airflow patterns. The spatial temperature distribution at seven thermocouple locations under 30A discharge for all five modes is presented in Supplementary Figure S2, providing the baseline data for future CFD validation.

4 Conclusion

This study developed and validated a lumped-parameter thermal model of a 13S8P lithium-ion pack (104 Panasonic NCR18650BD cells) and used it to evaluate five cooling modes, an adaptive PID control strategy, and four design parameters. The model reproduced published experimental data with R 2 = 0.9515 and MAPE = 6.2%.

• PCM-enhanced hybrid cooling gave the lowest peak temperature at every discharge rate, reducing 30A peaks by 49.8% (48.9 °C versus 97.5 °C) and temperature non-uniformity by 58.6% (8.87 K versus 21.45 K), because isothermal melting engages selectively at the hottest interior cells.
• Adaptive PID control cut cooling energy by 29.75% ± 0.09% relative to a threshold-triggered baseline while retaining a 10.8 °C margin below the 60 °C limit, by escalating actuation in graded steps rather than binary switching.
• PCM mass was the dominant design parameter (−14.1% sensitivity index), roughly three times more influential than internal resistance and an order of magnitude more than airflow speed or TEC current.
• Ageing to 2000 cycles (2 × R base) raised peak temperature from 42.4 °C to 69.0 °C at 25A, breaching the safety limit near 1500 cycles and requiring age-aware controller recalibration.
• Pack topology mattered through n p alone: 10S10P reached 33.3 °C ± 0.04 °C against 38.1 °C ± 0.01 °C for 13S8P, consistent with Q pack ∝ n s / n p.

Collectively these results indicate that latent-heat buffering and graded control address different limitations—peak temperature and uniformity in the first case, energy cost in the second—and are most effective when combined. Priorities for future work are experimental validation of the PCM-enhanced mode, formal PID gain optimisation, extension to −10 °C to 45 °C ambients, and CFD resolution of intra-pack airflow.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmech.2026.1886296/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Anatoliy Levin, Melentiev Energy Systems Institute (RAS), Russia

Reviewed by: Nur Irmawati Om, Universiti Tenaga Nasional, Malaysia

Deepak Kumar Singh, Dr DY Patil Institute of Technology Department of Mechanical Engineering, India