2D Heat Conduction with Python

Question

the output of the study///// I just started learning Python, so I am new with python. I have written a simple code for 2D Heat Conduction. I don't know what is the problem with my code. The result is so strange.I think the Temperature distribution is not shown correctly. I have searched about it a lot but unfortunately I could not find any answer for my problem . can anyone help me? Thank you

# Library

import numpy                           
from matplotlib import pyplot     



# Grid Generation

nx = 200
ny = 200                                    
dx = 2 / (nx-1)
dy = 2 / (ny-1)

# Time Step

nt = 50                                                                   
alpha = 1                                 
dt = 0.001                     

# Initial Condition (I.C) and Boundry Condition (B.C)

T = numpy.ones((nx, ny))                         # I.C (U = Velocity)
x = numpy.linspace(0,2,nx)                       # B.C
y = numpy.linspace(0,2,ny)                       # B.C

Tn = numpy.empty_like(T)                         #initialize a temporary array
X, Y = numpy.meshgrid(x,y)

T[0, :] = 20          #  B.C
T[-1,:] = -100        #  B.C
T[:, 0] = 150         #  B.C
T[:,-1] = 100         #  B.C
# Solver
###Run through nt timesteps
    
for n in range(nt + 1): 
    Tn = T.copy()
        
    T[1:-1, 1:-1] = (Tn[1:-1,1:-1] + 
                        ((alpha * dt / dx**2) * 
                        (Tn[1:-1, 2:] - 2 * Tn[1:-1, 1:-1] + Tn[1:-1, 0:-2])) +
                        ((alpha * dt / dy**2) * 
                        (Tn[2:,1: -1] - 2 * Tn[1:-1, 1:-1] + Tn[0:-2, 1:-1])))
        
    T[0, :] = 20          # From B.C
    T[-1,:] = -100        # From B.C
    T[:, 0] = 150         # From B.C
    T[:,-1] = 100         # From B.C

   
fig = pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.contourf(X, Y, T)
pyplot.colorbar()
pyplot.contour(X, Y, T)
pyplot.xlabel('X')
pyplot.ylabel('Y');

What is "strange"? Please clearly explain what the problem is. Show the output of your code and what the expected output is. Is the problem with the calculation? With the vizualization? — Diziet Asahi, Nov 11 '20 at 13:18
Thank you. I added the photo of the result. The problem is, I think the temperature distribution is not shown correctly. any help will be appreciated. — Mohammad, Nov 11 '20 at 14:08

score 2 · Accepted Answer · answered Nov 11 '20 at 14:06

You are using a Forward Time Centered Space discretisation scheme to solve your heat equation which is stable if and only if alpha*dt/dx**2 + alpha*dt/dy**2 < 0.5. With your values for dt, dx, dy, and alpha you get

alpha*dt/dx**2 + alpha*dt/dy**2 = 19.8 > 0.5

Which means your numerical solution will diverge very quickly. To get around this you need to make dt smaller and/or dx and dy larger. For example, for dt=2.5e-5 and the rest as before you get alpha*dt/dx**2 + alpha*dt/dy**2 = 0.495, and the solution will look like this after 1000 iterations: Alternatively, you could use a different discretisation scheme like for ex the API scheme which is unconditionally stable but will be harder to implement.

2D Heat Conduction with Python

1 Answers1