用Python（NumPy）解决热传导方程。

Question

用Python（NumPy）解决热传导方程。

8

我解决了金属棒的热方程，其中一端保持在100°C，另一端保持在0°C。

import numpy as np
import matplotlib.pyplot as plt

dt = 0.0005
dy = 0.0005
k = 10**(-4)
y_max = 0.04
t_max = 1
T0 = 100

def FTCS(dt,dy,t_max,y_max,k,T0):
    s = k*dt/dy**2
    y = np.arange(0,y_max+dy,dy) 
    t = np.arange(0,t_max+dt,dt)
    r = len(t)
    c = len(y)
    T = np.zeros([r,c])
    T[:,0] = T0
    for n in range(0,r-1):
        for j in range(1,c-1):
            T[n+1,j] = T[n,j] + s*(T[n,j-1] - 2*T[n,j] + T[n,j+1]) 
    return y,T,r,s

y,T,r,s = FTCS(dt,dy,t_max,y_max,k,T0)

plot_times = np.arange(0.01,1.0,0.01)
for t in plot_times:
    plt.plot(y,T[t/dt,:])

如果将纽曼边界条件更改为一端绝缘（非通量），

那么，计算术语如何处理

T[n+1,j] = T[n,j] + s*(T[n,j-1] - 2*T[n,j] + T[n,j+1])

应该被修改吗？

- Googlebot

1

这是一个编程网站，不是物理学网站，因此您需要提供相关的方程式，无论是您展示的代码还是您想要的代码。我们中的许多人不明白“一端绝缘（非通量）”是什么意思。 - Rory Daulton

1

你的热方程式不正确。左手边应该是∂u/∂t，而不是∂u/∂x。否则你的热量根本不会随时间变化，也没有通量。 - dROOOze

1

关键字是"Neumann边界条件"。 - Andras Deak -- Слава Україні

1

由于您正在使用有限差分逼近，请参阅此处。您需要做的就是找出有限差分逼近中的边界条件，然后在有限差分逼近达到这些条件时将表达式替换为0。这不是一个真正的python或实现问题，因为在尝试编码之前，您还没有弄清楚FD离散化（或者至少没有在问题中显示更改的BC的离散化）。 - dROOOze

2

请查看scipy.integrate教程中的示例，它展示了具有Neumann（“无流量”）边界条件的PDE系统的空间离散化。（时间演化使用scipy.integrate.odeint解决；该教程是“线性方法”的一个示例。） - Warren Weckesser

显示剩余4条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user6655984 · Accepted Answer

典型的 Neumann 边界条件处理方法是想象一个“ghost point”（幽灵点），它在定义域之外一步，然后使用边界条件计算其值；接下来正常地进行（使用 PDE）内部网格点计算，包括 Neumann 边界。

幽灵点允许我们在边界处使用对称的有限差分逼近导数，即如果 y 是空间变量，则为 (T[n, j+1] - T[n, j-1]) / (2*dy)。不对称的逼近 (T[n, j] - T[n, j-1]) / dy 不涉及幽灵点，精度远低于对称逼近：引入的误差比离散化 PDE 本身的误差大一个数量级。

因此，当 j 是 T 的最大可能索引时，边界条件指出 "T[n, j+1]" 应被理解为 T[n, j-1]，下面就是这样做的。

for j in range(1, c-1):
    T[n+1,j] = T[n,j] + s*(T[n,j-1] - 2*T[n,j] + T[n,j+1])  # as before
j = c-1 
T[n+1, j] = T[n,j] + s*(T[n,j-1] - 2*T[n,j] + T[n,j-1])  # note the last term here