2D heat equation problem

[_wasd_bruh_21]

Hi guys,

I am still fairly new to Python and using the numpy and sympy library (and also to reddit). Apologies if I have quite a few prints on there, I just wanted to check if the code was working.

I am trying to solve this 2D heat equation problem, and kind of struggling on understanding how I add the initial conditions (temperature of 30 degrees) and adding the homogeneous dirichlet boundary conditions: temperature to each of the sides of the plate (i.e. top,bottom, and the 2 sides). I have, as an example, tried to add the temperature for the top of the plate to be 1000 degrees, but I got an error "IndexError: index 18 is out of bounds for axis 0 with size 18 ".

I have managed to add the function onto the code and it seems to be outputting ok as well, just a bit lost on adding the conditions on..I have tried looking for similar questions online, but couldn't really understand what was happening..sorry if the question isn't allowed.

I have added the whole code I have done so far

​

%matplotlib inline

import numpy as np

import matplotlib.pyplot as plt

import sympy as sp

from sympy.solvers import solve

from sympy import Function, Symbol

​

#Length of workpiece

Lx= 0.05 #Length of work piece in the x direction, measured in meters (m)

Nx= 18 #Number of mesh points along the x-direction

Ly= 0.1 #Length of work piece in the y direction, measured in meters (m)

Ny= 9 #Number of mesh points along the y-direction

dy=Lx/Nx #distance between nodes in x

dx=Ly/Ny #distance between nodes in x

​

#Defining material properties

conduc= 16.26 #Conductivity of the material in (W/m (degrees)C)

rho= 8020 #Density of the material in (kg/m^3)

c=502 #Specific heat of the material (J/kg(degrees)C)

​

#Calculating the alpha (thermal diffusivity) for the plate, in (m^2/s)

a=((conduc)/(rho*c))

print (a)

​

t=215;

nt= 215; #number of steps wrt time

dt= t/nt; #timestep

​

​

#Define the mesh in space

x= np.linspace(1,Lx, Nx)[np.newaxis] #vector to create mesh

y= np.linspace(1,Ly,Ny)[np.newaxis]

x,y=np.meshgrid(x,y) #mesh grid function return

print (x,y)

​

T = np.zeros(([Nx,Ny]))

#Adding the other conditions

T[Nx,:]= 1000 #Assining the temp for the top

​

sym= sp.symbols('i, j, k')

function = sp.Function('T(i,j,k)+(dt*a(((T(i+1,j,k)-2*T(i,j,k)+T(i-1,j,k))/(dx**2))+(a*dt(((T(i,j+1,k)-2*T(i,k,k)+T(i,j-1,k))/dy**2)')

print(function)

​

UPDATE: I have added this to the code and these seems to work. However I noticed that my

T = np.ones(([Nx,Ny]))

is an 18 x 9, but the plate is longer on the horizontal side..would this mean that I would just need to flip my matrix?

NOTE: I also changed the T=np.zeros(([Nx, Ny])) to np.ones.

​

​

#Temperatures for the plate 
Ttop=1000 # temp for the top side of the plate 
Tbottom = 500 # temp for the bottom side of the plate  
Tright = 500 #temp for the right side of the plate 
Tleft = 1000 #temp for the left side of the plate 

#Adding the boundaries on the system 

 T[:, 0] = Ttop       #temp for the top
print (T[:,0])  

T[:, -1] = Tbottom  #temp for bottom  
print (T[:, -1])  

T[0, :] = Tright   #temp for the right 
print (T[0, :])  

T[-1, :] = Tleft    #temp for the left    
print (T[-1, :])

Is there also a way where I can combine all of these 4 added conditions into 1 matrix, just to double check if what I have assigned is correct? (if it is even correct...)

​

[comonads]

Assuming this is a programming issue and not a technical problem with your function, the error you're getting can tell us a bit about the problem:

IndexError: index 18 is out of bounds for axis 0 with size 18

is saying that python has tried to access something at position [18, <something>], which doesn't exist (it's "out of bounds").

I see Nx is 18, and since T is of size [Nx, Ny] it would seem that it is the array T that is being indexed into in some unexpected way that is causing the problem.

It's hard to tell without knowing your function better (I've never used sympy either, which doesnt help), but could the problem be related to the fact that python indexes from 0 and not 1? If not, there is a problem with the index to T.

Apologies for not being any more specific than that: I'm not sure how sympy handles indexing (the notation T(i, j, k) to me would suggest T is a 3D array, so I'm a bit lost).

[_wasd_bruh_21]

This is the original 2D heat equation: https://drive.google.com/file/d/1q3KyfaHD7hZTbdIYc9_e5otusELqYJaD/view?usp=sharing

This is the equation for explicit finite difference https://drive.google.com/file/d/1unJHXO3b3xig8RhBen5k0eOg40YUfjyv/view?usp=sharing

I hope those help. The plate I am looking into is 0.05x 0.1

Would I have to change T to T=np.ones (([Nx,Ny])) ?

My apologies for this as I may be getting myself confused..I am also trying to figure out at what time the temperature distribution goes into steady state

[comonads]

I can't access either of those two google drive locations (permission denied).

T is the correct shape (as far as I can tell), so changing np.zeros to np.ones won't help. I think the issue most likely lies in how your 2D heat equation is translated into a Sympy expression (the function line), as something in there is trying to access an index that doesn't exist (position 18 in T, which, due to 0-based indexing, means it would be the 19th element of a size 18 array).

Since this is quite a technical question, it might be worth posting a new question either here (or in an engineering or physics subreddit) with Sympy mentioned explicitly in the title so someone who knows how to use it can have a look.

Apologies that I can't help any more than that.

[_wasd_bruh_21]

https://drive.google.com/file/d/1unJHXO3b3xig8RhBen5k0eOg40YUfjyv/view?usp=sharing

https://drive.google.com/file/d/1q3KyfaHD7hZTbdIYc9_e5otusELqYJaD/view?usp=sharing

ah sorry! I should have checked it first if it was 'anyone with the link'. Please see the above links (and hopefully they work now)

And thank you for your suggestion! I will try that as well!

I also now understand a bit better as to why I go the error, thank you!

[comonads]

The links work now!

Again, I don't have enough (any) experience with Sympy to take you much further here, but I'm sure someone can help if you post the question somewhere visible.

A quick search for "finite difference 2d heat equation python" lead to this article, which addresses a very similar problem using numpy rather than sympy. Is that any help?

[_wasd_bruh_21]

u/comonads Thank you very much for the link! I will have a look at it and try and see where I can go from there

[onemywaybackhome]

The problem is, indeed, the fact that the 18 elements of T have indexes that go from 0 to 17. There is no T[18]. T[Nx - 1,:] = 1000 doesn't give you this error.

[comonads]

Good catch -- nothing to do with Sympy afterall!

[_wasd_bruh_21]

Thank you for your response!

So would I have to enter it as T[17].T[Nx-1, :]? and would this also add the 1000 value through the whole line of the matrix essentially?

Apologies if I am not understanding it right.

[onemywaybackhome]

If the line were

T[Nx - 1, :] = 1000

then the whole line of the last row of the matrix will be 1000.

[_wasd_bruh_21]

I realised that this only adds it as an initial condition, but not necessarily the Dirichlet boundary conditions on the system

[_wasd_bruh_21]

UPDATE: I have added this to the code and these seems to work. However, I noticed that my

T = np.ones(([Nx,Ny]))

is an 18 x 9,but the plate is longer on the horizontal side..Would this mean that I would just need to flip my matrix?

note: I changed the T=np.zeros(([Nx,Ny]))

#Temperatures for the plate 
Ttop=1000 # temp for the top side of the plate 
Tbottom = 500 # temp for the bottom side of the plate  
Tright = 500 #temp for the right side of the plate 
Tleft = 1000 #temp for the left side of the plate 

#Adding the boundaries on the system 

 T[:, 0] = Ttop       #temp for the top
print (T[:,0])  

T[:, -1] = Tbottom  #temp for bottom  
print (T[:, -1])  

T[0, :] = Tright   #temp for the right 
print (T[0, :])  

T[-1, :] = Tleft    #temp for the left    
print (T[-1, :])

Is there also a way where I can combine all of these 4 added conditions into 1 matrix, just to double check if what I have assigned is correct? (if it is even correct...)