Simulation 101: Conductive Heat Transfer

A gentle introduction to computational physics

Conduction, or heat transfer between objects, is something we experience everyday. Putting a pan on the stove or sitting on a hot park bench gives us an intuitive sense of conductive heat transfer but here we will formalize the process and build a basic computational framework to simulate it. Conduction is an excellent first simulation problem to tackle because it uses the basic tools found in many computational physics problems.

In this article we will:

Create a mesh grid to represent materials
Learn basic heat transfer equations and their computational equivalents
Update the values in our mesh grid based on the underlying Physics
Simulate conductive heat transfer

Creating a Mesh Grid

A mesh grid is a computational tool used to discretize a continuous space. That is, we can’t perform calculations on all time and space in our problem, so we chose a representative subset of points, usually at regular intervals, to look at instead.

In figure 1 below we can see an example of a mesh grid. Here a space is subdivided into evenly spaced cells which is common practice is physics Simulation. Instead of running calculations/simulation on the entire surface we can now work with only our grid points which makes our problem more feasible.

Figure 1: Example of a mesh grid. In a simulation we break up our space into such a grid and compute values at every dotted grid point.

The mesh grid above was created using Python’s numpy meshgrid function which can take in a set of 1 dimensional arrays and create a mesh grid for us. For our simulation, we want to model a 2 dimensional surface, so we are going to generate 2 arrays filled with the starting values we want with a length of how many intervals we want to evaluate our simulation on. See the code snippet below where we create a 100×100 mesh grid of zeros as the basis of our simulation.

import numpy as np

#Define how many intervals we want per axis
resolution = 100

#Create x and Y arrays of zeros of length 100
x = np.zeros(resolution)
y = np.zeros(resolution)

#Create a mesh grid from above arrays. Outputs are all x and Y values in grid
gridX, gridY = np.meshgrid(x,y)

Before moving on from mesh grids, it’s important to note that every point on our grid is typically not a single value. Usually at every grid point we want an array of values that represent the properties we will be working with. In our heat transfer simulation we will need to know material properties as well as the temperature at every point on our mesh grid. So all the points on our grid will be an array that contains information like temperature, material density, material conductivity and anything else we need to know.

An easy way to make a mesh grid that is an array of arrays is to create multiple mesh grids that contain a single value at each point and then stack then together as seen in figure 2.

Figure 2: Stacking 2 mesh grids to create a new mesh grid that contain multiple values per point.

Numpy comes in handy here again with the dstack function which will stack two arrays element wise for us. Below is a code snippet that will create 2 mesh grids and stack them.

import numpy as np

#Define how many intervals we want per axis
resolution = 100

#Create x and Y arrays of zeros of length 100
x1 = np.zeros(resolution)
y1 = np.zeros(resolution)

x2 = np.full(resolution, 1)
y2 = np.full(resolution, 1)

#Create a mesh grids from above arrays. Outputs are all x and Y values in grid
gridX, gridY = np.meshgrid(x1,y1)
gridX2, gridY2 = np.meshgrid(x2,y2)

#Stack both of the mesh grids we have created so every element is [0,1]
fullGridX = np.dstack([gridX,gridX2])
fullGridY = np.dstack([gridY,gridY2])

With the tools to create our own mesh grid to represent our simulation environment, we can move on to the physics.

Basics of Heat Transfer

We will first start with the 1 dimensional time dependent heat conduction equation.

Equation 1: 1-dimensional time dependent heat conduction equation.

The equation states the the change in temperature over time is based on the material’s thermal properties scaled by the change in temperature across the material (technically the change of the change of the temperature since it is a 2nd derivative). We don’t need to worry about the thermal properties of the material since they can be found in look up tables. What we need to do is convert the derivatives in the equation into something we can handle computationally. Luckily for us, this can be done with the finite difference method which will convert equation 1 into equation 2 below.

Equation 2: Computational 1-dimensional time dependent heat conduction equation.

Equation 2 is a computational equivalent of the time dependent heat conduction equation that we can use to update the temperature in our simulation. Going through each cell in our mesh grid, we can update the temperature of that cell by its current temperature plus the difference in temperature between that cell and its neighboring cells (here i represents the mesh grid index) scaled by the thermal properties of the material and a time step of our choosing.

Equation 1 & 2 are the 1-dimenstional heat equation but we want to run our simulation in 2-dimensions. Adding in a 2nd dimension is straight forward, especially having the computational version of the heat equation, we just need to add in all neighboring cells as seen below in equation 3.

Equation 3: Computational 2-dimenstional time dependent heat conduction equation.

Now that we have translated the time dependent heat conduction equation into a computational expression, we can use it to update our mesh grid.

Updating Our Mesh Grid

Let’s first initialize a mesh grid using the tools we have developed so far. We know we need at least 2 features to plug into our temperature equation which is a current temperature and a thermal diffusivity constant for a material of our choosing. We will start with our material being at room temperature (20°C), and our material will be copper which has a thermal diffusivity of 1.11×10^-2 cm²/s. We fill 2 mesh grids with these values and stack them to make our overall mesh.

import numpy as np

#Define how many intervals we want per axis
resolution = 100
startingTemperature = 20
thermalDiffusivity = 1.11*10**-2

#Create x and Y arrays of zeros of length 100
x1 = np.full(resolution, 20)
y1 = np.full(resolution, 20)

x2 = np.full(resolution, 1.11*10**-2)
y2 = np.full(resolution, 1.11*10**-2)

#Create a mesh grids from above arrays. Outputs are all x and Y values in grid
gridX, gridY = np.meshgrid(x1,y1)
gridX2, gridY2 = np.meshgrid(x2,y2)

#Stack both of the mesh grids we have created so every element is
#[20,1.11*10**-2]
fullGridX = np.dstack([gridX,gridX2])
fullGridY = np.dstack([gridY,gridY2])

In our specific case, our mesh grid is symmetric (the feature values are the same in the x and y direction), so defining our x and y grid separately is not needed but it is still good practice because we will not always have such symmetry when we simulate. Now that we have our material represented by a mesh grid we need to add heat to it (as we saw in our physics equations, there has to be a difference in temperature for heat to transfer). To add heat to our material, we need to select the cells in our mesh grid we want to heat up and add value to their temper index. We will continue to keep things basic and add heat, let’s say 1000°C, in a square from index 35–65 in both the X and Y direction. The code to do this and a resulting plot of this can be seen below.

import matplotlib.pyplot as plt

#fullGridX[:,:,0] gets the first index of every element in our mesh grid
#fullGridX[:,:,0][35:65,35:65] selects a box of values
fullGridX[:,:,0][35:65,35:65] += 1000

plt.imshow(fullGridX[:,:,0])

Figure 3: Visualization of added heat to our mesh grid.

We’ve added a patch of heat to our mesh grid, which is a good start but now we need to evolve our system through time to see how the heat transfers. We will need to develop a few more tools to do this. The first will be a find neighbors’ temperature function, which given a cell can get the temperature off all neighboring cells. This function is given below and loops over all 8 bordering cells of a given cell but ignores the corners to only retrieve cells that are up, down, left and right. Another thing we have to consider when looking at all of a cell’s neighbors is boundary conditions or what we want to do when we hit the end of our grid. If we are at a grid cell that is already at the edge of the grid, looking for all of its neighbors will throw an error since they don’t all exist. We can handle this with a try except statement in our loop which will give our boundary condition when we try to look for a cell that does not exist. Our boundary condition can be filled with the temperature we want to assume outside of our grid. In our simulation, we are going to say this is room temperature, so when we calculate the heat transfer at our mesh grid edges, room temperature will always be given as a neighboring value to an edge cell.

def getNeighborsTemperature(grid, point, boundaryTemp):
    #List to collect all neighboring temperatures
    neighbors = []
    #Loop over all neighboring cells
    for i in range(-1,2):
        for j in range(-1,2):
            try:
                #Ignore corner cells
                if abs(i) != abs(j):
                    neighbors.append(grid[point[0] + i][point[1]+ j])
            except:
                #Apply boundary condition
                neighbors.append(boundaryTemp)

    return neighbors

The next function we have to develop will implement the 2-dimensional time dependent heat conduction equation we found in the last section. We solved for the computational equivalent for this equation so the implementation is straight forward and given below. The function completes the heat conduction equation by taking in the temperature of a given cell, the temperature of the cells around it, a time step of our choosing, and the thermal diffusivity of the cell material.

def calculateHeat(cellTemp, neighborTemps, timeStep, thermalDiffusivity):
    #Converting equation 3 into code
    cellTemp = cellTemp + timeStep*thermalDiffusivity*((neighborTemps[0] -2*cellTemp + neighborTemps[-1]) + 
                                                       (neighborTemps[1] -2*cellTemp + neighborTemps[-2]))
    return cellTemp

The last function we need will combine the two previous functions to loop over every cell in our mesh grid and update its temperature. For every cell in our grid we first run "getNeighbors" to get all of the neighboring temperatures and then we pass the neighboring temperatures plus the current cell temperature and other arguments to update the cell’s temperature.

def heatTransfer(grid, timeStep, boundaryTemp):
    #Loop over all grid cells 
    for i in range(0,len(grid)):
        for j in range(0,len(grid)):
            #Get neighboring cell temperatures
            neighbors = getNeighborsTemperature(grid[:,:,0], (i,j), boundaryTemp)
            #Update current cell temperature
            grid[:,:,0][i][j] = calculateHeat(grid[:,:,0][i][j], neighbors, timeStep, grid[:,:,1][i][j])

    return grid

Finally, we can use our heat transfer equation to update our original mesh grid. The final code snippet will update our grid with a 30 second time step and compare the before and after.

#Applying our heat transfer equation on current grid with a 30 secound time step
nextGrid = heatTransfer(fullGridX.copy(),30, 20)

#Matplotlib subplots to show the original and updated mesh grid
plt.figure(figsize = (12,6))

plt.subplot(1,2,1)
plt.imshow(fullGridX[:,:,0])
plt.title("t = 0")

plt.subplot(1,2,2)
plt.imshow(nextGrid[:,:,0])
plt.title("t = 30")

Figure 4: Updating our mesh grid with our heat transfer heat with a 30 second time step.

We now have all the tools need to run conductive heat transfer simulations, we will run a few in the next section.

Simulation

With all of out simulation tools created, all we have to do is put our heat transfer equation inside of a loop and evolve the system for as many time steps as we want. The code to do this is given as well as a function to visualize our simulation as a gif.

import imageio
"""
This make gif function takes in a set of images and turns them into a gif!
Give the frames as an array of mesh grids (only the temperature)and a name 
for the gif as well as the temperature bounds for the heat map.
"""
def makeGif(frames,name,minTemp,maxTemp):
    !mkdir frames

    counter=0
    images = []
    for i in range(0,len(frames)):
        plt.figure()
        plt.imshow(frames[i], cmap = "inferno", vmin = minTemp, vmax = maxTemp)
        plt.savefig("frames/" + str(counter)+ ".png")
        images.append(imageio.imread("frames/" + str(counter)+ ".png"))
        counter += 1
        plt.close()

    imageio.mimsave(name, images)

    !rm -r frames

#RUNNING SIMULATION
#Keep track of our mesh grid frames to make into gif
frames = [fullGridX[:,:,0].copy()]

timeStep = 30
boundaryTemp = 20

#Run for 500 time steps
for t in range(0,500):
    fullGridX = heatTransfer(fullGridX.copy(),timeStep,boundaryTemp)
    frames.append(fullGridX[:,:,0].copy())

makeGif(frames,"simulation.gif",20,500)

Figure 5: Gif of our conductive heat transfer simulation

And there we have it. We have ran a simulation of conductive heat transfer and visualized it (with the cool "inferno" color map). With the generalized set of tools we developed, we can run many different simulations with different time steps, starting temperatures, materials and geometries. The next and final section will show off a couple more scenarios made with our tools.

Example Scenarios

Welding

The following is a simulation of a sheet of brass and steel being welded together. What is interesting is the brass (left side) has about 3x the thermal conductivity as steel (right side) and we can see the difference in the material properties as the heat dissipates.

Note that we can add heat into out simulation iteratively as well by adding to the temperature layer of our mesh grid inside of our simulating loop.

Figure 6: Simulation of welding brass and steel together.

Bruner

Here a stove top burner is being simulated by adding a circular patch of heat to the mesh grid. The heat is being added gradually (3°C per second) until the burner is turned off and allowed to cool.

Figure 7: Simulation of a stove top burner being heated up and turned off.

Full Code

import numpy as np
import matplotlib.pyplot as plt
import imageio

def getNeighborsTemperature(grid, point, boundaryTemp):
    neighbors = []
    for i in range(-1,2):
        for j in range(-1,2):
            try:
                if abs(i) != abs(j):
                    neighbors.append(grid[point[0] + i][point[1]+ j])
            except:
                neighbors.append(boundaryTemp)

    return neighbors

def calculateHeat(cellTemp, neighborTemps, timeStep, thermalDiffusivity):
    cellTemp = cellTemp + timeStep*thermalDiffusivity*((neighborTemps[0] -2*cellTemp + neighborTemps[-1]) + 
                                                       (neighborTemps[1] -2*cellTemp + neighborTemps[-2]))
    return cellTemp

def heatTransfer(grid, timeStep, boundaryTemp):

    for i in range(0,len(grid)):
        for j in range(0,len(grid)):
            neighbors = getNeighborsTemperature(grid[:,:,0], (i,j), boundaryTemp)
            grid[:,:,0][i][j] = calculateHeat(grid[:,:,0][i][j], neighbors, timeStep, grid[:,:,1][i][j])

    return grid

def makeGif(frames,name,minTemp,maxTemp):
    !mkdir frames

    counter=0
    images = []
    for i in range(0,len(frames)):
        plt.figure()
        plt.imshow(frames[i], cmap = "inferno", vmin = minTemp, vmax = maxTemp)
        plt.savefig("frames/" + str(counter)+ ".png")
        images.append(imageio.imread("frames/" + str(counter)+ ".png"))
        counter += 1
        plt.close()

    imageio.mimsave(name, images)

    !rm -r frames

#Make mesh grid
resolution = 100
startingTemperature = 20
thermalDiffusivity = 1.11*10**-2

x1 = np.full(resolution, 20)
y1 = np.full(resolution, 20)

x2 = np.full(resolution, 1.11*10**-2)
y2 = np.full(resolution, 1.11*10**-2)

gridX, gridY = np.meshgrid(x1,y1)
gridX2, gridY2 = np.meshgrid(x2,y2)

fullGridX = np.dstack([gridX,gridX2])
fullGridY = np.dstack([gridY,gridY2])

#Add heat
fullGridX[:,:,0][35:65,35:65] += 1000

#Run simulation
frames = [fullGridX[:,:,0].copy()]

timeStep = 30
boundaryTemp = 20

for t in range(0,500):
    fullGridX = heatTransfer(fullGridX.copy(),timeStep,boundaryTemp)
    frames.append(fullGridX[:,:,0].copy())

#Make Gif, saves as "simulation.gif"
makeGif(frames,"simulation.gif",20,500)