3D-Heat Equation Numerical Simulation

3D Heat Equation Numerical Simulation

Introduction

The heat equation is a common thermodynamics equation first introduced to undergraduate students. The heat equation is a parabolic differential equation that describes the variation in temperature in any given region over time. If material properties are known this can describe the temperature variation in almost any shape and substance. Although in the example below is only in 2D space, I use the term 3D because we are calculating the variation in time and space. Thus time is our third dimension, mathematically speaking.

I utilized the finite difference method to solve the partial differential equation. This is a more advanced numerical solving technique as compared to the previous Euler method. I used this method as its relatively intuitive to those with a calculus background. This method is computationally expensive compare to more implicit methods however I believe utilizing intuitive methods first, before moving onto more complex and abstract numerical solving techniques.

Mathematical Background

Before we get into the programming let me lay out the mathematical/physical ground work we will program into MATLAB. to provide a simplistic example we will determine the temperature variation in a 2D square plate with a heat source applied to all four sides of the plate. One can mentally visualize what would happen in this case. The edges of the plate will be the hottest in the beginning and over time the center of the plate will become warmer. Starting with a physically intuitive problem is a great way to check your numerical approximation of the governing equations.

First the 3D heat equation:

Equation (1) states that the temperature over time of your object changes based off the surrounding finite elements around your desired point. If for example you care how hot the middle of a plate will be, the temperature will solely be dependent on the area around the point. Utilizing the finite difference method we can numerically approximate the temperature of any point in an object. Next I will explain the finite difference method and then how to apply it to equation (1).

The finite difference method with a second order central estimate numerically approximates a value based off the surrounding conditions. for second order derivatives, the finite difference method estimates that the change in Temperature (T), its dependent on the temperature immediately next to it in the x and y directions.

Where i and j are the points immediately adjacent to the point we are trying to estimate (x and y). This can be shown in the graphic below:

Figure 1: x,y is the point we are trying to estimate, dependent on x-i and x+i

Figure 2: x,y is the point we are trying to estimate, dependent on y-j and x+j

Figure 1 is the visual representation of equation (2) and Figure 2 is the visual representation of equation (3). The point x,y is solely dependent on the temperature directly adjacent to it in the x and y directions. As the numerical approximation steps through time, we can now estimate how the temperature changes over time instead of just a snap shot in time. This is accomplished by plugging equation 2, 3 and equation 4 into equation 1 and moving dt to the right hand side of the equation. This can be seen below:

Equation 4 demonstrates how we plug the numerical approximation into the partial differential equation into the heat equation. Notice the final term. This is a by product of the fact that we are calculating the CHANGE in temperature thus we must incorporate this change into the overall temperature (moved the temperature to the right hand side of the equation).

Finally we add the change in temp (equation 5) to calculate the next temperature in time (T). Note that T when starting the calculation is the initial starting conditions (the temperature profile of the object).

Equation 6 utilizes the calculated differential change in temperature (equation 5) added to the current temperature to predict the "future" temperature. The finite difference method utilizes the current conditions to predict the future temperature. Then iterates over time to find a steady state solution. Next I will go into the python code example to simulate the temperature of a flat plate with 300 degrees Celsius applied to the out boundaries and how the entire plate changes temperature over time.

3D-Heat Equation Python Code

Equation 6 utilizes the calculated differential change in temperature (equation 5) added to the current temperature to predict the "future" temperature. The finite difference method utilizes the the current conditions to predict the future temperature. Then iterates over time to find a steady state solution. Next I will go into the python code example to simulate the temperature of a flat plate with 200 degrees Celsius applied to the outer boundaries and how the entire plate changes temperature over time.

First I import the functions I used to conduct the math functions and movie export functions:

import numpy as np
import matplotlib.pyplot as plot #import plotting library
from matplotlib import cm #import color for surface plot
import os, sys
import matplotlib.animation as ani #importing animation library

Now I set up the initial conditions. Alpha is the thermal diffusivity and then I set the starting temperatures and number of points in my grid.

alpha=97 # (mm^2/s) thermal diffusivity, alluminum
#length of sides of heated square plate
Lx=152 #(mm)
Ly=152 #(mm)
#number of points
N=50 #number of x and y points x*y = total number of points
#Discretize my space
Xvec=np.linspace(0,Lx,N)
Yvec=np.linspace(0,Ly,N)
dx=Xvec[2]-Xvec[1] #since equally spaced, i set a constant dx
dy=Yvec[2]-Yvec[1] #since equally spaced, I set a constant dy
#Discretize time
#dt=0.5*(dx**2)/(2*alpha) #dt needed for stability can do larger however good rule of thumb
dt=0.01
print(dt)
time=15 #sec, how long i want to run the simulations for
tvec=np.linspace(0,1,100) #this is how long i run my numerical approximation. 0 to 100 sec
#Inital Boundary Conditions
T=np.full((Yvec.size,Xvec.size),20.0) #entire plate is at 20 degrees C

T[:,0]=200.0 #200 degrees C applied to left side of plate
T[-1,:]=200.0 #200 degrees C applied to bottom of plate
T[0,:]=200.0 #200 degrees C applied to top of plate
T[:,-1]=200.0 #200 degrees C applied to right side of plate

Now I implement the numerical solver for the 3D heat equation as shown in eq 5 and 6.

or i in range (0,frn):
Told=T

for ty in range(1,Yvec.size-1):

for tx in range(1,Xvec.size-1):

du=(dt*(alpha*(Told[tx+1,ty]-2*Told[tx,ty]+Told[tx-1,ty])/dx**2+alpha*...(Told[tx,ty+1]-2*Told[tx,ty]+Told[tx,ty1])/dy**2)+Told[tx,ty])-T[tx,ty]
#print(du)

T[tx,ty]=T[tx,ty]+du

#img.append((plot.contour(X,Y,T,30)))
#Tout.append(T)
Tout[i]=T

Finally I set up the graph the way I desire and turn it into a video:

def frame(frame_number,Tout,plot):
plot[0].remove()
plot[0]=ax1.plot_surface(X,Y,Tout[frame_number,:,:],cmap=cm.jet)

fig = plot.figure()
ax1 = plot.axes(projection="3d")
ax1.set_xlabel('x axis (mm)')
ax1.set_ylabel('y axis (mm)')
ax1.set_zlabel('Temperature (C)')
ax1.set_xlim(0,Lx)
ax1.set_ylim(0,Ly)
ax1.set_zlim(0,200)
x = np.linspace(0, Lx, N)
y = np.linspace(0, Ly, N)
X, Y = np.meshgrid(x, y)
plot1=[ax1.plot_surface(X,Y,Tout[0,:,:])]

#plot=[ax1.plot_surface(X,Y,Tout[:,:,0])]
animation=ani.FuncAnimation(fig,frame,frn,fargs=(Tout,plot1),interval=60,blit=False)
#plot.show()

#save video
#path for ffmpeg using imagemagick
ff_path = os.path.join('C:/', 'ImageMagick-7.0.10-Q16', 'ffmpeg.exe')
plot.rcParams['animation.ffmpeg_path'] = ff_path
#path if i want to convert the .mp4 to a gif
#imgk_path = os.path.join('C:/', 'ImageMagick-7.0.10-Q16', 'convert.exe')
#plot.rcParams['animation.convert_path'] = imgk_path
FFwriter=ani.FFMpegWriter(fps=220) #declar my writer and frame rate
animation.save(os.path.join('C:/','Users','Michael','Desktop','video','tesst.mp4'),writer=FFwriter) #saving function

The frame function takes the plot and writes it to a frame then stores it and then using imageMagick and animation.save to convert and save the video to the file path that is defined the animation.save (final line in the code). I utilized imageMagick as it was the easiest for me to use however there are plenty of other packages out there that you can use. If I recieve enough requests i can post a tutorial in how to set up imageMagick but a quick google search yields enough information.

When all this code is put together it will out put the following video:

As you can see in the video, the solution starts out with the edges at 200 degrees Celsius. As time progresses, the center of the plate heats up with no odd anomalies. Now that the solver is confirmed to work, one can change the grid shape and apply it to the solver and change the boundry conditions to see how any 2D shape will heat up. One could apply a point heat source in the middle and see how it will propagate through the aluminum plate.

As you can see using a numerical solver to see how the physical behavior of any shape or material could be a huge benefit. Just by changing the thermal diffusivity one could design an optimal material to retain or dissipate heat. I chose aluminum because I use in my home an aluminum cookie sheet for cooking. I was curious to see how the heat moves throughout the pan and how long it takes to heat up over time. Although this is a simple case but this can be extrapolated to many other problems and situations just by implementing the finite central difference method.