python - 二维传热 - 溢出错误

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 类型的极限。

您可以通过在 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]) 在 sd=10 和 st=1000 的情况下。

这不是 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