Solving the Heat Diffusion Equation (1D PDE) in Python

[BelmanCinematography]

import numpy as np
import matplotlib.pyplot as plt


L = 0.1 #Grosor de la pared en metros
n = 10 #Numero de nodos utilizados
T0 = 20 #Temperatura inicial
T1s = 40 #Temperatura superficial en cara 1
T2s = 20 #Temperatura superficial en cara 2
dx = L/n
alpha = #Difusividad termica K/(Rho*C_p)
t_final = 1800 #Tiempo final en segundos
dt = 60


x = np.linspace(dx/2, L-dx/2, n)

T = np.ones(n)*T0
dTdt = np.empty(n)

t = np.arange(0, t_final, dt)


for j in range(1, len(t)):
plt.clf()
for i in range(1, n-1):
dTdt[i] =
dTdt[0] =
dTdt[n-1] =
T = T + dTdt*dt
plt.figure(1)
plt.plot(x, T)
plt.axis([0, L, 0, 50])
plt.xlabel('Distance (m)')
plt.ylabel('Temperature (C)')
plt.show()
plt.pause(0.01)

[riyazhudda]

Excellent explanation. Thank you!

I would replace all the dT/dt calculations as follows:

# define this once
discrete_calc = lambda Ti_1, Ti, Ti1 : alpha * (-(Ti-Ti_1)/dx**2 + (Ti1-Ti)/dx**2)

Now replace the definition of dT/dt vector everywhere

dTdt[i] = discrete_calc(T[i-1], T[i], T[i+1])
dtDt[0] = discrete_calc(T1s, T[0], T[1])
dtDt[n-1] = discrete_calc(T[n-2], T[n-1], T2s)

[manuelo]

Hello Mr. Powell,

I would like to ask you a question.
One the boundary conditions shound't the distance be dx/2?
Also, on the first for loop the end should be len(t)+1 to finish in the last time step.
Kind regards

[sidebo]

This is a great video, with one small "insidious" bug (meaning it's a bug that introduces an error, but won't cause program to crash, and in this case, would not introduce much error). It's with the 2 boundary conditions, where the distance from the node to the surface should be dx/2, not dx. It's an easy fix in those 2 lines of code.

[boudzzzzsssdddfffg]

For those struggeling with getting the animation right, it might help to add the command plt.ion() in the first for loop. That should stop it from opening the tab every run.

[TriThom]

Is there a way to extend this to be able to apply heat at different points along the x direction? Also, would this result still hold if this was a thin rod, meaning could we consider a thin rod using the 1D equation too? Because in one of the lectures I know you mentioned a wall.

[nabeelkubba]

I'm trying to run this, but the figure only shows one profile then stops, after I close it, it then shows the next profile in time, it is not flowing as a video! Any ideas what is going on? I am using Python 3.7 and attempting to run it from Visual Studio. Thanks!

[gluijk]

Very didactic example. It can be easily replicated in R using vector notation (hyper fast update since the 'for i' spatial loop becomes unnecesary):

T=matrix(T0, nrow=1, ncol=n)
T[1]=Tleft
T[n]=Tright
indices=which(col(T)!=1 & col(T)!=n)
x=seq(from=dx/2, to=L-dx/2, length.out=n) # position array

val=alpha*(dt/dx^2)
for (j in 0:N) {
plot(x, T, type='l', col='red', lwd=2,
xlim=c(0, L), ylim=c(min(Tleft, Tright), max(Tleft, Tright)),
ylab='T (ºC)',
main=paste0('Iteration: ', j, '/', N, ', t=', j*dt, 's'))
abline(v=c(x[1], x[n]), lty=2)
abline(h=0)

# Iterate T in vector notation
- 2*T[indices] + T[indices-1])
}

[AdityaRoyiseast]

hey, why did you use linspace in 14 and arange in 19? why does error pops when i interchange them?

[m]

Quick question, what numerical method technique should I use if I want to find heat flow through a series of flat plates? So imagine we take your picture and keep the heat flow the same direction but push the wall on it's side so its layers are stacked up like pancakes. Any help with that?

And awesome video.

[JamesVestal-dzqm]

The heat transfer equation is from advanced transport phenomena 2, and were coding in python in reaction engineering. Brandon Tatum messaged me what does chemical engineering have to do with my presidential campaign and honestly i dont know. Brandon, my goal is to get healthy as a chemical engineer. Im learning about chemical engineering while waiting for my presidency to happen.

[LateMax]

I wrote my own (very similar) program independently before discovering your guide. I've found that both our programs have the same bug.

This may be a conceptual error with how I understand the heat transfer equation, but as dx decreases, I would expect the simulation to get more accurate, and not do anything weird.
But, as I lower the thickness of the material or increase the number of nodes I find that the numbers begin blasting off by several orders of magnitudes immediately.

I'm not sure why this occurs, but I'm (marginally) happier knowing that this isn't a unique problem. My current method is by lowering the number of nodes as much as is reasonable, and if continuing to lower it would be unreasonable I assume that steady state heat transfer begins almost immediately.

I'd be interested to know what your opinion on the issue is.

[sarazahoor]

I tried to run your code and it is giving me the following error,
$python main.py
File "main.py", line 25
dTdt [0] =
^
SyntaxError: invalid syntax

[bryanvazquez]

thank's a lot for this video!, i've got a question...
how did you make to "animate" the profile temperature at 23:47?

[BelmanCinematography]

Instead of an animation it gives me a list of images at the end, I'm using Python via Jupyter 5.5.0, does anyone know how to turn it into an animation?

[KV-gymr]

This is so cool. Could you also show us how to do this by solving a system of equations using linalg i.e implicit finite difference method?

[AnupumPant]

Amazing! Works like a charm. I added a gaussian source in the middle too. But can you explain how can I add a zero flux boundary condition.

[alip]

sorry can you explain more about T=T+dTdt*dt, because my phthon says, operands could not be broadcast together with shapes (600, ) (10, )

[camimartinez]

thank you i love you you saved my life

[laraibquamar]

sir could you please explain ADI scheme for 2D heat equation in python.