Heat Diffusion Simulation

A physical 3D printed model of heat diffusion and 2D simulation with Numpy. Built as a foundational learning project.

Progress

• Updating simulations to be able to have options/sliders etc.
• Creating 3D model/physical testing

---

Goals

• Understand how continuous physics gets discretized into a grid a computer can simulate
• Internalize the NumPy vectorized slicing pattern - transfers to image processing, signal processing, IK solvers
• Build intuition for stability, convergence, and adaptive control

---

Files

File	What it is
heat_diffusion.py	Normalized 2D simulation — unitless, focused on behavior and pattern
heat_diffusion_copper.py	Real units 2D simulation — actual copper properties, real seconds, thermistor validation targets

---

The Physics

Heat diffusion is governed by:

$$\frac{\partial T}{\partial t} = \alpha \cdot \nabla^2 T$$

In plain terms: the rate at which a point's temperature changes is proportional to how different it is from its neighbors.

Variables

Symbol	Name	Meaning
$\partial T / \partial t$	Time derivative	How fast temperature is changing at this point right now
$\alpha$	Thermal diffusivity	How fast the material spreads heat (copper: $1.17 \times 10^{-4}$ m²/s)
$\nabla^2 T$	Laplacian	How different this point is from its surroundings

The Laplacian

$$\nabla^2 T = \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}$$

• $\partial^2 T / \partial x^2$ — curvature in the x direction (left/right neighbors)
• $\partial^2 T / \partial y^2$ — curvature in the y direction (top/bottom neighbors)

Three cases:

• $\nabla^2 T < 0$ — hotter than neighbors, will cool down
• $\nabla^2 T > 0$ — cooler than neighbors, will warm up
• $\nabla^2 T = 0$ — in equilibrium, steady state reached

Discretization

The continuous equation becomes a finite difference update rule on a grid. For each interior cell at position (i, j):

T_new[i,j] = T[i,j] + α · dt · (T[i+1,j] + T[i-1,j] + T[i,j+1] + T[i,j-1] - 4·T[i,j])

In NumPy this is written without any loops using array slicing — shifting the entire grid by one row or column and adding:

laplacian = (
    T[2:,   1:-1] +   # neighbor below  → ∂²T/∂y²
    T[:-2,  1:-1] +   # neighbor above  → ∂²T/∂y²
    T[1:-1, 2:]   +   # neighbor right  → ∂²T/∂x²
    T[1:-1, :-2]  -   # neighbor left   → ∂²T/∂x²
    4 * T[1:-1, 1:-1]
)

This is the most important pattern in the codebase. The math and the code are the same thing written two different ways.

CFL Stability Condition

The time step dt cannot be chosen freely — it must satisfy:

$$dt \leq \frac{dx^2}{4\alpha} \quad \text{(2D)}$$

If dt is too large the simulation diverges deterministically — small errors amplify each step until temperatures explode. This is not random, it is guaranteed by the math.

Normalized vs Real Units

	Normalized	Real (Copper)
Purpose	Understand behavior and shape	Validate against physical measurements
ALPHA	0.25 (chosen for stability)	1.17e-4 m²/s (real copper)
dt	Abstract steps	Real seconds
BOUNDARY_TEMP	0.0	22.0°C (room temperature)
Axes	Grid cells	Centimeters

The normalized version is a map without a scale bar — correct shape, no real units. The real version adds the scale bar. The simulation structure is identical between both.

---

Connection to Eragon & Future Projects

This project is a direct foundation for upcoming systems:

Eragon Control Loop

The Medulla layer on Eragon runs inverse kinematics — matrix math transforming target foot positions into joint angles. The NumPy vectorized operations practiced here (slicing, linear algebra, iterative solvers) are exactly what IK solvers use. The adaptive controller maps directly onto servo control: a PID loop managing joint position has the same stability threshold, the same overshoot risk, and the same need for asymmetric correction as dt management here.

The Medulla/Reflex split mirrors the architecture of this simulation:

• Reflex (ESP32) — executes the update rule each step, low latency
• Medulla (Pi) — runs the controller, monitors stability, adjusts parameters