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

问题描述

我对使用numpy和sympy库还是很陌生。抱歉，如果我在上面有很多打印件，我只是想检查代码是否正常工作。我正在尝试解决这个2D热方程问题，并且在理解如何添加初始条件（温度为30度）和添加齐次Dirichlet边界条件（例如，板的两边（即顶部，底部和2面）。作为示例，我尝试将板顶部的温度添加为1000度，但出现错误“ IndexError：索引18超出轴18尺寸18的范围”。我已经设法将函数添加到代码上，并且似乎也可以正常输出，但在添加条件上就有点丢失了。我试图在线查找类似的问题，但无法真正理解正在发生的事情。

以下是原始2D热方程的链接：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,k)+(dt*a(((T(i+1,k)-2*T(i,k)+T(i-1,k))/(dx**2))+(a*dt(((T(i,j+1,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,:])

解决方法

我认为我已经设法将其整理出来，并根据@chris的建议添加了每一侧的温度后检查了这些值。https://scipython.com/book/chapter-7-matplotlib/examples/the-two-dimensional-diffusion-equation/

dirichlet numpy pde python sympy