Assignment 09

Heat Conduction Simulation

Simulate transient heat conduction in a 2D rectangular domain until a steady state is reached. Visualize the final steady-state temperature distribution. Follow these steps:

1. Grid Initialization

Create a 2D NumPy array representing the temperature grid. Initially, set all temperatures to zero. (See Parameter section below for specifics on the grid.) Also create a 3D array that will store the 2D array at each time step.

2. Boundary Conditions

Bottom row: Constant temperature of 100 degrees.
Left and right boundary columns: Constant temperature of 50 degrees.
Top boundary row: Insulated (no heat loss), implemented by setting these cells equal to their neighboring cells. At every time step, the top row equals the second top row.

3. Simulation Loop

Iterate over time steps
At each time step, the new temperature matrix computes from the previous matrix using the discretized heat equation (see Background section below). You will have to iterate through all matrix cells in order to implement it.
Once you computed the new grid, you will have to update the boundary conditions of the top boundary cells.
At each time step, check for convergence, and break the loop if the rule is fulfilled (see Convergence Check below).
Store the grid state at each time step in a 3D NumPy array.

4. Convergence Check

Stop the simulation when the steady state is reached. The steady state is reached when the maximum temperature change between two consecutive time steps is below a certain threshold (see Parameters section below).

5. Visualization

Use Matplotlib to visualize the final temperature distribution. At which time step is steady state reached?

Parameters to Use

Matrix dimensions: 10x10
Maximum number of iterations: 1000
Convergence threshold: 0.01
Time step size (\(\Delta t\)): 0.1

Background

Heat conduction equation

The heat equation in two dimensions is given by:

\[ \frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) \]

where \(T\) is temperature, \(t\) is time, and \(\alpha\) is the thermal diffusivity of the material.

Numerical discretization

Using the finite difference method, the heat equation can be approximated as:

\[ T_{i,j}^{new} = T_{i,j} + \Delta t \cdot \left( T_{i+1,j} + T_{i-1,j} + T_{i,j+1} + T_{i,j-1} - 4T_{i,j} \right) \]

In [1]:

import numpy as np
import matplotlib.pyplot as plt

# Parameters
height, width = 10, 10  # Dimensions of the matrix
max_iter = 1000         # Maximum number of iterations
tolerance = 0.01        # Convergence tolerance
dt = 0.1                # Time step size

# Initialize the temperature matrix and the matrix for temperature history
temperature = np.zeros((height, width))
temp_history = np.zeros((height, width, max_iter))

# Boundary conditions
temperature[0, :] = 100  # Bottom edge
temperature[:, 0] = 50    # Left edge
temperature[:, -1] = 50   # Right edge
# Top edge is handled within the loop


# Heat conduction loop
for t in range(max_iter):
    temp_old = temperature.copy()

    # Update temperature (excluding boundaries)
    for i in range(1, height - 1):
        for j in range(1, width - 1):
            temperature[i, j] = temp_old[i, j] + dt * (
                temp_old[i+1, j] + temp_old[i-1, j] +
                temp_old[i, j+1] + temp_old[i, j-1] - 4 * temp_old[i, j]
            )

    # Neumann boundary condition for the top edge (i.e., insulated: temperature gradient is zero)
    temperature[-1, :] = temperature[-2, :]

    # Check for convergence
    if np.max(np.abs(temperature - temp_old)) < tolerance:
        print(f"Steady state reached at iteration {t}")
        break

    # Store the result
    temp_history[:, :, t] = temperature

# Truncate the history array to the number of iterations actually performed
temp_history = temp_history[:, :, :t+1]

# Visualization of the final steady state
fig, ax = plt.subplots()
im = ax.pcolormesh(temperature, cmap='hot')
fig.colorbar(im, label='Temperature')
ax.set_title('Steady State Temperature Distribution')
fig.savefig("heat_final.png")

Steady state reached at iteration 321

09_tasks_solved.ipynb