二维传热问题(python实现)

二维传热问题

如图所示为一个厚度为1cm的板。板材的导热系数为k=1000W/m.k。西侧边界施加500kW/$m^2$的稳定热通量，南侧边界和东侧边界为绝热，北侧边界的温度保持在100°。计算稳定状态下板的温度分布。
控制方程:
$$\frac{\partial }{\partial x}(k\frac{\partial T}{\partial x})+\frac{\partial }{\partial y}(k\frac{\partial T}{\partial y})=0$$
代码实现，使用的直接法求解(LU分解)：

import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from scipy import interpolate


# 绘制云图
def plot_coutour(Coordx, Coordy, T, minX, maxX, minY, maxY):
    X = np.linspace(minX, maxX, 1000)
    Y = np.linspace(minY, maxY, 1000)
    # 生成二维数据坐标点
    X1, Y1 = np.meshgrid(X, Y)
    # 通过griddata函数插值得到所有的(X1，Y1)处对应的值，原始数据为Coordx,Coordy（原始控制点）,T（原始控制点对应的T）
    Z = interpolate.griddata((Coordx, Coordy), T, (X1, Y1), method="cubic")
    fig, ax = plt.subplots(figsize=(6, 8))  # 创建一个6x8大小的子图

    # level设置云图颜色范围以及颜色梯度划分，每间隔5颜色变化，为显示完全，将上限提高10(也可以变小点，这样上限不需要提升太高)
    levels = range((int)(T.min()), (int)(T.max())+10, 5)

    # 设置cmap为jet，最小值为蓝色，最大值为红色
    # 绘制轮廓图
    cset1 = ax.contourf(X1, Y1, Z, levels, cmap=cm.jet)

    # 设置云图坐标范围以及坐标轴标签
    ax.set_xlim(minX, maxX)
    ax.set_ylim(minY, maxY)
    ax.set_xlabel("X(mm)", size=15)
    ax.set_ylabel("Y(mm)", size=15)

    # 设置colorbar的刻度，标签
    cbar = fig.colorbar(cset1)
    cbar.set_label('T', size=18)

if __name__ == '__main__':
    k = 1000      # 导热系数
    q = 500e3     # 西侧边界通量
    T_top = 100   # 北侧边界温度

    L = 0.3  # 长度
    H = 0.4  # 高度

    nx = 100  # x方向网格个数
    ny = 100  # y方向网格个数

    dx = L / nx  # 网格大小
    dy = H / ny

    A = np.zeros((nx * ny, nx * ny))  # 构造系数矩阵
    B = np.zeros((nx * ny))           # 构造源项矩阵

    gdiff_e = gdiff_w = dy / dx
    gdiff_n = gdiff_s = dx / dy

    for i in range(nx * ny): # 填充nx * ny行系数矩阵
        # 内部单元（去除底部、左侧、顶部以及右侧单元）
        if np.floor(i / nx) != 0 and i % nx != 0 and np.floor(i / nx) != ny - 1 and (i + 1) % nx != 0:
            A[i, i - 1] = -k * gdiff_w  # 西侧面
            A[i, i + 1] = -k * gdiff_e  # 东侧面
            A[i, i + nx] = -k * gdiff_n  # 北侧面
            A[i, i - nx] = -k * gdiff_s  # 南侧面
            A[i, i] = -(A[i, i-1] + A[i, i+1] + A[i, i+nx] + A[i, i-nx])

        # 底部单元，不考虑角点
        elif np.floor(i / nx) == 0 and i != 0 and i != nx - 1:
            A[i, i - 1] = -k * gdiff_w  # 西侧面
            A[i, i + 1] = -k * gdiff_e  # 东侧面
            A[i, i + nx] = -k * gdiff_n  # 北侧面
            A[i, i] = -(A[i, i-1] + A[i, i+1] + A[i, i+nx])  # 因为底部边界是绝热边界所以fluxCb=0 fluxV_b=0（第二类边界条件）

        # 右侧单元,不考虑角点
        elif (i+1) % nx == 0 and i != nx - 1 and i != nx * ny - 1:
            A[i, i - 1] = -k * gdiff_w  # 西侧面
            A[i, i + nx] = -k * gdiff_n  # 北侧面
            A[i, i - nx] = -k * gdiff_s  # 南侧面
            A[i, i] = -(A[i, i - 1] +A[i, i + nx] + A[i, i - nx])  # 因为右侧边界是绝热边界所以fluxCb=0 fluxVb=0（第二类边界条件）

        # 左侧单元， 不考虑角点
        elif i % nx == 0 and i != 0 and i != (ny - 1) * nx:
            A[i, i + 1] = -k * gdiff_e  # 东侧面
            A[i, i + nx] = -k * gdiff_n  # 北侧面
            A[i, i - nx] = -k * gdiff_s  # 南侧面
            A[i, i] = -(A[i, i + 1] + A[i, i + nx] + A[i, i - nx])
            # （第二类边界条件）源项多了fluxVb = q*dy
            B[i] = q * dy

        # 顶部，不考虑角点
        elif np.floor(i / nx) == ny - 1 and i != nx * ny -1 and i != (ny - 1) * nx:
            A[i, i - 1] = -k * gdiff_w  # 西侧面
            A[i, i + 1] = -k * gdiff_e  # 东侧面
            A[i, i - nx] = -k * gdiff_s  # 南侧面
            A[i, i] = -(A[i, i - 1] + A[i, i + 1] + A[i, i - nx]) + 2 * k * gdiff_n   # 第一类边界条件 边界项为k*dx/(dy/2)
            B[i] = 2 * k * gdiff_n * T_top

        # 左下角
        elif i == 0:
            A[i, i + 1] = -k * gdiff_e  # 东侧面
            A[i, i + nx] = -k * gdiff_n  # 北侧面
            A[i, i] = -(A[i, i + 1] + A[i, i + nx])
            B[i] = q * dy

        # 右下角
        elif i == nx - 1:
            A[i, i - 1] = -k * gdiff_w  # 西侧面
            A[i, i + nx] = -k * gdiff_n  # 北侧面
            A[i, i] = -(A[i, i - 1] + A[i, i + nx])

        # 左上角
        elif i == nx * (ny - 1):
            A[i, i + 1] = -k * gdiff_e  # 东侧面
            A[i, i - nx] = -k * gdiff_s  # 南侧面
            A[i, i] = -(A[i, i + 1] + A[i, i - nx]) + 2 * k * gdiff_n
            B[i] = q * dy + 2 * k * gdiff_n * T_top

        # 右上角
        elif i == nx * ny - 1:
            A[i, i - 1] = -k * gdiff_w  # 西侧面
            A[i, i - nx] = -k * gdiff_s  # 南侧面
            A[i, i] = -(A[i, i - 1] + A[i, i - nx]) + 2 * k * gdiff_n
            B[i] = 2 * k * gdiff_n * T_top

    T = np.linalg.solve(A,B)  # LU分解法
    print(T)

    # 创建数组用于存储各单元体心坐标
    cellC = np.ones((nx * ny, 2))
    for i in range(len(cellC)):
        cellC[i, 0] = dx / 2 + (i % nx) * dx
        cellC[i, 1] = dy / 2 + (np.floor(i / nx)) * dy
    x = cellC[:, 0]
    y = cellC[:, 1]
    plot_coutour(x, y, T, 0, 0.3, 0, 0.4)
    plt.show()

最终结果：