Overview

#

A transient thermal simulation tool for predicting battery pack temperatures under various heat loads and cooling configurations. The simulator models forced-air convection through a tube bank (cylindrical battery cells) with fan control, enabling design optimization of the cooling system for a solar racing vehicle.

Key Simulation Assumptions:

#

1. Incompressible flow
   • Proven valid by low pressure drops across system an low Mach numbers (<0.3)
2. Viscous forces are non-negligible
   • Adiabatic flow is assumed in pipe sections, with pressure losses due to friction only (Darcy-Weisbach) and pipe bends (K factors)
3. Steady-state flow at each time step
4. Constant Mass Flow Rate at each time step
   • Air can have two states: stagnant and flowing. When flowing, mass flow rate is constant throughout the system.
   • This permits simpler series component modeling and permits simple usage of a fan in the middle

Key Features

#

Language: Python 3.12

System Architecture

#

┌─────────┐   ┌─────────┐   ┌─────────────┐   ┌─────┐   ┌─────────┐   ┌─────────┐
│  Inlet  │ → │  Pipe1  │ → │ ThermalLoad │ → │ Fan │ → │  Pipe2  │ → │ Exhaust │
└─────────┘   └─────────┘   │  (Battery)  │   └─────┘   └─────────┘   └─────────┘
                            └─────────────┘

Each component implements a process() method with signature:

@abstractmethod
def process(self, air_flow, air_temp, velocity, pressure, density, **kwargs):
    pressure_out = pressure
    return air_flow, air_temp, velocity, pressure_out, density

Outputs are fed in series through the defined cooling system. This allows for compartmentalization and modular testing of each component, and together as a system. Code can be written for specific parts, and improved independently. Complete systems can be swapped out or reconfigured easily in code.

Physics & Equations

#

Governing Energy Balance

#

The battery temperature evolves according to:

Where:

• \(\dot{Q}_{gen}\) = Internal heat generation (I²R losses)
• \(\dot{Q}_{conv}\) = Convective heat transfer to air
• \(\dot{Q}_{rad}\) = Radiative heat transfer
• \(C_{batt}(T)\) = Temperature-dependent battery thermal capacitance

dT_dt = (self.heat_load_W - self.Q_to_air_W) / self.batt_thermal_cap_JK_lam_T(self.T_batt_K)
self.T_batt_K += dT_dt * dt_s

Convective Heat Transfer

#

Zukauskas Correlation for tube banks (cylindrical cells):

Heat Transfer Coefficient from Nusselt number:

From the definition of Nusselt number:

From the heat transfer, we can compute the overall conductance:

we assume \(U = h\), the overall heat transfer is equal to the convective heat transfer coefficient for the air side only

NTU-Effectiveness Method for flowing air:

Where:

• \(C\) = Constant defined by Zukauskas
• \(h\) = Convective heat transfer coefficient [W/(m²·K)]
• \(k\) = Thermal conductivity of air [W/(m·K)]
• \(D\) = Characteristic length (tube diameter) [m]
• \(A = N \cdot \pi D \cdot L\) = Total surface area of all tubes [m²]
• \(L\) = Tube length [m]
• \(N\) = Number of tubes in the bank

from ht import Nu_ESDU_73031 
def _Nu_ESDU(self, T_air_K, p_Pa, velocity, rho_air):
    n_rows = int(getattr(self.zakauskas_data, "n_rows"))
    st = getattr(self.zakauskas_data, "st")
    sl = getattr(self.zakauskas_data, "sl")

    prop = self.air_props(T_air_K, p_Pa=p_Pa)
    Pr, mu = prop["Pr"], prop["mu"]

    prop_surface = self.air_props(self.T_batt_K, p_Pa=p_Pa)
    Pr_s = prop_surface["Pr"]

    Re = self._tube_bank_reynolds(velocity, rho_air, mu)

    nu = Nu_ESDU_73031(
        Re=Re,
        Pr=Pr,
        tube_rows=n_rows,
        pitch_parallel=sl,
        pitch_normal=st,
        Pr_wall=Pr_s,
    )

    return nu

Log-Mean Temperature Difference

#

For heat exchanger calculations:

Pressure Drop

#

Darcy-Weisbach for tubes

#

With friction factor:

• Laminar \((Re < 2300): f = 64/Re\)
• Turbulent: \(f = 0.3164 \cdot Re^{-0.25}\) (Blasius)

Zukauskas Pressure Drop for tube banks

#

Where \(f_{bank}\) is obtained from Zukauskas correlations.

from ht import dP_Zukauskas
def _pressure_drop_ht(self, air_flow, air_temp, velocity, pressure, density):
    if (
        not air_flow
        or air_flow <= self.FLOW_THRESHOLD
        or self.pa_method != "Zukauskas"
    ):
        return 0.0

    n_rows = int(getattr(self.zakauskas_data, "n_rows"))
    st = getattr(self.zakauskas_data, "st")
    sl = getattr(self.zakauskas_data, "sl")
    D = self.cell_diameter_m

    prop = self.air_props(air_temp, p_Pa=pressure)
    mu = prop["mu"]

    V_max = self._tube_bank_max_velocity(velocity)
    Re = self._tube_bank_reynolds(velocity, density, mu)

    try:
        dp = dP_Zukauskas(
            Re=Re, n=n_rows, ST=st, SL=sl, D=D, rho=density, Vmax=V_max
        )
    except Exception as _e:
        logger.warning(f"{self.name}: dP_Zukauskas raised: {_e}")
        dp = 0.0

    return float(dp)

Note that it is assumed that there is no pressure drop when air is stagnant

Fan Control System

#

PWM-Based Pressure Curve

#

The fan pressure-flow relationship is modeled symbolically, as an approximation of manufacturer data:

Where \(Q \left[\frac{m³}{min}\right] \) is volumetric flow and \(p \in [0, 1]\) is the duty cycle.

The slope and the intercept coefficients \((m_i, b_i)\) are determined by curve-fitting manufacturer fan curve data at 0%, 50%, and 100% duty cycles.

Temperature-Based Duty Control

#

A Fan duty cycle control law is defined piecewise based on battery temperature \(t_K\) to optimize cooling. At lower temperatures, the fan is off or at minimum speed to save power. As temperature rises, the duty cycle ramps up linearly to maximum at a defined threshold.

from numpy import interp

T_LOW_K = 273.15 + 20  # Minimum temperature for fan operation
T_MAX_K = 273.15 + 45  # Maximum allowable temperature
T_MID_K = (T_LOW_K + T_MAX_K) / 2 + 5  # Slope change point
DUTY_MIN_REGS = 0.01  # Minimum duty cycle for regulations
DUTY_MIN_OPER = 0.20  # Minimum duty cycle for actualoperation
DUTY_MAX = 1.00

duty_curve_points = [
    (T_LOW_K, DUTY_MIN_REGS),
    (T_LOW_K, DUTY_MIN_OPER),  # Jump from regs minimum to operational minimum
    (T_MID_K, 0.5),  # slope change
    (T_MAX_K, DUTY_MAX),  # Maximum duty point
]

"""Functions merged for brevity."""

sorted_points = sorted(duty_curve_points)
T_points = [p[0] for p in sorted_points]
duty_points = [p[1] for p in sorted_points]

self._duty_interp = lambda x_1: interp(
    x_1, T_points, duty_points, left=duty_points[0], right=duty_points[-1]
)

Operating Point Determination

#

The system finds equilibrium where fan curve intersects system impedance:

Battery Heat Generation via Pre-computed Lookup Tables

#

The Challenge

#

Direct PyBaMM simulation of battery heat generation takes ~700ms per operating point. For an 8-hour solar car race with 1-second timesteps (28,800 points), this means 5+ hours just calculating heat—making design optimization infeasible. Using a pregenerated lookup table can reduce this to seconds added over an entire race simulation.

The Solution

#

Pre-compute heat generation across the operating space once, then query via fast interpolation during simulation:

1. Offline Phase: Generate a 3D lookup table across 1.4M operating points using PyBaMM (fine resolution):

   • SoC: 0.5–99.5% (80 points, non-uniform—denser at extremes where resistance peaks)
   • Temperature: 10–50°C (40 points, denser in thermal stress zone 35–50°C)
   • Power: −2.5 to +8 kW (450 points, extremely dense above 3 kW where I²R dominates)

2. Runtime Phase: Query via 3D linear interpolation (~90 microseconds)

Three heat components are stored at each point:

• Ohmic heat (I²R losses): Dominant at high discharge
• Entropic heat (reversible): Entropy-driven heating/cooling due to electrochemical reactions
• Total heat: Sum fed to the cooling system

Performance & Accuracy

#

Calculted over 43k random test points

Average error of 0.5% and max deviation of 0.01 W is negligible for thermal design, enabling rapid iteration without sacrificing meaningful accuracy.

Code Integration

#

from Objects.BatteryModel import BatteryModel

# Initialize once
battery = BatteryModel(series_cells=36, parallel_strings=8, initial_soc=0.85)

# Query each timestep (~80 microseconds)
heat = battery.calculate_heat_generation(power_W=3500, dt_s=1.0, temperature_K=310)
print(f"Heat: {heat['heat_total_W']:.0f} W (ohmic: {heat['heat_ohmic_W']:.0f} W)")

Grid Design Strategy

#

Non-uniform spacing targets thermally critical regions:

• SoC extremes: Resistance varies sharply near 0% and 100%
• High discharge (>3 kW): Every ~13 watts resolved (I²R quadratic sensitivity)
• Thermal stress (35–50°C): Every degree resolved (critical for fan control)

Four resolution options balance accuracy vs. generation time. Coarse and medium resolutions can be run locally for quick tests, while fine and ultra require HPC resources provided by the Purdue Aero/Astro Department.

• Coarse: 14K points, ~3 hours
• Medium: 360K points, ~70 hours
• Fine: 1.4M points, ~285 hours
• Ultra: 34M points, ~7000 hours

Times are given for single-core runs; parallel computation(multi-threading) reduces generation time significantly.

Integration with Thermal Simulation

#

The lookup table feeds directly into the cooling system:

Driving Cycle → Power Management
    ↓
Battery Model ← Lookup Table
    ↓
Heat Generation (W) = function(Power, SOC, Temperature)
    ↓
Cooling Simulation (fan control, temperature prediction)
    ↓
Results: Battery temperature, state of charge, accurate power losses

This enables real-time design optimization—sensitivity analysis on cooling configurations now takes seconds instead of hours, transforming thermal management from offline analysis into an interactive design tool.

Interactive Heat Generation Surface

#

The 3D surface below visualizes heat generation as a function of battery power, temperature, and state of charge at a medium resolution (360K points).

Battery Pack Parameters

#

Key Results

#

The simulation produces:

1. Battery temperature vs. time for various heat loads
2. Fan operating point (duty cycle, flow rate) vs. time
3. System impedance curves with fan curve overlays
4. Steady-state temperatures at different power levels

Code Quality Highlights

#

• Caching: LRU cache on CoolProp calls with temperature/pressure quantization to avoid redundant lookups, with over 99% of CoolProp calls hitting the cache during typical simulations
• Pregenerated lookup tables: Battery heat generation via 3D interpolation for rapid querying
• Validation warnings: Biot number checks for lumped capacitance and incompressible flow assumption validity
• Modular design: Easy to swap components or add new correlation methods

from functools import lru_cache
from thermo.mixture import Mixture

@staticmethod
def air_props(T_K, p_Pa=None, h=0):
    # Quantized + cached with lru_cache
    T_STEP = 0.1  # K
    P_STEP = 5.0  # Pa
    if p_Pa is None:
        p_Pa = CoolingComponent.P0 * math.exp(
            -CoolingComponent.g * h / (CoolingComponent.R * T_K)
        )
    # quantize inputs
    T_bin = int(T_K // T_STEP * T_STEP)  # floor to lower step
    P_bin = int(p_Pa // P_STEP * P_STEP)  # floor to lower step
    
    result = CoolingComponent._air_props_cached(T_bin, P_bin)
    
    return result

@staticmethod
@lru_cache(maxsize=4096)
def _air_props_cached(T_K_q, P_Pa_q):
    # Cached version of air_props using thermo
    air = Mixture('air', T=T_K_q, P=P_Pa_q)
    
    return {
        "rho": air.rho,           # density, kg/m³
        "mu": air.mu,             # dynamic viscosity, Pa·s
        "k": air.k,               # thermal conductivity, W/m·K
        "Cp": air.Cp,             # specific heat, J/kg·K
        "Pr": air.Pr,        # Prandtl number, dimensionless
        "p": P_Pa_q,              # pressure, Pa
        "gamma": air.Cp / air.Cvg # heat capacity ratio (using gas Cv)
    }

Data quantization and caching reduces expensive Thermo calls significantly during time-stepping

2206 ms - WARNING - ThermalLoad: Mild Internal Gradients! Bi = 0.203
2572 ms - INFO - ThermalLoad: Biot Number Summary over 5760 steps:
2572 ms - INFO -   Max Bi = 0.370
2572 ms - INFO -   Mild (0.1-0.5):         3145 steps (54.6%)
2580 ms - INFO -   Significant (0.5-1.0):     0 steps (0.0%)
2580 ms - INFO -   Severe (>=1.0):            0 steps (0.0%)
2580 ms - INFO - Completed run for 15W heat load!
2581 ms - INFO - Fan energy used: 98.94 Wh over 8.00 hours
2581 ms - INFO - Running simulation for 20W heat load...
2686 ms - WARNING - ThermalLoad: Mild Internal Gradients! Bi = 0.203
3289 ms - INFO - ThermalLoad: Biot Number Summary over 5760 steps:
3289 ms - INFO -   Max Bi = 0.371
3289 ms - INFO -   Mild (0.1-0.5):         3798 steps (65.9%)
3289 ms - INFO -   Significant (0.5-1.0):     0 steps (0.0%)
3289 ms - INFO -   Severe (>=1.0):            0 steps (0.0%)
3289 ms - INFO - Completed run for 20W heat load!

Biot number warnings confirm the validity of the lumped capacity assumption of battery thermal model

Future Work

#

Completed Improvements

#

• Transition off SymPy symbolic math for performance

• Now using SciPy and Numpy for numerical root finding and interpolation

• Migrate to a more robust Heat Exchanger library like HT, implementing ESDU 73031

• ESDU 73031 implemented using ht library

• Refine Battery thermal model with more accurate geometry and conduction paths

• Now using high fidelity thermal model from PyBaMM module

Planned Enhancements

#

• Integrate with vehicle-level simulation for co-simulation of electrical and thermal systems
• Create more robust web/local interface for design exploration

Results

#

• Numerically quantified choices of Fans and Cooling Configurations for optimal battery temperatures
• Created integral part of the vehicle simulation toolchain for Artemis 2026 and future solar cars

Key Skills

#

• Programming: Python, OOP, Numerical Methods, Data Analysis
• Fluids: Thermal Modeling, Heat Transfer, Fluid Dynamics
• Simulation: Vehicle-level co-simulation, Battery thermal modeling

Computational Tools & Libraries

#

This project leverages several open-source libraries for thermal modeling, data analysis, and visualization:

Scientific Computing:

• NumPy 2.3.5 – Numerical array operations and mathematical functions
• SciPy 1.16.3 – Scientific computing and optimization
• Pandas 2.3.3 – Data manipulation and analysis

Thermodynamic & Heat Transfer:

• Thermo 0.6.0 – Thermodynamic properties and calculations
• PyBaMM 25.10 – Battery and electrochemical modeling
• ht 1.2.0 – Heat transfer correlations and calculations

Visualization & Output:

• Matplotlib 3.10.7 – Data plotting and visualization
• Pillow 12.0.0 – Image processing and handling

Utilities:

• tqdm 4.67.1 – Progress bar visualization
• Colorama 0.4.6 – Terminal color output formatting