计算科学 - 求解具有 Neumann 边界的 3-D 热方程 - 吾爱随笔录

我想在给定尺寸的 3-D 立体立方体中求解热流的泊松 PDE $x$ , $y$ ，和 $z$ ：

ρ C \frac{\partial T}{\partial t} = k Δ T

$\rho C\frac{\partial T}{\partial t} = k \Delta T$

立方体被恒定的热通量照射 $I(x,y)$ 在 $z=0$ 表面。边界应满足由下式给出的辐射冷却的诺伊曼条件 $\frac{\partial T}{\partial x_i} = \pm\epsilon\sigma(T^4-300^4)$ , 符号取决于你看的地方 $x_i=$ 最大或 $x_i=0$ . 由于辐照 $z=0$ 表面是一个例外，条件应该是 $\frac{\partial T}{\partial x_i} = \epsilon\sigma(T^4-300^4) + I(x,y)$ . 最后一项表示辐射通量的强度分布。

我决定采用我在这个线程中使用 Python 找到的数值方法。结果，我想出了以下代码。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplotn3d

dt = 3
di = 0.01

#Thermal conductivity
k = 1.38
#Density
rho = 2202
#Specific Heat capacity
cp = 745
#Thermal diffusivity
alpha = k/(rho*cp)
#Emissivity
emissivity = 0.79

sigma = 5.670367e-8

y_max = 0.3
x_max = 0.3
z_max = 0.3
t_max = 60

P = 1
omega = 0.02

#Function that generates a gaussian intensity pattern I(x,y) with power P and beam waist omega

def get_beam(P, omega, di, x_max, ymax, z_max):
  x = np.arange(0,x_max+di,di)
  y = np.arange(0,y_max+di,di)
  z = np.arange(0,z_max+di,di)
  t = np.arange(0,t_max+dt,dt)
  r = len(t)
  cx = len(x)
  cy = len(y)
  cz = len(z)
  I = np.zeros([r, cx, cy, cz])
  for jx in range(0,cx-1):
      for jy in range(0,cy-1):
          I[:,jx,jy,0] = 2*P/(np.pi*omega**2)*np.exp(-2*(abs(di*(jx-cx/2))**2 + abs(di*(jy-cy/2))**2)/(omega**2))
  return I

def FTCS(dt,dy,t_max,x_max,y_max,z_max,k,T0, I):
  s = alpha*dt/di**2
  x = np.arange(0,x_max+di,di)
  y = np.arange(0,y_max+di,di)
  z = np.arange(0,z_max+di,di)
  t = np.arange(0,t_max+dt,dt)
  r = len(t)
  cx = len(x)
  cy = len(y)
  cz = len(z)
  #Initialize mesh
  T = np.ones([r,cx, cy, cz])*T0
  #time loop
  for n in range(0,r-1):  
      print(n)
      #y loop
      for jy in range(0,cy-1):   
          #x loop 
          for jx in range(0,cx-1):   
              #z loop 
              for jz in range(0,cz-1):  
                  #Compute T at next time step: Tn+1 = Tn + d^2T/dx^2 + d^2T/dy^2 + d^2T/dz^2  
                  T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 


      #Boundary Condition
      #Set boundary condition for x on both sides   
      T[n+1,cx-1,:,:] = T[n,cx-1,:,:] + s*(T[n,cx-2,:,:] - 2*T[n,cx-1,:,:] + T[n,cx-2,:,:] - 2*di*emissivity*sigma/k*(T[n,cx-1, :,:]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,cx-1,:,:]*300**3))
      T[n+1,0,:,:] = T[n,0,:,:] + s*(T[n,1,:,:] - 2*T[n,0,:,:] + T[n,1,:,:] - 2*di*emissivity*sigma/k*(T[n,0,:,:]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,0,:,:]*300**3))

      #Set boundary condition for y on both sides 
      T[n+1,:,cy-1,:] = T[n,:,cy-1,:] + s*(T[n,:,cy-2,:] - 2*T[n,:,cy-1,:] + T[n,:,cy-2,:] - 2*di*emissivity*sigma/k*(T[n,:,cy-1,:]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,:,cy-1,:]*300**3))
      T[n+1,:,0,:] = T[n,:,0,:] + s*(T[n,:,1,:] - 2*T[n,:,0,:] + T[n,:,1,:] - 2*di*emissivity*sigma/k*(T[n,:,0,:]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,:,0,:]*300**3))

      #Set boundary condition for z on both sides 
      T[n+1,:,:,cz-1] = T[n,:,:,cz-1] + s*(T[n,:,:,cz-2] - 2*T[n,:,:,cz-1] + T[n,:,:,cz-2] - 2*di*emissivity*sigma/k*(T[n,:,:,cz-1]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,:,:,cz-1]*300**3))
      T[n+1,:,:,0] = T[n,:,:,0] + s*(T[n,:,:,1] - 2*T[n,:,:,0] + T[n,:,:,1] - 2*di*emissivity*sigma/k*(T[n,:,:,0]**4-300**4) + 2*di/k*I[n,:,:,0])#- 2*di*4*emissivity*sigma/k*(T[n,:,:,0]*300**3) + 2*di*I[n,:,:,0])

  return x,y,z,t,T,r,s 

I = get_beam(P, omega, di, x_max, y_max, z_max)
x,y,z,t,T,r,s = FTCS(dt,di,t_max,x_max,y_max,z_max,k,300,I)


fig = plt.figure()
X, Y = np.meshgrid(x, y)
ax = plt.axes(projection='3d')
ax.plot_surface(X, Y, T[-1,:,:,0])
ax.set_xlabel('x-axis')
ax.set_ylabel('y-axis')

结果看起来很有希望，但不等于我找到的参考资料。此外，结果随着网格大小迅速扩展 $d_i$ 并因巨大的点密度而崩溃。我不知道我是否犯了错误以及是否犯了错误。你能帮助我吗？

//------------------------------------------------//

我更新了代码以包含正确的边界。循环现在看起来像这样。我知道它有点笨重，但我想我以后会照顾它。尽管如此，网格尺寸的依赖性仍然存在。恐怕我只是忘记了一个 $di$ 某处，但我找不到在哪里。

    for n in range(0,r-1):
    print(n)
    for jy in range(0,cy-1):
        for jx in range(0,cx-1):
            for jz in range(0,cz-1):
                T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 

                if jx == cx-1:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx-1,jy,jz] + 2*di*emissivity*sigma/k*(T[n,jx, jy, jz]**4-300**4)) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1])
                elif jx == 0:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx-1,jy,jz] - 2*di*emissivity*sigma/k*(T[n,jx, jy, jz]**4-300**4)) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 
                elif jy == cy-1:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy-1,jz] + 2*di*emissivity*sigma/k*(T[n,jx,jy,jz]**4-300**4)) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 
                elif jy == 0:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy-1,jz] - 2*di*emissivity*sigma/k*(T[n,jx,jy,jz]**4-300**4)) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 
                elif jz == cz-1:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz-1] + 2*di*emissivity*sigma/k*(T[n,jx,jy,jz]**4-300**4)) 
                elif jz == 0:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz-1] - 2*di*emissivity*sigma/k*(T[n,jx,jy,jz]**4-300**4) + 2*di*I[n,jx,jy,0])