2D 热方程 - 添加初始条件并检查 Dirichlet 边界条件是否正确

Question

我对使用 numpy 和 sympy 库还是很陌生。抱歉，如果我在那里打印了很多照片，我只是想检查代码是否正常工作。 我正在尝试解决这个 2D 热方程问题，并且有点难以理解我如何添加初始条件（30 度的温度）并添加均匀的狄利克雷边界条件：温度到板的每一侧（即顶部，底部和两侧）。例如，我尝试将板顶部的温度添加为 1000 度，但出现错误“IndexError: index 18 is out of bounds for axis 0 with size 18”。 我已经设法将函数添加到代码中，它似乎也可以输出，只是在添加条件时有点迷失......我尝试在网上寻找类似的问题，但无法真正理解发生了什么。 .

这里是原始二维热方程的链接：https://drive.google.com/file/d/1q3KyfaHD7hZTbdIYc9_e5otusELqYJaD/view?usp=sharing

这里是显式有限差分的最终方程的链接：https://drive.google.com/file/d/1unJHXO3b3xig8RhBen5k0eOg40YUfjyv/view?usp=sharing

i;j 是位置（节点编号），k= 时间（时间步长编号）

我也在尝试找出温度达到稳定状态的时间

我已经添加了到目前为止我所做的全部代码

%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)

---

更新：我已将此添加到代码中，这些似乎有效。但是，我注意到我的T = np.ones(([Nx,Ny])) 是 18 x 9，但水平边的盘子更长。这是否意味着我只需要翻转矩阵？

#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, :])

Answer 1

我想我已经设法对其进行了整理，并在根据@chris 的建议查看每一侧的温度后检查了值。https://scipython.com/book/chapter-7-matplotlib/examples/the-two-dimensional-diffusion-equation/