Modelling Two Dimensional Heat Conduction Problem using Python

In this tutorial, we will see how to model the 2D heat conduction equation using Python. A 2D, steady-state heat conduction equation with heat generation can be written in Cartesian coordinates as follows

$$\mathrm{\nabla^{2} T \: + \: \frac{q_{g}}{k} \: = \: \frac{\partial^{2}T}{\partial x^{2}} \: + \: \frac{\partial^{2}T}{\partial y^{2}} \: + \: \frac{q_{g}}{k} \: = \: 0 \:\:\dotso\dotso (1)}$$

This equation must be discretized to obtain a finite difference equation. Let us consider a rectangular grid as shown below.

The index i runs vertically (row-wise) whereas the index j runs horizontally (column-wise). Any interior node (i,j) is surrounded by four neighboring nodes, therefore information at this node can only be transferred from these four nodes. For boundary nodes, the information comes from the boundary conditions.

Finite Difference Discretization

To discretize Eq. 1, we use the central difference scheme which is 2nd order accurate. The spacing between grid points in both directions is kept the same: $\mathrm{\Delta x \: = \: \Delta y \: = \: \Delta}$. The central difference formula for the 2nd order derivative is

$$\mathrm{\frac{d^{2}f(x)}{dx^{2}} \: = \: \frac{f_{i-1} \: - \: 2f_{i} \: + \: f_{i+1}}{\Delta x^{2}}\:\:\dotso\dotso (2)}$$

Applying this to Eq. 1, the discretized form becomes

$$\mathrm{\frac{T_{i-1, \: j} \: - \: 2T_{i,\: j} \: + \: T_{i+1, \: j}}{\Delta^{2}} \: + \: \frac{T_{i, \: j-1} \: - \: 2T_{i,\: j} \: + \: T_{i, \: j+1}}{\Delta^{2}} \: + \: \frac{q_{g}}{k} \: = \: 0 \:\:\dotso\dotso (3)}$$

Simplifying Eq. 3, we get the iterative formula

$$\mathrm{T_{i,\: j} \: = \: \frac{1}{4}(T_{i-1, \: j}\: + \: T_{i+1, \: j} \: + \: T_{i,\: j-1} \: + \: T_{i, \: j+1} \: + \: \Delta^{2} \frac{q_{g}}{k}) \:\:\dotso\dotso (4)}$$

Boundary Conditions

Three types of boundary conditions are commonly used

• Dirichlet: Temperature value is specified at the boundary

  $$\mathrm{x \: = \: 0 \: \rightarrow \: T \: = \: 500 \: K }$$

• Neumann: Temperature gradient is specified at the boundary

  $$\mathrm{x \: = \: 0 \: \rightarrow \: \frac{\partial T}{\partial x} \: = \: 0 }$$

• Robin: Convective heat transfer boundary condition

  $$\mathrm{- \: k \: (\frac{\partial T}{\partial x})_{x=0} \: = \: h(T \: - \: T_{\infty}) }$$

Algorithm

1. Define the computational domain
2. Set number of grid points in x and y directions
3. Calculate grid spacing $\mathrm{\Delta x}$ and $\mathrm{\Delta y}$
4. Generate the grid
5. Provide initial temperature guess
6. Apply boundary conditions
7. Solve Eq. 4 iteratively for internal nodes
8. Check convergence criteria
9. Repeat until convergence is achieved

Python Implementation

We will solve two cases: with and without heat generation using Dirichlet boundary conditions ?

# Importing modules
from pylab import *
from numpy import *

# Defining thermal properties
# Case 1: qg = 0 (no heat generation)
# Case 2: qg = 1000000 (with heat generation)
qg = 1000000  # W/m³
k = 100       # W/m·K

# Define domain
Lx = 1.0  # Length in x-direction (m)
Ly = 1.0  # Length in y-direction (m)

# Number of grid points
nx = 20
ny = 20

# Grid spacing
dx = Lx / (nx - 1)
dy = Ly / (ny - 1)
delta = dx  # Assuming dx = dy

# Grid generation
x = linspace(0, Lx, nx)
y = linspace(0, Ly, ny)
X, Y = meshgrid(x, flipud(y))

# Initial temperature guess
T = zeros(shape=(nx, ny))

# Boundary conditions (Dirichlet)
T[:, 0] = 300   # Left boundary: 300 K
T[:, -1] = 300  # Right boundary: 300 K
T[0, :] = 800   # Top boundary: 800 K
T[-1, :] = 300  # Bottom boundary: 300 K

# Copy for comparison
Tg = T.copy()

# Convergence parameters
error = 1.0
tolerance = 1e-5
count = 1

# Iterative solution
while error > tolerance:
    # Sweep through internal nodes
    for i in range(1, nx - 1):
        for j in range(1, ny - 1):
            T[i, j] = (T[i, j-1] + T[i, j+1] + T[i-1, j] + T[i+1, j] + 
                      delta**2 * qg / k) / 4
    
    # Calculate error
    error = sqrt(sum(abs(T - Tg)**2))
    
    # Print progress every 100 iterations
    if count % 100 == 0:
        print(f"Iteration {count}: Error = {error:.2e}")
    
    # Update guess for next iteration
    Tg = T.copy()
    count += 1

print(f"Converged in {count-1} iterations")

# Plot temperature contours
figure(figsize=(8, 6))
contourf(X, Y, T, levels=20, cmap='jet')
colorbar(label='Temperature (K)')
contour(X, Y, T, levels=10, colors='black', linewidths=0.5)
title('2D Temperature Distribution')
xlabel('x (m)')
ylabel('y (m)')
show()

Results and Analysis

The solution shows distinct behavior for the two cases ?

Case 1: No Heat Generation (qg = 0)

Temperature distribution is dominated by boundary conditions, with heat flow from the hot top boundary (800 K) toward the cooler boundaries (300 K).

Case 2: With Heat Generation (qg = 1,000,000 W/m³)

Internal heat generation creates high core temperatures that exceed boundary values. Heat flows outward from the center while still satisfying boundary conditions.

Key Observations

• Heat generation significantly increases core temperatures
• Boundary conditions are always satisfied regardless of internal heat generation
• The finite difference method converges efficiently for both cases
• Higher heat generation changes the dominant heat flow direction

Conclusion

This tutorial demonstrated solving the 2D heat conduction equation using Python's finite difference method. The comparison between cases with and without heat generation shows how internal heat sources dramatically affect temperature distribution while boundary conditions remain enforced.