|
|
本帖最后由 丁小白 于 2023-5-21 18:21 编辑
import numpy as np
import matplotlib.pyplot as plt
# 设置初始条件
L = 1.0 # 区域长度
W = 1.0 # 区域宽度
Nx = 100 # x方向网格数
Ny = 100 # y方向网格数
dx = L / (Nx) # x方向网格间距
dy = W / (Ny) # y方向网格间距
nt = 22 # 迭代次数
sigma = 0.25 # 热传导系数
dt = sigma * dx * dy / 2 # 时间步长
x = np.linspace(0, L, Nx)
y = np.linspace(0, W, Ny)
# 生成网格点坐标矩阵
X, Y = np.meshgrid(x, y)
# 设置初始温度
T0 = np.zeros((Ny, Nx))
T0[:, :] = 0
# 设置边界条件
T0[:, 0] =0 # x=0边界
T0[:, -1] = 0 # x=L边界
T0[0, :] = 0 # y=0边界
T0[-1, :] = 100 # y=W边界
# 定义克兰克-尼科尔方法对应的系数矩阵
A = np.zeros((Nx-2, Ny-2))
np.fill_diagonal(A, -1 - 2*sigma*dt/(dx**2))
np.fill_diagonal(A[1:], sigma*dt/(dx**2))
np.fill_diagonal(A[:, 1:], sigma*dt/(dx**2))
At=np.matrix(A)
print(At)
for n in range(nt):
b = T0[1:-1,1:-1] + (2*sigma*dt/(dx**2))*(T0[2:,1:-1] - 2 * T0[1:-1,1:-1] + T0[:98,1:-1] + T0[1:-1, 2:] - 2 * T0[1:-1, 1:-1] + T0[1:-1, :98])
bt=np.matrix(b)
T0[1:-1,1:-1]=np.linalg.solve(A,b)
print(T0)
# 绘制温度分布
fig, ax = plt.subplots()
im = ax.imshow(T0, extent=(0, L, 0, W), cmap='coolwarm', origin='lower')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Temperature distribution')
plt.colorbar(im)
plt.show()
|
|
