平板导热系数计算代码循环仅统计首尾三个元素问题求助

核心问题原因 #

• 离散方程缺失关键项：一维非稳态导热控制方程的空间二阶偏导差分形式需要除以空间步长的平方h²，现有代码直接省略了该项，导致差分项的计算系数仅为9.16e-6量级，迭代过程中温度增量几乎可忽略，只有靠近边界的区域因边界温度与初始温度差足够大，能产生可观测的数值变化。
• 迭代逻辑错误：当前时间步的温度计算应该以上一个时间步的同位置温度为基础迭代，现有代码以当前时间步的初始值350为基础计算，没有继承上一步的迭代结果，本质上没有完成时间维度的递进计算。
• 稳定性风险：显式差分求解导热问题需要满足稳定性判据Fo = lam*tau/(rho*c*h²) ≤ 0.5，修正前两个问题后如果步长不满足该条件，会出现数值振荡，结果无物理意义。

修正后的参考代码 #

timend = 8
H = 0.10
# 调整空间、时间步长满足稳定性条件
n, m = 20, 20000
h = H / (n-1)
tau = timend / m
t0, th, ts = 550, 700, 350
rho, c, lam = 8800, 381, 384
x = [[350 for _ in range(n)] for _ in range(m)]

# 初始边界条件
x[0][0] = t0
x[0][n-1] = th
for i in range(1, m):
    x[i][0], x[i][n-1] = t0, th
    for j in range(1, n-1):
        # 补上除以h²的项，使用上一时间步的温度计算当前步温度
        x[i][j] = x[i-1][j] + tau * lam / rho / c / (h**2) * (x[i-1][j+1] - 2 * x[i-1][j] + x[i-1][j-1])

print(x[m-1])

内容的提问来源于stack exchange，提问作者One-eight greek