我想用有限差分法来解三维的扩散方程。我想我的主回路有问题。具体而言,离散方程为:
使用Neumann边界条件(仅以一个面为例):
现在代码是:
import numpy as np
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D ##library for 3d projection plots
%matplotlib inline
kx = 15 #Number of points
ky = 15
kz = 15
largx = 90 #Domain length.
largy = 90
largz = 90
dt4 = 1/2 #Time delta (arbitrary for the time).
dx4 = largx/(kx-1) #Position deltas.
dy4 = largy/(ky-1)
dz4 = largz/(kz-1)
Tin = 25 #Initial temperature
kapp = 0.23
Tamb3d = 150 #Ambient temperature
#Heat per unit of area. One for each face.
qq1=0 / (largx*largz)
qq2=0 / (largx*largz)
qq3=0 / (largz*largy)
qq4=0 / (largz*largy)
qq5=0 / (largx*largy)
qq6=1000 / (largx*largy)
x4 = np.linspace(0,largx,kx)
y4 = np.linspace(0,largy,ky)
z4 = np.linspace(0,largz,kz)
#Defining the function.
def diff3d(tt):
w2 = np.ones((kx,ky,kz))*Tin #Temperature array
wn2 = np.ones((kx,ky,kz))*Tin
for k in range(tt+2):
wn2 = w2.copy()
w2[1:-1,1:-1,1:-1] = (wn2[1:-1,1:-1,1:-1] +
kapp*dt4 / dy4**2 *
(wn2[1:-1, 2:,1:-1] - 2 * wn2[1:-1, 1:-1,1:-1] + wn2[1:-1, 0:-2,1:-1]) +
kapp*dt4 / dz4**2 *
(wn2[1:-1,1:-1,2:] - 2 * wn2[1:-1, 1:-1,1:-1] + wn2[1:-1, 1:-1,0:-2]) +
kapp*dt4 / dx4**2 *
(wn2[2:,1:-1,1:-1] - 2 * wn2[1:-1, 1:-1,1:-1] + wn2[0:-2, 1:-1,1:-1]))
#Neumann boundary (dx=dy=dz for the time)
w2[0,:,:] = w2[0,:,:] + 2*kapp* (dt4/(dx4**2)) * (w2[1,:,:] - w2[0,:,:] - qq1 * dx4/kapp)
w2[-1,:,:] = w2[-1,:,:] + 2* kapp*(dt4/(dx4**2)) * (w2[-2,:,:] - w2[-1,:,:] + qq2 * dx4/kapp)
w2[:,0,:] = w2[:,0,:] + 2*kapp* (dt4/(dx4**2)) * (w2[:,1,:] - w2[:,0,:] - qq3 * dx4/kapp)
w2[:,-1,:] = w2[:,-1,:] + 2*kapp* (dt4/(dx4**2)) * (w2[:,-2,:] - w2[:,-1,:] + qq4 * dx4/kapp)
w2[:,:,0] = w2[:,:,0] + 2 *kapp* (dt4/(dx4**2)) * (w2[:,:,-1] - w2[:,:,0] - qq5 * dx4/kapp)
w2[:,:,-1] = w2[:,:,-1] + 2 *kapp* (dt4/(dx4**2)) * (w2[:,:,-2] - w2[:,:,-1] + qq6 * dx4/kapp)
w2[1:,:-1,:-1] = np.nan #We'll only plot the "outside" points.
w2_uno = np.reshape(w2,-1)
#Plotting
fig = pyplot.figure()
X4, Y4, Z4 = np.meshgrid(x4, y4,z4)
ax = fig.add_subplot(111, projection='3d')
img = ax.scatter(X4, Y4, Z4, c=w2_uno, cmap=pyplot.jet())
fig.colorbar(img)
pyplot.show()
对于5000次迭代(qq6=1000/面积),我们得到:
我只在上面加热。不知何故,我最终的底部表面加热
除此之外,我还试图限制我加热的区域(只是一张脸的一小部分)。当我试着去做的时候,似乎热量只向一个方向传递,而忽略了所有其他方向。将(qq1=1000/面积)应用于一半正面(另一部分是绝热的,q=0),我们最终得出:
这很奇怪。我怀疑我在主循环(或者可能在边界条件)中遇到了一些我找不到的问题
也许您可以尝试一种更受测试驱动的方法。一个帮助是Pythonslice-命令。它在numpy数组访问中速度很快,并且在有限差分中使用时不太容易出错。试一试:
我不太清楚舒尔您代码中的图形显示了什么。因此,通过使用不同的T(或w2)初始化和不完全相同的
kx, ky, kz
值仔细检查它。注意np.meshgrid()
命令中的indexing="ij"
选项。使用np.nan
在某种程度上是危险的,因为它改变了变量globaly这是运行热方程求解器的主程序。请注意,我们不能使用在曲面上显示点的图形来查看任何结果!Neumann条件位于我们看不到的底部,面上的值是不变的定值边界条件
下面是一个切片,沿JY1=常数显示结果:
此代码有问题: w2[:,:,0]=w2[:,:,0]+2kapp(dt4/(dx4**2))*(w2[:,:,,-1]-w2[,:,:,0]-qq5*dx4/kapp) 请再检查一遍。 我最近在做一个类似的项目。你有时间和我分享一些经验吗
相关问题 更多 >
编程相关推荐