二维传热 - 溢出错误

Question

我正在尝试使用以下公式模拟 python 铝中的二维热传递：

dT/dt = K*(d^2T/d^2x + d^2T/d^2y)

来源： https://www.youtube.com/watch?v=V00p-TgQML0

python 代码适用于 dx = dy = 1 (mm)，但如果 dx 和 dy 变小，我会收到溢出错误，我不知道如何避免。

import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy
from matplotlib.animation import FuncAnimation
import time


x = 11
y = 11

sd = 1
nx = x*sd
ny = y*sd

dx = 1/float(sd)
dy = dx

#Initial Temperature
T0 = 25

# Source: https://en.wikipedia.org/wiki/Thermal_diffusivity  
# Aluminium Thermal diffusivity (mm**2/s) 
K = 97

#Time
t0 = 0
te = 1
st = 1000
dt = 1/float(st)


#Iterations
N = (te - t0)*st

T = np.zeros(shape=(nx, ny))

T[:,:] = T0

# Dirichlet Condition
T[nx/2,ny/2] = 1000


MM = []

for n in range(N):

    Tb = deepcopy(T)

    for i in range(nx):
        for j in range(ny):

            #Neumann Boundary Conditions
            #TOP
            if i == 0:
                T[i,j] = Tb[i+1,j]

            #RIGHT
            elif j == ny1-1 and i != 0 and i != nx1-1: 
                T[i,j] = Tb[i,j-1]

            #BOTTOM
            elif i == nx-1:  
                T[i,j] = Tb[i-1,j]

            #LEFT
            elif j==0 and i != 0 and i != nx-1:
                T[i,j] = Tb[i,j+1]

            else:

                T[i,j] = Tb[i,j] + K*(dt/dx**2)*(Tb[i+1,j]+Tb[i-1,j]+Tb[i,j+1]+Tb[i,j-1]-4*Tb[i,j])
            T[nx/2,ny/2] = 200

    MM.append(T.copy())

fig = plt.figure()
pcm = plt.pcolormesh(MM[0])
plt.colorbar()

# Function called to update the graphic
def step(i):
    if i >= len(MM): return 0
    pcm.set_array(MM[i].ravel())
    plt.title("Time: {0} s\n".format(i*dt))
    plt.draw()

anim = FuncAnimation(fig, step, interval=3)
plt.show()

我为了重现溢出错误将sd的值改为10（每个mm将被分成10个元素）。

Answer 1

基本上 T[i,j] 正在发散并达到非常高的 + 和 - 值，并且在某个点达到 double_scale 类型的极限。

在 sd=10 和 st=1000 的情况下，您可以通过在 T[i,j] = Tb[i,j] + K*(dt/dx**2)*(Tb[i+1,j]+Tb[i-1,j]+Tb[i,j+1]+Tb[i,j-1]-4*Tb[i,j]) 之后添加 print(T[i,j]) 来检查。

这不是 python 或 numpy 问题，而是在尝试使用您正在使用的方法对刚性微分方程进行数值求解时出现的数值数学问题。

当您决定以更高的空间分辨率求解时，您也应该以更高的时间分辨率求解。我测试了代码，它适用于：sd=2、st=5000 和 sd=4、st=10000。你看到了模式。

或者为您的微分方程使用更好的数值解。就像向后微分公式 (BDF) 一样，您可以在其中采取更大的时间步长而不会导致数值求解器发散。在这里寻找灵感：

https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.ode.html

Answer 2

算法的核心是

T[i,j] = Tb[i,j] + K*(dt/dx**2)*(Tb[i+1,j]+Tb[i-1,j]+Tb[i,j+1]+Tb[i,j-1]-4*Tb[i,j])

注意这里Tb[i,j]的系数：是1 - 4*K*(dt/dx**2)。为了使算法起作用，这必须是正数；否则你从火中创造冰（从正到负），解决方案没有物理意义并且数字爆炸。

因此，请确保 K*(dt/dx**2) 小于 1/4。这意味着当 dx 变小时，dt 可能也需要减小。例如，K=1 和 dx=0.001 将要求 dt < 2.5e-7.

二维传热 - 溢出错误

2D Heat Transfer - Overflow Error

python

simulation

numpy