传热学作业：python非稳态二维矩形导热计算

import matplotlib.pyplot as plt
import numpy as np

""" 声明变量 """
lx = 21 #x方向网格数量
ly = 21 #y方向网格数量
k, rho, cp = 237, 2700, 880  # 材料为金属铝:导热系数，密度。比热容
a = k / (rho * cp)  # 热扩散率
""" 空间微元 """
nx = 0.1 # x方向长0.1m
ny = 0.05 # y方向长0.05m

dx = nx /lx #x微元长度
dy = ny /ly #y微元长度
""" 时间微元 """
dt = 0.001  # 限制时间步长0.001s，防止计算发散
lt = 200 #每dt计算的时间步数

h = 1000000 #对流换热系数
tf = 20.0
""" 边界条件 """
T_top = 0.0
T_left = 10.0
T_bottom = 0.0
T_right = 0.0
Fo = (a * dt) / (dx ** 2)
""" 初始化 """
T_initial = 0.0
T = T_initial*np.ones((ly, lx))#初始温度为10℃
T[:1, :] = T_bottom
T[:, :1] = T_left
T[:, -1:] = T_right
T[-1:, :] = T_top
eps = 1e-5 #残差
iterStep = 20 #计算步数
def plot_contour(Tem,lx,ly,Name,levelsYmin=0, levelsYmax=0):#定义显示云图
    levelsminY = np.min(Tem)  # 找出数组的最小值
    levelsmaxY = np.max(Tem) # 找出数组的最大值
    fig, ax = plt.subplots(figsize=(8, 8)) #用来创建总画布
    x = np.linspace(0, 0.2, lx)#定义x轴
    y = np.linspace(0, 0.2, ly)#定义y轴
    X, Y = np.meshgrid(x, y) #生成网格
    levels = np.linspace(levelsminY, levelsmaxY, 21)#生成y轴从小到大的数列
    cset1 = plt.contourf(X, Y, Tem, levels, cmap=plt.cm.jet,zorder=1) #云图
    ##cset = plt.contour(X, Y, Tem, levels, cmap=plt.cm.jet,zorder=1) #等值线
    plt.title('Time = ' + str(Name) + 's', size=20)
    plt.xlabel('X(mm)',fontdict={'family' : 'Times New Roman', 'size':20},labelpad=5, loc='center', color='green')
    plt.ylabel('Y(mm)',fontdict={'family' : 'Times New Roman', 'size':20},labelpad=5, loc='center', color='red')
    cbarbar = fig.colorbar(cset1)  # 设置图例
    cbarbar.set_label('Temperature /℃', fontsize=12)# 设置图例文字和大小

    plt.show()

def thermal_conducitivity(T):
    for k in range(lt):#时间循环
        T_kn=T.copy() #储存上一个时间步数的数据
        for i in range(0, lx-1, 1):
            for j in range(0, ly-1, 1):
                T_k = T.copy()
             #中心节点
                T[i][j] = T_kn[i][j]+Fo * (T[i + 1][j] - 2*T[i][j] + T[i-1][j] + T[i][j-1]) + Fo * (T[i][j+1] - 2*T[i][j] + T[i][j-1])
              #下边缘为第三类边界条件
                T[0][j]= T_kn[0][j] + (h * dx * (tf - T[0][j]) - k* (T[0][j] - T[0][j+1]) * dy / (2 * dx) -k* (T[0][j]- T[0][j+1]) * dy / (2 * dx) -k* (T[0][j+1]- T[0][j-1]) * dx / dy) * dt / (rho * cp * 0.5 * dx * dy)
              # 上边缘为第三类边界条件
                T[lx-1][j] = T_kn[lx-1][j] + (h * dx * (tf - T[lx-1][j]) - k * (T[lx-1][j] - T[lx-1][j + 1]) * dy / (2 * dx) - k * (T[lx-1][j] - T[lx-1][j + 1]) * dy / (2 * dx) - k * (T[lx-1][j + 1] - T[0][j - 1]) * dx / dy) * dt / (rho * cp * 0.5 * dx * dy)
                # 残差
                error = np.max(np.abs(T - T_k))
                if error < eps:
                   break
                   print('时间： ', (k + 1) * dt,'s,第',k + 1,'次计算')
                # 计算达到设定计算步数
                if i == iterStep - 1:
                   print('时间： ', (k + 1) * dt,'s,第',k + 1,'次计算')

    #plot_contour(T, lx, ly, lt*dt,300, 900)#每一时刻输出一个云图
    return T





T=thermal_conducitivity(T)
plot_contour(T, lx, ly, lt*dt,300, 900)

print(T)
print('点温度：', T[15][20], ' K', '时间：', k* dt, ' s')#输出点[15][20]的温度和时间