3D 扩散/热方程的有限差分法

Question

我正在尝试使用有限差分来求解 3D 扩散方程。我认为我的主循环有问题。特别地，离散方程为：

使用诺依曼边界条件（仅以一个面为例）：

现在的代码：

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/面积），我们得到：

我只在顶部表面加热。不知何故，我最终导致底面变热。

除此之外，我试图限制加热的区域（只是一张脸的一小部分）。当我尝试这样做时，似乎热量只向一个方向传递，忽略了所有其他方向。对正面的一半（另一部分是绝热的，q = 0）应用（qq1 = 1000/面积），我们最终得到：

这很奇怪。我怀疑我在主循环（或者可能在边界条件）中遇到了一些我没有发现的问题。

Answer 1

也许您可以尝试一种更加测试驱动的方法。Python slice命令是一个帮助。它在 numpy 数组访问中速度很快，并且在有限差分中更不容易出错。试一试：

\n

import numpy as np\n\ndef get_slices(mx, my, mz):\n    jxw, jxp, jxe = slice(0,mx-2), slice(1,mx-1), slice(2,mx)  # west,  center, east\n    jys, jyp, jyn = slice(0,my-2), slice(1,my-1), slice(2,my)  # south, center, north\n    jzb, jzp, jzt = slice(0,mz-2), slice(1,mz-1), slice(2,mz)  # botoom,center, top\n    return jxw, jxp, jxe,   jys, jyp, jyn,    jzb, jzp, jzt\n\nnx, ny, nz = 15, 10, 15\njxw, jxp, jxe, \\\njys, jyp, jyn, \\\njzb, jzpc, jzt = get_slices(nx, ny, nz)\n\nax, ay, az = np.arange(nx), np.arange(ny), np.arange(nz)\n\nprint(\'axW\', ax[jxw])\nprint(\'axP\', ax[jxp])\nprint(\'axE\', ax[jxe])\n\naxW [ 0  1  2  3  4  5  6  7  8  9 10 11 12]\naxP [ 1  2  3  4  5  6  7  8  9 10 11 12 13]\naxE [ 2  3  4  5  6  7  8  9 10 11 12 13 14]\n

\n

我不太清楚代码中的图形向我们展示了什么。kx, ky, kz因此，请通过对 T（或 w2）使用不同的初始化以及不完全相同的值来仔细检查。请注意命令indexing="ij"中的选项np.meshgrid()。使用np.nan有点危险，因为它会全局更改变量。

\n

import numpy as np\nfrom matplotlib import pyplot as plt\nfrom matplotlib import cm\nfrom mpl_toolkits.mplot3d import Axes3D\n\ndef get_grid(mx, my, mz, Lx,Ly,Lz):\n    ix, iy, iz = Lx*np.linspace(0,1,mx), Ly*np.linspace(0,1,my), Lz*np.linspace(0,1,mz)\n    x, y, z = np.meshgrid(ix,iy,iz, indexing=\'ij\')\n    print(\'ix\', ix), print(\'iy\', iy), print(\'iz\', iz)\n    return x,y,z\n\ndef plot_grid(x,y,z,T):\n    def plot_boundary_only(x,y,z,T):\n        mx, my, mz = x.shape\n        x[1:-1, 1:-1, 1:-1],y[1:-1, 1:-1, 1:-1],z[1:-1, 1:-1, 1:-1],T[1:-1, 1:-1, 1:-1] = np.nan, np.nan, np.nan, np.nan     # remove interior\n        x[1:-1, 1:,0], y[1:-1, 1:,0], z[1:-1, 1:,0],  T[1:-1, 1:,0] = np.nan, np.nan, np.nan, np.nan                         # remove bottom\n        x[1:-1, my-1, 1:-1], y[1:-1, my-1, 1:-1], z[1:-1, my-1, 1:-1], T[1:-1, my-1, 1:-1] = np.nan, np.nan, np.nan, np.nan  # remove north face\n        x[0, 1:, :-1], y[0, 1:, :-1], z[0, 1:, :-1], T[0, 1:, :-1] = np.nan, np.nan, np.nan, np.nan                          # remove west face\n        return x,y,z,T\n    \n    x,y,z,T = plot_boundary_only(x,y,z,T)   \n    fig = plt.figure(figsize=(15,15))\n    ax = fig.add_subplot(111, projection=\'3d\')\n    img = ax.scatter(x,y,z, c=T.reshape(-1), s=150, cmap=plt.jet())\n    fig.colorbar(img)\n    \n    #ax.zaxis.set_rotate_label(False) # To disable automatic label rotation\n    ax.set_ylabel(\'Y\')\n    ax.set_xlabel(\'X\', rotation=0)\n    ax.set_zlabel(\'Z\', rotation=0)\n    plt.savefig("cube.png")\n    plt.show()\n    \ndef init_T(x,y,z):\n    T = np.zeros_like(x)\n    case = \'xz\'\n    if case == \'x\': T = x\n    if case == \'y\': T = y\n    if case == \'z\': T = z\n    if case == \'xz\': T = x+z\n    return T\n\nnx, ny, nz = 35, 20, 30\nLx, Ly, Lz = nx-1, ny-1, nz-1\nx,y,z = get_grid(nx, ny, nz, Lx,Ly,Lz)  # generate a grid with mesh size \xce\x94x = \xce\x94y = \xce\x94z = 1\nT = init_T(x,y,z)\nplot_grid(x,y,z,T)\n

\n

\n

这是运行热方程求解器的主程序。请注意，我们不能使用显示表面点的图形来查看任何结果！诺依曼条件位于我们看不到的底部，并且面处的值是固定值边界条件，不会改变。

\n

def set_GradBC(u):\n    #--- Neuman BC at bottom boundary ---\n    u[:, :, 0] = u[:, :, 1]\n    return u\n\n\ndef dT(u,DU):\n    DU[jxp,jyp,jzp] = (u[jxe,jyp,jzp] - 2.0*u[jxp,jyp,jzp] + u[jxw,jyp,jzp])*Dx + \\\n                      (u[jxp,jyn,jzp] - 2.0*u[jxp,jyp,jzp] + u[jxp,jys,jzp])*Dy + \\\n                      (u[jxp,jyp,jzt] - 2.0*u[jxp,jyp,jzp] + u[jxp,jyp,jzb])*Dz\n    return DU\n\n#mmmmmmmmmmmmm main mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm\n#---- set the parameters ------------------------\ndt = 0.01                                    # time step size\nnIterations = 5000                           # number of iterations\nnx, ny, nz = 40,20,30                       # number of grid points\nLx, Ly, Lz = nx-1, ny-1, nz-1                 # physical grid size\nKx, Ky, Kz = 1, 1, 1                          # diffusion coeficients\ndx, dy, dz = Lx/(nx-1), Ly/(ny-1), Lz/(nz-1)  # mesh size\nDx, Dy, Dz = Kx/dx**2, Ky/dy**2, Kz/dz**2     # equation coefficients\n\n#---- initialize the field variables -----------------\nx,y,z = get_grid(nx, ny, nz, Lx,Ly,Lz)       # generate the grid\nT = init_T(x,y,z)                            # initialize the temperature\nDT = np.zeros_like(T)                        # initialize the change of temperature\njxw, jxp, jxe, \\\njys, jyp, jyn, \\\njzb, jzp, jzt = get_slices(nx, ny, nz)       # slices for finite differencing\n\n#---- run the solving algorithm ----------------------\nfor jter in range(nIterations):\n    T = set_GradBC(T)              # set Neumann boundary condition\n    T = T + dt*dT(T,DT)            # update\n

\n

这是一个沿 JY1=constant 显示结果的切片：

\n

def plot_y_slice(JY1, x,y,z,T):\n    T2 = T[:,JY1,:]\n    X =  x[:,JY1,:]\n    Z =  z[:,JY1,:]\n\n    fig = plt.figure(figsize=(15,15))\n    ax = fig.add_subplot(111)\n    ax.set_aspect(\'equal\')\n    plt.contourf(X,Z,T2, 40,alpha=0.99)\n    plt.savefig("y_slice.png")\n    plt.show()\n    \nJY1 = 10\nplot_y_slice(JY1, x,y,z,T)\n

\n

\n