一维无热源传导代码实现

一维无热源传导的代码实现(均匀网格)

问题描述： 有一个均匀圆棒，导热系数为k=1000W/(m.k) 横截面积S=0.01$m^2$，长度为0.5m。

两端温度分别为100°和500°，求解稳态下棒内的温度分布

控制方程$\frac{\partial }{\partial x}(\frac{\partial T}{\partial x})=0$

求解出来解析解为T(x)=800x+100

import numpy as np

if __name__ == '__main__':
    nx = 5  # 定义离散的单元个数
    L = 0.5
    k = 1000
    S = 0.01
    Ta = 100
    Tb = 500

    dx = L / nx  # 网格间距
    A = np.zeros((nx, nx))  # 创建系数矩阵
    B = np.zeros(nx)        # 源项

    g = k * S / dx
    # 对于2 3 4单元其为内部网格，因此a_E = -kS/dx  a_W = -kS/dx   a_C = -(a_E+a_W)  b_c = 0

    # 对于1单元，a_W=0  fluxC_a = kS/(dx/2) = 2kS/dx  a_E = -kS/dx
    # 左侧还有一个源项fluxV_a =  -2kSTa/dx
    # a_C = -a_E + fluxC_a
    # b_c = -fluxV_a

    # 组装系数矩阵
    for i in range(nx):
        if i == 0:  # 第一个单元(左侧单元)
            A[i, i+1] = -g  # 东侧面
            A[i, i] = 3 * g # 当前网格C的系数  因为左侧是2*g 右侧为g
            B[i] = 2 * g * Ta
        elif i == nx-1:  # 最后一个单元(右侧单元)
            A[i, i-1] = -g # 西侧面
            A[i, i] = 3 * g
            B[i] = 2 * g * Tb
        else:  # 中间单元
            A[i, i-1] = A[i, i+1] = -g
            A[i, i] = 2 * g

    T = np.linalg.solve(A, B)  # 求解
    x = np.linspace(L/nx/2, L-L/nx/2, nx)
    y = 800 * x + 100
    plt.plot(x, y, 'bo-', label='Analytical solution')
    plt.plot(x, T, 'rs', label='Numerical solution')
    plt.legend(loc='best')
    plt.show()

结果比对：

一维无热源传导的代码实现(非均匀网格)

import numpy as np
import matplotlib.pyplot as plt

if __name__ == '__main__':
    L = 0.5
    k = 1000
    S = 0.01
    Ta = 100
    Tb = 500

    # 网格体心点的坐标
    cellx = np.array([0.07, 0.12, 0.18, 0.24, 0.26, 0.37, 0.48])
    nx = len(cellx)  # 网格数量
    A = np.zeros((nx, nx))  # 创建系数矩阵
    B = np.zeros(nx)  # 源项
    g = k * S
    # 组装系数矩阵
    for i in range(nx):
        if i == 0:  # 第一个单元(左侧单元)
            A[i, i + 1] = -g / (cellx[i+1] - cellx[i])  # 东侧面
            A[i, i] = g / (cellx[i+1] - cellx[i]) + g / cellx[i]  # 东侧面加上西侧边界面
            B[i] = g / cellx[i] * Ta
        elif i == nx - 1:  # 最后一个单元(右侧单元)
            A[i, i - 1] = -g / (cellx[i] - cellx[i-1])  # 西侧面
            A[i, i] = g / (cellx[i] - cellx[i-1]) + g / (L - cellx[i])  # 西侧面加上东侧面
            B[i] = g / (L - cellx[i]) * Tb
        else:  # 中间单元
            A[i, i - 1] = -g / (cellx[i] - cellx[i-1])
            A[i, i + 1] = -g / (cellx[i+1] - cellx[i])
            A[i, i] = -(A[i, i-1] + A[i, i+1])
    T = np.linalg.solve(A, B)  # 求解
    y = 800 * cellx + 100
    plt.plot(cellx, y, 'bo-', label='Analytical solution')
    plt.plot(cellx, T, 'rs', label='Numerical solution')
    plt.legend(loc='best')
    plt.show()

结果比对：