I am still fairly new to using the numpy and sympy library. 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..
Here is a link to the original 2D heat equation:https://drive.google.com/file/d/1q3KyfaHD7hZTbdIYc9_e5otusELqYJaD/view?usp=sharing
Here is the link to the final equation for explicit finite differencing: https://drive.google.com/file/d/1unJHXO3b3xig8RhBen5k0eOg40YUfjyv/view?usp=sharing
i;j are locations (node numbers) and k= time (time step number)
I am also trying to find the time to which the temperature reaches the steady state
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.ones(([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?
#Temperatures for the plate
Ttop=1000 # temp for the top side of the plate
Tbottom = 500 # temnp 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, :])
I think I have managed to sort it out and have checked the values after adding in the temperatures for each of the side from @chris's recommendation on looking at.https://scipython.com/book/chapter-7-matplotlib/examples/the-two-dimensional-diffusion-equation/